Replaced doubles by Reals, fixed two bugs
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 27 Oct 2009 13:00:56 +0000 (14:00 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 27 Oct 2009 13:00:56 +0000 (14:00 +0100)
* Replaced all doubles by Reals where necessary
* Replaced constructor
    TFactor<T>::TFactor( const VarSet& vars, TIterator begin )
  by two constructors
    TFactor<T>::TFactor( const VarSet& vars, const T* p )
    TFactor<T>::TFactor( const VarSet& vars, const std::vector<S> &x )
* Fixed bugs (min and max were reversed) in
    TFactor<T> min( const TFactor<T> &, const TFactor<T> & )
    TFactor<T> max( const TFactor<T> &, const TFactor<T> & )

39 files changed:
ChangeLog
examples/example.cpp
examples/example_sprinkler.cpp
include/dai/bbp.h
include/dai/bp.h
include/dai/cbp.h
include/dai/daialg.h
include/dai/exactinf.h
include/dai/factor.h
include/dai/gibbs.h
include/dai/hak.h
include/dai/jtree.h
include/dai/lc.h
include/dai/mf.h
include/dai/mr.h
include/dai/prob.h
include/dai/regiongraph.h
include/dai/treeep.h
include/dai/util.h
src/bbp.cpp
src/bp.cpp
src/cbp.cpp
src/daialg.cpp
src/exactinf.cpp
src/factorgraph.cpp
src/gibbs.cpp
src/hak.cpp
src/jtree.cpp
src/lc.cpp
src/mf.cpp
src/mr.cpp
src/properties.cpp
src/regiongraph.cpp
src/treeep.cpp
src/util.cpp
tests/testbbp.cpp
tests/testdai.cpp
utils/createfg.cpp
utils/fginfo.cpp

index 17dc97e..ccafce5 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,11 @@
+* Replaced constructor
+    TFactor<T>::TFactor( const VarSet& vars, TIterator begin )
+  by two constructors
+    TFactor<T>::TFactor( const VarSet& vars, const T* p )
+    TFactor<T>::TFactor( const VarSet& vars, const std::vector<S> &x )
+* Fixed bugs (min and max were reversed) in
+  TFactor<T> min( const TFactor<T> &, const TFactor<T> & )
+  TFactor<T> max( const TFactor<T> &, const TFactor<T> & )
 * Renamed RegionGraph::Check_Counting_Numbers into
   RegionGraph::checkCountingNumbers
   (and provided an alias for backwards compatibility)
index e8ed770..210a293 100644 (file)
@@ -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;
index 85df931..0d1e570 100644 (file)
@@ -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;
index 9dc9132..ddad9c5 100644 (file)
@@ -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<size_t>
  *  \param cfn Cost function to be used.
  *  \param h controls size of perturbation.
  */
-double numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h );
+Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h );
 
 
 } // end of namespace dai
index 2e985c4..c571b4c 100644 (file)
@@ -32,18 +32,18 @@ namespace dai {
 class BP : public DAIAlgFG {
     private:
         typedef std::vector<size_t> ind_t;
-        typedef std::multimap<double, std::pair<std::size_t, std::size_t> > LutType;
+        typedef std::multimap<Real, std::pair<std::size_t, std::size_t> > LutType;
         struct EdgeProp {
             ind_t  index;
             Prob   message;
             Prob   newMessage;
-            double residual;
+            Real   residual;
         };
         std::vector<std::vector<EdgeProp> > _edges;
         std::vector<std::vector<LutType::iterator> > _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;
index 8fb446b..d9032ef 100644 (file)
@@ -49,10 +49,10 @@ class CBP : public DAIAlgFG {
         /// Factor beliefs
         std::vector<Factor> _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<size_t> clamped_vars_list, size_t &num_leaves,
-                         size_t &choose_count, double &sum_level, Real &lz_out, std::vector<Factor> &beliefs_out );
+        void runRecurse( InfAlg *bp, Real orig_logZ, std::vector<size_t> clamped_vars_list, size_t &num_leaves,
+                         size_t &choose_count, Real &sum_level, Real &lz_out, std::vector<Factor> &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<Factor> &bs, double logZ );
+        void setBeliefs( const std::vector<Factor> &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; }
 };
 
 
index 067eb13..1833e5d 100644 (file)
@@ -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
index ae4aef7..7a0de7e 100644 (file)
@@ -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; }
         //@}
 
index e058bce..70042ba 100644 (file)
@@ -98,13 +98,21 @@ template <typename T> 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<typename TIterator>
-        TFactor( const VarSet& vars, TIterator begin ) : _vs(vars), _p(begin, begin + _vs.nrStates(), _vs.nrStates()) {}
+        template<typename S>
+        TFactor( const VarSet& vars, const std::vector<S> &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<T> &p ) : _vs(vars), _p(p) {
@@ -539,7 +547,7 @@ template<typename T, typename binaryOp> TFactor<T> pointwiseOp( const TFactor<T>
             result[i] = op( result[i], g[i] );
         return result;
     } else {
-        TFactor<T> result( f.vars() | g.vars(), 0.0 );
+        TFactor<T> result( f.vars() | g.vars(), (T)0 );
 
         IndexFor i1(f.vars(), result.vars());
         IndexFor i2(g.vars(), result.vars());
@@ -584,7 +592,7 @@ template<typename T> T dist( const TFactor<T> &f, const TFactor<T> &g, Prob::Dis
  */
 template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g ) {
     DAI_ASSERT( f._vs == g._vs );
-    return TFactor<T>( f._vs, min( f.p(), g.p() ) );
+    return TFactor<T>( f._vs, max( f.p(), g.p() ) );
 }
 
 
@@ -594,7 +602,7 @@ template<typename T> TFactor<T> max( const TFactor<T> &f, const TFactor<T> &g )
  */
 template<typename T> TFactor<T> min( const TFactor<T> &f, const TFactor<T> &g ) {
     DAI_ASSERT( f._vs == g._vs );
-    return TFactor<T>( f._vs, max( f.p(), g.p() ) );
+    return TFactor<T>( f._vs, min( f.p(), g.p() ) );
 }
 
 
index c1e94d3..d33f2b8 100644 (file)
@@ -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; }
         //@}
 
index 96a2f00..4d17923 100644 (file)
@@ -38,7 +38,7 @@ class HAK : public DAIAlgRG {
         std::vector<std::vector<Factor> >  _muab;
         std::vector<std::vector<Factor> >  _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:
index 8107757..3e07356 100644 (file)
@@ -37,7 +37,7 @@ namespace dai {
 class JTree : public DAIAlgRG {
     private:
         std::vector<std::vector<Factor> >  _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; }
         //@}
 
index f4cacda..8f642e4 100644 (file)
@@ -41,7 +41,7 @@ class LC : public DAIAlgFG {
         std::vector<Factor>      _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<Factor> &Q );
 
         Factor NewPancake (size_t i, size_t _I, bool & hasNaNs);
index 644ddd4..c5fca87 100644 (file)
@@ -32,7 +32,7 @@ class MF : public DAIAlgFG {
     private:
         std::vector<Factor>  _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; }
         //@}
 
index c772a7b..62da8b9 100644 (file)
@@ -39,20 +39,20 @@ class MR : public DAIAlgFG {
 
         std::vector<size_t>                             con;       // con[i] = connectivity of spin i
         std::vector<std::vector<size_t> >               nb;        // nb[i] are the neighbours of spin i
-        std::vector<std::vector<double> >               tJ;        // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
-        std::vector<double>                             theta;     // theta[i] is the local field on spin i
-        std::vector<std::vector<double> >               M;         // M[i][_j] is M^{(i)}_j
+        std::vector<std::vector<Real> >                 tJ;        // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
+        std::vector<Real>                               theta;     // theta[i] is the local field on spin i
+        std::vector<std::vector<Real> >                 M;         // M[i][_j] is M^{(i)}_j
         std::vector<std::vector<size_t> >               kindex;    // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
-        std::vector<std::vector<std::vector<double> > > cors;
+        std::vector<std::vector<std::vector<Real> > >   cors;
 
         static const size_t kmax = 31;
         typedef boost::dynamic_bitset<> sub_nb;
 
         size_t N;
 
-        std::vector<double> Mag;
+        std::vector<Real> 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;
index 0caa188..0018bcf 100644 (file)
@@ -301,10 +301,10 @@ template <typename T> 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;
         }
 
index 231e8a4..07ffbaf 100644 (file)
@@ -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; }
 };
 
 
index b479323..86e2b55 100644 (file)
@@ -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<Factor> &Qa, const std::vector<Factor> &Qb );
                 void HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb );
-                double logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const;
+                Real logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &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; }
         //@}
 
index e2b4f13..62db510 100644 (file)
@@ -79,7 +79,7 @@
     double log1p( double x );
 
     /// Define INFINITY
-    #define INFINITY (std::numeric_limits<double>::infinity())
+    #define INFINITY (std::numeric_limits<Real>::infinity())
 #endif
 
 
@@ -197,28 +197,28 @@ std::vector<T> concat( const std::vector<T>& u, const std::vector<T>& v ) {
 void tokenizeString( const std::string& s, std::vector<std::string>& outTokens, const std::string& delim="\t\n" );
 
 /// Used to keep track of the progress made by iterative algorithms
-class Diffs : public std::vector<double> {
+class Diffs : public std::vector<Real> {
     private:
         size_t _maxsize;
-        double _def;
-        std::vector<double>::iterator _pos;
-        std::vector<double>::iterator _maxpos;
+        Real _def;
+        std::vector<Real>::iterator _pos;
+        std::vector<Real>::iterator _maxpos;
     public:
         /// Constructor
-        Diffs(long maxsize, double def) : std::vector<double>(), _maxsize(maxsize), _def(def) {
+        Diffs(long maxsize, Real def) : std::vector<Real>(), _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();
index 96f8b19..973293e 100644 (file)
@@ -678,7 +678,7 @@ std::vector<Prob> 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<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) {
+Real numericBBPTest( const InfAlg &bp, const vector<size_t> *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<size_t> *state, const Prop
 
         // for each variable i
         for( size_t i = 0; i < fg.nrVars(); i++ ) {
-            vector<double> adj_est;
+            vector<Real> 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<size_t> *state, const Prop
         for(size_t i=0; i<bp_dual.nrVars(); i++) {
             // for each factor I ~ i
             foreach(size_t I, bp_dual.nbV(i)) {
-                vector<double> adj_n_est;
+                vector<Real> adj_n_est;
                 // for each value xi
                 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                     BP_dual bp_dual_prb(bp_dual);
@@ -842,7 +842,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
                     adj_n_est.push_back((cf_prb-cf0)/h);
                 }
 
-                vector<double> adj_m_est;
+                vector<Real> adj_m_est;
                 // for each value xi
                 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                     BP_dual bp_dual_prb(bp_dual);
@@ -876,7 +876,7 @@ double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Prop
     /*    if(0) {
         // verify bbp.adj_b_V
         for(size_t i=0; i<bp_dual.nrVars(); i++) {
-            vector<double> adj_b_V_est;
+            vector<Real> adj_b_V_est;
             // for each value xi
             for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                 BP_dual bp_dual_prb(bp_dual);
@@ -965,7 +965,7 @@ void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vec
                 size_t dim = fg.factor(I).states();
                 Prob p( dim, 0.0 );
                 for( size_t xI = 0; xI < dim; xI++ ) {
-                    double bIxI = ia.beliefF(I)[xI];
+                    Real bIxI = ia.beliefF(I)[xI];
                     if( bIxI < 1.0e-15 )
                         p[xI] = -1.0e10;
                     else
@@ -982,7 +982,7 @@ void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vec
                 size_t dim = fg.var(i).states();
                 Prob p( dim, 0.0 );
                 for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
-                    double bixi = ia.beliefV(i)[xi];
+                    Real bixi = ia.beliefV(i)[xi];
                     if( bixi < 1.0e-15 )
                         p[xi] = -1.0e10;
                     else
@@ -1009,7 +1009,7 @@ void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vec
                 size_t n = fg.var(i).states();
                 Prob delta( n, 0.0 );
                 DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
-                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:
                         delta[state[i]] = 1.0;
@@ -1047,7 +1047,7 @@ void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vec
                 size_t x_I = getFactorEntryForState( fg, I, state );
                 DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
 
-                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:
                         delta[x_I] = 1.0;
@@ -1072,7 +1072,7 @@ void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vec
 
 
 Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *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<size_t> *stat
             vector<size_t> 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<size_t> *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<size_t>("verbose");
     maxiter = opts.getStringAs<size_t>("maxiter");
-    tol = opts.getStringAs<double>("tol");
-    damping = opts.getStringAs<double>("damping");
+    tol = opts.getStringAs<Real>("tol");
+    damping = opts.getStringAs<Real>("damping");
     updates = opts.getStringAs<UpdateType>("updates");
 }
 PropertySet BBP::Properties::get() const {
index e3043cc..39c70ce 100644 (file)
@@ -38,7 +38,7 @@ void BP::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("logdomain") );
     DAI_ASSERT( opts.hasKey("updates") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.logdomain = opts.getStringAs<bool>("logdomain");
     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
@@ -48,7 +48,7 @@ void BP::setProperties( const PropertySet &opts ) {
     else
         props.verbose = 0;
     if( opts.hasKey("damping") )
-        props.damping = opts.getStringAs<double>("damping");
+        props.damping = opts.getStringAs<Real>("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<size_t> 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] )
index 28c10ee..d3e960f 100644 (file)
@@ -32,7 +32,7 @@ using boost::shared_ptr;
 const char *CBP::Name = "CBP";
 
 
-void CBP::setBeliefs( const std::vector<Factor> &bs, double logZ ) {
+void CBP::setBeliefs( const std::vector<Factor> &bs, Real logZ ) {
     size_t i = 0;
     _beliefsV.clear();
     _beliefsV.reserve( nrVars() );
@@ -87,7 +87,7 @@ static vector<Factor> mixBeliefs( Real p, const vector<Factor> &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<size_t> clamped_vars_list, size_t &num_leaves,
-                      size_t &choose_count, double &sum_level, Real &lz_out, vector<Factor>& beliefs_out) {
+void CBP::runRecurse( InfAlg *bp, Real orig_logZ, vector<size_t> clamped_vars_list, size_t &num_leaves,
+                      size_t &choose_count, Real &sum_level, Real &lz_out, vector<Factor>& beliefs_out) {
     // choose a variable/states to clamp:
     size_t i;
     vector<size_t> xis;
@@ -195,7 +195,7 @@ void CBP::runRecurse( InfAlg *bp, double orig_logZ, vector<size_t> 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<double>("tol");
+    tol = opts.getStringAs<Real>("tol");
     updates = opts.getStringAs<UpdateType>("updates");
     maxiter = opts.getStringAs<size_t>("maxiter");
-    rec_tol = opts.getStringAs<double>("rec_tol");
+    rec_tol = opts.getStringAs<Real>("rec_tol");
     if(opts.hasKey("max_levels")) {
         max_levels = opts.getStringAs<size_t>("max_levels");
     } else {
         max_levels = 10;
     }
-    min_max_adj = opts.getStringAs<double>("min_max_adj");
+    min_max_adj = opts.getStringAs<Real>("min_max_adj");
     choose = opts.getStringAs<ChooseMethodType>("choose");
     recursion = opts.getStringAs<RecurseType>("recursion");
     clamp = opts.getStringAs<ClampType>("clamp");
index 241dc43..b1cd4f8 100644 (file)
@@ -130,7 +130,7 @@ vector<Factor> 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<Factor> 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
index 12b1253..312b797 100644 (file)
@@ -67,7 +67,7 @@ void ExactInf::init() {
 }
 
 
-double ExactInf::run() {
+Real ExactInf::run() {
     if( props.verbose >= 1 )
         cerr << "Starting " << identify() << "...";
 
index f94586d..b737864 100644 (file)
@@ -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<size_t> 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<size_t,Factor> 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<VarSet> 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<size_t, Factor> 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<size_t> &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<size_t, Factor> newFacs;
@@ -332,7 +332,7 @@ void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
 
 void FactorGraph::clampFactor( size_t I, const vector<size_t> &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<Factor> clamped_facs;
     for( size_t I = 0; I < nrFactors(); I++ ) {
         VarSet v_I = factor(I).vars();
index c1bff37..cd29490 100644 (file)
@@ -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();
 }
 
 
index 87fc45b..dab801a 100644 (file)
@@ -51,7 +51,7 @@ void HAK::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("doubleloop") );
     DAI_ASSERT( opts.hasKey("clusters") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.doubleloop = opts.getStringAs<bool>("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<double>("damping");
+        props.damping = opts.getStringAs<Real>("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<FRegion> org_ORs = ORs;
 
     // Save original inner counting numbers and set negative counting numbers to zero
-    vector<double> org_IR_cs( nrIRs(), 0.0 );
+    vector<Real> 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
index 976562f..61d24e1 100644 (file)
@@ -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)
index d6ba935..5c2e1e3 100644 (file)
@@ -34,7 +34,7 @@ void LC::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("cavity") );
     DAI_ASSERT( opts.hasKey("updates") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
@@ -46,7 +46,7 @@ void LC::setProperties( const PropertySet &opts ) {
     if( opts.hasKey("reinit") )
         props.reinit = opts.getStringAs<bool>("reinit");
     if( opts.hasKey("damping") )
-        props.damping = opts.getStringAs<double>("damping");
+        props.damping = opts.getStringAs<Real>("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;
 
index 23ab0e2..3e96c82 100644 (file)
@@ -30,14 +30,14 @@ void MF::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("tol") );
     DAI_ASSERT( opts.hasKey("maxiter") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     if( opts.hasKey("verbose") )
         props.verbose = opts.getStringAs<size_t>("verbose");
     else
         props.verbose = 0U;
     if( opts.hasKey("damping") )
-        props.damping = opts.getStringAs<double>("damping");
+        props.damping = opts.getStringAs<Real>("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 )
index 0045a7a..5cb1b55 100644 (file)
@@ -35,7 +35,7 @@ void MR::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("updates") );
     DAI_ASSERT( opts.hasKey("inits") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.updates = opts.getStringAs<Properties::UpdateType>("updates");
     props.inits = opts.getStringAs<Properties::InitType>("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<vector<double> > tJ_org;
+    vector<vector<Real> > tJ_org;
     vector<vector<size_t> > nb_org;
     vector<size_t> con_org;
-    vector<double> theta_org;
+    vector<Real> theta_org;
 
-    vector<double> xfield(N*kmax,0.0);
-    vector<double> rfield(N*kmax,0.0);
-    vector<double> Hfield(N,0.0);
-    vector<double> devs(N*kmax,0.0);
-    vector<double> devs2(N*kmax,0.0);
-    vector<double> dev(N,0.0);
-    vector<double> avmag(N,0.0);
+    vector<Real> xfield(N*kmax,0.0);
+    vector<Real> rfield(N*kmax,0.0);
+    vector<Real> Hfield(N,0.0);
+    vector<Real> devs(N*kmax,0.0);
+    vector<Real> devs2(N*kmax,0.0);
+    vector<Real> dev(N,0.0);
+    vector<Real> 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<con[i]; _k++) if(_k != _j) {
                         sub_nb _nbi_min_jk(_nbi_min_j);
                         _nbi_min_jk.reset(_k);
@@ -431,7 +431,7 @@ void MR::solvemcav() {
                             newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
                 }
 
-                double dev = newM - M[i][_j];
+                Real dev = newM - M[i][_j];
 //              dev *= 0.02;
                 if( fabs(dev) >= 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<Factor> 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;
index d7a87eb..628ef98 100644 (file)
@@ -25,6 +25,8 @@ std::ostream& operator<< (std::ostream & os, const Property & p) {
         os << boost::any_cast<std::string>(p.second);
     else if( p.second.type() == typeid(double) )
         os << boost::any_cast<double>(p.second);
+    else if( p.second.type() == typeid(long double) )
+        os << boost::any_cast<long double>(p.second);
     else if( p.second.type() == typeid(bool) )
         os << boost::any_cast<bool>(p.second);
     else if( p.second.type() == typeid(PropertySet) )
index 1c20e1f..5b1ef7e 100644 (file)
@@ -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<size_t>::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<Var>::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();
index 90883ac..a914223 100644 (file)
@@ -32,7 +32,7 @@ void TreeEP::setProperties( const PropertySet &opts ) {
     DAI_ASSERT( opts.hasKey("verbose") );
     DAI_ASSERT( opts.hasKey("type") );
 
-    props.tol = opts.getStringAs<double>("tol");
+    props.tol = opts.getStringAs<Real>("tol");
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.type = opts.getStringAs<Properties::TypeType>("type");
@@ -176,8 +176,8 @@ void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector<Factor> &Qa, std::vector<F
 }
 
 
-double TreeEP::TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const {
-    double s = 0.0;
+Real TreeEP::TreeEPSubTree::logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &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<double> wg;
+            WeightedGraph<Real> 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++ )
index 73fdf8f..d2a67b7 100644 (file)
@@ -75,7 +75,7 @@ boost::uniform_real<Real> _uni_dist(0,1);
 boost::variate_generator<_rnd_gen_type&, boost::uniform_real<Real> > _uni_rnd(_rnd_gen, _uni_dist);
 
 /// Normal distribution with mean 0 and standard deviation 1.
-boost::normal_distribution<> _normal_dist;
+boost::normal_distribution<Real> _normal_dist;
 
 /// Global random number generator with standard normal distribution
 boost::variate_generator<_rnd_gen_type&, boost::normal_distribution<Real> > _normal_rnd(_rnd_gen, _normal_dist);
index 9a4a8b9..45cfbdd 100644 (file)
@@ -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;
         }
     }
index 6704b66..9b61d57 100644 (file)
@@ -30,12 +30,12 @@ class TestDAI {
     protected:
         InfAlg          *obj;
         string          name;
-        vector<double>  err;
+        vector<Real>    err;
 
     public:
         vector<Factor>  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<string, PropertySet> parseMethod( const string &_s, const map<string,string
 }
 
 
-double clipdouble( double x, double minabs ) {
-    if( fabs(x) < minabs )
+Real clipReal( Real x, Real minabs ) {
+    if( abs(x) < minabs )
         return minabs;
     else
         return x;
@@ -204,7 +205,7 @@ int main( int argc, char *argv[] ) {
     string filename;
     string aliases;
     vector<string> 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<Factor> 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 ) )
index f495424..5edf476 100644 (file)
 using namespace std;
 using namespace dai;
 namespace po = boost::program_options;
-typedef boost::numeric::ublas::compressed_matrix<double> matrix;
+typedef boost::numeric::ublas::compressed_matrix<Real> 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<Var> vars;
     vector<Factor> 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<double> &th, FactorGraph &fg ) {
+void WTh2FG( const matrix &w, const vector<Real> &th, FactorGraph &fg ) {
     vector<Var>    vars;
     vector<Factor> factors;
 
@@ -127,7 +127,7 @@ void WTh2FG( const matrix &w, const vector<double> &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<double> &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<double> th(N,0.0);
+    vector<Real> 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<Var> vars;
     vector<Factor> 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<double> th(N,0.0);
+    vector<Real> 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<Var>    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<double> th(N,0.0);
+    vector<Real> 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<Var>    vars;
     vector<Factor> 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<double> th(N,0.0);
+    vector<Real> 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<double> th(N,0.0);
+    vector<Real> 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<size_t>(&d),        "variable connectivity\n\t(only for type=='dreg')")
             ("j",        po::value<size_t>(&j),        "number of parity checks per bit\n\t(only for type=='ldpc_{random,group}')")
             ("prime",    po::value<size_t>(&prime),    "prime number for construction of LDPC code\n\t(only for type=='ldpc_group')")
-            ("beta",     po::value<double>(&beta),     "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
-            ("mean_w",   po::value<double>(&mean_w),   "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("mean_th",  po::value<double>(&mean_th),  "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("sigma_w",  po::value<double>(&sigma_w),  "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
-            ("sigma_th", po::value<double>(&sigma_th), "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
-            ("noise",    po::value<double>(&noise),    "bitflip probability for binary symmetric channel (only for type=='ldpc')")
+            ("beta",     po::value<Real>(&beta),       "stddev of log-factor entries\n\t(only for type=='hoi', 'potts3d', 'grid' if states>2)")
+            ("mean_w",   po::value<Real>(&mean_w),     "mean of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("mean_th",  po::value<Real>(&mean_th),    "mean of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("sigma_w",  po::value<Real>(&sigma_w),    "stddev of pairwise interactions w_{ij}\n\t(not for type=='hoi', 'ldpc_*', 'potts3d')")
+            ("sigma_th", po::value<Real>(&sigma_th),   "stddev of singleton interactions th_i\n\t(not for type=='hoi', 'ldpc_*', 'potts3d'")
+            ("noise",    po::value<Real>(&noise),      "bitflip probability for binary symmetric channel (only for type=='ldpc')")
             ("states",   po::value<size_t>(&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
index eaefc78..612b52c 100644 (file)
@@ -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<size_t,size_t> cavsizes;
         for( size_t i = 0; i < fg.nrVars(); i++ ) {