+* 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 )
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 );
Improve documentation:
- daialg.h
- alldai.h
exactinf.h
mf.h
bp.h
/// \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
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.
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
/// \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
/// 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;
/// 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
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
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 ); }
+ //@}
};
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
}
-/// \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;
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 {
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() );
}
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];
}
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();
# ({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))
# ({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))
projects / libdai.git / commitdiff