Improved documentation of include/dai/alldai.h and daialg.h, and
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Nov 2009 15:08:26 +0000 (16:08 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Nov 2009 15:08:26 +0000 (16:08 +0100)
declared calcPairBeliefsNew() and calcMarginal2ndO() as obsolete

ChangeLog
OBSOLETE
TODO
include/dai/alldai.h
include/dai/daialg.h
src/alldai.cpp
src/daialg.cpp
src/lc.cpp
src/mr.cpp
tests/testfast.out

index 63168a6..7a70370 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,6 @@
+* Declared calcPairBeliefsNew() as obsolete (its functionality is now part
+  of calcPairBeliefs())
+* Declared calcMarginal2ndO() as obsolete
 * Removed DAI_ACCMUT macro because it didn't yield satisfactory doxygen documentation
 * Replaced constructor
     TFactor<T>::TFactor( const VarSet& vars, TIterator begin )
index 7ca1197..df7e4d3 100644 (file)
--- a/OBSOLETE
+++ b/OBSOLETE
@@ -25,3 +25,5 @@ size_t Permute::convert_linear_index( size_t li ) const;
 TProb<T> TProb<T>::sgn() const;
 bool RegionGraph::Check_Counting_Numbers();
 void RegionGraph::Calc_Counting_Numbers();
+std::vector<Factor> calcPairBeliefsNew( const InfAlg& obj, const VarSet& vs, bool reInit );
+Factor calcMarginal2ndO( const InfAlg& obj, const VarSet& vs, bool reInit );
diff --git a/TODO b/TODO
index 50d6081..633d2d7 100644 (file)
--- a/TODO
+++ b/TODO
@@ -2,8 +2,6 @@ To do for the next release (0.2.3):
 
 Improve documentation:
 
-       daialg.h
-       alldai.h
        exactinf.h
        mf.h
        bp.h
index 54d4d12..1998a24 100644 (file)
@@ -10,8 +10,7 @@
 
 
 /// \file
-/// \brief Main libDAI header file
-/// \todo Improve documentation
+/// \brief Main libDAI header file. It \#includes all other libDAI headers.
 
 
 #ifndef __defined_libdai_alldai_h
@@ -57,8 +56,8 @@
 namespace dai {
 
 
-/// Constructs a new approximate inference algorithm.
-/** \param name The name of the approximate inference algorithm (should be one of the names in DAINames).
+/// Constructs a new inference algorithm.
+/** \param name The name of the inference algorithm (should be one of the names in DAINames).
  *  \param fg The FactorGraph that the algorithm should be applied to.
  *  \param opts A PropertySet specifying the options for the algorithm.
  *  \return Returns a pointer to the new InfAlg object; it is the responsibility of the caller to delete it later.
@@ -66,15 +65,16 @@ namespace dai {
 InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const PropertySet &opts );
 
 
-/// Constructs a new approximate inference algorithm.
-/** \param nameOpts The name and options of the approximate inference algorithm (should be in the format "name[opts]").
+/// Constructs a new inference algorithm.
+/** \param nameOpts The name and options of the inference algorithm (should be in the format "name[key1=val1,key2=val2,...,keyn=valn]").
  *  \param fg The FactorGraph that the algorithm should be applied to.
  *  \return Returns a pointer to the new InfAlg object; it is the responsibility of the caller to delete it later.
+ *  \todo Support aliases like in testdai
  */
 InfAlg *newInfAlgFromString( const std::string &nameOpts, const FactorGraph &fg );
 
 
-/// Contains the names of all approximate inference algorithms compiled into libDAI.
+/// Contains the names of all inference algorithms compiled into libDAI.
 static const char* DAINames[] = {
     ExactInf::Name,
 #ifdef DAI_WITH_BP
index 1833e5d..bd0a8e0 100644 (file)
@@ -10,8 +10,7 @@
 
 
 /// \file
-/// \brief Defines abstract base class InfAlg, its descendants DAIAlg<>, the specializations DAIAlgFG and DAIAlgRG and some generic inference methods.
-/// \todo Improve documentation
+/// \brief Defines abstract base class InfAlg, its descendant DAIAlg<>, the specializations DAIAlgFG and DAIAlgRG and some generic inference methods.
 
 
 #ifndef __defined_libdai_daialg_h
@@ -29,66 +28,100 @@ namespace dai {
 
 
 /// InfAlg is an abstract base class, defining the common interface of all inference algorithms in libDAI.
-/** \todo General marginalization functions like calcMarginal now copy a complete InfAlg object. Instead,
+/** \todo General marginalization functions like calcMarginal() now copy a complete InfAlg object. Instead,
  *  it would make more sense that they construct a new object without copying the FactorGraph or RegionGraph.
  *  Or they can simply be made methods of the general InfAlg class.
  *  \idea Use a PropertySet as output of an InfAlg, instead of functions like maxDiff() and Iterations().
  */
 class InfAlg {
     public:
-        /// Virtual desctructor (needed because this class contains virtual functions)
+    /// \name Constructors/destructors
+    //@{
+        /// Virtual destructor (needed because this class contains virtual functions)
         virtual ~InfAlg() {}
 
-    public:
-        /// Returns a pointer to a new, cloned copy of *this (i.e., virtual copy constructor)
+        /// Returns a pointer to a new, cloned copy of \c *this (i.e., virtual copy constructor)
         virtual InfAlg* clone() const = 0;
+    //@}
 
+    /// \name Queries
+    //@{
         /// Identifies itself for logging purposes
         virtual std::string identify() const = 0;
 
-        /// Returns the "belief" (i.e., approximate marginal probability distribution) of a variable
-        virtual Factor belief( const Var &n ) const = 0;
+        /// Returns reference to underlying FactorGraph.
+        virtual FactorGraph &fg() = 0;
+
+        /// Returns constant reference to underlying FactorGraph.
+        virtual const FactorGraph &fg() const = 0;
+    //@}
+
+    /// \name Inference interface
+    //@{
+        /// Initializes all data structures of the approximate inference algorithm.
+        /** \note This method should be called at least once before run() is called.
+         */
+        virtual void init() = 0;
+
+        /// Initializes all data structures corresponding to some set of variables.
+        /** This method can be used to do a partial initialization after a part of the factor graph has changed.
+         *  Instead of initializing all data structures, it only initializes those involving the variables in \a vs.
+         */
+        virtual void init( const VarSet &vs ) = 0;
 
-        /// Returns the "belief" (i.e., approximate marginal probability distribution) of a set of variables
-        virtual Factor belief( const VarSet &n ) const = 0;
+        /// Runs the approximate inference algorithm.
+        /** \note Before run() is called the first time, init() should have been called.
+         */
+        virtual Real run() = 0;
+
+        /// Returns the (approximate) marginal probability distribution of a variable.
+        /** \note Before this method is called, run() should have been called.
+         */
+        virtual Factor belief( const Var &v ) const = 0;
 
-        /// Returns marginal for a variable.
-        /** Sometimes preferred to belief() for performance reasons.
-          * Faster implementations exist in e.g. BP.
-          */
+        /// Returns the (approximate) marginal probability distribution of a set of variables.
+        /** \note Before this method is called, run() should have been called.
+         */
+        virtual Factor belief( const VarSet &vs ) const = 0;
+
+        /// Returns the (approximate) marginal probability distribution of the variable with index \a i.
+        /** For some approximate inference algorithms, using beliefV() is preferred to belief() for performance reasons.
+         *  \note Before this method is called, run() should have been called.
+         */
         virtual Factor beliefV( size_t i ) const { return belief( fg().var(i) ); }
 
-        /// Returns marginal for a factor.
-        /** Sometimes preferred to belief() for performance reasons.
-          * Faster implementations exist in e.g. BP.
-          */
+        /// Returns the (approximate) marginal probability distribution of the variables on which factor \a I depends.
+        /** For some approximate inference algorithms, using beliefF() is preferred to belief() for performance reasons.
+         *  \note Before this method is called, run() should have been called.
+         */
         virtual Factor beliefF( size_t I ) const { return belief( fg().factor(I).vars() ); }
 
-        /// Returns all "beliefs" (i.e., approximate marginal probability distribution) calculated by the algorithm
+        /// Returns all beliefs (approximate marginal probability distributions) calculated by the algorithm.
+        /** \note Before this method is called, run() should have been called.
+         */
         virtual std::vector<Factor> beliefs() const = 0;
 
-        /// Returns the logarithm of the (approximated) partition sum (normalizing constant of the factor graph)
-        virtual Real logZ() const = 0;
-
-        /// Initializes all data structures of the approximate inference algorithm
-        /** This method should be called at least once before run() is called
+        /// Returns the logarithm of the (approximated) partition sum (normalizing constant of the factor graph).
+        /** \note Before this method is called, run() should have been called.
+         *  \throw NOT_IMPLEMENTED if not implemented/supported
          */
-        virtual void init() = 0;
+        virtual Real logZ() const = 0;
 
-        /// Initializes all data structures corresponding to some set of variables
-        /** This method can be used to do a partial initialization after a part of the factor graph has changed.
-         *  Instead of initializing all data structures, it only initializes those involving the variables in ns.
+        /// Returns maximum difference between single variable beliefs in the last iteration.
+        /** \throw NOT_IMPLEMENTED if not implemented/supported
          */
-        virtual void init( const VarSet &ns ) = 0;
+        virtual Real maxDiff() const = 0;
 
-        /// Runs the approximate inference algorithm
-        /*  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().
+        /// Returns number of iterations done (one iteration passes over the complete factorgraph).
+        /** \throw NOT_IMPLEMENTED if not implemented/supported
          */
-        virtual Real run() = 0;
+        virtual size_t Iterations() const = 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
+    /// \name Changing the factor graph
+    //@{
+        /// Clamp variable with index \a i to value \a x (i.e. multiply with a Kronecker delta \f$\delta_{x_i, x}\f$)
+        /** If \a backup == \c true, make a backup of all factors that are changed.
          */
         virtual void clamp( size_t i, size_t x, bool backup = false ) = 0;
 
@@ -96,38 +129,32 @@ class InfAlg {
         /// Only for backwards compatibility (to be removed soon)
         virtual void clamp( const Var &v, size_t x, bool backup = false ) = 0;
 
-        /// Set all factors interacting with var(i) to 1
+        /// Sets all factors interacting with variable with index \a i to one.
+        /** If \a backup == \c true, make a backup of all factors that are changed.
+         */
         virtual void makeCavity( size_t i, bool backup = false ) = 0;
+    //@}
 
-        /// Return maximum difference between single node beliefs in the last pass
-        /// \throw Exception if not implemented/supported
-        virtual Real maxDiff() const = 0;
-
-        /// Return number of passes over the factorgraph
-        /// \throw Exception if not implemented/supported
-        virtual size_t Iterations() const = 0;
-
-
-        /// Get reference to underlying FactorGraph
-        virtual FactorGraph &fg() = 0;
-
-        /// Get const reference to underlying FactorGraph
-        virtual const FactorGraph &fg() const = 0;
-
-        /// Save factor I
+    /// \name Backup/restore mechanism for factors
+    //@{
+        /// Make a backup copy of factor \a I
         virtual void backupFactor( size_t I ) = 0;
-        /// Save Factors involving ns
-        virtual void backupFactors( const VarSet &ns ) = 0;
+        /// Make backup copies of all factors involving the variables in \a vs
+        virtual void backupFactors( const VarSet &vs ) = 0;
 
-        /// Restore factor I
+        /// Restore factor \a I from its backup copy
         virtual void restoreFactor( size_t I ) = 0;
-        /// Restore Factors involving ns
-        virtual void restoreFactors( const VarSet &ns ) = 0;
+        /// Restore the factors involving the variables in \a vs from their backup copies
+        virtual void restoreFactors( const VarSet &vs ) = 0;
+    //@}
 };
 
 
-/// Combines an InfAlg and a graphical model, e.g., a FactorGraph or RegionGraph
-/** \tparam GRM Should be castable to FactorGraph
+/// Combines the abstract base class InfAlg with a graphical model (e.g., a FactorGraph or RegionGraph).
+/** Inference algorithms in libDAI directly inherit from a DAIAlg, currently either
+ *  from a DAIAlg<FactorGraph> or from a DAIAlg<RegionGraph>.
+ *
+ *  \tparam GRM Should be castable to FactorGraph
  *  \todo A DAIAlg should not inherit from a FactorGraph or RegionGraph, but should
  *  store a reference to the graphical model object. This prevents needless copying
  *  of (possibly large) data structures. Disadvantage: the caller must not change
@@ -137,23 +164,29 @@ class InfAlg {
 template <class GRM>
 class DAIAlg : public InfAlg, public GRM {
     public:
+    /// \name Constructors/destructors
+    //@{
         /// Default constructor
         DAIAlg() : InfAlg(), GRM() {}
 
         /// Construct from GRM
         DAIAlg( const GRM &grm ) : InfAlg(), GRM(grm) {}
+    //@}
 
-        /// Save factor I
-        void backupFactor( size_t I ) { GRM::backupFactor( I ); }
-        /// Save Factors involving ns
-        void backupFactors( const VarSet &ns ) { GRM::backupFactors( ns ); }
+    /// \name Queries
+    //@{
+        /// Returns reference to underlying FactorGraph.
+        FactorGraph &fg() { return (FactorGraph &)(*this); }
 
-        /// Restore factor I
-        void restoreFactor( size_t I ) { GRM::restoreFactor( I ); }
-        /// Restore Factors involving ns
-        void restoreFactors( const VarSet &ns ) { GRM::restoreFactors( ns ); }
+        /// Returns constant reference to underlying FactorGraph.
+        const FactorGraph &fg() const { return (const FactorGraph &)(*this); }
+    //@}
 
-        /// Clamp variable with index i to value x (i.e. multiply with a Kronecker delta \f$\delta_{x_i, x}\f$)
+    /// \name Changing the factor graph
+    //@{
+        /// Clamp variable with index \a i to value \a x (i.e. multiply with a Kronecker delta \f$\delta_{x_i, x}\f$)
+        /** If \a backup == \c true, make a backup of all factors that are changed.
+         */
         void clamp( size_t i, size_t x, bool backup = false ) { GRM::clamp( i, x, backup ); }
 
         // OBSOLETE
@@ -163,14 +196,24 @@ class DAIAlg : public InfAlg, public GRM {
             std::cerr << "Warning: this DAIAlg<...>::clamp(const Var&,...) interface is obsolete!" << std::endl;
         }
 
-        /// Set all factors interacting with var(i) to 1
+        /// Sets all factors interacting with variable with index \a i to one.
+        /** If \a backup == \c true, make a backup of all factors that are changed.
+         */
         void makeCavity( size_t i, bool backup = false ) { GRM::makeCavity( i, backup ); }
+    //@}
 
-        /// Get reference to underlying FactorGraph
-        FactorGraph &fg() { return (FactorGraph &)(*this); }
+    /// \name Backup/restore mechanism for factors
+    //@{
+        /// Make a backup copy of factor \a I
+        void backupFactor( size_t I ) { GRM::backupFactor( I ); }
+        /// Make backup copies of all factors involving the variables in \a vs
+        void backupFactors( const VarSet &vs ) { GRM::backupFactors( vs ); }
 
-        /// Get const reference to underlying FactorGraph
-        const FactorGraph &fg() const { return (const FactorGraph &)(*this); }
+        /// Restore factor \a I from its backup copy
+        void restoreFactor( size_t I ) { GRM::restoreFactor( I ); }
+        /// Restore the factors involving the variables in \a vs from their backup copies
+        void restoreFactors( const VarSet &vs ) { GRM::restoreFactors( vs ); }
+    //@}
 };
 
 
@@ -181,10 +224,36 @@ typedef DAIAlg<FactorGraph> DAIAlgFG;
 typedef DAIAlg<RegionGraph> DAIAlgRG;
 
 
-Factor calcMarginal( const InfAlg & obj, const VarSet & ns, bool reInit );
-std::vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit );
-std::vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit );
-Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit );
+/// Calculates the marginal probability distribution for \a vs using inference algorithm \a obj.
+/** calcMarginal() works by clamping all variables in \a vs and calculating the partition sum for each clamped state.
+ *  Therefore, it can be used in combination with any inference algorithm that can calculate/approximate partition sums.
+ *  \param obj instance of inference algorithm to be used 
+ *  \param vs variables for which the marginal should be calculated
+ *  \param reInit should be set to \c true if at least one of the possible clamped states would be invalid (leading to a factor graph with zero partition sum).
+ */
+Factor calcMarginal( const InfAlg& obj, const VarSet& vs, bool reInit );
+
+/// Calculates beliefs for all pairs of variables in \a vs using inference algorithm \a obj.
+/** calcPairBeliefs() works by 
+ *  - clamping single variables in \a vs and calculating the partition sum and the single variable beliefs for each clamped state, if \a accurate == \c false;
+ *  - clamping pairs of variables in \a vs and calculating the partition sum for each clamped state, if \a accurate == \c true.
+ *
+ *  Therefore, it can be used in combination with any inference algorithm that can calculate/approximate partition sums (and single variable beliefs, if
+ *  \a accurate == \c true).
+ *  \param obj instance of inference algorithm to be used 
+ *  \param vs variables for which the pair beliefs should be calculated
+ *  \param reInit should be set to \c true if at least one of the possible clamped states would be invalid (leading to a factor graph with zero partition sum).
+ *  \param accurate if \c true, uses a slower but more accurate approximation algorithm
+ */
+std::vector<Factor> calcPairBeliefs( const InfAlg& obj, const VarSet& vs, bool reInit, bool accurate=false );
+
+// OBSOLETE
+/// Only for backwards compatibility (to be removed soon)
+std::vector<Factor> calcPairBeliefsNew( const InfAlg& obj, const VarSet& vs, bool reInit );
+
+// OBSOLETE
+/// Only for backwards compatibility (to be removed soon)
+Factor calcMarginal2ndO( const InfAlg& obj, const VarSet& vs, bool reInit );
 
 
 } // end of namespace dai
index e1c4119..fe12e3a 100644 (file)
@@ -64,7 +64,6 @@ InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const Propert
 }
 
 
-/// \todo Make alias file non-testdai-specific, and use it in newInfAlgFromString
 InfAlg *newInfAlgFromString( const std::string &nameOpts, const FactorGraph &fg ) {
     string::size_type pos = nameOpts.find_first_of('[');
     string name;
index b1cd4f8..a04bcc3 100644 (file)
@@ -19,34 +19,31 @@ namespace dai {
 using namespace std;
 
 
-/// Calculates the marginal of obj on ns by clamping all variables in ns and calculating logZ for each joined state.
-/*  reInit should be set to true if at least one of the possible clamped states would be invalid (leading to a factor graph with zero partition sum).
- */
-Factor calcMarginal( const InfAlg &obj, const VarSet &ns, bool reInit ) {
-    Factor Pns (ns);
+Factor calcMarginal( const InfAlg &obj, const VarSet &vs, bool reInit ) {
+    Factor Pvs (vs);
 
     InfAlg *clamped = obj.clone();
     if( !reInit )
         clamped->init();
 
     map<Var,size_t> varindices;
-    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
+    for( VarSet::const_iterator n = vs.begin(); n != vs.end(); n++ )
         varindices[*n] = obj.fg().findVar( *n );
 
     Real logZ0 = -INFINITY;
-    for( State s(ns); s.valid(); s++ ) {
-        // save unclamped factors connected to ns
-        clamped->backupFactors( ns );
+    for( State s(vs); s.valid(); s++ ) {
+        // save unclamped factors connected to vs
+        clamped->backupFactors( vs );
 
         // set clamping Factors to delta functions
-        for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
+        for( VarSet::const_iterator n = vs.begin(); n != vs.end(); n++ )
             clamped->clamp( varindices[*n], s(*n) );
 
-        // run DAIAlg, calc logZ, store in Pns
+        // run DAIAlg, calc logZ, store in Pvs
         if( reInit )
             clamped->init();
         else
-            clamped->init(ns);
+            clamped->init(vs);
 
         Real logZ;
         try {
@@ -64,195 +61,168 @@ Factor calcMarginal( const InfAlg &obj, const VarSet &ns, bool reInit ) {
                 logZ0 = logZ;
 
         if( logZ == -INFINITY )
-            Pns[s] = 0;
+            Pvs[s] = 0;
         else
-            Pns[s] = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
+            Pvs[s] = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
 
         // restore clamped factors
-        clamped->restoreFactors( ns );
+        clamped->restoreFactors( vs );
     }
 
     delete clamped;
 
-    return( Pns.normalized() );
+    return( Pvs.normalized() );
 }
 
 
-/// Calculates beliefs of all pairs in ns (by clamping nodes in ns and calculating logZ and the beliefs for each state).
-/*  reInit should be set to true if at least one of the possible clamped states would be invalid (leading to a factor graph with zero partition sum).
- */
-vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& ns, bool reInit ) {
-    // convert ns to vector<VarSet>
-    size_t N = ns.size();
-    vector<Var> vns;
-    vns.reserve( N );
-    map<Var,size_t> varindices;
-    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) {
-        vns.push_back( *n );
-        varindices[*n] = obj.fg().findVar( *n );
-    }
-
-    vector<Factor> pairbeliefs;
-    pairbeliefs.reserve( N * N );
-    for( size_t j = 0; j < N; j++ )
-        for( size_t k = 0; k < N; k++ )
-            if( j == k )
-                pairbeliefs.push_back( Factor() );
-            else
-                pairbeliefs.push_back( Factor( VarSet(vns[j], vns[k]) ) );
+vector<Factor> calcPairBeliefs( const InfAlg & obj, const VarSet& vs, bool reInit, bool accurate ) {
+    vector<Factor> result;
+    size_t N = vs.size();
+    result.reserve( N * (N - 1) / 2 );
 
     InfAlg *clamped = obj.clone();
     if( !reInit )
         clamped->init();
 
-    Real logZ0 = -INFINITY;
-    for( size_t j = 0; j < N; j++ ) {
-        // clamp Var j to its possible values
-        for( size_t j_val = 0; j_val < vns[j].states(); j_val++ ) {
-            clamped->clamp( varindices[vns[j]], j_val, true );
-            if( reInit )
-                clamped->init();
-            else
-                clamped->init(ns);
-
-            Real logZ;
-            try {
-                clamped->run();
-                logZ = clamped->logZ();
-            } catch( Exception &e ) {
-                if( e.code() == Exception::NOT_NORMALIZABLE )
-                    logZ = -INFINITY;
-                else
-                    throw;
-            }
+    map<Var,size_t> varindices;
+    for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ )
+        varindices[*v] = obj.fg().findVar( *v );
+
+    if( accurate ) {
+        Real logZ0 = 0.0;
+        VarSet::const_iterator nj = vs.begin();
+        for( long j = 0; j < (long)N - 1; j++, nj++ ) {
+            size_t k = 0;
+            for( VarSet::const_iterator nk = nj; (++nk) != vs.end(); k++ ) {
+                Factor pairbelief( VarSet(*nj, *nk) );
+
+                // clamp Vars j and k to their possible values
+                for( size_t j_val = 0; j_val < nj->states(); j_val++ )
+                    for( size_t k_val = 0; k_val < nk->states(); k_val++ ) {
+                        // save unclamped factors connected to vs
+                        clamped->backupFactors( vs );
+
+                        clamped->clamp( varindices[*nj], j_val );
+                        clamped->clamp( varindices[*nk], k_val );
+                        if( reInit )
+                            clamped->init();
+                        else
+                            clamped->init(vs);
+
+                        Real logZ;
+                        try {
+                            clamped->run();
+                            logZ = clamped->logZ();
+                        } catch( Exception &e ) {
+                            if( e.code() == Exception::NOT_NORMALIZABLE )
+                                logZ = -INFINITY;
+                            else
+                                throw;
+                        }
+
+                        if( logZ0 == -INFINITY )
+                            if( logZ != -INFINITY )
+                                logZ0 = logZ;
+
+                        Real Z_xj;
+                        if( logZ == -INFINITY )
+                            Z_xj = 0;
+                        else
+                            Z_xj = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
 
-            if( logZ0 == -INFINITY )
-                if( logZ != -INFINITY )
-                    logZ0 = logZ;
+                        // we assume that j.label() < k.label()
+                        // i.e. we make an assumption here about the indexing
+                        pairbelief[j_val + (k_val * nj->states())] = Z_xj;
 
-            Real Z_xj;
-            if( logZ == -INFINITY )
-                Z_xj = 0;
-            else
-                Z_xj = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
+                        // restore clamped factors
+                        clamped->restoreFactors( vs );
+                    }
+
+                result.push_back( pairbelief.normalized() );
+            }
+        }
+    } else {
+        // convert vs to vector<VarSet>
+        vector<Var> vvs( vs.begin(), vs.end() );
 
+        vector<Factor> pairbeliefs;
+        pairbeliefs.reserve( N * N );
+        for( size_t j = 0; j < N; j++ )
             for( size_t k = 0; k < N; k++ )
-                if( k != j ) {
-                    Factor b_k = clamped->belief(vns[k]);
-                    for( size_t k_val = 0; k_val < vns[k].states(); k_val++ )
-                        if( vns[j].label() < vns[k].label() )
-                            pairbeliefs[j * N + k][j_val + (k_val * vns[j].states())] = Z_xj * b_k[k_val];
-                        else
-                            pairbeliefs[j * N + k][k_val + (j_val * vns[k].states())] = Z_xj * b_k[k_val];
+                if( j == k )
+                    pairbeliefs.push_back( Factor() );
+                else
+                    pairbeliefs.push_back( Factor( VarSet(vvs[j], vvs[k]) ) );
+
+        Real logZ0 = -INFINITY;
+        for( size_t j = 0; j < N; j++ ) {
+            // clamp Var j to its possible values
+            for( size_t j_val = 0; j_val < vvs[j].states(); j_val++ ) {
+                clamped->clamp( varindices[vvs[j]], j_val, true );
+                if( reInit )
+                    clamped->init();
+                else
+                    clamped->init(vs);
+
+                Real logZ;
+                try {
+                    clamped->run();
+                    logZ = clamped->logZ();
+                } catch( Exception &e ) {
+                    if( e.code() == Exception::NOT_NORMALIZABLE )
+                        logZ = -INFINITY;
+                    else
+                        throw;
                 }
 
-            // restore clamped factors
-            clamped->restoreFactors( ns );
-        }
-    }
+                if( logZ0 == -INFINITY )
+                    if( logZ != -INFINITY )
+                        logZ0 = logZ;
 
-    delete clamped;
+                Real Z_xj;
+                if( logZ == -INFINITY )
+                    Z_xj = 0;
+                else
+                    Z_xj = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
+
+                for( size_t k = 0; k < N; k++ )
+                    if( k != j ) {
+                        Factor b_k = clamped->belief(vvs[k]);
+                        for( size_t k_val = 0; k_val < vvs[k].states(); k_val++ )
+                            if( vvs[j].label() < vvs[k].label() )
+                                pairbeliefs[j * N + k][j_val + (k_val * vvs[j].states())] = Z_xj * b_k[k_val];
+                            else
+                                pairbeliefs[j * N + k][k_val + (j_val * vvs[k].states())] = Z_xj * b_k[k_val];
+                    }
 
-    // Calculate result by taking the geometric average
-    vector<Factor> result;
-    result.reserve( N * (N - 1) / 2 );
-    for( size_t j = 0; j < N; j++ )
-        for( size_t k = j+1; k < N; k++ )
-            result.push_back( ((pairbeliefs[j * N + k] * pairbeliefs[k * N + j]) ^ 0.5).normalized() );
+                // restore clamped factors
+                clamped->restoreFactors( vs );
+            }
+        }
 
+        // Calculate result by taking the geometric average
+        for( size_t j = 0; j < N; j++ )
+            for( size_t k = j+1; k < N; k++ )
+                result.push_back( ((pairbeliefs[j * N + k] * pairbeliefs[k * N + j]) ^ 0.5).normalized() );
+    }
+    delete clamped;
     return result;
 }
 
 
-/// Calculates beliefs of all pairs in ns (by clamping pairs in ns and calculating logZ for each joined state).
-/*  reInit should be set to true if at least one of the possible clamped states would be invalid (leading to a factor graph with zero partition sum).
- */
-Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& ns, bool reInit ) {
-    // returns a a probability distribution whose 1st order interactions
-    // are unspecified, whose 2nd order interactions approximate those of
-    // the marginal on ns, and whose higher order interactions are absent.
-
-    vector<Factor> pairbeliefs = calcPairBeliefs( obj, ns, reInit );
-
-    Factor Pns (ns);
-    for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
-        Pns *= pairbeliefs[ij];
-
-    return( Pns.normalized() );
+std::vector<Factor> calcPairBeliefsNew( const InfAlg& obj, const VarSet& vs, bool reInit ) { 
+    return calcPairBeliefs( obj, vs, reInit, true );
 }
 
 
-/// Calculates 2nd order interactions of the marginal of obj on ns.
-vector<Factor> calcPairBeliefsNew( const InfAlg & obj, const VarSet& ns, bool reInit ) {
-    vector<Factor> result;
-    result.reserve( ns.size() * (ns.size() - 1) / 2 );
+Factor calcMarginal2ndO( const InfAlg & obj, const VarSet& vs, bool reInit ) {
+    vector<Factor> pairbeliefs = calcPairBeliefs( obj, vs, reInit );
 
-    InfAlg *clamped = obj.clone();
-    if( !reInit )
-        clamped->init();
-
-    map<Var,size_t> varindices;
-    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
-        varindices[*n] = obj.fg().findVar( *n );
-
-    Real logZ0 = 0.0;
-    VarSet::const_iterator nj = ns.begin();
-    for( long j = 0; j < (long)ns.size() - 1; j++, nj++ ) {
-        size_t k = 0;
-        for( VarSet::const_iterator nk = nj; (++nk) != ns.end(); k++ ) {
-            Factor pairbelief( VarSet(*nj, *nk) );
-
-            // clamp Vars j and k to their possible values
-            for( size_t j_val = 0; j_val < nj->states(); j_val++ )
-                for( size_t k_val = 0; k_val < nk->states(); k_val++ ) {
-                    // save unclamped factors connected to ns
-                    clamped->backupFactors( ns );
-
-                    clamped->clamp( varindices[*nj], j_val );
-                    clamped->clamp( varindices[*nk], k_val );
-                    if( reInit )
-                        clamped->init();
-                    else
-                        clamped->init(ns);
-
-                    Real logZ;
-                    try {
-                        clamped->run();
-                        logZ = clamped->logZ();
-                    } catch( Exception &e ) {
-                        if( e.code() == Exception::NOT_NORMALIZABLE )
-                            logZ = -INFINITY;
-                        else
-                            throw;
-                    }
-
-                    if( logZ0 == -INFINITY )
-                        if( logZ != -INFINITY )
-                            logZ0 = logZ;
-
-                    Real Z_xj;
-                    if( logZ == -INFINITY )
-                        Z_xj = 0;
-                    else
-                        Z_xj = exp(logZ - logZ0); // subtract logZ0 to avoid very large numbers
-
-                    // we assume that j.label() < k.label()
-                    // i.e. we make an assumption here about the indexing
-                    pairbelief[j_val + (k_val * nj->states())] = Z_xj;
-
-                    // restore clamped factors
-                    clamped->restoreFactors( ns );
-                }
-
-            result.push_back( pairbelief.normalized() );
-        }
-    }
-
-    delete clamped;
-
-    DAI_ASSERT( result.size() == (ns.size() * (ns.size() - 1) / 2) );
+    Factor Pvs (vs);
+    for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
+        Pvs *= pairbeliefs[ij];
 
-    return result;
+    return( Pvs.normalized() );
 }
 
 
index 5c2e1e3..4595c8f 100644 (file)
@@ -134,10 +134,12 @@ Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &o
 
         if( props.cavity == Properties::CavityType::FULL )
             Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
-        else if( props.cavity == Properties::CavityType::PAIR )
-            Bi = calcMarginal2ndO( *cav, cav->fg().delta(i), props.reinit );
-        else if( props.cavity == Properties::CavityType::PAIR2 ) {
-            vector<Factor> pairbeliefs = calcPairBeliefsNew( *cav, cav->fg().delta(i), props.reinit );
+        else if( props.cavity == Properties::CavityType::PAIR ) {
+            vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, false );
+            for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
+                Bi *= pairbeliefs[ij];
+        } else if( props.cavity == Properties::CavityType::PAIR2 ) {
+            vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, true );
             for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
                 Bi *= pairbeliefs[ij];
         }
index 5cb1b55..fde5f0f 100644 (file)
@@ -488,11 +488,11 @@ void MR::init_cor() {
         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 );
+            pairq = calcPairBeliefs( bpcav, delta(i), false, true );
         } else if( props.inits == Properties::InitType::EXACT ) {
             JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
             jtcav.makeCavity( i );
-            pairq = calcPairBeliefs( jtcav, delta(i), false );
+            pairq = calcPairBeliefs( jtcav, delta(i), false, true );
         }
         for( size_t jk = 0; jk < pairq.size(); jk++ ) {
             VarSet::const_iterator kit = pairq[jk].vars().begin();
index 74c39f1..ae8f3bb 100644 (file)
@@ -459,23 +459,23 @@ MR_RESPPROP_FULL                          1.676e-02       4.933e-03       N/A             1.000e-09
 # ({x13}, (5.407e-01, 4.593e-01))
 # ({x14}, (6.252e-01, 3.748e-01))
 # ({x15}, (1.346e-01, 8.654e-01))
-MR_CLAMPING_FULL                               8.359e-03       2.508e-03       N/A             1.000e-09       
-# ({x0}, (3.900e-01, 6.100e-01))
-# ({x1}, (5.560e-01, 4.440e-01))
-# ({x2}, (4.589e-01, 5.411e-01))
-# ({x3}, (5.480e-01, 4.520e-01))
-# ({x4}, (6.626e-01, 3.374e-01))
-# ({x5}, (2.105e-01, 7.895e-01))
-# ({x6}, (8.168e-01, 1.832e-01))
-# ({x7}, (2.332e-01, 7.668e-01))
-# ({x8}, (2.242e-01, 7.758e-01))
-# ({x9}, (2.136e-01, 7.864e-01))
-# ({x10}, (7.621e-01, 2.379e-01))
-# ({x11}, (1.206e-01, 8.794e-01))
-# ({x12}, (4.180e-01, 5.820e-01))
-# ({x13}, (5.404e-01, 4.596e-01))
-# ({x14}, (6.273e-01, 3.727e-01))
-# ({x15}, (1.344e-01, 8.656e-01))
+MR_CLAMPING_FULL                               6.060e-03       1.754e-03       N/A             1.000e-09       
+# ({x0}, (3.872e-01, 6.128e-01))
+# ({x1}, (5.585e-01, 4.415e-01))
+# ({x2}, (4.575e-01, 5.425e-01))
+# ({x3}, (5.486e-01, 4.514e-01))
+# ({x4}, (6.672e-01, 3.328e-01))
+# ({x5}, (2.118e-01, 7.882e-01))
+# ({x6}, (8.174e-01, 1.826e-01))
+# ({x7}, (2.329e-01, 7.671e-01))
+# ({x8}, (2.189e-01, 7.811e-01))
+# ({x9}, (2.075e-01, 7.925e-01))
+# ({x10}, (7.651e-01, 2.349e-01))
+# ({x11}, (1.201e-01, 8.799e-01))
+# ({x12}, (4.174e-01, 5.826e-01))
+# ({x13}, (5.409e-01, 4.591e-01))
+# ({x14}, (6.292e-01, 3.708e-01))
+# ({x15}, (1.338e-01, 8.662e-01))
 MR_EXACT_FULL                                  3.527e-03       1.038e-03       N/A             1.000e-09       
 # ({x0}, (3.867e-01, 6.133e-01))
 # ({x1}, (5.572e-01, 4.428e-01))
@@ -510,23 +510,23 @@ MR_RESPPROP_LINEAR                        1.932e-02       5.506e-03       N/A             1.000e-09
 # ({x13}, (5.388e-01, 4.612e-01))
 # ({x14}, (6.239e-01, 3.761e-01))
 # ({x15}, (1.369e-01, 8.631e-01))
-MR_CLAMPING_LINEAR                             1.089e-02       3.076e-03       N/A             1.000e-09       
-# ({x0}, (3.884e-01, 6.116e-01))
-# ({x1}, (5.551e-01, 4.449e-01))
-# ({x2}, (4.598e-01, 5.402e-01))
-# ({x3}, (5.469e-01, 4.531e-01))
-# ({x4}, (6.610e-01, 3.390e-01))
-# ({x5}, (2.089e-01, 7.911e-01))
-# ({x6}, (8.182e-01, 1.818e-01))
-# ({x7}, (2.325e-01, 7.675e-01))
-# ({x8}, (2.280e-01, 7.720e-01))
-# ({x9}, (2.148e-01, 7.852e-01))
-# ({x10}, (7.613e-01, 2.387e-01))
-# ({x11}, (1.228e-01, 8.772e-01))
-# ({x12}, (4.254e-01, 5.746e-01))
-# ({x13}, (5.383e-01, 4.617e-01))
-# ({x14}, (6.257e-01, 3.743e-01))
-# ({x15}, (1.370e-01, 8.630e-01))
+MR_CLAMPING_LINEAR                             5.992e-03       1.960e-03       N/A             1.000e-09       
+# ({x0}, (3.856e-01, 6.144e-01))
+# ({x1}, (5.573e-01, 4.427e-01))
+# ({x2}, (4.586e-01, 5.414e-01))
+# ({x3}, (5.474e-01, 4.526e-01))
+# ({x4}, (6.653e-01, 3.347e-01))
+# ({x5}, (2.101e-01, 7.899e-01))
+# ({x6}, (8.189e-01, 1.811e-01))
+# ({x7}, (2.321e-01, 7.679e-01))
+# ({x8}, (2.231e-01, 7.769e-01))
+# ({x9}, (2.089e-01, 7.911e-01))
+# ({x10}, (7.643e-01, 2.357e-01))
+# ({x11}, (1.224e-01, 8.776e-01))
+# ({x12}, (4.252e-01, 5.748e-01))
+# ({x13}, (5.388e-01, 4.612e-01))
+# ({x14}, (6.275e-01, 3.725e-01))
+# ({x15}, (1.366e-01, 8.634e-01))
 MR_EXACT_LINEAR                                5.617e-03       1.742e-03       N/A             1.000e-09       
 # ({x0}, (3.853e-01, 6.147e-01))
 # ({x1}, (5.560e-01, 4.440e-01))