Even more cleanup of BBP code
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 13:36:37 +0000 (15:36 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 4 Aug 2009 13:36:37 +0000 (15:36 +0200)
include/dai/bbp.h
src/bbp.cpp
tests/aliases.conf
tests/testbbp.cpp

index bf68121..2b41568 100644 (file)
@@ -21,7 +21,6 @@
 /// \file
 /// \brief Defines class BBP [\ref EaG09]
 /// \todo Improve documentation
-/// \todo Clean up
 /// \todo Debug clean_updates
 
 
 namespace dai {
 
 
-/// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg
-std::vector<Prob> get_zero_adj_F( const FactorGraph& fg );
+/// 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 );
 
-/// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph fg
-std::vector<Prob> get_zero_adj_V( const FactorGraph& fg );
+/// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state
+std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters );
 
 
 /// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
@@ -178,6 +177,19 @@ class BBP {
         void doParUpdate();
         //@}
 
+        /// @name Sequential algorithm
+        //@{
+        /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
+        void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
+        void updateSeqMsgM( size_t i, size_t _I );
+        /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint
+        void setSeqMsgM( size_t i, size_t _I, const Prob &p ); 
+        /// Implements routine Send-n in Figure 5 in [\ref EaG09]
+        void sendSeqMsgN( size_t i, size_t _I, const Prob &f );
+        /// Implements routine Send-m in Figure 5 in [\ref EaG09]
+        void sendSeqMsgM( size_t i, size_t _I );
+        //@}
+
         /// Calculates averaged L-1 norm of unnormalized message adjoints
         Real getUnMsgMag();
         /// Calculates averaged L-1 norms of current and new normalized message adjoints
@@ -191,22 +203,8 @@ class BBP {
                 _adj_b_F.push_back( Prob( _fg->factor(I).states(), Real( 0.0 ) ) );
         }
 
-        /// @name Sequential algorithm
-        //@{
-        /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
-        void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
-        void updateSeqMsgM( size_t i, size_t _I );
-        /// Implements routine Send-n in Figure 5 in [\ref EaG09]
-        void sendSeqMsgN( size_t i, size_t _I, const Prob &f );
-        /// Implements routine Send-m in Figure 5 in [\ref EaG09]
-        void sendSeqMsgM( size_t i, size_t _I );
-        /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint
-        void setSeqMsgM( size_t i, size_t _I, const Prob &p ); 
-        //@}
-
         /// 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
@@ -225,6 +223,11 @@ class BBP {
             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
         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;
@@ -236,12 +239,12 @@ class BBP {
 
         /// Initializes belief adjoints and with zero initial factor adjoints and regenerates
         void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
-            init( adj_b_V, adj_b_F, get_zero_adj_V(*_fg), get_zero_adj_F(*_fg) );
+            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, get_zero_adj_F(*_fg));
+            init(adj_b_V, getZeroAdjF(*_fg));
         }
 
         /// Run until change is less than given tolerance
@@ -317,43 +320,32 @@ class BBP {
 /* }}} END OF GENERATED CODE */
 };
 
-/// Cost functions. Not used by BBP class, only used by following functions.
-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);
+
+/// 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);
+void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
 
-/// Answers question: Does given cost function depend on having a Gibbs state?
-bool needGibbsState(bbp_cfn_t cfn);
+/// 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, use numerical derivatives
- *  (multiplying a factor containing a variable by psi_1 adjustments)
- *  to verify accuracy of _adj_psi_V.
- *  'h' controls size of perturbation.
- *  'bbpTol' controls tolerance of BBP run.
- */
-double numericBBPTest(const InfAlg& bp, const std::vector<size_t> *state, const PropertySet& bbp_props, bbp_cfn_t cfn, double h);
-
-// ----------------------------------------------------------------
-// Utility functions, some of which are used elsewhere
+Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
 
-/// 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 );
-
-/// Runs Gibbs sampling for 'iters' iterations on ia.fg(), and returns state
-std::vector<size_t> getGibbsState(const InfAlg& ia, size_t iters);
+/// 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 h controls size of perturbation.
+ */
+double numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h );
 
 
 } // end of namespace dai
index 06750db..1e55c2e 100644 (file)
@@ -47,6 +47,16 @@ Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_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();
+}
+
+
 size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
     size_t f_entry = 0;
     for( int _j = fg.nbF(I).size() - 1; _j >= 0; _j-- ) {
@@ -60,269 +70,43 @@ size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<siz
 }
 
 
-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++ ) {
-                double 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(get_zero_adj_V(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++ ) {
-                double 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;
-        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);
-            assert(/*0<=state[i] &&*/ state[i] < n);
-            double 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(INTERNAL_ERROR);
-            }
-            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;
-        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);
-            assert(/*0<=x_I &&*/ x_I < n);
-
-            double 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(INTERNAL_ERROR);
-            }
-            b2_adj.push_back(delta);
-        }
-        bbp.init(get_zero_adj_V(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 ) {
-    double 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: {
-        assert(stateP != NULL);
-        vector<size_t> state = *stateP;
-        assert(state.size()==fg.nrVars());
-        for( size_t i = 0; i < fg.nrVars(); i++ ) {
-            double 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; break;
-            case bbp_cfn_t::cfn_gibbs_exp: cf += exp(b); break;
-            default: DAI_THROW(INTERNAL_ERROR);
-            }
-        }
-        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: {
-        assert(stateP != NULL);
-        vector<size_t> state = *stateP;
-        assert(state.size()==fg.nrVars());
-        for( size_t I = 0; I < fg.nrFactors(); I++ ) {
-            size_t x_I = getFactorEntryForState(fg, I, state);
-            double b = ia.beliefF(I)[x_I];
-            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; break;
-            case bbp_cfn_t::cfn_gibbs_exp_factor: cf += exp(b); break;
-            default: DAI_THROW(INTERNAL_ERROR);
-            }
-        }
-        break;
-    }
-    default: 
-        cerr << "Unknown cost function " << cfn_type << endl;
-        DAI_THROW(UNKNOWN_ENUM_VALUE);
-    }
-    return cf;
-}
-
-
-vector<Prob> get_zero_adj_F( const FactorGraph &fg ) {
-    vector<Prob> adj_2;
-    adj_2.reserve( fg.nrFactors() );
-    for( size_t I = 0; I < fg.nrFactors(); I++ )
-        adj_2.push_back( Prob( fg.factor(I).states(), Real(0.0) ) );
-    return adj_2;
-}
-
-
-vector<Prob> get_zero_adj_V( const FactorGraph &fg ) {
-    vector<Prob> adj_1;
-    adj_1.reserve( fg.nrVars() );
-    for( size_t i = 0; i < fg.nrVars(); i++ )
-        adj_1.push_back( Prob( fg.var(i).states(), Real(0.0) ) );
-    return adj_1;
+#define LOOP_ij(body) {                             \
+    size_t i_states = _fg->var(i).states();         \
+    size_t j_states = _fg->var(j).states();         \
+    if(_fg->var(i) > _fg->var(j)) {                 \
+        size_t xij=0;                               \
+        for(size_t xi=0; xi<i_states; xi++) {       \
+            for(size_t xj=0; xj<j_states; xj++) {   \
+                body;                               \
+                xij++;                              \
+            }                                       \
+        }                                           \
+    } else {                                        \
+        size_t xij=0;                               \
+        for(size_t xj=0; xj<j_states; xj++) {       \
+            for(size_t xi=0; xi<i_states; xi++) {   \
+                body;                               \
+                xij++;                              \
+            }                                       \
+        }                                           \
+    }                                               \
 }
 
 
-#define LOOP_ij(body) {                                 \
-        size_t i_states = _fg->var(i).states();         \
-        size_t j_states = _fg->var(j).states();         \
-        if(_fg->var(i) > _fg->var(j)) {                 \
-            size_t xij=0;                               \
-            for(size_t xi=0; xi<i_states; xi++) {       \
-                for(size_t xj=0; xj<j_states; xj++) {   \
-                    body;                               \
-                    xij++;                              \
-                }                                       \
-            }                                           \
-        } else {                                        \
-            size_t xij=0;                               \
-            for(size_t xj=0; xj<j_states; xj++) {       \
-                for(size_t xi=0; xi<i_states; xi++) {   \
-                    body;                               \
-                    xij++;                              \
-                }                                       \
-            }                                           \
-        }                                               \
-    }
-
-
 void BBP::RegenerateInds() {
     // initialise _indices
-    //     typedef std::vector<size_t>  _ind_t;
+    //     typedef std::vector<size_t>        _ind_t;
     //     std::vector<std::vector<_ind_t> >  _indices; 
-    _indices.resize(_fg->nrVars());
+    _indices.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _indices[i].reserve(_fg->nbV(i).size());
-        foreach(const Neighbor &I, _fg->nbV(i)) {
+        _indices[i].clear();
+        _indices[i].reserve( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
             _ind_t index;
-            index.reserve(_fg->factor(I).states());
-            for(IndexFor k( _fg->var(i), _fg->factor(I).vars() ); k >= 0; ++k) {
-                index.push_back(k);
-            }
-            _indices[i].push_back(index);
+            index.reserve( _fg->factor(I).states() );
+            for( IndexFor k(_fg->var(i), _fg->factor(I).vars()); k >= 0; ++k )
+                index.push_back( k );
+            _indices[i].push_back( index );
         }
     }
 }
@@ -330,7 +114,6 @@ void BBP::RegenerateInds() {
 
 void BBP::RegenerateT() {
     // _T[i][_I]
-    _T.clear();
     _T.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
         _T[i].resize( _fg->nbV(i).size() );
@@ -347,7 +130,6 @@ void BBP::RegenerateT() {
 
 void BBP::RegenerateU() {
     // _U[I][_i]
-    _U.clear();
     _U.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
         _U[I].resize( _fg->nbF(I).size() );
@@ -369,7 +151,6 @@ void BBP::RegenerateU() {
 
 void BBP::RegenerateS() {
     // _S[i][_I][_j]
-    _S.clear();
     _S.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
         _S[i].resize( _fg->nbV(i).size() );
@@ -387,7 +168,6 @@ void BBP::RegenerateS() {
                         }
                     }
                     // "Marginalize" onto i|j (unnormalized)
-                    // XXX improve efficiency?
                     Prob marg;
                     marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p();
                     _S[i][I.iter][j.iter] = marg;
@@ -399,7 +179,6 @@ void BBP::RegenerateS() {
 
 void BBP::RegenerateR() {
     // _R[I][_i][_J]
-    _R.clear();
     _R.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
         _R[I].resize( _fg->nbF(I).size() );
@@ -410,7 +189,7 @@ void BBP::RegenerateR() {
                     Prob prod( _fg->var(i).states(), 1.0 );
                     foreach( const Neighbor &K, _fg->nbV(i) )
                         if( K.node != I && K.node != J.node ) 
-                            prod *= _bp_dual.msgM(i,K.iter);
+                            prod *= _bp_dual.msgM( i, K.iter );
                     _R[I][i.iter][J.iter] = prod;
                 }
             }
@@ -421,109 +200,216 @@ void BBP::RegenerateR() {
 
 void BBP::RegenerateInputs() {
     _adj_b_V_unnorm.clear();
-    _adj_b_V_unnorm.reserve(_fg->nrVars());
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _adj_b_V_unnorm.push_back(
-            unnormAdjoint(_bp_dual.beliefV(i).p(),
-                _bp_dual.beliefVZ(i), _adj_b_V[i]));
-    }
+    _adj_b_V_unnorm.reserve( _fg->nrVars() );
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        _adj_b_V_unnorm.push_back( unnormAdjoint( _bp_dual.beliefV(i).p(), _bp_dual.beliefVZ(i), _adj_b_V[i] ) );
     _adj_b_F_unnorm.clear();
-    _adj_b_F_unnorm.reserve(_fg->nrFactors());
-    for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        _adj_b_F_unnorm.push_back(
-            unnormAdjoint(_bp_dual.beliefF(I).p(),
-                _bp_dual.beliefFZ(I), _adj_b_F[I]));
-    }
+    _adj_b_F_unnorm.reserve( _fg->nrFactors() );
+    for( size_t I = 0; I < _fg->nrFactors(); I++ )
+        _adj_b_F_unnorm.push_back( unnormAdjoint( _bp_dual.beliefF(I).p(), _bp_dual.beliefFZ(I), _adj_b_F[I] ) );
 }
 
 
 void BBP::RegeneratePsiAdjoints() {
     _adj_psi_V.clear();
-    _adj_psi_V.reserve(_fg->nrVars());
+    _adj_psi_V.reserve( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        Prob p(_adj_b_V_unnorm[i]);
-        assert(p.size()==_fg->var(i).states());
-        foreach(const Neighbor& I, _fg->nbV(i)) {
-            p *= _bp_dual.msgM(i,I.iter);
-        }
+        Prob p( _adj_b_V_unnorm[i] );
+        assert( p.size() == _fg->var(i).states() );
+        foreach( const Neighbor &I, _fg->nbV(i) )
+            p *= _bp_dual.msgM( i, I.iter );
         p += _init_adj_psi_V[i];
-        _adj_psi_V.push_back(p);
+        _adj_psi_V.push_back( p );
     }
     _adj_psi_F.clear();
-    _adj_psi_F.reserve(_fg->nrFactors());
+    _adj_psi_F.reserve( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        Prob p(_adj_b_F_unnorm[I]);
-        assert(p.size()==_fg->factor(I).states());
-        foreach(const Neighbor& i, _fg->nbF(I)) {
-            Prob n_iI(_bp_dual.msgN(i,i.dual));
-            const _ind_t& ind = _index(i,i.dual);
+        Prob p( _adj_b_F_unnorm[I] );
+        assert( p.size() == _fg->factor(I).states() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            Prob n_iI( _bp_dual.msgN( i, i.dual ) );
+            const _ind_t& ind = _index( i, i.dual );
             // multiply prod with n_jI
-            for(size_t x_I = 0; x_I < p.size(); x_I++)
+            for( size_t x_I = 0; x_I < p.size(); x_I++ )
                 p[x_I] *= n_iI[ind[x_I]];
         }
         p += _init_adj_psi_F[I];
-        _adj_psi_F.push_back(p);
+        _adj_psi_F.push_back( p );
     }
 }
 
 
 void BBP::RegenerateParMessageAdjoints() {
     size_t nv = _fg->nrVars();
-    _adj_n.resize(nv);
-    _adj_m.resize(nv);
-    _adj_n_unnorm.resize(nv);
-    _adj_m_unnorm.resize(nv);
-    _new_adj_n.clear();
-    _new_adj_n.resize(nv);
-    _new_adj_m.clear();
-    _new_adj_m.resize(nv);
+    _adj_n.resize( nv );
+    _adj_m.resize( nv );
+    _adj_n_unnorm.resize( nv );
+    _adj_m_unnorm.resize( nv );
+    _new_adj_n.resize( nv );
+    _new_adj_m.resize( nv );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
         size_t n_i = _fg->nbV(i).size();
-        _adj_n[i].resize(n_i);
-        _adj_m[i].resize(n_i);
-        _adj_n_unnorm[i].resize(n_i);
-        _adj_m_unnorm[i].resize(n_i);
-        _new_adj_n[i].resize(n_i);
-        _new_adj_m[i].resize(n_i);
-        foreach(const Neighbor& I, _fg->nbV(i)) {
-            // calculate adj_n
-            {
-                Prob prod(_fg->factor(I).p());
+        _adj_n[i].resize( n_i );
+        _adj_m[i].resize( n_i );
+        _adj_n_unnorm[i].resize( n_i );
+        _adj_m_unnorm[i].resize( n_i );
+        _new_adj_n[i].resize( n_i );
+        _new_adj_m[i].resize( n_i );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            { // calculate adj_n
+                Prob prod( _fg->factor(I).p() );
                 prod *= _adj_b_F_unnorm[I];
-                foreach(const Neighbor& j, _fg->nbF(I)) {
-                    if(i != j) {
-                        Prob n_jI(_bp_dual.msgN(j,j.dual));
-                        const _ind_t& ind = _index(j,j.dual);
+                foreach( const Neighbor &j, _fg->nbF(I) )
+                    if( i != j ) {
+                        Prob n_jI( _bp_dual.msgN( j, j.dual ) );
+                        const _ind_t &ind = _index( j, j.dual );
                         // multiply prod with n_jI
                         for( size_t x_I = 0; x_I < prod.size(); x_I++ )
                             prod[x_I] *= n_jI[ind[x_I]];
                     }
-                }
-                Prob marg(_fg->var(i).states(), 0.0);
-                const _ind_t &ind = _index(i,I.iter);
+                Prob marg( _fg->var(i).states(), 0.0 );
+                const _ind_t &ind = _index( i, I.iter );
                 for( size_t r = 0; r < prod.size(); r++ )
                     marg[ind[r]] += prod[r];
                 _new_adj_n[i][I.iter] = marg;
-                upMsgN(i,I.iter);
+                upMsgN( i, I.iter );
             }
 
-            // calculate adj_m
-            {
-                Prob prod(_adj_b_V_unnorm[i]);
-                assert(prod.size()==_fg->var(i).states());
-                foreach(const Neighbor& J, _fg->nbV(i)) {
-                    if(size_t(J) != size_t(I)) {
+            { // calculate adj_m
+                Prob prod( _adj_b_V_unnorm[i] );
+                assert( prod.size() == _fg->var(i).states() );
+                foreach( const Neighbor &J, _fg->nbV(i) )
+                    if( J.node != I.node )
                         prod *= _bp_dual.msgM(i,J.iter);
-                    }
-                }
                 _new_adj_m[i][I.iter] = prod;
-                upMsgM(i,I.iter);
+                upMsgM( i, I.iter );
             }
         }
     }
 }
 
 
+void BBP::RegenerateSeqMessageAdjoints() {
+    size_t nv = _fg->nrVars();
+    _adj_m.resize( nv );
+    _adj_m_unnorm.resize( nv );
+    _new_adj_m.resize( nv );
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        size_t n_i = _fg->nbV(i).size();
+        _adj_m[i].resize( n_i );
+        _adj_m_unnorm[i].resize( n_i );
+        _new_adj_m[i].resize( n_i );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            // calculate adj_m
+            Prob prod( _adj_b_V_unnorm[i] );
+            assert( prod.size() == _fg->var(i).states() );
+            foreach( const Neighbor &J, _fg->nbV(i) )
+                if( J.node != I.node )
+                    prod *= _bp_dual.msgM( i, J.iter );
+            _adj_m[i][I.iter] = prod;
+            calcUnnormMsgM( i, I.iter );
+            _new_adj_m[i][I.iter] = Prob( _fg->var(i).states(), 0.0 );
+        }
+    }
+    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            // calculate adj_n
+            Prob prod( _fg->factor(I).p() );
+            prod *= _adj_b_F_unnorm[I];
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i != j ) {
+                    Prob n_jI( _bp_dual.msgN( j, j.dual) );
+                    const _ind_t& ind = _index( j, j.dual );
+                    // multiply prod with n_jI
+                    for( size_t x_I = 0; x_I < prod.size(); x_I++ )
+                        prod[x_I] *= n_jI[ind[x_I]];
+                }
+            Prob marg( _fg->var(i).states(), 0.0 );
+            const _ind_t &ind = _index( i, I.iter );
+            for( size_t r = 0; r < prod.size(); r++ )
+                marg[ind[r]] += prod[r];
+            sendSeqMsgN( i, I.iter,marg );
+        }
+    }
+}
+
+
+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];
+    new_adj_n_iI = Prob( _fg->var(i).states(), 0.0 );
+    size_t I = _fg->nbV(i)[_I];
+    foreach( const Neighbor &j, _fg->nbF(I) )
+        if( j != i ) {
+            const Prob &p = _S[i][_I][j.iter];
+            const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][j.dual];
+            LOOP_ij(
+                new_adj_n_iI[xi] += p[xij] * _adj_m_unnorm_jI[xj];
+            );
+            /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
+            Var vi = _fg->var(i);
+            Var vj = _fg->var(j);
+            new_adj_n_iI = (Factor(VarSet(vi, vj), p) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p();
+            */
+        }
+}
+
+
+void BBP::calcNewM( size_t i, size_t _I ) {
+    const Neighbor &I = _fg->nbV(i)[_I];
+    Prob p( U(I, I.dual) );
+    const Prob &adj = _adj_m_unnorm[i][_I];
+    const _ind_t &ind = _index(i,_I);
+    for( size_t x_I = 0; x_I < p.size(); x_I++ )
+        p[x_I] *= adj[ind[x_I]];
+    _adj_psi_F[I] += p;
+    /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
+    _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
+    */
+
+    _new_adj_m[i][_I] = Prob( _fg->var(i).states(), 0.0 );
+    foreach( const Neighbor &J, _fg->nbV(i) )
+        if( J != I )
+            _new_adj_m[i][_I] += _R[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
+}
+
+
+void BBP::calcUnnormMsgN( size_t i, size_t _I ) {
+    _adj_n_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), _adj_n[i][_I] );
+}
+
+
+void BBP::calcUnnormMsgM( size_t i, size_t _I ) {
+    _adj_m_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgM(i,_I), _bp_dual.zM(i,_I), _adj_m[i][_I] );
+}
+
+
+void BBP::upMsgN( size_t i, size_t _I ) {
+    _adj_n[i][_I] = _new_adj_n[i][_I];
+    calcUnnormMsgN( i, _I );
+}
+
+
+void BBP::upMsgM( size_t i, size_t _I ) {
+    _adj_m[i][_I] = _new_adj_m[i][_I];
+    calcUnnormMsgM( i, _I );
+}
+
+
+void BBP::doParUpdate() {
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            calcNewM( i, I.iter );
+            calcNewN( i, I.iter );
+        }
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            upMsgM( i, I.iter );
+            upMsgN( i, I.iter );
+        }
+}
+
+
 void BBP::incrSeqMsgM( size_t i, size_t _I, const Prob &p ) {
     if( props.clean_updates )
         _new_adj_m[i][_I] += p;
@@ -561,7 +447,7 @@ void BBP::updateSeqMsgM( size_t i, size_t _I ) {
 
 void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) {
     _adj_m[i][_I] = p;
-    calcUnnormMsgM(i, _I);
+    calcUnnormMsgM( i, _I );
 }
 
 
@@ -569,7 +455,7 @@ void BBP::sendSeqMsgN( size_t i, size_t _I, const Prob &f ) {
     Prob f_unnorm = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), f );
     const Neighbor &I = _fg->nbV(i)[_I];
     assert( I.iter == _I );
-    _adj_psi_V[i] += f_unnorm * T(i,_I);
+    _adj_psi_V[i] += f_unnorm * T( i, _I );
 #if 0
     if(f_unnorm.sumAbs() > pv_thresh) {
         DAI_DMSG("in sendSeqMsgN");
@@ -665,146 +551,6 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) {
 }
 
 
-void BBP::RegenerateSeqMessageAdjoints() {
-    size_t nv = _fg->nrVars();
-    _adj_m.resize(nv);
-    _adj_m_unnorm.resize(nv);
-    _new_adj_m.resize(nv);
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        size_t n_i = _fg->nbV(i).size();
-        _adj_m[i].resize(n_i);
-        _adj_m_unnorm[i].resize(n_i);
-        _new_adj_m[i].resize(n_i);
-        foreach(const Neighbor& I, _fg->nbV(i)) {
-            // calculate adj_m
-            Prob prod(_adj_b_V_unnorm[i]);
-            assert(prod.size()==_fg->var(i).states());
-            foreach(const Neighbor& J, _fg->nbV(i)) {
-                if(size_t(J) != size_t(I)) {
-                    prod *= _bp_dual.msgM(i,J.iter);
-                }
-            }
-            _adj_m[i][I.iter] = prod;
-            calcUnnormMsgM(i, I.iter);
-            _new_adj_m[i][I.iter] = Prob(_fg->var(i).states(),0.0);
-        }
-    }
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor& I, _fg->nbV(i)) {
-            // calculate adj_n
-            Prob prod(_fg->factor(I).p());
-            prod *= _adj_b_F_unnorm[I];
-            foreach(const Neighbor& j, _fg->nbF(I)) {
-                if(i != j) {
-                    Prob n_jI(_bp_dual.msgN(j,j.dual));
-                    const _ind_t& ind = _index(j,j.dual);
-                    // multiply prod with n_jI
-                    for( size_t x_I = 0; x_I < prod.size(); x_I++ )
-                        prod[x_I] *= n_jI[ind[x_I]];
-                }
-            }
-            Prob marg(_fg->var(i).states(), 0.0);
-            const _ind_t &ind = _index(i,I.iter);
-            for( size_t r = 0; r < prod.size(); r++ )
-                marg[ind[r]] += prod[r];
-            sendSeqMsgN(i,I.iter,marg);
-        }
-    }
-}
-
-
-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];
-    new_adj_n_iI = Prob( _fg->var(i).states(), 0.0 );
-    size_t I = _fg->nbV(i)[_I];
-    foreach( const Neighbor &j, _fg->nbF(I) )
-        if( j != i ) {
-            const Prob &p = _S[i][_I][j.iter];
-            const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][j.dual];
-            LOOP_ij(
-                new_adj_n_iI[xi] += p[xij] * _adj_m_unnorm_jI[xj];
-            );
-            /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
-            Var vi = _fg->var(i);
-            Var vj = _fg->var(j);
-            new_adj_n_iI = (Factor(VarSet(vi, vj), p) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p();
-            */
-        }
-}
-
-
-void BBP::calcNewM( size_t i, size_t _I ) {
-    const Neighbor &I = _fg->nbV(i)[_I];
-    Prob p( U(I, I.dual) );
-    const Prob &adj = _adj_m_unnorm[i][_I];
-    const _ind_t &ind = _index(i,_I);
-    for( size_t x_I = 0; x_I < p.size(); x_I++ )
-        p[x_I] *= adj[ind[x_I]];
-    _adj_psi_F[I] += p;
-    /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
-    _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
-    */
-
-    _new_adj_m[i][_I] = Prob( _fg->var(i).states(), 0.0 );
-    foreach( const Neighbor &J, _fg->nbV(i) )
-        if( J != I )
-            _new_adj_m[i][_I] += _R[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
-}
-
-
-void BBP::calcUnnormMsgN( size_t i, size_t _I ) {
-    _adj_n_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), _adj_n[i][_I] );
-}
-
-
-void BBP::calcUnnormMsgM( size_t i, size_t _I ) {
-    _adj_m_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgM(i,_I), _bp_dual.zM(i,_I), _adj_m[i][_I] );
-}
-
-
-void BBP::upMsgN( size_t i, size_t _I ) {
-    _adj_n[i][_I] = _new_adj_n[i][_I];
-    calcUnnormMsgN( i, _I );
-}
-
-
-void BBP::upMsgM( size_t i, size_t _I ) {
-    _adj_m[i][_I] = _new_adj_m[i][_I];
-    calcUnnormMsgM( i, _I );
-}
-
-
-void BBP::doParUpdate() {
-    for( size_t i = 0; i < _fg->nrVars(); i++ )
-        foreach( const Neighbor &I, _fg->nbV(i) ) {
-            calcNewM( i, I.iter );
-            calcNewN( i, I.iter );
-        }
-    for( size_t i = 0; i < _fg->nrVars(); i++ )
-        foreach( const Neighbor &I, _fg->nbV(i) ) {
-            upMsgM( i, I.iter );
-            upMsgN( i, I.iter );
-        }
-}
-
-
 Real BBP::getUnMsgMag() {
     Real s = 0.0;
     size_t e = 0;
@@ -903,183 +649,178 @@ Real BBP::getTotalMsgN() {
 }
 
 
-void BBP::run() {
-    typedef BBP::Properties::UpdateType UT;
-    Real tol = props.tol;
-    UT updates = props.updates;
-    Real tic=toc();
-    switch((size_t)updates) {
-    case UT::SEQ_MAX: {
-        size_t i, _I;
-        Real mag;
-        do {
-            _iters++;
-            getArgmaxMsgM(i,_I,mag);
-            sendSeqMsgM(i,_I);
-        } while(mag > tol && _iters < props.maxiter);
-
-        if(_iters >= props.maxiter) {
-            cerr << "Warning: BBP didn't converge in " << _iters
-                 << " iterations (greatest message magnitude = " << mag << ")"
-                 << endl;
-        }
-        break;
-    }
-    case UT::SEQ_FIX: {
-        Real mag;
-        do {
-            _iters++;
-            mag = getTotalMsgM();
-            if(mag < tol) break;
-
-            for( size_t i = 0; i < _fg->nrVars(); i++ )
-                foreach(const Neighbor &I, _fg->nbV(i))
-                    sendSeqMsgM(i, I.iter);
-            for( size_t i = 0; i < _fg->nrVars(); i++ )
-                foreach(const Neighbor &I, _fg->nbV(i))
-                    updateSeqMsgM(i, I.iter);
-        } while(mag > tol && _iters < props.maxiter);
-
-        if(_iters >= props.maxiter) {
-            cerr << "Warning: BBP didn't converge in " << _iters
-                 << " iterations (greatest message magnitude = " << mag << ")"
-                 << endl;
-            break;
-        }
-        break;
-    }
-    case UT::SEQ_BP_REV:
-    case UT::SEQ_BP_FWD: {
-        const BP *bp = static_cast<const BP*>(_ia);
-        vector<pair<size_t, size_t> > sentMessages = bp->getSentMessages();
-        size_t totalMessages = sentMessages.size();
-        if(totalMessages==0) {
-            cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl;
-            DAI_THROW(INTERNAL_ERROR);
-        }
-        if(updates==UT::SEQ_BP_FWD) {
-            reverse(sentMessages.begin(), sentMessages.end());
-        }
-//         DAI_PV(sentMessages.size());
-//         DAI_PV(_iters);
-//         DAI_PV(props.maxiter);
-        while(sentMessages.size()>0 && _iters < props.maxiter) {
-//             DAI_PV(sentMessages.size());
-//             DAI_PV(_iters);
-            _iters++;
-            pair<size_t, size_t> e = sentMessages.back();
-            sentMessages.pop_back();
-            size_t i = e.first, _I = e.second;
-            sendSeqMsgM(i,_I);
-        }
-        if(_iters >= props.maxiter) {
-            cerr << "Warning: BBP updates limited to " << props.maxiter
-                 << " iterations, but using UpdateType " << updates
-                 << " with " << totalMessages << " messages"
-                 << endl;
-        }
-        break;
-    }
-    case UT::PAR: {
-        do {
-            _iters++;
-            doParUpdate();
-        } while((_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter);
-        if(_iters==props.maxiter) {
-            Real s, new_s;
-            getMsgMags(s,new_s);
-            cerr << "Warning: BBP didn't converge in " << _iters
-                 << " iterations (unnorm message magnitude = " << getUnMsgMag()
-                 << ", norm message mags = " << s << " -> " << new_s
-                 << ")" << endl;
-        }
-        break;
-    }
-    }
-    if(props.verbose >= 3) {
-        cerr << "BBP::run() took " << toc()-tic << " seconds "
-             << doneIters() << " iterations" << endl;
-    }
+void BBP::Regenerate() {
+    RegenerateInds();
+    RegenerateT();
+    RegenerateU();
+    RegenerateS();
+    RegenerateR();
+    RegenerateInputs();
+    RegeneratePsiAdjoints();
+    if( props.updates == Properties::UpdateType::PAR )
+        RegenerateParMessageAdjoints();
+    else
+        RegenerateSeqMessageAdjoints();
+    _iters = 0;
 }
 
 
-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;
-    }
+std::vector<Prob> BBP::getZeroAdjF( const FactorGraph &fg ) {
+    vector<Prob> adj_2;
+    adj_2.reserve( fg.nrFactors() );
+    for( size_t I = 0; I < fg.nrFactors(); I++ )
+        adj_2.push_back( Prob( fg.factor(I).states(), 0.0 ) );
+    return adj_2;
 }
 
 
-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();
+std::vector<Prob> BBP::getZeroAdjV( const FactorGraph &fg ) {
+    vector<Prob> adj_1;
+    adj_1.reserve( fg.nrVars() );
+    for( size_t i = 0; i < fg.nrVars(); i++ )
+        adj_1.push_back( Prob( fg.var(i).states(), 0.0 ) );
+    return adj_1;
 }
 
 
-double numericBBPTest(const InfAlg& bp, const vector<size_t> *state, const PropertySet& bbp_props, bbp_cfn_t cfn, double h) {
+void BBP::run() {
+    typedef BBP::Properties::UpdateType UT;
+    Real &tol = props.tol;
+    UT &updates = props.updates;
+
+    Real tic = toc();
+    switch( (size_t)updates ) {
+        case UT::SEQ_MAX: {
+            size_t i, _I;
+            Real mag;
+            do {
+                _iters++;
+                getArgmaxMsgM( i, _I, mag );
+                sendSeqMsgM( i, _I );
+            } while( mag > tol && _iters < props.maxiter );
+
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl;
+            break;
+        } case UT::SEQ_FIX: {
+            Real mag;
+            do {
+                _iters++;
+                mag = getTotalMsgM();
+                if( mag < tol ) 
+                    break;
+
+                for( size_t i = 0; i < _fg->nrVars(); i++ )
+                    foreach( const Neighbor &I, _fg->nbV(i) )
+                        sendSeqMsgM( i, I.iter );
+                for( size_t i = 0; i < _fg->nrVars(); i++ )
+                    foreach( const Neighbor &I, _fg->nbV(i) )
+                        updateSeqMsgM( i, I.iter );
+            } while( mag > tol && _iters < props.maxiter );
+
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl;
+            break;
+        } case UT::SEQ_BP_REV:
+          case UT::SEQ_BP_FWD: {
+            const BP *bp = static_cast<const BP*>(_ia);
+            vector<pair<size_t, size_t> > sentMessages = bp->getSentMessages();
+            size_t totalMessages = sentMessages.size();
+            if( totalMessages == 0 ) {
+                cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl;
+                DAI_THROW(INTERNAL_ERROR);
+            }
+            if( updates==UT::SEQ_BP_FWD )
+                reverse( sentMessages.begin(), sentMessages.end() );
+//          DAI_PV(sentMessages.size());
+//          DAI_PV(_iters);
+//          DAI_PV(props.maxiter);
+            while( sentMessages.size() > 0 && _iters < props.maxiter ) {
+//              DAI_PV(sentMessages.size());
+//              DAI_PV(_iters);
+                _iters++;
+                pair<size_t, size_t> e = sentMessages.back();
+                sentMessages.pop_back();
+                size_t i = e.first, _I = e.second;
+                sendSeqMsgM( i, _I );
+            }
+            if( _iters >= props.maxiter )
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP updates limited to " << props.maxiter << " iterations, but using UpdateType " << updates << " with " << totalMessages << " messages" << endl;
+            break;
+        } case UT::PAR: {
+            do {
+                _iters++;
+                doParUpdate();
+            } while( (_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter );
+            if( _iters == props.maxiter ) {
+                Real s, new_s;
+                getMsgMags( s, new_s );
+                if( props.verbose >= 1 )
+                    cerr << "Warning: BBP didn't converge in " << _iters << " iterations (unnorm message magnitude = " << getUnMsgMag() << ", norm message mags = " << s << " -> " << new_s << ")" << endl;
+            }
+            break;
+        }
+    }
+    if( props.verbose >= 3 )
+        cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl;
+}
+
+
+double numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) {
     // calculate the value of the unperturbed cost function
-    Real cf0 = getCostFn(bp, cfn, state);
+    Real cf0 = getCostFn( bp, cfn, state );
 
     // run BBP to estimate adjoints
-    BBP bbp(&bp,bbp_props);
-    initBBPCostFnAdj(bbp, bp, cfn, state);
+    BBP bbp( &bp, bbp_props );
+    initBBPCostFnAdj( bbp, bp, cfn, state );
     bbp.run();
 
-    Real d=0;
+    Real d = 0;
     const FactorGraph& fg = bp.fg();
 
-    if(1) {
+    if( 1 ) {
         // verify bbp.adj_psi_V
 
         // for each variable i
-        for(size_t i=0; i<fg.nrVars(); i++) {
+        for( size_t i = 0; i < fg.nrVars(); i++ ) {
             vector<double> adj_est;
             // for each value xi
-            for(size_t xi=0; xi<fg.var(i).states(); xi++) {
+            for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
                 // Clone 'bp' (which may be any InfAlg)
                 InfAlg *bp_prb = bp.clone();
 
                 // perturb it
                 size_t n = bp_prb->fg().var(i).states();
-                Prob psi_1_prb(n,1.0);
+                Prob psi_1_prb( n, 1.0 );
                 psi_1_prb[xi] += h;
 //                 psi_1_prb.normalize();
                 size_t I = bp_prb->fg().nbV(i)[0]; // use first factor in list of neighbors of i
-                bp_prb->fg().factor(I) *= Factor(bp_prb->fg().var(i), psi_1_prb);
+                bp_prb->fg().factor(I) *= Factor( bp_prb->fg().var(i), psi_1_prb );
                 
                 // call 'init' on the perturbed variables
-                bp_prb->init(bp_prb->fg().var(i));
+                bp_prb->init( bp_prb->fg().var(i) );
                 
                 // run copy to convergence
                 bp_prb->run();
 
                 // calculate new value of cost function
-                Real cf_prb = getCostFn(*bp_prb, cfn, state);
+                Real cf_prb = getCostFn( *bp_prb, cfn, state );
 
                 // use to estimate adjoint for i
-                adj_est.push_back((cf_prb-cf0)/h);
+                adj_est.push_back( (cf_prb - cf0) / h );
                 
                 // free cloned InfAlg
                 delete bp_prb;
             }
-            Prob p_adj_est(adj_est.begin(), adj_est.end());
+            Prob p_adj_est( adj_est.begin(), adj_est.end() );
             // compare this numerical estimate to the BBP estimate; sum the distances
-            cerr << "i: " << i
+            cout << "i: " << i
                  << ", p_adj_est: " << p_adj_est
                  << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
-            d += dist(p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1);
+            d += dist( p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1 );
         }
     }
     /*    if(1) {
@@ -1173,6 +914,243 @@ double numericBBPTest(const InfAlg& bp, const vector<size_t> *state, const Prope
 }
 
 
+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++ ) {
+                    double 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++ ) {
+                    double 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;
+            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 );
+                assert(/*0<=state[i] &&*/ state[i] < n);
+                double 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(INTERNAL_ERROR);
+                }
+                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;
+            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 );
+                assert(/*0<=x_I &&*/ x_I < n);
+
+                double 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(INTERNAL_ERROR);
+                }
+                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 ) {
+    double 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: {
+            assert( stateP != NULL );
+            vector<size_t> state = *stateP;
+            assert( state.size() == fg.nrVars() );
+            for( size_t i = 0; i < fg.nrVars(); i++ ) {
+                double 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(INTERNAL_ERROR);
+                }
+            }
+            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: {
+            assert( stateP != NULL );
+            vector<size_t> state = *stateP;
+            assert( state.size() == fg.nrVars() );
+            for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+                size_t x_I = getFactorEntryForState( fg, I, state );
+                double b = ia.beliefF(I)[x_I];
+                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(INTERNAL_ERROR);
+                }
+            }
+            break;
+        } default: 
+            cerr << "Unknown cost function " << cfn_type << endl;
+            DAI_THROW(UNKNOWN_ENUM_VALUE);
+    }
+    return cf;
+}
+
+
 } // end of namespace dai
 
 
index dfa4802..dc9a71b 100644 (file)
@@ -104,5 +104,5 @@ GIBBS_1e9:                      GIBBS[iters=1000000000]
 
 # --- CBP ---------------------
 
-CBP:                            CBP[max_levels=12,updates=SEQMAX,tol=1e-9,rec_tol=1e-9,maxiter=500,choose=CHOOSE_RANDOM,recursion=REC_FIXED,clamp=CLAMP_VAR,min_max_adj=1.0e-9,bbp_cfn=cfn_factor_ent,verbose=0,rand_seed=0,bbp_props=[verbose=0,tol=1.0e-9,maxiter=10000,damping=0,updates=SEQ_BP_REV,clean_updates=0],clamp_outfile=]
+CBP:                            CBP[max_levels=12,updates=SEQMAX,tol=1e-9,rec_tol=1e-9,maxiter=500,choose=CHOOSE_RANDOM,recursion=REC_FIXED,clamp=CLAMP_VAR,min_max_adj=1.0e-9,bbp_cfn=CFN_FACTOR_ENT,verbose=0,rand_seed=0,bbp_props=[verbose=0,tol=1.0e-9,maxiter=10000,damping=0,updates=SEQ_BP_REV,clean_updates=0],clamp_outfile=]
 BBP:                            CBP[choose=CHOOSE_BBP]
index 502d869..8c061cb 100644 (file)
@@ -85,31 +85,31 @@ int main( int argc, char *argv[] ) {
             bbp_cfn_t cfn;
             switch( (t / 10) % 9 ) {
                 case 0:
-                    cfn = bbp_cfn_t::cfn_gibbs_b;
+                    cfn = bbp_cfn_t::CFN_GIBBS_B;
                     break;
                 case 1:
-                    cfn = bbp_cfn_t::cfn_gibbs_b2;
+                    cfn = bbp_cfn_t::CFN_GIBBS_B2;
                     break;
                 case 2:
-                    cfn = bbp_cfn_t::cfn_gibbs_exp;
+                    cfn = bbp_cfn_t::CFN_GIBBS_EXP;
                     break;
                 case 3:
-                    cfn = bbp_cfn_t::cfn_gibbs_b_factor;
+                    cfn = bbp_cfn_t::CFN_GIBBS_B_FACTOR;
                     break;
                 case 4:
-                    cfn = bbp_cfn_t::cfn_gibbs_b2_factor;
+                    cfn = bbp_cfn_t::CFN_GIBBS_B2_FACTOR;
                     break;
                 case 5:
-                    cfn = bbp_cfn_t::cfn_gibbs_exp_factor;
+                    cfn = bbp_cfn_t::CFN_GIBBS_EXP_FACTOR;
                     break;
                 case 6:
-                    cfn = bbp_cfn_t::cfn_var_ent;
+                    cfn = bbp_cfn_t::CFN_VAR_ENT;
                     break;
                 case 7:
-                    cfn = bbp_cfn_t::cfn_factor_ent;
+                    cfn = bbp_cfn_t::CFN_FACTOR_ENT;
                     break;
                 case 8:
-                    cfn = bbp_cfn_t::cfn_bethe_ent;
+                    cfn = bbp_cfn_t::CFN_BETHE_ENT;
                     break;
             }