Miscellaneous smaller improvements
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Thu, 14 Jan 2010 12:59:47 +0000 (13:59 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Thu, 14 Jan 2010 12:59:47 +0000 (13:59 +0100)
include/dai/bbp.h
include/dai/doc.h
include/dai/jtree.h
include/dai/mr.h
include/dai/treeep.h
src/jtree.cpp
src/mr.cpp
src/treeep.cpp

index fda1036..18fd489 100644 (file)
@@ -9,7 +9,8 @@
 
 
 /// \file
-/// \brief  Defines class BBP, which implements Back-Belief-Propagation
+/// \brief Defines class BBP, which implements Back-Belief-Propagation
+/// \todo Clean up code
 
 
 #ifndef ___defined_libdai_bbp_h
index b15aa7c..dbc7904 100644 (file)
@@ -22,9 +22,9 @@
  *
  *  \todo Add FAQ
  *
- *  \todo Adapt (part of the) guidelines in http://www.boost.org/development/requirements.html#Design_and_Programming
+ *  \idea Adapt (part of the) guidelines in http://www.boost.org/development/requirements.html#Design_and_Programming
  *
- *  \todo Use "gcc -MM" to generate dependencies for targets: http://make.paulandlesley.org/autodep.html
+ *  \idea Use "gcc -MM" to generate dependencies for targets: http://make.paulandlesley.org/autodep.html
  *
  *  \todo Replace VarSets by SmallSet<size_t> where appropriate, in order to minimize the use of FactorGraph::findVar().
  *
index 37f6981..fad066b 100644 (file)
@@ -46,7 +46,7 @@ class JTree : public DAIAlgRG {
         /// Stores the messages
         std::vector<std::vector<Factor> >  _mes;
 
-        /// Stores the logarithm of the partition sum
+        /// Stores the logarithm of the partition sum (not used - why not?)
         Real _logZ;
 
     public:
@@ -150,8 +150,8 @@ class JTree : public DAIAlgRG {
         /// Returns reference to the message from outer region \a alpha to its \a _beta 'th neighboring inner region
         Factor & message( size_t alpha, size_t _beta ) { return _mes[alpha][_beta]; }
 
-        /// Runs junction tree algorithm using HUGIN updates
-        /** \note The initial messages may be arbitrary.
+        /// Runs junction tree algorithm using HUGIN (message-free) updates
+        /** \note The initial messages may be arbitrary; actually they are not used at all.
          */
         void runHUGIN();
 
index ce8a7c6..38b9377 100644 (file)
@@ -142,7 +142,7 @@ class MR : public DAIAlgFG {
         void makekindex();
         
         /// Initialize cors
-        void init_cor();
+        Real init_cor();
         
         /// Calculate cors using response propagation
         Real init_cor_resp();
index 5a3e1ac..585613e 100644 (file)
@@ -11,7 +11,8 @@
 
 /// \file
 /// \brief Defines class TreeEP, which implements Tree Expectation Propagation
-/// \todo Clean up the TreeEP code
+/// \todo Clean up the TreeEP code (by making JTree more powerful, e.g., by
+/// adding Pearl's cutset algorithm and propagation with an arbitrary root)
 
 
 #ifndef __defined_libdai_treeep_h
@@ -65,6 +66,9 @@ class TreeEP : public JTree {
 
             /// How to choose the tree
             TypeType type;
+
+            /// Optimize within-loop propagation? (currently buggy, so disabled by default)
+            bool optimize;
         } props;
 
         /// Name of this inference method
@@ -72,24 +76,37 @@ class TreeEP : public JTree {
 
     private:
         /// Stores the data structures needed to efficiently update the approximation of an off-tree factor
-        /// \todo Write documentation
         class TreeEPSubTree {
             private:
+                /// Outer region pseudomarginals (corresponding with the \f$\tilde f_i(x_j,x_k)\f$ in [\ref MiQ04])
                 std::vector<Factor>  _Qa;
+                /// Inner region pseudomarginals (corresponding with the \f$\tilde f_i(x_s)\f$ in [\ref MiQ04])
                 std::vector<Factor>  _Qb;
+                /// The junction tree (stored as a rooted tree)
                 RootedTree           _RTree;
-                std::vector<size_t>  _a;        // _Qa[alpha]  <->  superTree.Qa[_a[alpha]]
-                std::vector<size_t>  _b;        // _Qb[beta]   <->  superTree.Qb[_b[beta]]
-                                                // _Qb[beta]   <->  _RTree[beta]
+                /// Index conversion table for outer region indices (_Qa[alpha] corresponds with Qa[_a[alpha]] of the supertree)
+                std::vector<size_t>  _a;        
+                /// Index conversion table for inner region indices (_Qb[beta] corresponds with Qb[_b[beta]] of the supertree)
+                std::vector<size_t>  _b;
+                /// Pointer to off-tree factor
                 const Factor *       _I;
+                /// Variables in off-tree factor
                 VarSet               _ns;
+                /// Variables in off-tree factor which are not in the root of this subtree
                 VarSet               _nsrem;
+                /// Used for calculating the free energy
                 Real                 _logZ;
 
-
             public:
+            /// \name Constructors/destructors
+            //@{
+                /// Default constructor
                 TreeEPSubTree() : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(NULL), _ns(), _nsrem(), _logZ(0.0) {}
-                TreeEPSubTree( const TreeEPSubTree &x) : _Qa(x._Qa), _Qb(x._Qb), _RTree(x._RTree), _a(x._a), _b(x._b), _I(x._I), _ns(x._ns), _nsrem(x._nsrem), _logZ(x._logZ) {}
+
+                /// Copy constructor
+                TreeEPSubTree( const TreeEPSubTree &x ) : _Qa(x._Qa), _Qb(x._Qb), _RTree(x._RTree), _a(x._a), _b(x._b), _I(x._I), _ns(x._ns), _nsrem(x._nsrem), _logZ(x._logZ) {}
+
+                /// Assignment operator
                 TreeEPSubTree & operator=( const TreeEPSubTree& x ) {
                     if( this != &x ) {
                         _Qa         = x._Qa;
@@ -105,11 +122,23 @@ class TreeEP : public JTree {
                     return *this;
                 }
 
+                /// Construct from super tree
                 TreeEPSubTree( const RootedTree &subRTree, const RootedTree &jt_RTree, const std::vector<Factor> &jt_Qa, const std::vector<Factor> &jt_Qb, const Factor *I );
+            //@}
+
+                /// Initializes beliefs of this subtree
                 void init();
+
+                /// Inverts this approximation and multiplies it by the (super) junction tree marginals \a Qa and \a Qb
                 void InvertAndMultiply( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb );
+
+                /// Runs junction tree algorithm (including off-tree factor I) storing the results in the (super) junction tree \a Qa and \a Qb
                 void HUGIN_with_I( std::vector<Factor> &Qa, std::vector<Factor> &Qb );
+
+                /// Returns energy (?) of this subtree
                 Real logZ( const std::vector<Factor> &Qa, const std::vector<Factor> &Qb ) const;
+
+                /// Returns constant reference to the pointer to the off-tree factor
                 const Factor *& I() { return _I; }
         };
 
index bcdbefc..037e868 100644 (file)
@@ -333,14 +333,16 @@ Real JTree::run() {
 
 
 Real JTree::logZ() const {
-    Real s = 0.0;
+/*    Real s = 0.0;
     for( size_t beta = 0; beta < nrIRs(); beta++ )
         s += IR(beta).c() * Qb[beta].entropy();
     for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
         s += OR(alpha).c() * Qa[alpha].entropy();
         s += (OR(alpha).log(true) * Qa[alpha]).sum();
     }
-    return s;
+    DAI_ASSERT( abs( _logZ - s ) < 1e-8 );
+    return s;*/
+    return _logZ;
 }
 
 
index 1a74352..29ade12 100644 (file)
@@ -457,13 +457,15 @@ void MR::solveM() {
 }
 
 
-void MR::init_cor() {
+Real MR::init_cor() {
+    Real md = 0.0;
     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", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
             bpcav.makeCavity( i );
             pairq = calcPairBeliefs( bpcav, delta(i), false, true );
+            md = std::max( md, bpcav.maxDiff() );
         } else if( props.inits == Properties::InitType::EXACT ) {
             JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
             jtcav.makeCavity( i );
@@ -482,6 +484,7 @@ void MR::init_cor() {
                 }
         }
     }
+    return md;
 }
 
 
@@ -516,10 +519,16 @@ Real MR::run() {
             Real md = init_cor_resp();
             if( md > _maxdiff )
                 _maxdiff = md;
-        } else if( props.inits == Properties::InitType::EXACT )
-            init_cor(); // FIXME no MaxDiff() calculation
-        else if( props.inits == Properties::InitType::CLAMPING )
-            init_cor(); // FIXME no MaxDiff() calculation
+        } else if( props.inits == Properties::InitType::EXACT ) {
+            Real md = init_cor();
+            if( md > _maxdiff )
+                _maxdiff = md;
+        }
+        else if( props.inits == Properties::InitType::CLAMPING ) {
+            Real md = init_cor();
+            if( md > _maxdiff )
+                _maxdiff = md;
+        }
 
         solvemcav();
 
@@ -529,7 +538,7 @@ Real MR::run() {
         if( props.verbose >= 1 )
             cerr << Name << " needed " << toc() - tic << " seconds." << endl;
 
-        return 0.0;
+        return _maxdiff;
     } else
         return 1.0;
 }
index 000ce59..5e3e8db 100644 (file)
@@ -36,6 +36,11 @@ void TreeEP::setProperties( const PropertySet &opts ) {
     props.maxiter = opts.getStringAs<size_t>("maxiter");
     props.verbose = opts.getStringAs<size_t>("verbose");
     props.type = opts.getStringAs<Properties::TypeType>("type");
+
+    if( opts.hasKey("optimize") )
+        props.optimize = opts.getStringAs<bool>("optimize");
+    else
+        props.optimize = false;
 }
 
 
@@ -45,6 +50,7 @@ PropertySet TreeEP::getProperties() const {
     opts.Set( "maxiter", props.maxiter );
     opts.Set( "verbose", props.verbose );
     opts.Set( "type", props.type );
+    opts.Set( "optimize", props.optimize );
     return opts;
 }
 
@@ -55,7 +61,8 @@ string TreeEP::printProperties() const {
     s << "tol=" << props.tol << ",";
     s << "maxiter=" << props.maxiter << ",";
     s << "verbose=" << props.verbose << ",";
-    s << "type=" << props.type << "]";
+    s << "type=" << props.type << ",";
+    s << "optimize=" << props.optimize << "]";
     return s.str();
 }
 
@@ -134,10 +141,12 @@ void TreeEP::construct( const RootedTree &tree ) {
             if( offtree(I) ) {
                 // find efficient subtree
                 RootedTree subTree;
-                /*size_t subTreeSize =*/ findEfficientTree( factor(I).vars(), subTree, PreviousRoot );
+                size_t subTreeSize = findEfficientTree( factor(I).vars(), subTree, PreviousRoot );
                 PreviousRoot = subTree[0].n1;
-                //subTree.resize( subTreeSize );  // FIXME
-                //cerr << "subtree " << I << " has size " << subTreeSize << endl;
+                if( props.optimize ) {
+                    subTree.resize( subTreeSize );  // FIXME
+                    cerr << "subtree " << I << " has size " << subTreeSize << endl;
+                }
                 _Q[I] = TreeEPSubTree( subTree, RTree, Qa, Qb, &factor(I) );
                 if( repeats == 1 )
                     break;