From 587a545494f2e9538688acc78c7fe900665b7f15 Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Tue, 27 Oct 2009 14:00:56 +0100 Subject: [PATCH] Replaced doubles by Reals, fixed two bugs * Replaced all doubles by Reals where necessary * Replaced constructor TFactor::TFactor( const VarSet& vars, TIterator begin ) by two constructors TFactor::TFactor( const VarSet& vars, const T* p ) TFactor::TFactor( const VarSet& vars, const std::vector &x ) * Fixed bugs (min and max were reversed) in TFactor min( const TFactor &, const TFactor & ) TFactor max( const TFactor &, const TFactor & ) --- ChangeLog | 8 +++ examples/example.cpp | 6 +-- examples/example_sprinkler.cpp | 5 +- include/dai/bbp.h | 10 ++-- include/dai/bp.h | 20 +++---- include/dai/cbp.h | 32 ++++++------ include/dai/daialg.h | 4 +- include/dai/exactinf.h | 4 +- include/dai/factor.h | 24 ++++++--- include/dai/gibbs.h | 4 +- include/dai/hak.h | 14 ++--- include/dai/jtree.h | 6 +-- include/dai/lc.h | 14 ++--- include/dai/mf.h | 10 ++-- include/dai/mr.h | 40 +++++++------- include/dai/prob.h | 4 +- include/dai/regiongraph.h | 16 +++--- include/dai/treeep.h | 14 ++--- include/dai/util.h | 16 +++--- src/bbp.cpp | 30 +++++------ src/bp.cpp | 14 ++--- src/cbp.cpp | 18 +++---- src/daialg.cpp | 4 +- src/exactinf.cpp | 2 +- src/factorgraph.cpp | 22 ++++---- src/gibbs.cpp | 6 +-- src/hak.cpp | 16 +++--- src/jtree.cpp | 4 +- src/lc.cpp | 18 +++---- src/mf.cpp | 6 +-- src/mr.cpp | 96 +++++++++++++++++----------------- src/properties.cpp | 2 + src/regiongraph.cpp | 4 +- src/treeep.cpp | 12 ++--- src/util.cpp | 2 +- tests/testbbp.cpp | 8 +-- tests/testdai.cpp | 31 +++++------ utils/createfg.cpp | 66 +++++++++++------------ utils/fginfo.cpp | 6 +-- 39 files changed, 318 insertions(+), 300 deletions(-) diff --git a/ChangeLog b/ChangeLog index 17dc97e..ccafce5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +* Replaced constructor + TFactor::TFactor( const VarSet& vars, TIterator begin ) + by two constructors + TFactor::TFactor( const VarSet& vars, const T* p ) + TFactor::TFactor( const VarSet& vars, const std::vector &x ) +* Fixed bugs (min and max were reversed) in + TFactor min( const TFactor &, const TFactor & ) + TFactor max( const TFactor &, const TFactor & ) * Renamed RegionGraph::Check_Counting_Numbers into RegionGraph::checkCountingNumbers (and provided an alias for backwards compatibility) diff --git a/examples/example.cpp b/examples/example.cpp index e8ed770..210a293 100644 --- a/examples/example.cpp +++ b/examples/example.cpp @@ -29,9 +29,9 @@ int main( int argc, char *argv[] ) { fg.ReadFromFile(argv[1]); // Set some constants - size_t maxiter = 10000; - double tol = 1e-9; - size_t verb = 1; + size_t maxiter = 10000; + Real tol = 1e-9; + size_t verb = 1; // Store the constants in a PropertySet object PropertySet opts; diff --git a/examples/example_sprinkler.cpp b/examples/example_sprinkler.cpp index 85df931..0d1e570 100644 --- a/examples/example_sprinkler.cpp +++ b/examples/example_sprinkler.cpp @@ -74,13 +74,12 @@ int main() { // Calculate joint probability of all four variables Factor P; - for( size_t I = 0; I < SprinklerNetwork.nrFactors(); I++ ) { + for( size_t I = 0; I < SprinklerNetwork.nrFactors(); I++ ) P *= SprinklerNetwork.factor( I ); - } // P.normalize(); // Not necessary: a Bayesian network is already normalized by definition // Calculate some probabilities - double denom = P.marginal( W )[1]; + Real denom = P.marginal( W )[1]; cout << "P(W=1) = " << denom << endl; cout << "P(S=1 | W=1) = " << P.marginal( VarSet( S, W ) )[3] / denom << endl; cout << "P(R=1 | W=1) = " << P.marginal( VarSet( R, W ) )[3] / denom << endl; diff --git a/include/dai/bbp.h b/include/dai/bbp.h index 9dc9132..ddad9c5 100644 --- a/include/dai/bbp.h +++ b/include/dai/bbp.h @@ -271,10 +271,10 @@ class BBP { size_t maxiter; /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD) - double tol; + Real tol; /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09] - double damping; + Real damping; /// Update schedule UpdateType updates; @@ -294,9 +294,9 @@ class BBP { /// Maximum number of iterations size_t maxiter; /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD) - double tol; + Real tol; /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09] - double damping; + Real damping; /// Update schedule UpdateType updates; @@ -339,7 +339,7 @@ Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector * \param cfn Cost function to be used. * \param h controls size of perturbation. */ -double numericBBPTest( const InfAlg &bp, const std::vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ); +Real numericBBPTest( const InfAlg &bp, const std::vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h ); } // end of namespace dai diff --git a/include/dai/bp.h b/include/dai/bp.h index 2e985c4..c571b4c 100644 --- a/include/dai/bp.h +++ b/include/dai/bp.h @@ -32,18 +32,18 @@ namespace dai { class BP : public DAIAlgFG { private: typedef std::vector ind_t; - typedef std::multimap > LutType; + typedef std::multimap > LutType; struct EdgeProp { ind_t index; Prob message; Prob newMessage; - double residual; + Real residual; }; std::vector > _edges; std::vector > _edge2lut; LutType _lut; /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Number of iterations needed size_t _iters; /// The history of message updates (only recorded if recordSentMessages is true) @@ -65,13 +65,13 @@ class BP : public DAIAlgFG { size_t maxiter; /// Tolerance - double tol; + Real tol; /// Do updates in logarithmic domain? bool logdomain; /// Damping constant - double damping; + Real damping; /// Update schedule UpdateType updates; @@ -133,8 +133,8 @@ class BP : public DAIAlgFG { virtual Real logZ() const; virtual void init(); virtual void init( const VarSet &ns ); - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} @@ -167,12 +167,12 @@ class BP : public DAIAlgFG { const Prob & newMessage(size_t i, size_t _I) const { return _edges[i][_I].newMessage; } ind_t & index(size_t i, size_t _I) { return _edges[i][_I].index; } const ind_t & index(size_t i, size_t _I) const { return _edges[i][_I].index; } - double & residual(size_t i, size_t _I) { return _edges[i][_I].residual; } - const double & residual(size_t i, size_t _I) const { return _edges[i][_I].residual; } + Real & residual(size_t i, size_t _I) { return _edges[i][_I].residual; } + const Real & residual(size_t i, size_t _I) const { return _edges[i][_I].residual; } void calcNewMessage( size_t i, size_t _I ); void updateMessage( size_t i, size_t _I ); - void updateResidual( size_t i, size_t _I, double r ); + void updateResidual( size_t i, size_t _I, Real r ); void findMaxResidual( size_t &i, size_t &_I ); /// Calculates unnormalized belief of variable void calcBeliefV( size_t i, Prob &p ) const; diff --git a/include/dai/cbp.h b/include/dai/cbp.h index 8fb446b..d9032ef 100644 --- a/include/dai/cbp.h +++ b/include/dai/cbp.h @@ -49,10 +49,10 @@ class CBP : public DAIAlgFG { /// Factor beliefs std::vector _beliefsF; /// Log-partition sum - double _logZ; + Real _logZ; /// Counts number of clampings at each leaf node - double _sum_level; + Real _sum_level; /// Number of leaves of recursion tree size_t _num_leaves; @@ -71,8 +71,8 @@ class CBP : public DAIAlgFG { * beliefs estimates of the children, and returns the improved * estimates in \a lz_out and \a beliefs_out to its parent */ - void runRecurse( InfAlg *bp, double orig_logZ, std::vector clamped_vars_list, size_t &num_leaves, - size_t &choose_count, double &sum_level, Real &lz_out, std::vector &beliefs_out ); + void runRecurse( InfAlg *bp, Real orig_logZ, std::vector clamped_vars_list, size_t &num_leaves, + size_t &choose_count, Real &sum_level, Real &lz_out, std::vector &beliefs_out ); /// Choose the next variable to clamp /** Choose the next variable to clamp, given a converged InfAlg (\a bp), @@ -91,12 +91,12 @@ class CBP : public DAIAlgFG { /// Numer of iterations needed size_t _iters; /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Sets variable beliefs, factor beliefs and logZ /** \param bs should be a concatenation of the variable beliefs followed by the factor beliefs */ - void setBeliefs( const std::vector &bs, double logZ ); + void setBeliefs( const std::vector &bs, Real logZ ); /// Constructor helper function void construct(); @@ -121,8 +121,8 @@ class CBP : public DAIAlgFG { virtual Real logZ() const { return _logZ; } virtual void init() {}; virtual void init( const VarSet & ) {}; - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} @@ -146,18 +146,18 @@ class CBP : public DAIAlgFG { size_t verbose = 0; /// Tolerance to use in BP - double tol; + Real tol; /// Update style for BP UpdateType updates; /// Maximum number of iterations for BP size_t maxiter; /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF) - double rec_tol; + Real rec_tol; /// Maximum number of levels of recursion (\a recurse is REC_FIXED) size_t max_levels = 10; /// If choose==CHOOSE_BBP and maximum adjoint is less than this value, don't recurse - double min_max_adj; + Real min_max_adj; /// Heuristic for choosing clamping variable ChooseMethodType choose; /// Method for deciding when to stop recursing @@ -190,17 +190,17 @@ class CBP : public DAIAlgFG { /// Verbosity size_t verbose; /// Tolerance to use in BP - double tol; + Real tol; /// Update style for BP UpdateType updates; /// Maximum number of iterations for BP size_t maxiter; /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF) - double rec_tol; + Real rec_tol; /// Maximum number of levels of recursion (\a recurse is REC_FIXED) size_t max_levels; /// If choose==CHOOSE_BBP and maximum adjoint is less than this value, don't recurse - double min_max_adj; + Real min_max_adj; /// Heuristic for choosing clamping variable ChooseMethodType choose; /// Method for deciding when to stop recursing @@ -234,9 +234,9 @@ class CBP : public DAIAlgFG { /// Returns maximum number of levels of recursion size_t maxClampLevel() { return props.max_levels; } /// Returns props.min_max_adj @see CBP::Properties::min_max_adj - double minMaxAdj() { return props.min_max_adj; } + Real minMaxAdj() { return props.min_max_adj; } /// Returns tolerance used for controlling recursion depth - double recTol() { return props.rec_tol; } + Real recTol() { return props.rec_tol; } }; diff --git a/include/dai/daialg.h b/include/dai/daialg.h index 067eb13..1833e5d 100644 --- a/include/dai/daialg.h +++ b/include/dai/daialg.h @@ -85,7 +85,7 @@ class InfAlg { /* Before run() is called the first time, init() should be called. * If run() returns successfully, the results can be queried using the methods belief(), beliefs() and logZ(). */ - virtual double run() = 0; + virtual Real run() = 0; /// Clamp variable with index i to value x (i.e. multiply with a Kronecker delta \f$\delta_{x_i, x}\f$) /** If backup == true, make a backup of all factors that are changed @@ -101,7 +101,7 @@ class InfAlg { /// Return maximum difference between single node beliefs in the last pass /// \throw Exception if not implemented/supported - virtual double maxDiff() const = 0; + virtual Real maxDiff() const = 0; /// Return number of passes over the factorgraph /// \throw Exception if not implemented/supported diff --git a/include/dai/exactinf.h b/include/dai/exactinf.h index ae4aef7..7a0de7e 100644 --- a/include/dai/exactinf.h +++ b/include/dai/exactinf.h @@ -65,8 +65,8 @@ class ExactInf : public DAIAlgFG { virtual Real logZ() const { return _logZ; } virtual void init(); virtual void init( const VarSet &/*ns*/ ) { DAI_THROW(NOT_IMPLEMENTED); } - virtual double run(); - virtual double maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } + virtual Real run(); + virtual Real maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } virtual size_t Iterations() const { DAI_THROW(NOT_IMPLEMENTED); return 0; } //@} diff --git a/include/dai/factor.h b/include/dai/factor.h index e058bce..70042ba 100644 --- a/include/dai/factor.h +++ b/include/dai/factor.h @@ -98,13 +98,21 @@ template class TFactor { /// Constructs factor depending on variables in \a vars with all values set to \a p TFactor( const VarSet& vars, T p ) : _vs(vars), _p(_vs.nrStates(),p) {} - /// Constructs factor depending on variables in \a vars, copying the values from a range - /** \tparam Iterator Iterates over instances of type \a T; should support addition of \c size_t. + /// Constructs factor depending on variables in \a vars, copying the values from a std::vector<> + /** \tparam S Type of values of \a x * \param vars contains the variables that the new factor should depend on. - * \param begin Points to first value to be added. + * \param x Vector with values to be copied. */ - template - TFactor( const VarSet& vars, TIterator begin ) : _vs(vars), _p(begin, begin + _vs.nrStates(), _vs.nrStates()) {} + template + TFactor( const VarSet& vars, const std::vector &x ) : _vs(vars), _p(x.begin(), x.begin() + _vs.nrStates(), _vs.nrStates()) { + DAI_ASSERT( x.size() == vars.nrStates() ); + } + + /// Constructs factor depending on variables in \a vars, copying the values from an array + /** \param vars contains the variables that the new factor should depend on. + * \param p Points to array of values to be added. + */ + TFactor( const VarSet& vars, const T* p ) : _vs(vars), _p(p, p + _vs.nrStates(), _vs.nrStates()) {} /// Constructs factor depending on variables in \a vars, copying the values from \a p TFactor( const VarSet& vars, const TProb &p ) : _vs(vars), _p(p) { @@ -539,7 +547,7 @@ template TFactor pointwiseOp( const TFactor result[i] = op( result[i], g[i] ); return result; } else { - TFactor result( f.vars() | g.vars(), 0.0 ); + TFactor result( f.vars() | g.vars(), (T)0 ); IndexFor i1(f.vars(), result.vars()); IndexFor i2(g.vars(), result.vars()); @@ -584,7 +592,7 @@ template T dist( const TFactor &f, const TFactor &g, Prob::Dis */ template TFactor max( const TFactor &f, const TFactor &g ) { DAI_ASSERT( f._vs == g._vs ); - return TFactor( f._vs, min( f.p(), g.p() ) ); + return TFactor( f._vs, max( f.p(), g.p() ) ); } @@ -594,7 +602,7 @@ template TFactor max( const TFactor &f, const TFactor &g ) */ template TFactor min( const TFactor &f, const TFactor &g ) { DAI_ASSERT( f._vs == g._vs ); - return TFactor( f._vs, max( f.p(), g.p() ) ); + return TFactor( f._vs, min( f.p(), g.p() ) ); } diff --git a/include/dai/gibbs.h b/include/dai/gibbs.h index c1e94d3..d33f2b8 100644 --- a/include/dai/gibbs.h +++ b/include/dai/gibbs.h @@ -70,8 +70,8 @@ class Gibbs : public DAIAlgFG { virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } virtual void init(); virtual void init( const VarSet &/*ns*/ ) { init(); } - virtual double run(); - virtual double maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } + virtual Real run(); + virtual Real maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } virtual size_t Iterations() const { return props.iters; } //@} diff --git a/include/dai/hak.h b/include/dai/hak.h index 96a2f00..4d17923 100644 --- a/include/dai/hak.h +++ b/include/dai/hak.h @@ -38,7 +38,7 @@ class HAK : public DAIAlgRG { std::vector > _muab; std::vector > _muba; /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Number of iterations needed size_t _iters; @@ -58,10 +58,10 @@ class HAK : public DAIAlgRG { size_t maxiter; /// Tolerance - double tol; + Real tol; /// Damping constant - double damping; + Real damping; /// How to choose the clusters ClustersType clusters; @@ -100,8 +100,8 @@ class HAK : public DAIAlgRG { virtual Real logZ() const; virtual void init(); virtual void init( const VarSet &ns ); - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} @@ -113,8 +113,8 @@ class HAK : public DAIAlgRG { const Factor& Qa( size_t alpha ) const { return _Qa[alpha]; }; const Factor& Qb( size_t beta ) const { return _Qb[beta]; }; - double doGBP(); - double doDoubleLoop(); + Real doGBP(); + Real doDoubleLoop(); //@} private: diff --git a/include/dai/jtree.h b/include/dai/jtree.h index 8107757..3e07356 100644 --- a/include/dai/jtree.h +++ b/include/dai/jtree.h @@ -37,7 +37,7 @@ namespace dai { class JTree : public DAIAlgRG { private: std::vector > _mes; - double _logZ; + Real _logZ; public: /// Rooted tree @@ -88,8 +88,8 @@ class JTree : public DAIAlgRG { virtual Real logZ() const; virtual void init() {} virtual void init( const VarSet &/*ns*/ ) {} - virtual double run(); - virtual double maxDiff() const { return 0.0; } + virtual Real run(); + virtual Real maxDiff() const { return 0.0; } virtual size_t Iterations() const { return 1UL; } //@} diff --git a/include/dai/lc.h b/include/dai/lc.h index f4cacda..8f642e4 100644 --- a/include/dai/lc.h +++ b/include/dai/lc.h @@ -41,7 +41,7 @@ class LC : public DAIAlgFG { std::vector _beliefs; /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Number of iterations needed size_t _iters; @@ -61,13 +61,13 @@ class LC : public DAIAlgFG { size_t maxiter; /// Tolerance - double tol; + Real tol; /// Complete or partial reinit of cavity graphs? bool reinit; /// Damping constant - double damping; + Real damping; /// How to initialize the cavities CavityType cavity; @@ -103,8 +103,8 @@ class LC : public DAIAlgFG { virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } virtual void init(); virtual void init( const VarSet &/*ns*/ ) { init(); } - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} @@ -112,8 +112,8 @@ class LC : public DAIAlgFG { /// @name Additional interface specific for LC //@{ - double CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts ); - double InitCavityDists( const std::string &name, const PropertySet &opts ); + Real CalcCavityDist( size_t i, const std::string &name, const PropertySet &opts ); + Real InitCavityDists( const std::string &name, const PropertySet &opts ); long SetCavityDists( std::vector &Q ); Factor NewPancake (size_t i, size_t _I, bool & hasNaNs); diff --git a/include/dai/mf.h b/include/dai/mf.h index 644ddd4..c5fca87 100644 --- a/include/dai/mf.h +++ b/include/dai/mf.h @@ -32,7 +32,7 @@ class MF : public DAIAlgFG { private: std::vector _beliefs; /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Number of iterations needed size_t _iters; @@ -46,10 +46,10 @@ class MF : public DAIAlgFG { size_t maxiter; /// Tolerance - double tol; + Real tol; /// Damping constant - double damping; + Real damping; } props; /// Name of this inference algorithm @@ -76,8 +76,8 @@ class MF : public DAIAlgFG { virtual Real logZ() const; virtual void init(); virtual void init( const VarSet &ns ); - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} diff --git a/include/dai/mr.h b/include/dai/mr.h index c772a7b..62da8b9 100644 --- a/include/dai/mr.h +++ b/include/dai/mr.h @@ -39,20 +39,20 @@ class MR : public DAIAlgFG { std::vector con; // con[i] = connectivity of spin i std::vector > nb; // nb[i] are the neighbours of spin i - std::vector > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j] - std::vector theta; // theta[i] is the local field on spin i - std::vector > M; // M[i][_j] is M^{(i)}_j + std::vector > tJ; // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j] + std::vector theta; // theta[i] is the local field on spin i + std::vector > M; // M[i][_j] is M^{(i)}_j std::vector > kindex; // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour - std::vector > > cors; + std::vector > > cors; static const size_t kmax = 31; typedef boost::dynamic_bitset<> sub_nb; size_t N; - std::vector Mag; + std::vector Mag; - double _maxdiff; + Real _maxdiff; size_t _iters; public: @@ -68,7 +68,7 @@ class MR : public DAIAlgFG { size_t verbose; /// Tolerance - double tol; + Real tol; /// Update equations UpdateType updates; @@ -98,8 +98,8 @@ class MR : public DAIAlgFG { virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; } virtual void init() {} virtual void init( const VarSet &/*ns*/ ) { DAI_THROW(NOT_IMPLEMENTED); } - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} @@ -109,25 +109,25 @@ class MR : public DAIAlgFG { //@} private: - void init(size_t Nin, double *_w, double *_th); + void init(size_t Nin, Real *_w, Real *_th); void makekindex(); void init_cor(); - double init_cor_resp(); + Real init_cor_resp(); void solvemcav(); void solveM(); - double _tJ(size_t i, sub_nb A); + Real _tJ(size_t i, sub_nb A); - double Omega(size_t i, size_t _j, size_t _l); - double T(size_t i, sub_nb A); - double T(size_t i, size_t _j); - double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2); - double Gamma(size_t i, size_t _l1, size_t _l2); + Real Omega(size_t i, size_t _j, size_t _l); + Real T(size_t i, sub_nb A); + Real T(size_t i, size_t _j); + Real Gamma(size_t i, size_t _j, size_t _l1, size_t _l2); + Real Gamma(size_t i, size_t _l1, size_t _l2); - double appM(size_t i, sub_nb A); - void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd); + Real appM(size_t i, sub_nb A); + void sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd); - double sign(double a) { return (a >= 0) ? 1.0 : -1.0; } + Real sign(Real a) { return (a >= 0) ? 1.0 : -1.0; } void setProperties( const PropertySet &opts ); PropertySet getProperties() const; diff --git a/include/dai/prob.h b/include/dai/prob.h index 0caa188..0018bcf 100644 --- a/include/dai/prob.h +++ b/include/dai/prob.h @@ -301,10 +301,10 @@ template class TProb { inv._p.reserve( size() ); if( zero ) for( size_t i = 0; i < size(); i++ ) - inv._p.push_back( _p[i] == 0.0 ? 0.0 : 1.0 / _p[i] ); + inv._p.push_back( _p[i] == (T)0 ? (T)0 : (T)1 / _p[i] ); else for( size_t i = 0; i < size(); i++ ) - inv._p.push_back( 1.0 / _p[i] ); + inv._p.push_back( (T)1 / _p[i] ); return inv; } diff --git a/include/dai/regiongraph.h b/include/dai/regiongraph.h index 231e8a4..07ffbaf 100644 --- a/include/dai/regiongraph.h +++ b/include/dai/regiongraph.h @@ -30,20 +30,20 @@ namespace dai { class Region : public VarSet { private: /// Counting number - double _c; + Real _c; public: /// Default constructor Region() : VarSet(), _c(1.0) {} /// Construct from a set of variables and a counting number - Region(const VarSet &x, double c) : VarSet(x), _c(c) {} + Region(const VarSet &x, Real c) : VarSet(x), _c(c) {} /// Returns constant reference to counting number - const double & c() const { return _c; } + const Real & c() const { return _c; } /// Returns reference to counting number - double & c() { return _c; } + Real & c() { return _c; } }; @@ -51,20 +51,20 @@ class Region : public VarSet { class FRegion : public Factor { private: /// Counting number - double _c; + Real _c; public: /// Default constructor FRegion() : Factor(), _c(1.0) {} /// Constructs from a factor and a counting number - FRegion( const Factor & x, double c ) : Factor(x), _c(c) {} + FRegion( const Factor & x, Real c ) : Factor(x), _c(c) {} /// Returns constant reference to counting number - const double & c() const { return _c; } + const Real & c() const { return _c; } /// Returns reference to counting number - double & c() { return _c; } + Real & c() { return _c; } }; diff --git a/include/dai/treeep.h b/include/dai/treeep.h index b479323..86e2b55 100644 --- a/include/dai/treeep.h +++ b/include/dai/treeep.h @@ -38,9 +38,9 @@ namespace dai { class TreeEP : public JTree { private: /// Maximum difference encountered so far - double _maxdiff; + Real _maxdiff; /// Number of iterations needed - size_t _iters; + size_t _iters; public: /// Parameters of this inference algorithm @@ -55,7 +55,7 @@ class TreeEP : public JTree { size_t maxiter; /// Tolerance - double tol; + Real tol; /// How to choose the tree TypeType type; @@ -76,7 +76,7 @@ class TreeEP : public JTree { const Factor * _I; VarSet _ns; VarSet _nsrem; - double _logZ; + Real _logZ; public: @@ -101,7 +101,7 @@ class TreeEP : public JTree { void init(); void InvertAndMultiply( const std::vector &Qa, const std::vector &Qb ); void HUGIN_with_I( std::vector &Qa, std::vector &Qb ); - double logZ( const std::vector &Qa, const std::vector &Qb ) const; + Real logZ( const std::vector &Qa, const std::vector &Qb ) const; const Factor *& I() { return _I; } }; @@ -144,8 +144,8 @@ class TreeEP : public JTree { virtual Real logZ() const; virtual void init(); virtual void init( const VarSet &/*ns*/ ) { init(); } - virtual double run(); - virtual double maxDiff() const { return _maxdiff; } + virtual Real run(); + virtual Real maxDiff() const { return _maxdiff; } virtual size_t Iterations() const { return _iters; } //@} diff --git a/include/dai/util.h b/include/dai/util.h index e2b4f13..62db510 100644 --- a/include/dai/util.h +++ b/include/dai/util.h @@ -79,7 +79,7 @@ double log1p( double x ); /// Define INFINITY - #define INFINITY (std::numeric_limits::infinity()) + #define INFINITY (std::numeric_limits::infinity()) #endif @@ -197,28 +197,28 @@ std::vector concat( const std::vector& u, const std::vector& v ) { void tokenizeString( const std::string& s, std::vector& outTokens, const std::string& delim="\t\n" ); /// Used to keep track of the progress made by iterative algorithms -class Diffs : public std::vector { +class Diffs : public std::vector { private: size_t _maxsize; - double _def; - std::vector::iterator _pos; - std::vector::iterator _maxpos; + Real _def; + std::vector::iterator _pos; + std::vector::iterator _maxpos; public: /// Constructor - Diffs(long maxsize, double def) : std::vector(), _maxsize(maxsize), _def(def) { + Diffs(long maxsize, Real def) : std::vector(), _maxsize(maxsize), _def(def) { this->reserve(_maxsize); _pos = begin(); _maxpos = begin(); } /// Returns maximum difference encountered - double maxDiff() { + Real maxDiff() { if( size() < _maxsize ) return _def; else return( *_maxpos ); } /// Register new difference x - void push(double x) { + void push(Real x) { if( size() < _maxsize ) { push_back(x); _pos = end(); diff --git a/src/bbp.cpp b/src/bbp.cpp index 96f8b19..973293e 100644 --- a/src/bbp.cpp +++ b/src/bbp.cpp @@ -678,7 +678,7 @@ std::vector BBP::getZeroAdjV( const FactorGraph &fg ) { void BBP::run() { typedef BBP::Properties::UpdateType UT; - Real &tol = props.tol; + Real tol = props.tol; UT &updates = props.updates; Real tic = toc(); @@ -760,7 +760,7 @@ void BBP::run() { } -double numericBBPTest( const InfAlg &bp, const vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) { +Real numericBBPTest( const InfAlg &bp, const vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h ) { // calculate the value of the unperturbed cost function Real cf0 = getCostFn( bp, cfn, state ); @@ -777,7 +777,7 @@ double numericBBPTest( const InfAlg &bp, const vector *state, const Prop // for each variable i for( size_t i = 0; i < fg.nrVars(); i++ ) { - vector adj_est; + vector adj_est; // for each value xi for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) { // Clone 'bp' (which may be any InfAlg) @@ -828,7 +828,7 @@ double numericBBPTest( const InfAlg &bp, const vector *state, const Prop for(size_t i=0; i adj_n_est; + vector adj_n_est; // for each value xi for(size_t xi=0; xi *state, const Prop adj_n_est.push_back((cf_prb-cf0)/h); } - vector adj_m_est; + vector adj_m_est; // for each value xi for(size_t xi=0; xi *state, const Prop /* if(0) { // verify bbp.adj_b_V for(size_t i=0; i adj_b_V_est; + vector adj_b_V_est; // for each value xi for(size_t xi=0; xi *stateP ) { - double cf = 0.0; + Real cf = 0.0; const FactorGraph &fg = ia.fg(); switch( (size_t)cfn_type ) { @@ -1094,7 +1094,7 @@ Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector *stat vector state = *stateP; DAI_ASSERT( state.size() == fg.nrVars() ); for( size_t i = 0; i < fg.nrVars(); i++ ) { - double b = ia.beliefV(i)[state[i]]; + Real b = ia.beliefV(i)[state[i]]; switch( (size_t)cfn_type ) { case bbp_cfn_t::CFN_GIBBS_B: cf += b; @@ -1118,7 +1118,7 @@ Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector *stat DAI_ASSERT( state.size() == fg.nrVars() ); for( size_t I = 0; I < fg.nrFactors(); I++ ) { size_t x_I = getFactorEntryForState( fg, I, state ); - double b = ia.beliefF(I)[x_I]; + Real b = ia.beliefF(I)[x_I]; switch( (size_t)cfn_type ) { case bbp_cfn_t::CFN_GIBBS_B_FACTOR: cf += b; @@ -1173,8 +1173,8 @@ void BBP::Properties::set(const PropertySet &opts) DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED,"BBP: Missing property \"updates\" for method \"BBP\""); verbose = opts.getStringAs("verbose"); maxiter = opts.getStringAs("maxiter"); - tol = opts.getStringAs("tol"); - damping = opts.getStringAs("damping"); + tol = opts.getStringAs("tol"); + damping = opts.getStringAs("damping"); updates = opts.getStringAs("updates"); } PropertySet BBP::Properties::get() const { diff --git a/src/bp.cpp b/src/bp.cpp index e3043cc..39c70ce 100644 --- a/src/bp.cpp +++ b/src/bp.cpp @@ -38,7 +38,7 @@ void BP::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("logdomain") ); DAI_ASSERT( opts.hasKey("updates") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.maxiter = opts.getStringAs("maxiter"); props.logdomain = opts.getStringAs("logdomain"); props.updates = opts.getStringAs("updates"); @@ -48,7 +48,7 @@ void BP::setProperties( const PropertySet &opts ) { else props.verbose = 0; if( opts.hasKey("damping") ) - props.damping = opts.getStringAs("damping"); + props.damping = opts.getStringAs("damping"); else props.damping = 0.0; if( opts.hasKey("inference") ) @@ -120,7 +120,7 @@ void BP::construct() { void BP::init() { - double c = props.logdomain ? 0.0 : 1.0; + Real c = props.logdomain ? 0.0 : 1.0; for( size_t i = 0; i < nrVars(); ++i ) { foreach( const Neighbor &I, nbV(i) ) { message( i, I.iter ).fill( c ); @@ -225,7 +225,7 @@ void BP::calcNewMessage( size_t i, size_t _I ) { // BP::run does not check for NANs for performance reasons // Somehow NaNs do not often occur in BP... -double BP::run() { +Real BP::run() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 3) @@ -460,7 +460,7 @@ void BP::init( const VarSet &ns ) { for( VarSet::const_iterator n = ns.begin(); n != ns.end(); ++n ) { size_t ni = findVar( *n ); foreach( const Neighbor &I, nbV( ni ) ) { - double val = props.logdomain ? 0.0 : 1.0; + Real val = props.logdomain ? 0.0 : 1.0; message( ni, I.iter ).fill( val ); newMessage( ni, I.iter ).fill( val ); if( props.updates == Properties::UpdateType::SEQMAX ) @@ -485,7 +485,7 @@ void BP::updateMessage( size_t i, size_t _I ) { } -void BP::updateResidual( size_t i, size_t _I, double r ) { +void BP::updateResidual( size_t i, size_t _I, Real r ) { EdgeProp* pEdge = &_edges[i][_I]; pEdge->residual = r; @@ -508,7 +508,7 @@ std::vector BP::findMaximum() const { // Maximise with respect to variable i Prob prod; calcBeliefV( i, prod ); - maximum[i] = max_element( prod.begin(), prod.end() ) - prod.begin(); + maximum[i] = prod.argmax().first; foreach( const Neighbor &I, nbV(i) ) if( !visitedFactors[I] ) diff --git a/src/cbp.cpp b/src/cbp.cpp index 28c10ee..d3e960f 100644 --- a/src/cbp.cpp +++ b/src/cbp.cpp @@ -32,7 +32,7 @@ using boost::shared_ptr; const char *CBP::Name = "CBP"; -void CBP::setBeliefs( const std::vector &bs, double logZ ) { +void CBP::setBeliefs( const std::vector &bs, Real logZ ) { size_t i = 0; _beliefsV.clear(); _beliefsV.reserve( nrVars() ); @@ -87,7 +87,7 @@ static vector mixBeliefs( Real p, const vector &b, const vector< } -double CBP::run() { +Real CBP::run() { size_t seed = props.rand_seed; if( seed > 0) rnd_seed( seed ); @@ -116,7 +116,7 @@ InfAlg* CBP::getInfAlg() { bpProps.Set("maxiter", props.maxiter); bpProps.Set("verbose", props.verbose); bpProps.Set("logdomain", false); - bpProps.Set("damping", 0.0); + bpProps.Set("damping", (Real)0.0); BP *bp = new BP( *this, bpProps ); bp->recordSentMessages = true; bp->init(); @@ -124,8 +124,8 @@ InfAlg* CBP::getInfAlg() { } -void CBP::runRecurse( InfAlg *bp, double orig_logZ, vector clamped_vars_list, size_t &num_leaves, - size_t &choose_count, double &sum_level, Real &lz_out, vector& beliefs_out) { +void CBP::runRecurse( InfAlg *bp, Real orig_logZ, vector clamped_vars_list, size_t &num_leaves, + size_t &choose_count, Real &sum_level, Real &lz_out, vector& beliefs_out) { // choose a variable/states to clamp: size_t i; vector xis; @@ -195,7 +195,7 @@ void CBP::runRecurse( InfAlg *bp, double orig_logZ, vector clamped_vars_ cmp_lz = cmp_bp_c->logZ(); cmp_b = cmp_bp_c->beliefs(); - double p = unSoftMax( lz, cmp_lz ); + Real p = unSoftMax( lz, cmp_lz ); Real bp__d = 0.0; if( Recursion() == Properties::RecurseType::REC_BDIFF && recTol() > 0 ) { @@ -567,16 +567,16 @@ void CBP::Properties::set(const PropertySet &opts) } else { verbose = 0; } - tol = opts.getStringAs("tol"); + tol = opts.getStringAs("tol"); updates = opts.getStringAs("updates"); maxiter = opts.getStringAs("maxiter"); - rec_tol = opts.getStringAs("rec_tol"); + rec_tol = opts.getStringAs("rec_tol"); if(opts.hasKey("max_levels")) { max_levels = opts.getStringAs("max_levels"); } else { max_levels = 10; } - min_max_adj = opts.getStringAs("min_max_adj"); + min_max_adj = opts.getStringAs("min_max_adj"); choose = opts.getStringAs("choose"); recursion = opts.getStringAs("recursion"); clamp = opts.getStringAs("clamp"); diff --git a/src/daialg.cpp b/src/daialg.cpp index 241dc43..b1cd4f8 100644 --- a/src/daialg.cpp +++ b/src/daialg.cpp @@ -130,7 +130,7 @@ vector calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reIni if( logZ != -INFINITY ) logZ0 = logZ; - double Z_xj; + Real Z_xj; if( logZ == -INFINITY ) Z_xj = 0; else @@ -230,7 +230,7 @@ vector calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool re if( logZ != -INFINITY ) logZ0 = logZ; - double Z_xj; + Real Z_xj; if( logZ == -INFINITY ) Z_xj = 0; else diff --git a/src/exactinf.cpp b/src/exactinf.cpp index 12b1253..312b797 100644 --- a/src/exactinf.cpp +++ b/src/exactinf.cpp @@ -67,7 +67,7 @@ void ExactInf::init() { } -double ExactInf::run() { +Real ExactInf::run() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; diff --git a/src/factorgraph.cpp b/src/factorgraph.cpp index f94586d..b737864 100644 --- a/src/factorgraph.cpp +++ b/src/factorgraph.cpp @@ -87,11 +87,11 @@ std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) { os << endl; size_t nr_nonzeros = 0; for( size_t k = 0; k < fg.factor(I).states(); k++ ) - if( fg.factor(I)[k] != 0.0 ) + if( fg.factor(I)[k] != (Real)0 ) nr_nonzeros++; os << nr_nonzeros << endl; for( size_t k = 0; k < fg.factor(I).states(); k++ ) - if( fg.factor(I)[k] != 0.0 ) + if( fg.factor(I)[k] != (Real)0 ) os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl; } @@ -164,7 +164,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) { vardims[labels[mi]] = dims[mi]; I_vars |= Var(labels[mi], dims[mi]); } - facs.push_back( Factor( I_vars, 0.0 ) ); + facs.push_back( Factor( I_vars, (Real)0 ) ); // calculate permutation sigma (internally, members are sorted) vector sigma(nr_members,0); @@ -189,7 +189,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) { cerr << " nonzeroes: " << nr_nonzeros << endl; for( size_t k = 0; k < nr_nonzeros; k++ ) { size_t li; - double val; + Real val; while( (is.peek()) == '#' ) getline(is,line); is >> li; @@ -240,7 +240,7 @@ void FactorGraph::makeCavity( size_t i, bool backup ) { // fills all Factors that include var(i) with ones map newFacs; foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i - newFacs[I] = Factor(factor(I).vars(), 1.0); + newFacs[I] = Factor( factor(I).vars(), (Real)1 ); setFactors( newFacs, backup ); } @@ -302,8 +302,8 @@ vector FactorGraph::Cliques() const { void FactorGraph::clamp( size_t i, size_t x, bool backup ) { DAI_ASSERT( x <= var(i).states() ); - Factor mask( var(i), 0.0 ); - mask[x] = 1.0; + Factor mask( var(i), (Real)0 ); + mask[x] = (Real)1; map newFacs; foreach( const Neighbor &I, nbV(i) ) @@ -316,11 +316,11 @@ void FactorGraph::clamp( size_t i, size_t x, bool backup ) { void FactorGraph::clampVar( size_t i, const vector &is, bool backup ) { Var n = var(i); - Factor mask_n( n, 0.0 ); + Factor mask_n( n, (Real)0 ); foreach( size_t i, is ) { DAI_ASSERT( i <= n.states() ); - mask_n[i] = 1.0; + mask_n[i] = (Real)1; } map newFacs; @@ -332,7 +332,7 @@ void FactorGraph::clampVar( size_t i, const vector &is, bool backup ) { void FactorGraph::clampFactor( size_t I, const vector &is, bool backup ) { size_t st = factor(I).states(); - Factor newF( factor(I).vars(), 0.0 ); + Factor newF( factor(I).vars(), (Real)0 ); foreach( size_t i, is ) { DAI_ASSERT( i <= st ); @@ -412,7 +412,7 @@ bool FactorGraph::isBinary() const { FactorGraph FactorGraph::clamped( size_t i, size_t state ) const { Var v = var( i ); - Real zeroth_order = 1.0; + Real zeroth_order = (Real)1; vector clamped_facs; for( size_t I = 0; I < nrFactors(); I++ ) { VarSet v_I = factor(I).vars(); diff --git a/src/gibbs.cpp b/src/gibbs.cpp index c1bff37..cd29490 100644 --- a/src/gibbs.cpp +++ b/src/gibbs.cpp @@ -157,7 +157,7 @@ void Gibbs::init() { } -double Gibbs::run() { +Real Gibbs::run() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 3 ) @@ -188,12 +188,12 @@ double Gibbs::run() { Factor Gibbs::beliefV( size_t i ) const { - return Factor( var(i), _var_counts[i].begin() ).normalized(); + return Factor( var(i), _var_counts[i] ).normalized(); } Factor Gibbs::beliefF( size_t I ) const { - return Factor( factor(I).vars(), _factor_counts[I].begin() ).normalized(); + return Factor( factor(I).vars(), _factor_counts[I] ).normalized(); } diff --git a/src/hak.cpp b/src/hak.cpp index 87fc45b..dab801a 100644 --- a/src/hak.cpp +++ b/src/hak.cpp @@ -51,7 +51,7 @@ void HAK::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("doubleloop") ); DAI_ASSERT( opts.hasKey("clusters") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.maxiter = opts.getStringAs("maxiter"); props.verbose = opts.getStringAs("verbose"); props.doubleloop = opts.getStringAs("doubleloop"); @@ -62,7 +62,7 @@ void HAK::setProperties( const PropertySet &opts ) { else DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP ); if( opts.hasKey("damping") ) - props.damping = opts.getStringAs("damping"); + props.damping = opts.getStringAs("damping"); else props.damping = 0.0; if( opts.hasKey("init") ) @@ -248,7 +248,7 @@ void HAK::init() { } -double HAK::doGBP() { +Real HAK::doGBP() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 3) @@ -373,7 +373,7 @@ double HAK::doGBP() { } -double HAK::doDoubleLoop() { +Real HAK::doDoubleLoop() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 3) @@ -385,7 +385,7 @@ double HAK::doDoubleLoop() { vector org_ORs = ORs; // Save original inner counting numbers and set negative counting numbers to zero - vector org_IR_cs( nrIRs(), 0.0 ); + vector org_IR_cs( nrIRs(), 0.0 ); for( size_t beta = 0; beta < nrIRs(); beta++ ) { org_IR_cs[beta] = IR(beta).c(); if( IR(beta).c() < 0.0 ) @@ -402,9 +402,9 @@ double HAK::doDoubleLoop() { Diffs diffs(nrVars(), 1.0); size_t outer_maxiter = props.maxiter; - double outer_tol = props.tol; + Real outer_tol = props.tol; size_t outer_verbose = props.verbose; - double org_maxdiff = _maxdiff; + Real org_maxdiff = _maxdiff; // Set parameters for inner loop props.maxiter = 5; @@ -469,7 +469,7 @@ double HAK::doDoubleLoop() { } -double HAK::run() { +Real HAK::run() { if( props.doubleloop ) return doDoubleLoop(); else diff --git a/src/jtree.cpp b/src/jtree.cpp index 976562f..61d24e1 100644 --- a/src/jtree.cpp +++ b/src/jtree.cpp @@ -314,7 +314,7 @@ void JTree::runShaferShenoy() { } -double JTree::run() { +Real JTree::run() { if( props.updates == Properties::UpdateType::HUGIN ) runHUGIN(); else if( props.updates == Properties::UpdateType::SHSH ) @@ -484,7 +484,7 @@ Factor JTree::calcMarginal( const VarSet& ns ) { // For all states of nsrem for( State s(nsrem); s.valid(); s++ ) { // CollectEvidence - double logZ = 0.0; + Real logZ = 0.0; for( size_t i = Tsize; (i--) != 0; ) { // Make outer region T[i].n1 consistent with outer region T[i].n2 // IR(i) = seperator OR(T[i].n1) && OR(T[i].n2) diff --git a/src/lc.cpp b/src/lc.cpp index d6ba935..5c2e1e3 100644 --- a/src/lc.cpp +++ b/src/lc.cpp @@ -34,7 +34,7 @@ void LC::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("cavity") ); DAI_ASSERT( opts.hasKey("updates") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.maxiter = opts.getStringAs("maxiter"); props.verbose = opts.getStringAs("verbose"); props.cavity = opts.getStringAs("cavity"); @@ -46,7 +46,7 @@ void LC::setProperties( const PropertySet &opts ) { if( opts.hasKey("reinit") ) props.reinit = opts.getStringAs("reinit"); if( opts.hasKey("damping") ) - props.damping = opts.getStringAs("damping"); + props.damping = opts.getStringAs("damping"); else props.damping = 0.0; } @@ -119,9 +119,9 @@ void LC::CalcBelief (size_t i) { } -double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) { +Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) { Factor Bi; - double maxdiff = 0; + Real maxdiff = 0; if( props.verbose >= 2 ) cerr << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl; @@ -151,7 +151,7 @@ double LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet } -double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) { +Real LC::InitCavityDists( const std::string &name, const PropertySet &opts ) { double tic = toc(); if( props.verbose >= 1 ) { @@ -166,9 +166,9 @@ double LC::InitCavityDists( const std::string &name, const PropertySet &opts ) { cerr << "Using pairwise(new) " << name << opts << "..."; } - double maxdiff = 0.0; + Real maxdiff = 0.0; for( size_t i = 0; i < nrVars(); i++ ) { - double md = CalcCavityDist(i, name, opts); + Real md = CalcCavityDist(i, name, opts); if( md > maxdiff ) maxdiff = md; } @@ -237,7 +237,7 @@ Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) { } -double LC::run() { +Real LC::run() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 2 ) @@ -246,7 +246,7 @@ double LC::run() { double tic = toc(); Diffs diffs(nrVars(), 1.0); - double md = InitCavityDists( props.cavainame, props.cavaiopts ); + Real md = InitCavityDists( props.cavainame, props.cavaiopts ); if( md > _maxdiff ) _maxdiff = md; diff --git a/src/mf.cpp b/src/mf.cpp index 23ab0e2..3e96c82 100644 --- a/src/mf.cpp +++ b/src/mf.cpp @@ -30,14 +30,14 @@ void MF::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("tol") ); DAI_ASSERT( opts.hasKey("maxiter") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.maxiter = opts.getStringAs("maxiter"); if( opts.hasKey("verbose") ) props.verbose = opts.getStringAs("verbose"); else props.verbose = 0U; if( opts.hasKey("damping") ) - props.damping = opts.getStringAs("damping"); + props.damping = opts.getStringAs("damping"); else props.damping = 0.0; } @@ -84,7 +84,7 @@ void MF::init() { } -double MF::run() { +Real MF::run() { double tic = toc(); if( props.verbose >= 1 ) diff --git a/src/mr.cpp b/src/mr.cpp index 0045a7a..5cb1b55 100644 --- a/src/mr.cpp +++ b/src/mr.cpp @@ -35,7 +35,7 @@ void MR::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("updates") ); DAI_ASSERT( opts.hasKey("inits") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.verbose = opts.getStringAs("verbose"); props.updates = opts.getStringAs("updates"); props.inits = opts.getStringAs("inits"); @@ -64,7 +64,7 @@ string MR::printProperties() const { // init N, con, nb, tJ, theta -void MR::init(size_t Nin, double *_w, double *_th) { +void MR::init(size_t Nin, Real *_w, Real *_th) { size_t i,j; N = Nin; @@ -91,33 +91,33 @@ void MR::init(size_t Nin, double *_w, double *_th) { // calculate cors -double MR::init_cor_resp() { +Real MR::init_cor_resp() { size_t j,k,l, runx,i2; - double variab1, variab2; - double md, maxdev; - double thbJsite[kmax]; - double xinter; - double rinter; - double res[kmax]; + Real variab1, variab2; + Real md, maxdev; + Real thbJsite[kmax]; + Real xinter; + Real rinter; + Real res[kmax]; size_t s2; size_t flag; size_t concav; size_t runs = 3000; - double eps = 0.2; + Real eps = 0.2; size_t cavity; - vector > tJ_org; + vector > tJ_org; vector > nb_org; vector con_org; - vector theta_org; + vector theta_org; - vector xfield(N*kmax,0.0); - vector rfield(N*kmax,0.0); - vector Hfield(N,0.0); - vector devs(N*kmax,0.0); - vector devs2(N*kmax,0.0); - vector dev(N,0.0); - vector avmag(N,0.0); + vector xfield(N*kmax,0.0); + vector rfield(N*kmax,0.0); + vector Hfield(N,0.0); + vector devs(N*kmax,0.0); + vector devs2(N*kmax,0.0); + vector dev(N,0.0); + vector avmag(N,0.0); // save original tJ, nb nb_org = nb; @@ -248,7 +248,7 @@ double MR::init_cor_resp() { } -double MR::T(size_t i, sub_nb A) { +Real MR::T(size_t i, sub_nb A) { // i is a variable index // A is a subset of nb[i] // @@ -258,7 +258,7 @@ double MR::T(size_t i, sub_nb A) { _nbi_min_A.set(); _nbi_min_A &= ~A; - double res = theta[i]; + Real res = theta[i]; for( size_t _j = 0; _j < _nbi_min_A.size(); _j++ ) if( _nbi_min_A.test(_j) ) res += atanh(tJ[i][_j] * M[i][_j]); @@ -266,46 +266,46 @@ double MR::T(size_t i, sub_nb A) { } -double MR::T(size_t i, size_t _j) { +Real MR::T(size_t i, size_t _j) { sub_nb j(con[i]); j.set(_j); return T(i,j); } -double MR::Omega(size_t i, size_t _j, size_t _l) { +Real MR::Omega(size_t i, size_t _j, size_t _l) { sub_nb jl(con[i]); jl.set(_j); jl.set(_l); - double Tijl = T(i,jl); + Real Tijl = T(i,jl); return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl); } -double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) { +Real MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) { sub_nb jll(con[i]); jll.set(_j); - double Tij = T(i,jll); + Real Tij = T(i,jll); jll.set(_l1); jll.set(_l2); - double Tijll = T(i,jll); + Real Tijll = T(i,jll); return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll); } -double MR::Gamma(size_t i, size_t _l1, size_t _l2) { +Real MR::Gamma(size_t i, size_t _l1, size_t _l2) { sub_nb ll(con[i]); - double Ti = T(i,ll); + Real Ti = T(i,ll); ll.set(_l1); ll.set(_l2); - double Till = T(i,ll); + Real Till = T(i,ll); return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till); } -double MR::_tJ(size_t i, sub_nb A) { +Real MR::_tJ(size_t i, sub_nb A) { // i is a variable index // A is a subset of nb[i] // @@ -319,7 +319,7 @@ double MR::_tJ(size_t i, sub_nb A) { } -double MR::appM(size_t i, sub_nb A) { +Real MR::appM(size_t i, sub_nb A) { // i is a variable index // A is a subset of nb[i] // @@ -334,7 +334,7 @@ double MR::appM(size_t i, sub_nb A) { else { sub_nb A_j(A); A_j.reset(_j); - double result = M[i][_j] * appM(i, A_j); + Real result = M[i][_j] * appM(i, A_j); for( size_t _k = 0; _k < A_j.size(); _k++ ) if( A_j.test(_k) ) { sub_nb A_jk(A_j); A_jk.reset(_k); @@ -346,7 +346,7 @@ double MR::appM(size_t i, sub_nb A) { } -void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) { +void MR::sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd) { // j is a variable index // A is a subset of nb[j] @@ -378,8 +378,8 @@ void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) { void MR::solvemcav() { - double sum_even, sum_odd; - double maxdev; + Real sum_even, sum_odd; + Real maxdev; size_t maxruns = 1000; makekindex(); @@ -397,7 +397,7 @@ void MR::solvemcav() { size_t j = nb[i][_j]; DAI_ASSERT( nb[j][_i] == i ); - double newM = 0.0; + Real newM = 0.0; if( props.updates == Properties::UpdateType::FULL ) { // find indices in nb[j] that do not correspond with i sub_nb _nbj_min_i(con[j]); @@ -413,8 +413,8 @@ void MR::solvemcav() { newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd); sum_subs(i, _nbi_min_j, &sum_even, &sum_odd); - double denom = sum_even + tanh(theta[i]) * sum_odd; - double numer = 0.0; + Real denom = sum_even + tanh(theta[i]) * sum_odd; + Real numer = 0.0; for(size_t _k=0; _k= maxdev ) maxdev = fabs(dev); @@ -463,7 +463,7 @@ void MR::solveM() { _nbi.set(); // calc numerator1 and denominator1 - double sum_even, sum_odd; + Real sum_even, sum_odd; sum_subs(i, _nbi, &sum_even, &sum_odd); Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd); @@ -486,7 +486,7 @@ void MR::init_cor() { for( size_t i = 0; i < nrVars(); i++ ) { vector pairq; if( props.inits == Properties::InitType::CLAMPING ) { - BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", 1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false)); + BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false)); bpcav.makeCavity( i ); pairq = calcPairBeliefs( bpcav, delta(i), false ); } else if( props.inits == Properties::InitType::EXACT ) { @@ -499,7 +499,7 @@ void MR::init_cor() { size_t j = findVar( *(kit) ); size_t k = findVar( *(++kit) ); pairq[jk].normalize(); - double cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]); + Real cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]); for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j ) for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) { cors[i][_j][_k] = cor; @@ -515,7 +515,7 @@ string MR::identify() const { } -double MR::run() { +Real MR::run() { if( supported ) { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; @@ -539,7 +539,7 @@ double MR::run() { kindex[i].resize(kmax); if( props.inits == Properties::InitType::RESPPROP ) { - double md = init_cor_resp(); + Real md = init_cor_resp(); if( md > _maxdiff ) _maxdiff = md; } else if( props.inits == Properties::InitType::EXACT ) @@ -577,7 +577,7 @@ Factor MR::belief( const Var &n ) const { if( supported ) { size_t i = findVar( n ); - Real x[2]; + Prob x(2); x[0] = 0.5 - Mag[i] / 2.0; x[1] = 0.5 + Mag[i] / 2.0; @@ -623,8 +623,8 @@ MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), support // create w and th size_t Nin = fg.nrVars(); - double *w = new double[Nin*Nin]; - double *th = new double[Nin]; + Real *w = new Real[Nin*Nin]; + Real *th = new Real[Nin]; for( size_t i = 0; i < Nin; i++ ) { th[i] = 0.0; diff --git a/src/properties.cpp b/src/properties.cpp index d7a87eb..628ef98 100644 --- a/src/properties.cpp +++ b/src/properties.cpp @@ -25,6 +25,8 @@ std::ostream& operator<< (std::ostream & os, const Property & p) { os << boost::any_cast(p.second); else if( p.second.type() == typeid(double) ) os << boost::any_cast(p.second); + else if( p.second.type() == typeid(long double) ) + os << boost::any_cast(p.second); else if( p.second.type() == typeid(bool) ) os << boost::any_cast(p.second); else if( p.second.type() == typeid(PropertySet) ) diff --git a/src/regiongraph.cpp b/src/regiongraph.cpp index 1c20e1f..5b1ef7e 100644 --- a/src/regiongraph.cpp +++ b/src/regiongraph.cpp @@ -149,7 +149,7 @@ void RegionGraph::calcCountingNumbers() { if( !assigned[*beta2] ) has_unassigned_ancestor = true; if( !has_unassigned_ancestor ) { - double c = 1.0; + Real c = 1.0; foreach( const Neighbor &alpha, nbIR(beta) ) c -= OR(alpha).c(); for( vector::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ ) @@ -169,7 +169,7 @@ bool RegionGraph::checkCountingNumbers() const { bool all_valid = true; for( vector::const_iterator n = vars().begin(); n != vars().end(); n++ ) { - double c_n = 0.0; + Real c_n = 0.0; for( size_t alpha = 0; alpha < nrORs(); alpha++ ) if( OR(alpha).vars().contains( *n ) ) c_n += OR(alpha).c(); diff --git a/src/treeep.cpp b/src/treeep.cpp index 90883ac..a914223 100644 --- a/src/treeep.cpp +++ b/src/treeep.cpp @@ -32,7 +32,7 @@ void TreeEP::setProperties( const PropertySet &opts ) { DAI_ASSERT( opts.hasKey("verbose") ); DAI_ASSERT( opts.hasKey("type") ); - props.tol = opts.getStringAs("tol"); + props.tol = opts.getStringAs("tol"); props.maxiter = opts.getStringAs("maxiter"); props.verbose = opts.getStringAs("verbose"); props.type = opts.getStringAs("type"); @@ -176,8 +176,8 @@ void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector &Qa, std::vector &Qa, const std::vector &Qb ) const { - double s = 0.0; +Real TreeEP::TreeEPSubTree::logZ( const std::vector &Qa, const std::vector &Qb ) const { + Real s = 0.0; for( size_t alpha = 0; alpha < _Qa.size(); alpha++ ) s += (Qa[_a[alpha]] * _Qa[alpha].log(true)).sum(); for( size_t beta = 0; beta < _Qb.size(); beta++ ) @@ -200,7 +200,7 @@ TreeEP::TreeEP( const FactorGraph &fg, const PropertySet &opts ) : JTree(fg, opt // ALT: construct weighted graph with as weights an upper bound on the // effective interaction strength between pairs of nodes - WeightedGraph wg; + WeightedGraph wg; for( size_t i = 0; i < nrVars(); ++i ) { Var v_i = var(i); VarSet di = delta(i); @@ -365,7 +365,7 @@ void TreeEP::init() { } -double TreeEP::run() { +Real TreeEP::run() { if( props.verbose >= 1 ) cerr << "Starting " << identify() << "..."; if( props.verbose >= 3) @@ -420,7 +420,7 @@ double TreeEP::run() { Real TreeEP::logZ() const { - double s = 0.0; + Real s = 0.0; // entropy of the tree for( size_t beta = 0; beta < nrIRs(); beta++ ) diff --git a/src/util.cpp b/src/util.cpp index 73fdf8f..d2a67b7 100644 --- a/src/util.cpp +++ b/src/util.cpp @@ -75,7 +75,7 @@ boost::uniform_real _uni_dist(0,1); boost::variate_generator<_rnd_gen_type&, boost::uniform_real > _uni_rnd(_rnd_gen, _uni_dist); /// Normal distribution with mean 0 and standard deviation 1. -boost::normal_distribution<> _normal_dist; +boost::normal_distribution _normal_dist; /// Global random number generator with standard normal distribution boost::variate_generator<_rnd_gen_type&, boost::normal_distribution > _normal_rnd(_rnd_gen, _normal_dist); diff --git a/tests/testbbp.cpp b/tests/testbbp.cpp index 9a4a8b9..45cfbdd 100644 --- a/tests/testbbp.cpp +++ b/tests/testbbp.cpp @@ -30,9 +30,9 @@ int main( int argc, char *argv[] ) { // Set some constants size_t verbose = 0; - double tol = 1.0e-9; + Real tol = 1.0e-9; size_t maxiter = 10000; - double damping = 0.0; + Real damping = 0.0; BBP::Properties::UpdateType updates = BBP::Properties::UpdateType::PAR; // Store the constants in a PropertySet object @@ -100,8 +100,8 @@ int main( int argc, char *argv[] ) { break; } - double h = 1e-4; - double result = numericBBPTest( bp, &state, opts("updates",updates), cfn, h ); + Real h = 1e-4; + Real result = numericBBPTest( bp, &state, opts("updates",updates), cfn, h ); cout << "result: " << result << ",\tupdates=" << updates << ", cfn=" << cfn << endl; } } diff --git a/tests/testdai.cpp b/tests/testdai.cpp index 6704b66..9b61d57 100644 --- a/tests/testdai.cpp +++ b/tests/testdai.cpp @@ -30,12 +30,12 @@ class TestDAI { protected: InfAlg *obj; string name; - vector err; + vector err; public: vector q; - double logZ; - double maxdiff; + Real logZ; + Real maxdiff; double time; size_t iters; bool has_logZ; @@ -45,7 +45,8 @@ class TestDAI { TestDAI( const FactorGraph &fg, const string &_name, const PropertySet &opts ) : obj(NULL), name(_name), err(), q(), logZ(0.0), maxdiff(0.0), time(0), iters(0U), has_logZ(false), has_maxdiff(false), has_iters(false) { double tic = toc(); if( name == "LDPC" ) { - double zero[2] = {1.0, 0.0}; + Prob zero(2,0.0); + zero[0] = 1.0; q.clear(); for( size_t i = 0; i < fg.nrVars(); i++ ) q.push_back( Factor(Var(i,2), zero) ); @@ -134,11 +135,11 @@ class TestDAI { err.push_back( dist( q[i], x[i], Prob::DISTTV ) ); } - double maxErr() { + Real maxErr() { return( *max_element( err.begin(), err.end() ) ); } - double avgErr() { + Real avgErr() { return( accumulate( err.begin(), err.end(), 0.0 ) / err.size() ); } }; @@ -192,8 +193,8 @@ pair parseMethod( const string &_s, const map methods; - double tol; + Real tol; size_t maxiter; size_t verbose; bool marginals = false; @@ -221,7 +222,7 @@ int main( int argc, char *argv[] ) { opts_optional.add_options() ("help", "produce help message") ("aliases", po::value< string >(&aliases), "Filename for aliases") - ("tol", po::value< double >(&tol), "Override tolerance") + ("tol", po::value< Real >(&tol), "Override tolerance") ("maxiter", po::value< size_t >(&maxiter), "Override maximum number of iterations") ("verbose", po::value< size_t >(&verbose), "Override verbosity") ("marginals", po::value< bool >(&marginals), "Output single node marginals?") @@ -282,7 +283,7 @@ int main( int argc, char *argv[] ) { fg.ReadFromFile( filename.c_str() ); vector q0; - double logZ0 = 0.0; + Real logZ0 = 0.0; cout.setf( ios_base::scientific ); cout.precision( 3 ); @@ -333,22 +334,22 @@ int main( int argc, char *argv[] ) { cout.setf( ios_base::scientific ); cout.precision( 3 ); - double me = clipdouble( piet.maxErr(), 1e-9 ); + Real me = clipReal( piet.maxErr(), 1e-9 ); cout << me << "\t"; - double ae = clipdouble( piet.avgErr(), 1e-9 ); + Real ae = clipReal( piet.avgErr(), 1e-9 ); cout << ae << "\t"; if( piet.has_logZ ) { cout.setf( ios::showpos ); - double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 ); + Real le = clipReal( piet.logZ / logZ0 - 1.0, 1e-9 ); cout << le << "\t"; cout.unsetf( ios::showpos ); } else cout << "N/A \t"; if( piet.has_maxdiff ) { - double md = clipdouble( piet.maxdiff, 1e-9 ); + Real md = clipReal( piet.maxdiff, 1e-9 ); if( isnan( me ) ) md = me; if( isnan( ae ) ) diff --git a/utils/createfg.cpp b/utils/createfg.cpp index f495424..5edf476 100644 --- a/utils/createfg.cpp +++ b/utils/createfg.cpp @@ -36,32 +36,32 @@ using namespace std; using namespace dai; namespace po = boost::program_options; -typedef boost::numeric::ublas::compressed_matrix matrix; +typedef boost::numeric::ublas::compressed_matrix matrix; typedef matrix::value_array_type::const_iterator matrix_vcit; typedef matrix::index_array_type::const_iterator matrix_icit; -Factor BinaryFactor( const Var &n, double field ) { +Factor BinaryFactor( const Var &n, Real field ) { DAI_ASSERT( n.states() == 2 ); - double buf[2]; + Real buf[2]; buf[0] = std::exp(-field); buf[1] = std::exp(field); return Factor(n, &buf[0]); } -Factor BinaryFactor( const Var &n1, const Var &n2, double coupling ) { +Factor BinaryFactor( const Var &n1, const Var &n2, Real coupling ) { DAI_ASSERT( n1.states() == 2 ); DAI_ASSERT( n2.states() == 2 ); DAI_ASSERT( n1 != n2 ); - double buf[4]; + Real buf[4]; buf[0] = (buf[3] = std::exp(coupling)); buf[1] = (buf[2] = std::exp(-coupling)); return Factor( VarSet(n1, n2), &buf[0] ); } -Factor RandomFactor( const VarSet &ns, double beta ) { +Factor RandomFactor( const VarSet &ns, Real beta ) { Factor fac( ns ); for( size_t t = 0; t < fac.states(); t++ ) fac[t] = std::exp(rnd_stdnormal() * beta); @@ -69,7 +69,7 @@ Factor RandomFactor( const VarSet &ns, double beta ) { } -Factor PottsFactor( const Var &n1, const Var &n2, double beta ) { +Factor PottsFactor( const Var &n1, const Var &n2, Real beta ) { Factor fac( VarSet( n1, n2 ), 1.0 ); DAI_ASSERT( n1.states() == n2.states() ); for( size_t s = 0; s < n1.states(); s++ ) @@ -78,7 +78,7 @@ Factor PottsFactor( const Var &n1, const Var &n2, double beta ) { } -void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) { +void MakeHOIFG( size_t N, size_t M, size_t k, Real sigma, FactorGraph &fg ) { vector vars; vector factors; @@ -106,7 +106,7 @@ void MakeHOIFG( size_t N, size_t M, size_t k, double sigma, FactorGraph &fg ) { // w should be upper triangular or lower triangular -void WTh2FG( const matrix &w, const vector &th, FactorGraph &fg ) { +void WTh2FG( const matrix &w, const vector &th, FactorGraph &fg ) { vector vars; vector factors; @@ -127,7 +127,7 @@ void WTh2FG( const matrix &w, const vector &th, FactorGraph &fg ) { while( pos == w.index1_data()[i+1] ) i++; size_t j = w.index2_data()[pos]; - double w_ij = w.value_data()[pos]; + Real w_ij = w.value_data()[pos]; factors.push_back( BinaryFactor( vars[i], vars[j], w_ij ) ); } for( size_t i = 0; i < N; i++ ) @@ -137,9 +137,9 @@ void WTh2FG( const matrix &w, const vector &th, FactorGraph &fg ) { } -void MakeFullFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { +void MakeFullFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) { matrix w(N,N,N*(N-1)/2); - vector th(N,0.0); + vector th(N,0.0); for( size_t i = 0; i < N; i++ ) { for( size_t j = i+1; j < N; j++ ) @@ -151,7 +151,7 @@ void MakeFullFG( size_t N, double mean_w, double mean_th, double sigma_w, double } -void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, double beta, FactorGraph &fg ) { +void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, Real beta, FactorGraph &fg ) { vector vars; vector factors; @@ -171,11 +171,11 @@ void Make3DPotts( size_t n1, size_t n2, size_t n3, size_t states, double beta, F } -void MakeGridFG( long periodic, size_t n, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { +void MakeGridFG( long periodic, size_t n, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) { size_t N = n*n; matrix w(N,N,2*N); - vector th(N,0.0); + vector th(N,0.0); for( size_t i = 0; i < n; i++ ) for( size_t j = 0; j < n; j++ ) { @@ -190,7 +190,7 @@ void MakeGridFG( long periodic, size_t n, double mean_w, double mean_th, double } -void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, double beta, FactorGraph &fg ) { +void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, Real beta, FactorGraph &fg ) { size_t N = n*n; vector vars; @@ -214,9 +214,9 @@ void MakeGridNonbinaryFG( bool periodic, size_t n, size_t states, double beta, F } -void MakeLoopFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { +void MakeLoopFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) { matrix w(N,N,N); - vector th(N,0.0); + vector th(N,0.0); for( size_t i = 0; i < N; i++ ) { w(i, (i+1)%N) = rnd_stdnormal() * sigma_w + mean_w; @@ -227,7 +227,7 @@ void MakeLoopFG( size_t N, double mean_w, double mean_th, double sigma_w, double } -void MakeLoopNonbinaryFG( size_t N, size_t states, double beta, FactorGraph &fg ) { +void MakeLoopNonbinaryFG( size_t N, size_t states, Real beta, FactorGraph &fg ) { vector vars; vector factors; @@ -244,9 +244,9 @@ void MakeLoopNonbinaryFG( size_t N, size_t states, double beta, FactorGraph &fg } -void MakeTreeFG( size_t N, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { +void MakeTreeFG( size_t N, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) { matrix w(N,N,N-1); - vector th(N,0.0); + vector th(N,0.0); for( size_t i = 0; i < N; i++ ) { th[i] = rnd_stdnormal() * sigma_th + mean_th; @@ -260,9 +260,9 @@ void MakeTreeFG( size_t N, double mean_w, double mean_th, double sigma_w, double } -void MakeDRegFG( size_t N, size_t d, double mean_w, double mean_th, double sigma_w, double sigma_th, FactorGraph &fg ) { +void MakeDRegFG( size_t N, size_t d, Real mean_w, Real mean_th, Real sigma_w, Real sigma_th, FactorGraph &fg ) { matrix w(N,N,(d*N)/2); - vector th(N,0.0); + vector th(N,0.0); UEdgeVec g = RandomDRegularGraph( N, d ); for( size_t i = 0; i < g.size(); i++ ) @@ -413,7 +413,7 @@ BipartiteGraph CreateGroupStructuredLDPCGraph( size_t p, size_t j, size_t k ) { // Make parity check table -void MakeParityCheck( double *result, size_t n, double eps ) { +void MakeParityCheck( Real *result, size_t n, Real eps ) { size_t N = 1 << n; for( size_t i = 0; i < N; i++ ) { size_t c = 0; @@ -447,7 +447,7 @@ int main( int argc, char *argv[] ) { size_t N, K, k, d, j, n1, n2, n3; size_t prime; size_t seed; - double beta, sigma_w, sigma_th, noise, mean_w, mean_th; + Real beta, sigma_w, sigma_th, noise, mean_w, mean_th; string type; size_t states = 2; @@ -466,12 +466,12 @@ int main( int argc, char *argv[] ) { ("d", po::value(&d), "variable connectivity\n\t(only for type=='dreg')") ("j", po::value(&j), "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')") ("prime", po::value(&prime), "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')") - ("beta", po::value(&beta), "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)") - ("mean_w", po::value(&mean_w), "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") - ("mean_th", po::value(&mean_th), "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") - ("sigma_w", po::value(&sigma_w), "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") - ("sigma_th", po::value(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'") - ("noise", po::value(&noise), "bitflip probability for binary symmetric channel (only for type=='ldpc')") + ("beta", po::value(&beta), "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)") + ("mean_w", po::value(&mean_w), "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("mean_th", po::value(&mean_th), "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("sigma_w", po::value(&sigma_w), "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')") + ("sigma_th", po::value(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'") + ("noise", po::value(&noise), "bitflip probability for binary symmetric channel (only for type=='ldpc')") ("states", po::value(&states), "number of states of each variable (should be 2 for all but type=='grid', 'grid_torus', 'loop', 'potts3d')") ; @@ -703,8 +703,8 @@ int main( int argc, char *argv[] ) { // p = 29, j = 2, k = 4 // Construct likelihood and paritycheck factors - double likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise}; - double *paritycheck = new double[1 << k]; + Real likelihood[4] = {1.0 - noise, noise, noise, 1.0 - noise}; + Real *paritycheck = new Real[1 << k]; MakeParityCheck(paritycheck, k, 0.0); // Create LDPC structure diff --git a/utils/fginfo.cpp b/utils/fginfo.cpp index eaefc78..612b52c 100644 --- a/utils/fginfo.cpp +++ b/utils/fginfo.cpp @@ -108,13 +108,13 @@ int main( int argc, char *argv[] ) { cout << "Treewidth: " << tw.first << endl; cout << "Largest cluster for JTree has " << tw.second << " states " << endl; } - double stsp = 1.0; + long double stsp = 1.0; for( size_t i = 0; i < fg.nrVars(); i++ ) stsp *= fg.var(i).states(); cout << "Total state space: " << stsp << endl; - double cavsum_lcbp = 0.0; - double cavsum_lcbp2 = 0.0; + long double cavsum_lcbp = 0.0; + long double cavsum_lcbp2 = 0.0; size_t max_Delta_size = 0; map cavsizes; for( size_t i = 0; i < fg.nrVars(); i++ ) { -- 2.20.1