Cleaned up BBP and improved documentation of include/dai/bbp.h
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 10 Nov 2009 15:51:03 +0000 (16:51 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 10 Nov 2009 15:51:03 +0000 (16:51 +0100)
TODO
doxygen.conf
include/dai/bbp.h
include/dai/cbp.h
include/dai/enum.h
include/dai/gibbs.h
include/dai/treeep.h
src/bbp.cpp
src/cbp.cpp
src/gibbs.cpp
tests/testbbp.cpp

diff --git a/TODO b/TODO
index b17456b..be66603 100644 (file)
--- a/TODO
+++ b/TODO
@@ -2,7 +2,6 @@ To do for the next release (0.2.3):
 
 Improve documentation:
 
-       bbp.h
        cbp.h
        evidence.h
        emalg.h
index 398cb83..4514775 100644 (file)
@@ -1265,7 +1265,7 @@ PREDEFINED             = DAI_WITH_BP \
 # The macro definition that is found in the sources will be used. 
 # Use the PREDEFINED tag if you want to use a different macro definition.
 
-EXPAND_AS_DEFINED      =
+EXPAND_AS_DEFINED      = DAI_ENUM
 
 # If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then 
 # doxygen's preprocessor will remove all function-like macros that are alone 
index 866193c..c361afe 100644 (file)
@@ -9,10 +9,7 @@
 
 
 /// \file
-/// \brief Defines class BBP [\ref EaG09]
-/// \todo Fit more closely into libDAI framework
-/// \todo Improve documentation
-/// \author Frederik Eaton
+/// \brief  Defines class BBP, which implements Back-Belief-Propagation
 
 
 #ifndef ___defined_libdai_bbp_h
 namespace dai {
 
 
-/// Computes the adjoint of the unnormed probability vector from the normalizer and the adjoint of the normalized probability vector @see eqn. (13) in [\ref EaG09]
-Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w );
+/// Enumeration of several cost functions that can be used with BBP
+/** \note This class is meant as a base class for BBPCostFunction, which provides additional functionality.
+ */
+DAI_ENUM(BBPCostFunctionBase,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
 
-/// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state
-std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters );
+
+/// Predefined cost functions that can be used with BBP
+class BBPCostFunction : public BBPCostFunctionBase {
+    public:
+        /// Returns whether this cost function depends on having a Gibbs state
+        bool needGibbsState() const;
+
+        /// Evaluates cost function in state \a stateP using the information in inference algorithm \a ia
+        Real evaluate( const InfAlg &ia, const std::vector<size_t> *stateP ) const;
+
+        /// Assignment operator
+        BBPCostFunction& operator=( const BBPCostFunctionBase &x ) {
+            if( this != &x ) {
+                (BBPCostFunctionBase)*this = x;
+            }
+            return *this;
+        }
+};
 
 
 /// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
 /** \author Frederik Eaton
  */
 class BBP {
-    protected:
-    /// \name Inputs
+    private:
+    /// \name Input variables
     //@{
+        /// Stores a BP_dual helper object
         BP_dual _bp_dual;
+        /// Pointer to the factor graph
         const FactorGraph *_fg;
+        /// Pointer to the approximate inference algorithm
         const InfAlg *_ia;
     //@}
 
-        /// Number of iterations done
-        size_t _iters;
-
-    /// \name Outputs
+    /// \name Output variables
     //@{
         /// Variable factor adjoints
         std::vector<Prob> _adj_psi_V;
@@ -70,23 +85,8 @@ class BBP {
         std::vector<Prob> _adj_b_F;
     //@}
 
-    /// \name Helper quantities computed from the BP messages
+    /// \name Internal state variables
     //@{
-        /// _T[i][_I] (see eqn. (41) in [\ref EaG09])
-        std::vector<std::vector<Prob > > _T;
-        /// _U[I][_i] (see eqn. (42) in [\ref EaG09])
-        std::vector<std::vector<Prob > > _U;
-        /// _S[i][_I][_j] (see eqn. (43) in [\ref EaG09])
-        std::vector<std::vector<std::vector<Prob > > > _S;
-        /// _R[I][_i][_J] (see eqn. (44) in [\ref EaG09])
-        std::vector<std::vector<std::vector<Prob > > > _R;
-    //@}
-
-        /// Unnormalized variable belief adjoints
-        std::vector<Prob> _adj_b_V_unnorm;
-        /// Unnormalized factor belief adjoints
-        std::vector<Prob> _adj_b_F_unnorm;
-
         /// Initial variable factor adjoints
         std::vector<Prob> _init_adj_psi_V;
         /// Initial factor adjoints
@@ -100,21 +100,39 @@ class BBP {
         std::vector<std::vector<Prob> > _new_adj_n;
         /// Updated normalized factor->variable message adjoint (indexed [i][_I])
         std::vector<std::vector<Prob> > _new_adj_m;
+        /// Unnormalized variable belief adjoints
+        std::vector<Prob> _adj_b_V_unnorm;
+        /// Unnormalized factor belief adjoints
+        std::vector<Prob> _adj_b_F_unnorm;
 
-    /// \name Optimized indexing (for performance)
-    //@{
-        /// Calculates _indices, which is a cache of IndexFor @see bp.cpp
-        void RegenerateInds();
+        /// _T[i][_I] (see eqn. (41) in [\ref EaG09])
+        std::vector<std::vector<Prob > > _T;
+        /// _U[I][_i] (see eqn. (42) in [\ref EaG09])
+        std::vector<std::vector<Prob > > _U;
+        /// _S[i][_I][_j] (see eqn. (43) in [\ref EaG09])
+        std::vector<std::vector<std::vector<Prob > > > _S;
+        /// _R[I][_i][_J] (see eqn. (44) in [\ref EaG09])
+        std::vector<std::vector<std::vector<Prob > > > _R;
+
+        /// Number of iterations done
+        size_t _iters;
+    //@}
 
+    /// \name Index cache management (for performance)
+    //@{
         /// Index type
         typedef std::vector<size_t>  _ind_t;
         /// Cached indices (indexed [i][_I])
         std::vector<std::vector<_ind_t> >  _indices;
+        /// Prepares index cache _indices
+        /** \see bp.cpp
+         */
+        void RegenerateInds();
         /// Returns an index from the cache
         const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
     //@}
 
-    /// \name Initialization
+    /// \name Initialization helper functions
     //@{
         /// Calculate T values; see eqn. (41) in [\ref EaG09]
         void RegenerateT();
@@ -126,17 +144,26 @@ class BBP {
         void RegenerateR();
         /// Calculate _adj_b_V_unnorm and _adj_b_F_unnorm from _adj_b_V and _adj_b_F
         void RegenerateInputs();
-        /// Initialise members for factor adjoints (call after RegenerateInputs)
+        /// Initialise members for factor adjoints
+        /** \pre RegenerateInputs() should be called first
+         */
         void RegeneratePsiAdjoints();
-        /// Initialise members for message adjoints (call after RegenerateInputs) for parallel algorithm
+        /// Initialise members for message adjoints for parallel algorithm
+        /** \pre RegenerateInputs() should be called first
+         */
         void RegenerateParMessageAdjoints();
-        /// Initialise members for message adjoints (call after RegenerateInputs) for sequential algorithm
+        /// Initialise members for message adjoints for sequential algorithm
         /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
          *  than updating _adj_n (and friends) which are unused in the sequential algorithm.
+         *  \pre RegenerateInputs() should be called first
          */
         void RegenerateSeqMessageAdjoints();
+        /// Called by \a init, recalculates intermediate values
+        void Regenerate();
     //@}
 
+    /// \name Accessors/mutators
+    //@{
         /// Returns reference to T value; see eqn. (41) in [\ref EaG09]
         Prob & T(size_t i, size_t _I) { return _T[i][_I]; }
         /// Returns constant reference to T value; see eqn. (41) in [\ref EaG09]
@@ -154,6 +181,16 @@ class BBP {
         /// Returns constant reference to R value; see eqn. (44) in [\ref EaG09]
         const Prob & R(size_t I, size_t _i, size_t _J) const { return _R[I][_i][_J]; }
 
+        /// Returns reference to variable->factor message adjoint
+        Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
+        /// Returns constant reference to variable->factor message adjoint
+        const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
+        /// Returns reference to factor->variable message adjoint
+        Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
+        /// Returns constant reference to factor->variable message adjoint
+        const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
+    //@}
+
     /// \name Parallel algorithm
     //@{
         /// Calculates new variable->factor message adjoint
@@ -180,7 +217,7 @@ class BBP {
 
     /// \name Sequential algorithm
     //@{
-        /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
+        /// Helper function for sendSeqMsgM(): increases factor->variable message adjoint by \a p and calculates the corresponding unnormalized adjoint
         void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
         //  DISABLED BECAUSE IT IS BUGGY:
         //  void updateSeqMsgM( size_t i, size_t _I );
@@ -192,23 +229,20 @@ class BBP {
         void sendSeqMsgM( size_t i, size_t _I );
     //@}
 
-        /// Calculates averaged L-1 norm of unnormalized message adjoints
+        /// Computes the adjoint of the unnormed probability vector from the normalizer and the adjoint of the normalized probability vector
+        /** \see eqn. (13) in [\ref EaG09]
+         */
+        Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w );
+
+        /// Calculates averaged L1 norm of unnormalized message adjoints
         Real getUnMsgMag();
-        /// Calculates averaged L-1 norms of current and new normalized message adjoints
+        /// Calculates averaged L1 norms of current and new normalized message adjoints
         void getMsgMags( Real &s, Real &new_s );
-
-        /// Sets all vectors _adj_b_F to zero
-        void zero_adj_b_F() {
-            _adj_b_F.clear();
-            _adj_b_F.reserve( _fg->nrFactors() );
-            for( size_t I = 0; I < _fg->nrFactors(); I++ )
-                _adj_b_F.push_back( Prob( _fg->factor(I).states(), Real( 0.0 ) ) );
-        }
-
         /// Returns indices and magnitude of the largest normalized factor->variable message adjoint
         void getArgmaxMsgM( size_t &i, size_t &_I, Real &mag );
         /// Returns magnitude of the largest (in L1-norm) normalized factor->variable message adjoint
         Real getMaxMsgM();
+
         /// Calculates sum of L1 norms of all normalized factor->variable message adjoints
         Real getTotalMsgM();
         /// Calculates sum of L1 norms of all updated normalized factor->variable message adjoints
@@ -216,21 +250,23 @@ class BBP {
         /// Calculates sum of L1 norms of all normalized variable->factor message adjoints
         Real getTotalMsgN();
 
-    public:
-        /// Called by \a init, recalculates intermediate values
-        void Regenerate();
+        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph \a fg
+        std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
+        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph \a fg
+        std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
 
-        /// Constructor
+    public:
+    /// \name Constructors/destructors
+    //@{
+        /// Construct from a InfAlg \a ia and a PropertySet \a opts
         BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
             props.set(opts);
         }
+    //@}
 
-        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg
-        std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
-        /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph fg
-        std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
-
-        /// Initializes belief adjoints and initial factor adjoints and regenerates
+    /// \name Initialization
+    //@{
+        /// Initializes from given belief adjoints \a adj_b_V, \a adj_b_F and initial factor adjoints \a adj_psi_V, \a adj_psi_F
         void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F, const std::vector<Prob> &adj_psi_V, const std::vector<Prob> &adj_psi_F ) {
             _adj_b_V = adj_b_V;
             _adj_b_F = adj_b_F;
@@ -239,22 +275,37 @@ class BBP {
             Regenerate();
         }
 
-        /// Initializes belief adjoints and with zero initial factor adjoints and regenerates
+        /// Initializes from given belief adjoints \a adj_b_V and \a adj_b_F (setting initial factor adjoints to zero)
         void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
             init( adj_b_V, adj_b_F, getZeroAdjV(*_fg), getZeroAdjF(*_fg) );
         }
 
-        /// Initializes variable belief adjoints (and sets factor belief adjoints to zero) and with zero initial factor adjoints and regenerates
-        void init( const std::vector<Prob> &adj_b_V ) {
-            init(adj_b_V, getZeroAdjF(*_fg));
+        /// Initializes variable belief adjoints \a adj_b_V (and sets factor belief adjoints and initial factor adjoints to zero)
+        void init_V( const std::vector<Prob> &adj_b_V ) {
+            init( adj_b_V, getZeroAdjF(*_fg) );
         }
 
-        /// Run until change is less than given tolerance
-        void run();
+        /// Initializes factor belief adjoints \a adj_b_F (and sets variable belief adjoints and initial factor adjoints to zero)
+        void init_F( const std::vector<Prob> &adj_b_F ) {
+            init( getZeroAdjV(*_fg), adj_b_F );
+        }
 
-        /// Return number of iterations done so far
-        size_t doneIters() { return _iters; }
+        /// Initializes with adjoints calculated from cost function \a cfn, and state \a stateP
+        /** Uses the internal pointer to an inference algorithm in combination with the cost function and state for initialization.
+         *  \param cfn Cost function used for initialization;
+         *  \param stateP is a Gibbs state and can be NULL; it will be initialised using a Gibbs run.
+         */
+        void initCostFnAdj( const BBPCostFunction &cfn, const std::vector<size_t> *stateP );
+    //@}
 
+    /// \name BBP Algorithm
+    //@{
+        /// Perform iterative updates until change is less than given tolerance
+        void run();
+    //@}
+
+    /// \name Query results
+    //@{
         /// Returns reference to variable factor adjoint
         Prob& adj_psi_V(size_t i) { return _adj_psi_V[i]; }
         /// Returns constant reference to variable factor adjoint
@@ -271,18 +322,11 @@ class BBP {
         Prob& adj_b_F(size_t I) { return _adj_b_F[I]; }
         /// Returns constant reference to factor belief adjoint
         const Prob& adj_b_F(size_t I) const { return _adj_b_F[I]; }
+        /// Return number of iterations done so far
+        size_t Iterations() { return _iters; }
+    //@}
 
-     protected:
-        /// Returns reference to variable->factor message adjoint
-        Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
-        /// Returns constant reference to variable->factor message adjoint
-        const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
-        /// Returns reference to factor->variable message adjoint
-        Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
-        /// Returns constant reference to factor->variable message adjoint
-        const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
-
-     public:
+    public:
         /// Parameters of this algorithm
         /* PROPERTIES(props,BBP) {
            /// Enumeration of possible update schedules
@@ -335,35 +379,16 @@ class BBP {
 };
 
 
-/// Enumeration of several cost functions that can be used with BBP.
-DAI_ENUM(bbp_cfn_t,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
-
-/// Initialise BBP using InfAlg, cost function, and stateP
-/** Calls bbp.init with adjoints calculated from ia.beliefV and
- *  ia.beliefF. stateP is a Gibbs state and can be NULL, it will be
- *  initialised using a Gibbs run of 2*fg.Iterations() iterations.
- */
-void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
-
-/// Answers question: does the given cost function depend on having a Gibbs state?
-bool needGibbsState( bbp_cfn_t cfn );
-
-/// Calculate actual value of cost function (cfn_type, stateP)
-/** This function returns the actual value of the cost function whose
- *  gradient with respect to singleton beliefs is given by
- *  gibbsToB1Adj on the same arguments
- */
-Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
-
-/// Function to test the validity of adjoints computed by BBP given a state for each variable using numerical derivatives.
-/** Factors containing a variable are multiplied by psi_1 adjustments to verify accuracy of _adj_psi_V.
- *  \param bp BP object.
- *  \param state Global state of all variables.
- *  \param bbp_props BBP Properties.
- *  \param cfn Cost function to be used.
- *  \param h controls size of perturbation.
+/// Function to verify the validity of adjoints computed by BBP using numerical differentiation.
+/** Factors containing a variable are multiplied by small adjustments to verify accuracy of calculated variable factor adjoints.
+ *  \param bp BP object;
+ *  \param state Global state of all variables;
+ *  \param bbp_props BBP parameters;
+ *  \param cfn Cost function to be used;
+ *  \param h Size of perturbation.
+ *  \relates BBP
  */
-Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h );
+Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h );
 
 
 } // end of namespace dai
index 253baf8..370228d 100644 (file)
@@ -30,7 +30,7 @@ namespace dai {
 
 /// Find a variable to clamp using BBP (goes with maximum adjoint)
 /// \see BBP
-std::pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut );
+std::pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real *maxVarOut );
 
 
 /// Class for CBP (Clamped Belief Propagation)
@@ -63,7 +63,7 @@ class CBP : public DAIAlgFG {
         boost::shared_ptr<std::ofstream> _clamp_ofstream;
 
         /// Returns BBP cost function used
-        bbp_cfn_t BBP_cost_function() { return props.bbp_cfn; }
+        BBPCostFunction BBP_cost_function() { return props.bbp_cfn; }
 
         /// Prints beliefs, variables and partition sum, in case of a debugging build
         void printDebugInfo();
@@ -171,7 +171,7 @@ class CBP : public DAIAlgFG {
             /// Properties to pass to BBP
             PropertySet bbp_props;
             /// Cost function to use for BBP
-            bbp_cfn_t bbp_cfn;
+            BBPCostFunction bbp_cfn;
             /// Random seed
             size_t rand_seed = 0;
 
@@ -214,7 +214,7 @@ class CBP : public DAIAlgFG {
             /// Properties to pass to BBP
             PropertySet bbp_props;
             /// Cost function to use for BBP
-            bbp_cfn_t bbp_cfn;
+            BBPCostFunction bbp_cfn;
             /// Random seed
             size_t rand_seed;
             /// If non-empty, write clamping choices to this file
index 07c5785..02dfb3f 100644 (file)
@@ -95,7 +95,7 @@
             return os;\
         }\
 \
-    private:\
+    protected:\
         value v;\
 };
 
index c4166a6..451941e 100644 (file)
@@ -116,6 +116,12 @@ class Gibbs : public DAIAlgFG {
 };
 
 
+/// Runs Gibbs sampling for \a iters iterations (of which \a burnin for burn-in) on FactorGraph \a fg, and returns the resulting state
+/** \relates Gibbs
+ */
+std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t iters, size_t burnin=0 );
+
+
 } // end of namespace dai
 
 
index 5cc269e..07a7897 100644 (file)
@@ -11,6 +11,7 @@
 
 /// \file
 /// \brief Defines class TreeEP, which implements Tree Expectation Propagation
+/// \todo Clean up the TreeEP code
 
 
 #ifndef __defined_libdai_treeep_h
index 7720881..1055c35 100644 (file)
@@ -25,29 +25,6 @@ using namespace std;
 typedef BipartiteGraph::Neighbor Neighbor;
 
 
-Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
-    DAI_ASSERT( w.size() == adj_w.size() );
-    Prob adj_w_unnorm( w.size(), 0.0 );
-    Real s = 0.0;
-    for( size_t i = 0; i < w.size(); i++ )
-        s += w[i] * adj_w[i];
-    for( size_t i = 0; i < w.size(); i++ )
-        adj_w_unnorm[i] = (adj_w[i] - s) / Z_w;
-    return adj_w_unnorm;
-//  THIS WOULD BE ABOUT 50% SLOWER:  return (adj_w - (w * adj_w).sum()) / Z_w;
-}
-
-
-std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters ) {
-    PropertySet gibbsProps;
-    gibbsProps.Set("iters", iters);
-    gibbsProps.Set("verbose", size_t(0));
-    Gibbs gibbs( ia.fg(), gibbsProps );
-    gibbs.run();
-    return gibbs.state();
-}
-
-
 /// Returns the entry of the I'th factor corresponding to a global state
 size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
     size_t f_entry = 0;
@@ -62,6 +39,91 @@ size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<siz
 }
 
 
+bool BBPCostFunction::needGibbsState() const {
+    switch( (size_t)(*this) ) {
+        case CFN_GIBBS_B:
+        case CFN_GIBBS_B2:
+        case CFN_GIBBS_EXP:
+        case CFN_GIBBS_B_FACTOR:
+        case CFN_GIBBS_B2_FACTOR:
+        case CFN_GIBBS_EXP_FACTOR:
+            return true;
+        default:
+            return false;
+    }
+}
+
+
+Real BBPCostFunction::evaluate( const InfAlg &ia, const vector<size_t> *stateP ) const {
+    Real cf = 0.0;
+    const FactorGraph &fg = ia.fg();
+
+    switch( (size_t)(*this) ) {
+        case CFN_BETHE_ENT: // ignores state
+            cf = -ia.logZ();
+            break;
+        case CFN_VAR_ENT: // ignores state
+            for( size_t i = 0; i < fg.nrVars(); i++ )
+                cf += -ia.beliefV(i).entropy();
+            break;
+        case CFN_FACTOR_ENT: // ignores state
+            for( size_t I = 0; I < fg.nrFactors(); I++ )
+                cf += -ia.beliefF(I).entropy();
+            break;
+        case CFN_GIBBS_B:
+        case CFN_GIBBS_B2:
+        case CFN_GIBBS_EXP: {
+            DAI_ASSERT( stateP != NULL );
+            vector<size_t> state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                Real b = ia.beliefV(i)[state[i]];
+                switch( (size_t)(*this) ) {
+                    case CFN_GIBBS_B:
+                        cf += b;
+                        break;
+                    case CFN_GIBBS_B2:
+                        cf += b * b / 2.0;
+                        break;
+                    case CFN_GIBBS_EXP:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+            }
+            break;
+        } case CFN_GIBBS_B_FACTOR:
+          case CFN_GIBBS_B2_FACTOR:
+          case CFN_GIBBS_EXP_FACTOR: {
+            DAI_ASSERT( stateP != NULL );
+            vector<size_t> state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                Real b = ia.beliefF(I)[x_I];
+                switch( (size_t)(*this) ) {
+                    case CFN_GIBBS_B_FACTOR:
+                        cf += b;
+                        break;
+                    case CFN_GIBBS_B2_FACTOR:
+                        cf += b * b / 2.0;
+                        break;
+                    case CFN_GIBBS_EXP_FACTOR:
+                        cf += exp( b );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+            }
+            break;
+        } default:
+            DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(*this));
+    }
+    return cf;
+}
+
+
 #define LOOP_ij(body) {                             \
     size_t i_states = _fg->var(i).states();         \
     size_t j_states = _fg->var(j).states();         \
@@ -326,6 +388,22 @@ void BBP::RegenerateSeqMessageAdjoints() {
 }
 
 
+void BBP::Regenerate() {
+    RegenerateInds();
+    RegenerateT();
+    RegenerateU();
+    RegenerateS();
+    RegenerateR();
+    RegenerateInputs();
+    RegeneratePsiAdjoints();
+    if( props.updates == Properties::UpdateType::PAR )
+        RegenerateParMessageAdjoints();
+    else
+        RegenerateSeqMessageAdjoints();
+    _iters = 0;
+}
+
+
 void BBP::calcNewN( size_t i, size_t _I ) {
     _adj_psi_V[i] += T(i,_I) * _adj_n_unnorm[i][_I];
     Prob &new_adj_n_iI = _new_adj_n[i][_I];
@@ -544,6 +622,19 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
 }
 
 
+Prob BBP::unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
+    DAI_ASSERT( w.size() == adj_w.size() );
+    Prob adj_w_unnorm( w.size(), 0.0 );
+    Real s = 0.0;
+    for( size_t i = 0; i < w.size(); i++ )
+        s += w[i] * adj_w[i];
+    for( size_t i = 0; i < w.size(); i++ )
+        adj_w_unnorm[i] = (adj_w[i] - s) / Z_w;
+    return adj_w_unnorm;
+//  THIS WOULD BE ABOUT 50% SLOWER:  return (adj_w - (w * adj_w).sum()) / Z_w;
+}
+
+
 Real BBP::getUnMsgMag() {
     Real s = 0.0;
     size_t e = 0;
@@ -642,22 +733,6 @@ Real BBP::getTotalMsgN() {
 }
 
 
-void BBP::Regenerate() {
-    RegenerateInds();
-    RegenerateT();
-    RegenerateU();
-    RegenerateS();
-    RegenerateR();
-    RegenerateInputs();
-    RegeneratePsiAdjoints();
-    if( props.updates == Properties::UpdateType::PAR )
-        RegenerateParMessageAdjoints();
-    else
-        RegenerateSeqMessageAdjoints();
-    _iters = 0;
-}
-
-
 std::vector<Prob> BBP::getZeroAdjF( const FactorGraph &fg ) {
     vector<Prob> adj_2;
     adj_2.reserve( fg.nrFactors() );
@@ -676,6 +751,157 @@ std::vector<Prob> BBP::getZeroAdjV( const FactorGraph &fg ) {
 }
 
 
+void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *stateP ) {
+    const FactorGraph &fg = _ia->fg();
+
+    switch( (size_t)cfn ) {
+        case BBPCostFunction::CFN_BETHE_ENT: {
+            vector<Prob> b1_adj;
+            vector<Prob> b2_adj;
+            vector<Prob> psi1_adj;
+            vector<Prob> psi2_adj;
+            b1_adj.reserve( fg.nrVars() );
+            psi1_adj.reserve( fg.nrVars() );
+            b2_adj.reserve( fg.nrFactors() );
+            psi2_adj.reserve( fg.nrFactors() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                int c = fg.nbV(i).size();
+                Prob p(dim,0.0);
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p[xi] = (1 - c) * (1 + log( _ia->beliefV(i)[xi] ));
+                b1_adj.push_back( p );
+
+                for( size_t xi = 0; xi < dim; xi++ )
+                    p[xi] = -_ia->beliefV(i)[xi];
+                psi1_adj.push_back( p );
+            }
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t dim = fg.factor(I).states();
+                Prob p( dim, 0.0 );
+                for( size_t xI = 0; xI < dim; xI++ )
+                    p[xI] = 1 + log( _ia->beliefF(I)[xI] / fg.factor(I).p()[xI] );
+                b2_adj.push_back( p );
+
+                for( size_t xI = 0; xI < dim; xI++ )
+                    p[xI] = -_ia->beliefF(I)[xI] / fg.factor(I).p()[xI];
+                psi2_adj.push_back( p );
+            }
+            init( b1_adj, b2_adj, psi1_adj, psi2_adj );
+            break;
+        } case BBPCostFunction::CFN_FACTOR_ENT: {
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrFactors() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t dim = fg.factor(I).states();
+                Prob p( dim, 0.0 );
+                for( size_t xI = 0; xI < dim; xI++ ) {
+                    Real bIxI = _ia->beliefF(I)[xI];
+                    if( bIxI < 1.0e-15 )
+                        p[xI] = -1.0e10;
+                    else
+                        p[xI] = 1 + log( bIxI );
+                }
+                b2_adj.push_back(p);
+            }
+            init_F( b2_adj );
+            break;
+        } case BBPCostFunction::CFN_VAR_ENT: {
+            vector<Prob> b1_adj;
+            b1_adj.reserve( fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                size_t dim = fg.var(i).states();
+                Prob p( dim, 0.0 );
+                for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
+                    Real bixi = _ia->beliefV(i)[xi];
+                    if( bixi < 1.0e-15 )
+                        p[xi] = -1.0e10;
+                    else
+                        p[xi] = 1 + log( bixi );
+                }
+                b1_adj.push_back( p );
+            }
+            init_V( b1_adj );
+            break;
+        } case BBPCostFunction::CFN_GIBBS_B:
+          case BBPCostFunction::CFN_GIBBS_B2:
+          case BBPCostFunction::CFN_GIBBS_EXP: {
+            // cost functions that use Gibbs sample, summing over variable marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+            else
+                state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+
+            vector<Prob> b1_adj;
+            b1_adj.reserve(fg.nrVars());
+            for( size_t i = 0; i < state.size(); i++ ) {
+                size_t n = fg.var(i).states();
+                Prob delta( n, 0.0 );
+                DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
+                Real b = _ia->beliefV(i)[state[i]];
+                switch( (size_t)cfn ) {
+                    case BBPCostFunction::CFN_GIBBS_B:
+                        delta[state[i]] = 1.0;
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_B2:
+                        delta[state[i]] = b;
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_EXP:
+                        delta[state[i]] = exp(b);
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+                b1_adj.push_back( delta );
+            }
+            init_V( b1_adj );
+            break;
+        } case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+          case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+          case BBPCostFunction::CFN_GIBBS_EXP_FACTOR: {
+            // cost functions that use Gibbs sample, summing over factor marginals
+            vector<size_t> state;
+            if( stateP == NULL )
+                state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+            else
+                state = *stateP;
+            DAI_ASSERT( state.size() == fg.nrVars() );
+
+            vector<Prob> b2_adj;
+            b2_adj.reserve( fg.nrVars() );
+            for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
+                size_t n = fg.factor(I).states();
+                Prob delta( n, 0.0 );
+
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
+
+                Real b = _ia->beliefF(I)[x_I];
+                switch( (size_t)cfn ) {
+                    case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+                        delta[x_I] = 1.0;
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+                        delta[x_I] = b;
+                        break;
+                    case BBPCostFunction::CFN_GIBBS_EXP_FACTOR:
+                        delta[x_I] = exp( b );
+                        break;
+                    default:
+                        DAI_THROW(UNKNOWN_ENUM_VALUE);
+                }
+                b2_adj.push_back( delta );
+            }
+            init_F( b2_adj );
+            break;
+        } default:
+            DAI_THROW(UNKNOWN_ENUM_VALUE);
+    }
+}
+
+
 void BBP::run() {
     typedef BBP::Properties::UpdateType UT;
     Real tol = props.tol;
@@ -756,17 +982,16 @@ void BBP::run() {
         }
     }
     if( props.verbose >= 3 )
-        cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl;
+        cerr << "BBP::run() took " << toc()-tic << " seconds " << Iterations() << " iterations" << endl;
 }
 
 
-Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h ) {
+Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h ) {
+    BBP bbp( &bp, bbp_props );
     // calculate the value of the unperturbed cost function
-    Real cf0 = getCostFn( bp, cfn, state );
-
+    Real cf0 = cfn.evaluate( bp, state );
     // run BBP to estimate adjoints
-    BBP bbp( &bp, bbp_props );
-    initBBPCostFnAdj( bbp, bp, cfn, state );
+    bbp.initCostFnAdj( cfn, state );
     bbp.run();
 
     Real d = 0;
@@ -798,7 +1023,7 @@ Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Proper
                 bp_prb->run();
 
                 // calculate new value of cost function
-                Real cf_prb = getCostFn( *bp_prb, cfn, state );
+                Real cf_prb = cfn.evaluate( *bp_prb, state );
 
                 // use to estimate adjoint for i
                 adj_est.push_back( (cf_prb - cf0) / h );
@@ -905,242 +1130,6 @@ Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const Proper
 }
 
 
-bool needGibbsState( bbp_cfn_t cfn ) {
-    switch( (size_t)cfn ) {
-        case bbp_cfn_t::CFN_GIBBS_B:
-        case bbp_cfn_t::CFN_GIBBS_B2:
-        case bbp_cfn_t::CFN_GIBBS_EXP:
-        case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-        case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-        case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-            return true;
-        default:
-            return false;
-    }
-}
-
-
-void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
-    const FactorGraph &fg = ia.fg();
-
-    switch( (size_t)cfn_type ) {
-        case bbp_cfn_t::CFN_BETHE_ENT: {
-            vector<Prob> b1_adj;
-            vector<Prob> b2_adj;
-            vector<Prob> psi1_adj;
-            vector<Prob> psi2_adj;
-            b1_adj.reserve( fg.nrVars() );
-            psi1_adj.reserve( fg.nrVars() );
-            b2_adj.reserve( fg.nrFactors() );
-            psi2_adj.reserve( fg.nrFactors() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                size_t dim = fg.var(i).states();
-                int c = fg.nbV(i).size();
-                Prob p(dim,0.0);
-                for( size_t xi = 0; xi < dim; xi++ )
-                    p[xi] = (1 - c) * (1 + log( ia.beliefV(i)[xi] ));
-                b1_adj.push_back( p );
-
-                for( size_t xi = 0; xi < dim; xi++ )
-                    p[xi] = -ia.beliefV(i)[xi];
-                psi1_adj.push_back( p );
-            }
-            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t dim = fg.factor(I).states();
-                Prob p( dim, 0.0 );
-                for( size_t xI = 0; xI < dim; xI++ )
-                    p[xI] = 1 + log( ia.beliefF(I)[xI] / fg.factor(I).p()[xI] );
-                b2_adj.push_back( p );
-
-                for( size_t xI = 0; xI < dim; xI++ )
-                    p[xI] = -ia.beliefF(I)[xI] / fg.factor(I).p()[xI];
-                psi2_adj.push_back( p );
-            }
-            bbp.init( b1_adj, b2_adj, psi1_adj, psi2_adj );
-            break;
-        } case bbp_cfn_t::CFN_FACTOR_ENT: {
-            vector<Prob> b2_adj;
-            b2_adj.reserve( fg.nrFactors() );
-            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t dim = fg.factor(I).states();
-                Prob p( dim, 0.0 );
-                for( size_t xI = 0; xI < dim; xI++ ) {
-                    Real bIxI = ia.beliefF(I)[xI];
-                    if( bIxI < 1.0e-15 )
-                        p[xI] = -1.0e10;
-                    else
-                        p[xI] = 1 + log( bIxI );
-                }
-                b2_adj.push_back(p);
-            }
-            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
-            break;
-        } case bbp_cfn_t::CFN_VAR_ENT: {
-            vector<Prob> b1_adj;
-            b1_adj.reserve( fg.nrVars() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                size_t dim = fg.var(i).states();
-                Prob p( dim, 0.0 );
-                for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
-                    Real bixi = ia.beliefV(i)[xi];
-                    if( bixi < 1.0e-15 )
-                        p[xi] = -1.0e10;
-                    else
-                        p[xi] = 1 + log( bixi );
-                }
-                b1_adj.push_back( p );
-            }
-            bbp.init( b1_adj );
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B:
-          case bbp_cfn_t::CFN_GIBBS_B2:
-          case bbp_cfn_t::CFN_GIBBS_EXP: {
-            // cost functions that use Gibbs sample, summing over variable marginals
-            vector<size_t> state;
-            if( stateP == NULL )
-                state = getGibbsState( ia, 2*ia.Iterations() );
-            else
-                state = *stateP;
-            DAI_ASSERT( state.size() == fg.nrVars() );
-
-            vector<Prob> b1_adj;
-            b1_adj.reserve(fg.nrVars());
-            for( size_t i = 0; i < state.size(); i++ ) {
-                size_t n = fg.var(i).states();
-                Prob delta( n, 0.0 );
-                DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
-                Real b = ia.beliefV(i)[state[i]];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B:
-                        delta[state[i]] = 1.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2:
-                        delta[state[i]] = b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP:
-                        delta[state[i]] = exp(b);
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-                b1_adj.push_back( delta );
-            }
-            bbp.init( b1_adj );
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
-            // cost functions that use Gibbs sample, summing over factor marginals
-            vector<size_t> state;
-            if( stateP == NULL )
-                state = getGibbsState( ia, 2*ia.Iterations() );
-            else
-                state = *stateP;
-            DAI_ASSERT( state.size() == fg.nrVars() );
-
-            vector<Prob> b2_adj;
-            b2_adj.reserve( fg.nrVars() );
-            for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
-                size_t n = fg.factor(I).states();
-                Prob delta( n, 0.0 );
-
-                size_t x_I = getFactorEntryForState( fg, I, state );
-                DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
-
-                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;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-                        delta[x_I] = b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-                        delta[x_I] = exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-                b2_adj.push_back( delta );
-            }
-            bbp.init( bbp.getZeroAdjV(fg), b2_adj );
-            break;
-        } default:
-            DAI_THROW(UNKNOWN_ENUM_VALUE);
-    }
-}
-
-
-Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
-    Real cf = 0.0;
-    const FactorGraph &fg = ia.fg();
-
-    switch( (size_t)cfn_type ) {
-        case bbp_cfn_t::CFN_BETHE_ENT: // ignores state
-            cf = -ia.logZ();
-            break;
-        case bbp_cfn_t::CFN_VAR_ENT: // ignores state
-            for( size_t i = 0; i < fg.nrVars(); i++ )
-                cf += -ia.beliefV(i).entropy();
-            break;
-        case bbp_cfn_t::CFN_FACTOR_ENT: // ignores state
-            for( size_t I = 0; I < fg.nrFactors(); I++ )
-                cf += -ia.beliefF(I).entropy();
-            break;
-        case bbp_cfn_t::CFN_GIBBS_B:
-        case bbp_cfn_t::CFN_GIBBS_B2:
-        case bbp_cfn_t::CFN_GIBBS_EXP: {
-            DAI_ASSERT( stateP != NULL );
-            vector<size_t> state = *stateP;
-            DAI_ASSERT( state.size() == fg.nrVars() );
-            for( size_t i = 0; i < fg.nrVars(); i++ ) {
-                Real b = ia.beliefV(i)[state[i]];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B:
-                        cf += b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2:
-                        cf += b * b / 2.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP:
-                        cf += exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-            }
-            break;
-        } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-          case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
-            DAI_ASSERT( stateP != NULL );
-            vector<size_t> state = *stateP;
-            DAI_ASSERT( state.size() == fg.nrVars() );
-            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-                size_t x_I = getFactorEntryForState( fg, I, state );
-                Real b = ia.beliefF(I)[x_I];
-                switch( (size_t)cfn_type ) {
-                    case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
-                        cf += b;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
-                        cf += b * b / 2.0;
-                        break;
-                    case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
-                        cf += exp( b );
-                        break;
-                    default:
-                        DAI_THROW(UNKNOWN_ENUM_VALUE);
-                }
-            }
-            break;
-        } default:
-            DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(cfn_type));
-    }
-    return cf;
-}
-
-
 } // end of namespace dai
 
 
index b33bb9e..d77a90f 100644 (file)
@@ -16,7 +16,7 @@
 
 #include <dai/util.h>
 #include <dai/properties.h>
-
+#include <dai/gibbs.h>
 #include <dai/bp.h>
 #include <dai/cbp.h>
 #include <dai/bbp.h>
@@ -328,8 +328,8 @@ bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, siz
                ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_CFN ) {
         bool doL1 = (ChooseMethod() == Properties::ChooseMethodType::CHOOSE_BP_L1);
         vector<size_t> state;
-        if( !doL1 && needGibbsState( BBP_cost_function() ) )
-            state = getGibbsState( *bp, 2*bp->Iterations() );
+        if( !doL1 && BBP_cost_function().needGibbsState() )
+            state = getGibbsState( bp->fg(), 2*bp->Iterations() );
         // try clamping each variable manually
         DAI_ASSERT( Clamping() == Properties::ClampType::CLAMP_VAR );
         Real max_cost = 0.0;
@@ -347,7 +347,7 @@ bool CBP::chooseNextClampVar( InfAlg *bp, vector<size_t> &clamped_vars_list, siz
                     for( size_t j = 0; j < nrVars(); j++ )
                         cost += dist( bp->beliefV(j), bp1->beliefV(j), Prob::DISTL1 );
                 else
-                    cost = getCostFn( *bp1, BBP_cost_function(), &state );
+                    cost = BBP_cost_function().evaluate( *bp1, &state );
                 if( cost > max_cost || win_k == -1 ) {
                     max_cost = cost;
                     win_k = k;
@@ -415,9 +415,9 @@ void CBP::printDebugInfo() {
  *  creates and runs a BBP object, finds best variable, returns
  *  (variable,state) pair for clamping
  */
-pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut ) {
+pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real *maxVarOut ) {
     BBP bbp( &in_bp, bbp_props );
-    initBBPCostFnAdj( bbp, in_bp, cfn, NULL );
+    bbp.initCostFnAdj( cfn, NULL );
     bbp.run();
 
     // find and return the (variable,state) with the largest adj_psi_V
@@ -581,7 +581,7 @@ void CBP::Properties::set(const PropertySet &opts)
     recursion = opts.getStringAs<RecurseType>("recursion");
     clamp = opts.getStringAs<ClampType>("clamp");
     bbp_props = opts.getStringAs<PropertySet>("bbp_props");
-    bbp_cfn = opts.getStringAs<bbp_cfn_t>("bbp_cfn");
+    bbp_cfn = opts.getStringAs<BBPCostFunction>("bbp_cfn");
     if(opts.hasKey("rand_seed")) {
         rand_seed = opts.getStringAs<size_t>("rand_seed");
     } else {
index 5655448..e4bf112 100644 (file)
@@ -234,4 +234,15 @@ Factor Gibbs::belief( const VarSet &ns ) const {
 }
 
 
+std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t iters, size_t burnin ) {
+    PropertySet gibbsProps;
+    gibbsProps.Set("iters", iters);
+    gibbsProps.Set("burnin", burnin);
+    gibbsProps.Set("verbose", size_t(0));
+    Gibbs gibbs( fg, gibbsProps );
+    gibbs.run();
+    return gibbs.state();
+}
+
+
 } // end of namespace dai
index 45cfbdd..91952dc 100644 (file)
@@ -69,34 +69,34 @@ int main( int argc, char *argv[] ) {
                     updates = BBP::Properties::UpdateType::PAR;
                     break;
             }
-            bbp_cfn_t cfn;
+            BBPCostFunction cfn;
             switch( (t / 5) % 9 ) {
                 case 0:
-                    cfn = bbp_cfn_t::CFN_GIBBS_B;
+                    cfn = BBPCostFunction::CFN_GIBBS_B;
                     break;
                 case 1:
-                    cfn = bbp_cfn_t::CFN_GIBBS_B2;
+                    cfn = BBPCostFunction::CFN_GIBBS_B2;
                     break;
                 case 2:
-                    cfn = bbp_cfn_t::CFN_GIBBS_EXP;
+                    cfn = BBPCostFunction::CFN_GIBBS_EXP;
                     break;
                 case 3:
-                    cfn = bbp_cfn_t::CFN_GIBBS_B_FACTOR;
+                    cfn = BBPCostFunction::CFN_GIBBS_B_FACTOR;
                     break;
                 case 4:
-                    cfn = bbp_cfn_t::CFN_GIBBS_B2_FACTOR;
+                    cfn = BBPCostFunction::CFN_GIBBS_B2_FACTOR;
                     break;
                 case 5:
-                    cfn = bbp_cfn_t::CFN_GIBBS_EXP_FACTOR;
+                    cfn = BBPCostFunction::CFN_GIBBS_EXP_FACTOR;
                     break;
                 case 6:
-                    cfn = bbp_cfn_t::CFN_VAR_ENT;
+                    cfn = BBPCostFunction::CFN_VAR_ENT;
                     break;
                 case 7:
-                    cfn = bbp_cfn_t::CFN_FACTOR_ENT;
+                    cfn = BBPCostFunction::CFN_FACTOR_ENT;
                     break;
                 case 8:
-                    cfn = bbp_cfn_t::CFN_BETHE_ENT;
+                    cfn = BBPCostFunction::CFN_BETHE_ENT;
                     break;
             }