Cleanup of BBP code
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 3 Aug 2009 15:44:14 +0000 (17:44 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 3 Aug 2009 15:44:14 +0000 (17:44 +0200)
12 files changed:
Makefile
include/dai/bbp.h
include/dai/bp.h
include/dai/cbp.h
include/dai/lc.h
include/dai/mr.h
scripts/regenerate-properties
src/bbp.cpp
tests/aliases.conf
tests/testall
tests/testbbp.cpp [new file with mode: 0644]
tests/testfast.out

index 19ac0fa..d47ce9e 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -109,7 +109,7 @@ examples : examples/example$(EE) examples/example_bipgraph$(EE) examples/example
 
 matlabs : matlab/dai$(ME) matlab/dai_readfg$(ME) matlab/dai_writefg$(ME) matlab/dai_potstrength$(ME)
 
-tests : tests/testdai$(EE)
+tests : tests/testdai$(EE) tests/testbbp$(EE)
 
 utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
 
@@ -208,6 +208,9 @@ examples/example_sprinkler$(EE) : examples/example_sprinkler.cpp $(HEADERS) $(LI
 tests/testdai$(EE) : tests/testdai.cpp $(HEADERS) $(LIB)/libdai$(LE)
        $(CC) $(CCO)tests/testdai$(EE) tests/testdai.cpp $(LIBS) $(BOOSTLIBS)
 
+tests/testbbp$(EE) : tests/testbbp.cpp $(HEADERS) $(LIB)/libdai$(LE)
+       $(CC) $(CCO)tests/testbbp$(EE) tests/testbbp.cpp $(LIBS)
+
 
 # MATLAB INTERFACE
 ###################
@@ -289,7 +292,7 @@ clean :
        -rm *$(OE)
        -rm matlab/*$(ME)
        -rm examples/example$(EE) examples/example_bipgraph$(EE) examples/example_varset$(EE) examples/example_sprinkler$(EE)
-       -rm tests/testdai$(EE)
+       -rm tests/testdai$(EE) tests/testbbp$(EE)
        -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE)
        -rm -R doc
        -rm -R lib
index a71bb7a..bf68121 100644 (file)
@@ -22,6 +22,7 @@
 /// \brief Defines class BBP [\ref EaG09]
 /// \todo Improve documentation
 /// \todo Clean up
+/// \todo Debug clean_updates
 
 
 #ifndef ___defined_libdai_bbp_h
 #include <dai/daialg.h>
 #include <dai/factorgraph.h>
 #include <dai/enum.h>
-
 #include <dai/bp_dual.h>
 
 
 namespace dai {
 
 
-std::vector<Prob> get_zero_adj_F(const FactorGraph&);
-std::vector<Prob> get_zero_adj_V(const FactorGraph&);
+/// 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 );
+
+/// 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 );
 
 
-/// Implements BBP (Back-belief-propagation) [\ref EaG09]
+/// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
 class BBP {
-protected:
-    // ----------------------------------------------------------------
-    // inputs
-    BP_dual _bp_dual;
-    const FactorGraph* _fg;
-    const InfAlg *_ia;
-
-    // iterations
-    size_t _iters;
-
-    // ----------------------------------------------------------------
-    // Outputs
-    std::vector<Prob> _adj_psi_V, _adj_psi_F;
-    // The following vectors are indexed [i][_I]
-    std::vector<std::vector<Prob> > _adj_n, _adj_m; 
-    std::vector<Prob> _adj_b_V, _adj_b_F;
-
-    // Helper quantities computed from the BP messages:
-    // _T[i][_I]
-    std::vector<std::vector<Prob > > _T;
-    // _U[I][_i]
-    std::vector<std::vector<Prob > > _U;
-    // _S[i][_I][_j]
-    std::vector<std::vector<std::vector<Prob > > > _S;
-    // _R[I][_i][_J]
-    std::vector<std::vector<std::vector<Prob > > > _R;
-
-    std::vector<Prob> _adj_b_V_unnorm, _adj_b_F_unnorm;
-    std::vector<Prob> _init_adj_psi_V;
-    std::vector<Prob> _init_adj_psi_F;
-
-    std::vector<std::vector<Prob> > _adj_n_unnorm, _adj_m_unnorm;
-    std::vector<std::vector<Prob> > _new_adj_n, _new_adj_m;
-
-    // ----------------------------------------------------------------
-    // Indexing for performance
-  
-    /// Calculates _indices, which is a cache of IndexFor (see bp.cpp)
-    void RegenerateInds();
-    
-    typedef std::vector<size_t>  _ind_t;
-    std::vector<std::vector<_ind_t> >  _indices; 
-    const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
-
-    // ----------------------------------------------------------------
-    // Initialization
-
-    /// Calculate T values (see paper)
-    void RegenerateT();
-    /// Calculate U values (see paper)
-    void RegenerateU();
-    /// Calculate S values (see paper)
-    void RegenerateS();
-    /// Calculate R values (see paper)
-    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)
-    void RegeneratePsiAdjoints();
-    /// Initialise members for messages adjoints (call after RegenerateInputs)
-    void RegenerateParMessageAdjoints();
-    /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
-     *  than updating _adj_n (and friends) which are unused in sequential algorithm
-     */
-    void RegenerateSeqMessageAdjoints(); 
-
-    DAI_ACCMUT(Prob & T(size_t i, size_t _I), { return _T[i][_I]; });
-    DAI_ACCMUT(Prob & U(size_t I, size_t _i), { return _U[I][_i]; });
-    DAI_ACCMUT(Prob & S(size_t i, size_t _I, size_t _j), { return _S[i][_I][_j]; });
-    DAI_ACCMUT(Prob & R(size_t I, size_t _i, size_t _J), { return _R[I][_i][_J]; });
-
-    void calcNewN(size_t i, size_t _I);
-    void calcNewM(size_t i, size_t _I);
-    void calcUnnormMsgM(size_t i, size_t _I);
-    void calcUnnormMsgN(size_t i, size_t _I);
-    void upMsgM(size_t i, size_t _I);
-    void upMsgN(size_t i, size_t _I);
-    void doParUpdate();
-    Real getUnMsgMag();
-    void getMsgMags(Real &s, Real &new_s);
-
-    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)));
+    protected:
+        /// @name Inputs
+        //@{
+        BP_dual _bp_dual;
+        const FactorGraph *_fg;
+        const InfAlg *_ia;
+        //@}
+
+        /// Number of iterations done
+        size_t _iters;
+
+        /// @name Outputs
+        //@{
+        /// Variable factor adjoints
+        std::vector<Prob> _adj_psi_V;
+        /// Factor adjoints
+        std::vector<Prob> _adj_psi_F;
+        /// Variable->factor message adjoints (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_n;
+        /// Factor->variable message adjoints (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_m;
+        /// Normalized variable belief adjoints
+        std::vector<Prob> _adj_b_V;
+        /// Normalized factor belief adjoints
+        std::vector<Prob> _adj_b_F;
+        //@}
+
+        /// @name Helper quantities computed from the BP messages
+        //@{
+        /// _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
+        std::vector<Prob> _init_adj_psi_F;
+
+        /// Unnormalized variable->factor message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_n_unnorm;
+        /// Unnormalized factor->variable message adjoint (indexed [i][_I])
+        std::vector<std::vector<Prob> > _adj_m_unnorm;
+        /// Updated normalized variable->factor message adjoint (indexed [i][_I])
+        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;
+
+        /// @name Optimized indexing (for performance)
+        //@{
+        /// Calculates _indices, which is a cache of IndexFor @see bp.cpp
+        void RegenerateInds();
+        
+        /// Index type
+        typedef std::vector<size_t>  _ind_t;
+        /// Cached indices (indexed [i][_I])
+        std::vector<std::vector<_ind_t> >  _indices; 
+        /// Returns an index from the cache
+        const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
+        //@}
+
+        /// @name Initialization
+        //@{
+        /// Calculate T values; see eqn. (41) in [\ref EaG09]
+        void RegenerateT();
+        /// Calculate U values; see eqn. (42) in [\ref EaG09]
+        void RegenerateU();
+        /// Calculate S values; see eqn. (43) in [\ref EaG09]
+        void RegenerateS();
+        /// Calculate R values; see eqn. (44) in [\ref EaG09]
+        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)
+        void RegeneratePsiAdjoints();
+        /// Initialise members for message adjoints (call after RegenerateInputs) for parallel algorithm
+        void RegenerateParMessageAdjoints();
+        /// Initialise members for message adjoints (call after RegenerateInputs) for sequential algorithm
+        /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
+         *  than updating _adj_n (and friends) which are unused in the sequential algorithm.
+         */
+        void RegenerateSeqMessageAdjoints();
+        //@}
+
+        /// Returns T value; see eqn. (41) in [\ref EaG09]
+        DAI_ACCMUT(Prob & T(size_t i, size_t _I), { return _T[i][_I]; });
+        /// Retunrs U value; see eqn. (42) in [\ref EaG09]
+        DAI_ACCMUT(Prob & U(size_t I, size_t _i), { return _U[I][_i]; });
+        /// Returns S value; see eqn. (43) in [\ref EaG09]
+        DAI_ACCMUT(Prob & S(size_t i, size_t _I, size_t _j), { return _S[i][_I][_j]; });
+        /// Returns R value; see eqn. (44) in [\ref EaG09]
+        DAI_ACCMUT(Prob & R(size_t I, size_t _i, size_t _J), { return _R[I][_i][_J]; });
+
+        /// @name Parallel algorithm
+        //@{
+        /// Calculates new variable->factor message adjoint
+        /** Increases variable factor adjoint according to eqn. (27) in [\ref EaG09] and
+         *  calculates the new variable->factor message adjoint according to eqn. (29) in [\ref EaG09].
+         */
+        void calcNewN( size_t i, size_t _I );
+        /// Calculates new factor->variable message adjoint
+        /** Increases factor adjoint according to eqn. (28) in [\ref EaG09] and
+         *  calculates the new factor->variable message adjoint according to the r.h.s. of eqn. (30) in [\ref EaG09].
+         */
+        void calcNewM( size_t i, size_t _I );
+        /// Calculates unnormalized variable->factor message adjoint from the normalized one
+        void calcUnnormMsgN( size_t i, size_t _I );
+        /// Calculates unnormalized factor->variable message adjoint from the normalized one
+        void calcUnnormMsgM( size_t i, size_t _I );
+        /// Updates (un)normalized variable->factor message adjoints
+        void upMsgN( size_t i, size_t _I );
+        /// Updates (un)normalized factor->variable message adjoints
+        void upMsgM( size_t i, size_t _I );
+        /// Do one parallel update of all message adjoints
+        void doParUpdate();
+        //@}
+
+        /// Calculates averaged L-1 norm of unnormalized message adjoints
+        Real getUnMsgMag();
+        /// Calculates averaged L-1 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 ) ) );
         }
-    }
-
-    //----------------------------------------------------------------
-    // new interface
-
-    void incrSeqMsgM(size_t i, size_t _I, const Prob& p);
-    void updateSeqMsgM(size_t i, size_t _I);
-    void sendSeqMsgN(size_t i, size_t _I, const Prob &f);
-    void sendSeqMsgM(size_t i, size_t _I);
-    /// used instead of upMsgM / calcNewM, calculates adj_m_unnorm as well
-    void setSeqMsgM(size_t i, size_t _I, const Prob &p); 
-
-    Real getMaxMsgM();
-    Real getTotalMsgM();
-    Real getTotalNewMsgM();
-    Real getTotalMsgN();
-
-    void getArgmaxMsgM(size_t &i, size_t &_I, Real &mag);
-
-public:
-    /// Called by 'init', recalculates intermediate values
-    void Regenerate();
-
-    BBP(const InfAlg *ia, const PropertySet &opts) :
-        _bp_dual(ia), _fg(&(ia->fg())), _ia(ia)
-    {
-        props.set(opts);
-    }
-
-    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;
-        _init_adj_psi_V = adj_psi_V;
-        _init_adj_psi_F = adj_psi_F;
-        Regenerate(); 
-    }
-    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));
-    }
-    void init(const std::vector<Prob> &adj_b_V) {
-        init(adj_b_V, get_zero_adj_F(*_fg));
-    }
-
-    /// run until change is less than given tolerance
-    void run();
-
-    size_t doneIters() { return _iters; }
-
-    DAI_ACCMUT(Prob& adj_psi_V(size_t i), { return _adj_psi_V[i]; });
-    DAI_ACCMUT(Prob& adj_psi_F(size_t I), { return _adj_psi_F[I]; });
-    DAI_ACCMUT(Prob& adj_b_V(size_t i), { return _adj_b_V[i]; });
-    DAI_ACCMUT(Prob& adj_b_F(size_t I), { return _adj_b_F[I]; });
- protected:
-    DAI_ACCMUT(Prob& adj_n(size_t i, size_t _I), { return _adj_n[i][_I]; });
-    DAI_ACCMUT(Prob& adj_m(size_t i, size_t _I), { return _adj_m[i][_I]; });
- public: 
-
-    /// Parameters of this inference algorithm
-/* PROPERTIES(props,BBP) {
-       DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
-       size_t verbose;
-       /// tolerance (not used for updates=SEQ_BP_{REV,FWD})
-       double tol;
-       size_t maxiter;
-       /// damping (0 for none)
-       double damping;
-       UpdateType updates;
-       bool clean_updates;
-    }
-*/
+
+        /// @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
+        Real getTotalMsgM();
+        /// Calculates sum of L1 norms of all updated normalized factor->variable message adjoints
+        Real getTotalNewMsgM();
+        /// Calculates sum of L1 norms of all normalized variable->factor message adjoints
+        Real getTotalMsgN();
+
+    public:
+        /// Called by \a init, recalculates intermediate values
+        void Regenerate();
+
+        /// Constructor
+        BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
+            props.set(opts);
+        }
+
+        /// 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;
+            _adj_b_F = adj_b_F;
+            _init_adj_psi_V = adj_psi_V;
+            _init_adj_psi_F = adj_psi_F;
+            Regenerate(); 
+        }
+
+        /// 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) );
+        }
+
+        /// 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));
+        }
+
+        /// Run until change is less than given tolerance
+        void run();
+
+        /// Return number of iterations done so far
+        size_t doneIters() { return _iters; }
+
+        /// Returns variable factor adjoint
+        DAI_ACCMUT(Prob& adj_psi_V(size_t i), { return _adj_psi_V[i]; });
+        /// Returns factor adjoint
+        DAI_ACCMUT(Prob& adj_psi_F(size_t I), { return _adj_psi_F[I]; });
+        /// Returns variable belief adjoint
+        DAI_ACCMUT(Prob& adj_b_V(size_t i), { return _adj_b_V[i]; });
+        /// Returns factor belief adjoint
+        DAI_ACCMUT(Prob& adj_b_F(size_t I), { return _adj_b_F[I]; });
+
+     protected:
+        /// Returns variable->factor message adjoint
+        DAI_ACCMUT(Prob& adj_n(size_t i, size_t _I), { return _adj_n[i][_I]; });
+        /// Returns factor->variable message adjoint
+        DAI_ACCMUT(Prob& adj_m(size_t i, size_t _I), { return _adj_m[i][_I]; });
+
+     public: 
+        /// Parameters of this algorithm
+        /* PROPERTIES(props,BBP) {
+           /// Enumeration of possible update schedules
+           DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
+
+           /// Verbosity
+           size_t verbose;
+
+           /// Maximum number of iterations
+           size_t maxiter;
+
+           /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
+           double tol;
+
+           /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
+           double damping;
+
+           /// Update schedule
+           UpdateType updates;
+
+           bool clean_updates;
+        } 
+        */
 /* {{{ GENERATED CODE: DO NOT EDIT. Created by 
     ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp 
 */
-    struct Properties {
-        DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
-        size_t verbose;
-        /// tolerance (not used for updates=SEQ_BP_{REV,FWD})
-        double tol;
-        size_t maxiter;
-        /// damping (0 for none)
-        double damping;
-        UpdateType updates;
-        bool clean_updates;
-
-        /// Set members from PropertySet
-        void set(const PropertySet &opts);
-        /// Get members into PropertySet
-        PropertySet get() const;
-        /// Convert to a string which can be parsed as a PropertySet
-        std::string toString() const;
-    } props;
+        struct Properties {
+            /// Enumeration of possible update schedules
+            DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
+            /// Verbosity
+            size_t verbose;
+            /// Maximum number of iterations
+            size_t maxiter;
+            /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
+            double tol;
+            /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
+            double damping;
+            /// Update schedule
+            UpdateType updates;
+            bool clean_updates;
+
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } props;
 /* }}} END OF GENERATED CODE */
 };
 
@@ -262,11 +349,8 @@ double numericBBPTest(const InfAlg& bp, const std::vector<size_t> *state, const
 // ----------------------------------------------------------------
 // Utility functions, some of which are used elsewhere
 
-/// Subtract 1 from a size_t, or return 0 if the argument is 0
-inline size_t oneLess(size_t v) { return v==0?v:v-1; }
-
-/// function to compute adj_w_unnorm from w, Z_w, adj_w
-Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w);
+/// 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);
index 73c2556..2d0a1f8 100644 (file)
@@ -64,10 +64,10 @@ class BP : public DAIAlgFG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible update schedules
-            DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL)
+            DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL);
 
             /// Enumeration of inference variants
-            DAI_ENUM(InfType,SUMPROD,MAXPROD)
+            DAI_ENUM(InfType,SUMPROD,MAXPROD);
         
             /// Verbosity
             size_t verbose;
index e138bcf..58b67f3 100644 (file)
@@ -58,192 +58,209 @@ std::pair<size_t, size_t> bbpFindClampVar(const InfAlg &in_bp, bool clampingVar,
  *  combined using logZ estimates.
  */
 class CBP : public DAIAlgFG {
-    std::vector<Factor> _beliefsV;
-    std::vector<Factor> _beliefsF;
-    double _logZ;
-    double _est_logZ;
-    double _old_est_logZ;
-
-    double _sum_level;
-    size_t _num_leaves;
-
-    boost::shared_ptr<std::ofstream> _clamp_ofstream;
-
-    bbp_cfn_t BBP_cost_function() {
-      return props.bbp_cfn;
-    }
-
-    void printDebugInfo();
-
-    /// Called by 'run', and by itself. Implements the main algorithm. 
-    /** Chooses a variable to clamp, recurses, combines the logZ and
-     *  beliefs estimates of the children, and returns the improved
-     *  estimates in \a lz_out and \a beliefs_out to its parent
-     */
-    void runRecurse(InfAlg *bp,
-                    double orig_logZ,
-                    std::vector<size_t> clamped_vars_list,
-                    size_t &num_leaves,
-                    size_t &choose_count,
-                    double &sum_level,
-                    Real &lz_out,
-                    std::vector<Factor>& beliefs_out);
-
-    /// Choose the next variable to clamp
-    /** Choose the next variable to clamp, given a converged InfAlg (\a bp),
-     *  and a vector of variables that are already clamped (\a
-     *  clamped_vars_list). Returns the chosen variable in \a i, and
-     *  the set of states in \a xis. If \a maxVarOut is non-NULL and
-     *  props.choose==CHOOSE_BBP then it is used to store the
-     *  adjoint of the chosen variable
-     */
-    virtual bool chooseNextClampVar(InfAlg* bp,
-                                    std::vector<size_t> &clamped_vars_list,
-                                    size_t &i, std::vector<size_t> &xis, Real *maxVarOut);
-
-    /// Return the InfAlg to use at each step of the recursion. 
-    /// \todo At present, only returns a BP instance
-    InfAlg* getInfAlg();
-
-    size_t _iters;
-    double _maxdiff;
-
-    void setBeliefs(const std::vector<Factor>& bs, double logZ) {
-        size_t i=0;
-        _beliefsV.clear(); _beliefsV.reserve(nrVars());
-        _beliefsF.clear(); _beliefsF.reserve(nrFactors());
-        for(i=0; i<nrVars(); i++) _beliefsV.push_back(bs[i]);
-        for(; i<nrVars()+nrFactors(); i++) _beliefsF.push_back(bs[i]);
-        _logZ = logZ;
-    }
-
-    void construct();
-
-  public:
-
-    //----------------------------------------------------------------
-
-    /// construct CBP object from FactorGraph
-    CBP(const FactorGraph &fg, const PropertySet &opts) : DAIAlgFG(fg) {
-        props.set(opts);
-        construct();
-    }
-
-    static const char *Name;
-
-    /// @name General InfAlg interface
-    //@{
-    virtual CBP* clone() const { return new CBP(*this); }
-    virtual CBP* create() const { DAI_THROW(NOT_IMPLEMENTED); }
-    virtual std::string identify() const { return std::string(Name) + props.toString(); }
-    virtual Factor belief (const Var &n) const { return _beliefsV[findVar(n)]; }
-    virtual Factor belief (const VarSet &) const { DAI_THROW(NOT_IMPLEMENTED); }
-    virtual std::vector<Factor> beliefs() const { return concat(_beliefsV, _beliefsF); }
-    virtual Real logZ() const { return _logZ; }
-    virtual void init() {};
-    virtual void init( const VarSet & ) {};
-    virtual double run();
-    virtual double maxDiff() const { return _maxdiff; }
-    virtual size_t Iterations() const { return _iters; }
-    //@}
-
-    Factor beliefV (size_t i) const { return _beliefsV[i]; }
-    Factor beliefF (size_t I) const { return _beliefsF[I]; }
-
-    //----------------------------------------------------------------
-
- public:
-    /// Parameters of this inference algorithm
-/* PROPERTIES(props,CBP) {
-        typedef BP::Properties::UpdateType UpdateType;
-        DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
-        DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
-        DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
-        
-        size_t verbose = 0;
-
-        /// Tolerance to use in BP
-        double tol;
-        /// Update style for BP
-        UpdateType updates;
-        /// Maximum number of iterations for BP
-        size_t maxiter;
-
-        /// Tolerance to use for controlling recursion depth (\a recurse
-        /// is REC_LOGZ or REC_BDIFF)
-        double rec_tol;
-        /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
-        size_t max_levels = 10;
-        /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
-        double min_max_adj;
-        /// Heuristic for choosing clamping variable
-        ChooseMethodType choose;
-        /// Method for deciding when to stop recursing
-        RecurseType recursion;
-        /// Whether to clamp variables or factors
-        ClampType clamp;
-        /// Properties to pass to BBP
-        PropertySet bbp_props;
-        /// Cost function to use for BBP
-        bbp_cfn_t bbp_cfn;
-        /// Random seed
-        size_t rand_seed = 0;
-
-        /// If non-empty, write clamping choices to this file
-        std::string clamp_outfile = "";
-   }
-*/
+    private:
+        std::vector<Factor> _beliefsV;
+        std::vector<Factor> _beliefsF;
+        double _logZ;
+        double _est_logZ;
+        double _old_est_logZ;
+
+        double _sum_level;
+        size_t _num_leaves;
+
+        boost::shared_ptr<std::ofstream> _clamp_ofstream;
+
+        bbp_cfn_t BBP_cost_function() {
+          return props.bbp_cfn;
+        }
+
+        void printDebugInfo();
+
+        /// Called by 'run', and by itself. Implements the main algorithm. 
+        /** Chooses a variable to clamp, recurses, combines the logZ and
+         *  beliefs estimates of the children, and returns the improved
+         *  estimates in \a lz_out and \a beliefs_out to its parent
+         */
+        void runRecurse(InfAlg *bp,
+                        double orig_logZ,
+                        std::vector<size_t> clamped_vars_list,
+                        size_t &num_leaves,
+                        size_t &choose_count,
+                        double &sum_level,
+                        Real &lz_out,
+                        std::vector<Factor>& beliefs_out);
+
+        /// Choose the next variable to clamp
+        /** Choose the next variable to clamp, given a converged InfAlg (\a bp),
+         *  and a vector of variables that are already clamped (\a
+         *  clamped_vars_list). Returns the chosen variable in \a i, and
+         *  the set of states in \a xis. If \a maxVarOut is non-NULL and
+         *  props.choose==CHOOSE_BBP then it is used to store the
+         *  adjoint of the chosen variable
+         */
+        virtual bool chooseNextClampVar(InfAlg* bp,
+                                        std::vector<size_t> &clamped_vars_list,
+                                        size_t &i, std::vector<size_t> &xis, Real *maxVarOut);
+
+        /// Return the InfAlg to use at each step of the recursion. 
+        /// \todo At present, only returns a BP instance
+        InfAlg* getInfAlg();
+
+        size_t _iters;
+        double _maxdiff;
+
+        void setBeliefs(const std::vector<Factor>& bs, double logZ) {
+            size_t i=0;
+            _beliefsV.clear(); _beliefsV.reserve(nrVars());
+            _beliefsF.clear(); _beliefsF.reserve(nrFactors());
+            for(i=0; i<nrVars(); i++) _beliefsV.push_back(bs[i]);
+            for(; i<nrVars()+nrFactors(); i++) _beliefsF.push_back(bs[i]);
+            _logZ = logZ;
+        }
+
+        void construct();
+
+    public:
+
+        //----------------------------------------------------------------
+
+        /// construct CBP object from FactorGraph
+        CBP(const FactorGraph &fg, const PropertySet &opts) : DAIAlgFG(fg) {
+            props.set(opts);
+            construct();
+        }
+
+        /// Name of this inference algorithm
+        static const char *Name;
+
+        /// @name General InfAlg interface
+        //@{
+        virtual CBP* clone() const { return new CBP(*this); }
+        virtual CBP* create() const { DAI_THROW(NOT_IMPLEMENTED); }
+        virtual std::string identify() const { return std::string(Name) + props.toString(); }
+        virtual Factor belief (const Var &n) const { return _beliefsV[findVar(n)]; }
+        virtual Factor belief (const VarSet &) const { DAI_THROW(NOT_IMPLEMENTED); }
+        virtual std::vector<Factor> beliefs() const { return concat(_beliefsV, _beliefsF); }
+        virtual Real logZ() const { return _logZ; }
+        virtual void init() {};
+        virtual void init( const VarSet & ) {};
+        virtual double run();
+        virtual double maxDiff() const { return _maxdiff; }
+        virtual size_t Iterations() const { return _iters; }
+        //@}
+
+        Factor beliefV (size_t i) const { return _beliefsV[i]; }
+        Factor beliefF (size_t I) const { return _beliefsF[I]; }
+
+        //----------------------------------------------------------------
+
+        /// Parameters of this inference algorithm
+        /* PROPERTIES(props,CBP) {
+            /// Enumeration of possible update schedules
+            typedef BP::Properties::UpdateType UpdateType;
+            /// Enumeration of possible methods for deciding when to stop recursing
+            DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
+            /// Enumeration of possible heuristics for choosing clamping variable
+            DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
+            /// Enumeration of possible clampings: variables or factors
+            DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
+            
+            /// Verbosity
+            size_t verbose = 0;
+
+            /// Tolerance to use in BP
+            double tol;
+            /// Update style for BP
+            UpdateType updates;
+            /// Maximum number of iterations for BP
+            size_t maxiter;
+
+            /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF)
+            double rec_tol;
+            /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
+            size_t max_levels = 10;
+            /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
+            double min_max_adj;
+            /// Heuristic for choosing clamping variable
+            ChooseMethodType choose;
+            /// Method for deciding when to stop recursing
+            RecurseType recursion;
+            /// Whether to clamp variables or factors
+            ClampType clamp;
+            /// Properties to pass to BBP
+            PropertySet bbp_props;
+            /// Cost function to use for BBP
+            bbp_cfn_t bbp_cfn;
+            /// Random seed
+            size_t rand_seed = 0;
+
+            /// If non-empty, write clamping choices to this file
+            std::string clamp_outfile = "";
+        }
+        */
 /* {{{ GENERATED CODE: DO NOT EDIT. Created by 
     ./scripts/regenerate-properties include/dai/cbp.h src/cbp.cpp 
 */
-    struct Properties {
-        typedef BP::Properties::UpdateType UpdateType;
-        DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
-        DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
-        DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
-        size_t verbose;
-        /// Tolerance to use in BP
-        double tol;
-        /// Update style for BP
-        UpdateType updates;
-        /// Maximum number of iterations for BP
-        size_t maxiter;
-        /// Tolerance to use for controlling recursion depth (\a recurse
-        /// is REC_LOGZ or REC_BDIFF)
-        double rec_tol;
-        /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
-        size_t max_levels;
-        /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
-        double min_max_adj;
-        /// Heuristic for choosing clamping variable
-        ChooseMethodType choose;
-        /// Method for deciding when to stop recursing
-        RecurseType recursion;
-        /// Whether to clamp variables or factors
-        ClampType clamp;
-        /// Properties to pass to BBP
-        PropertySet bbp_props;
-        /// Cost function to use for BBP
-        bbp_cfn_t bbp_cfn;
-        /// Random seed
-        size_t rand_seed;
-        /// If non-empty, write clamping choices to this file
-        std::string clamp_outfile;
-
-        /// Set members from PropertySet
-        void set(const PropertySet &opts);
-        /// Get members into PropertySet
-        PropertySet get() const;
-        /// Convert to a string which can be parsed as a PropertySet
-        std::string toString() const;
-    } props;
+        struct Properties {
+            /// Enumeration of possible update schedules
+            typedef BP::Properties::UpdateType UpdateType;
+            /// Enumeration of possible methods for deciding when to stop recursing
+            DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
+            /// Enumeration of possible heuristics for choosing clamping variable
+            DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_MAXENT,CHOOSE_BBP,CHOOSE_BP_L1,CHOOSE_BP_CFN);
+            /// Enumeration of possible clampings: variables or factors
+            DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
+            /// Verbosity
+            size_t verbose;
+            /// Tolerance to use in BP
+            double tol;
+            /// Update style for BP
+            UpdateType updates;
+            /// Maximum number of iterations for BP
+            size_t maxiter;
+            /// Tolerance to use for controlling recursion depth (\a recurse is REC_LOGZ or REC_BDIFF)
+            double rec_tol;
+            /// Maximum number of levels of recursion (\a recurse is REC_FIXED)
+            size_t max_levels;
+            /// If choose=CHOOSE_BBP and maximum adjoint is less than this value, don't recurse
+            double min_max_adj;
+            /// Heuristic for choosing clamping variable
+            ChooseMethodType choose;
+            /// Method for deciding when to stop recursing
+            RecurseType recursion;
+            /// Whether to clamp variables or factors
+            ClampType clamp;
+            /// Properties to pass to BBP
+            PropertySet bbp_props;
+            /// Cost function to use for BBP
+            bbp_cfn_t bbp_cfn;
+            /// Random seed
+            size_t rand_seed;
+            /// If non-empty, write clamping choices to this file
+            std::string clamp_outfile;
+
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } props;
 /* }}} END OF GENERATED CODE */
 
+    /// Returns heuristic used for clamping variable
     Properties::ChooseMethodType ChooseMethod() { return props.choose; }
+
+    /// Returns method used for deciding when to stop recursing
     Properties::RecurseType Recursion() { return props.recursion; }
+
+    /// Returns clamping type used
     Properties::ClampType Clamping() { return props.clamp; }
+    /// Returns maximum number of levels of recursion
     size_t maxClampLevel() { return props.max_levels; }
+    /// Returns props.min_max_adj @see CBP::Properties::min_max_adj
     double minMaxAdj() { return props.min_max_adj; }
+    /// Returns tolerance used for controlling recursion depth
     double recTol() { return props.rec_tol; }
 };
 
index 3eca1a6..3bfb51a 100644 (file)
@@ -60,10 +60,10 @@ class LC : public DAIAlgFG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of possible ways to initialize the cavities
-            DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM)
+            DAI_ENUM(CavityType,FULL,PAIR,PAIR2,UNIFORM);
 
             /// Enumeration of different update schedules
-            DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE)
+            DAI_ENUM(UpdateType,SEQFIX,SEQRND,NONE);
 
             /// Verbosity
             size_t verbose;
index 1ff061f..39bf995 100644 (file)
@@ -72,10 +72,10 @@ class MR : public DAIAlgFG {
         /// Parameters of this inference algorithm
         struct Properties {
             /// Enumeration of different types of update equations
-            DAI_ENUM(UpdateType,FULL,LINEAR)
+            DAI_ENUM(UpdateType,FULL,LINEAR);
 
             /// Enumeration of different ways of initializing the cavity correlations
-            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT)
+            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
 
             /// Verbosity
             size_t verbose;
index 565bf4b..1c3b45a 100755 (executable)
@@ -158,9 +158,9 @@ sub process_properties ($$$) {
   
   my ($stext) = "";
   my ($text) = <<EOF;
-    struct Properties {
+        struct Properties {
 EOF
-  my $indent = (' 'x8);
+  my $indent = (' 'x12);
   for my $d (@typedecls) {
     my ($decl,$cmt) = @$d;
     if($cmt ne '') {
@@ -178,13 +178,13 @@ EOF
 
   $text .= <<EOF;
 
-        /// Set members from PropertySet
-        void set(const PropertySet &opts);
-        /// Get members into PropertySet
-        PropertySet get() const;
-        /// Convert to a string which can be parsed as a PropertySet
-        std::string toString() const;
-    } $struct_name;
+            /// Set members from PropertySet
+            void set(const PropertySet &opts);
+            /// Get members into PropertySet
+            PropertySet get() const;
+            /// Convert to a string which can be parsed as a PropertySet
+            std::string toString() const;
+        } $struct_name;
 EOF
 
   $stext .= <<EOF;
index 73dfac7..06750db 100644 (file)
 #include <dai/util.h>
 #include <dai/bipgraph.h>
 
-// for makeBBPPlot: {
-#include <sys/stat.h>
-#include <sys/types.h>
-#include <iostream>
-#include <fstream>
-// }
-
 
 namespace dai {
 
@@ -41,21 +34,20 @@ using namespace std;
 typedef BipartiteGraph::Neighbor Neighbor;
 
 
-Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w) {
-    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;
-    }
+Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
+    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;
 }
 
 
-size_t getFactorEntryForState(const FactorGraph &fg,  size_t I, const vector<size_t>& 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-- ) {
         // note that iterating over nbF(I) yields the same ordering
@@ -68,10 +60,10 @@ 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) {
+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) {
+    switch( (size_t)cfn_type ) {
     case bbp_cfn_t::cfn_bethe_ent: {
         vector<Prob> b1_adj;
         vector<Prob> b2_adj;
@@ -81,31 +73,27 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
         psi1_adj.reserve(fg.nrVars());
         b2_adj.reserve(fg.nrFactors());
         psi2_adj.reserve(fg.nrFactors());
-        for(size_t i=0; i<fg.nrVars(); i++) {
+        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++) {
+            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++) {
+            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++) {
+        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++) {
+            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++) {
+            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);
@@ -114,16 +102,15 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
     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++) {
+        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++) {
+            for( size_t xI = 0; xI < dim; xI++ ) {
                 double bIxI = ia.beliefF(I)[xI];
-                if(bIxI<1.0e-15) {
+                if( bIxI < 1.0e-15 )
                     p[xI] = -1.0e10;
-                } else {
+                else
                     p[xI] = 1+log(bIxI);
-                }
             }
             b2_adj.push_back(p);
         }
@@ -136,13 +123,12 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
         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++) {
+            for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
                 double bixi = ia.beliefV(i)[xi];
-                if(bixi<1.0e-15) {
+                if( bixi < 1.0e-15 )
                     p[xi] = -1.0e10;
-                } else {
+                else
                     p[xi] = 1+log(bixi);
-                }
             }
             b1_adj.push_back(p);
         }
@@ -155,28 +141,27 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
         // cost functions that use Gibbs sample, summing over variable
         // marginals
         vector<size_t> state;
-        if(stateP==NULL) {
+        if(stateP==NULL)
             state = getGibbsState(ia,2*ia.Iterations());
-        } else {
+        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++) {
+        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) {
+            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: abort();
+            default: DAI_THROW(INTERNAL_ERROR);
             }
             b1_adj.push_back(delta);
         }
@@ -189,16 +174,15 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
         // cost functions that use Gibbs sample, summing over factor
         // marginals
         vector<size_t> state;
-        if(stateP==NULL) {
+        if(stateP==NULL)
             state = getGibbsState(ia,2*ia.Iterations());
-        } else {
+        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++) {
+        for( size_t I = 0; I <  fg.nrFactors(); I++ ) {
             size_t n = fg.factor(I).states();
             Prob delta(n,0.0);
 
@@ -206,42 +190,40 @@ void initBBPCostFnAdj(BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vect
             assert(/*0<=x_I &&*/ x_I < n);
 
             double b = ia.beliefF(I)[x_I];
-            switch((size_t)cfn_type) {
+            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: abort();
+            default: DAI_THROW(INTERNAL_ERROR);
             }
             b2_adj.push_back(delta);
         }
         bbp.init(get_zero_adj_V(fg), b2_adj);
         break;
     }
-    default: abort();
+    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;
+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) {
+    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++) {
+        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++) {
+        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:
@@ -249,13 +231,13 @@ Real getCostFn(const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *state
         assert(stateP != NULL);
         vector<size_t> state = *stateP;
         assert(state.size()==fg.nrVars());
-        for(size_t i=0; i<fg.nrVars(); i++) {
+        for( size_t i = 0; i < fg.nrVars(); i++ ) {
             double b = ia.beliefV(i)[state[i]];
-            switch((size_t)cfn_type) {
+            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: abort();
+            default: DAI_THROW(INTERNAL_ERROR);
             }
         }
         break;
@@ -266,42 +248,40 @@ Real getCostFn(const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *state
         assert(stateP != NULL);
         vector<size_t> state = *stateP;
         assert(state.size()==fg.nrVars());
-        for(size_t I=0; I<fg.nrFactors(); I++) {
+        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: abort();
+            default: DAI_THROW(INTERNAL_ERROR);
             }
         }
         break;
     }
     default: 
         cerr << "Unknown cost function " << cfn_type << endl;
-        abort();
+        DAI_THROW(UNKNOWN_ENUM_VALUE);
     }
     return cf;
 }
 
 
-vector<Prob> get_zero_adj_F(const FactorGraph& fg) {
+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)));
-    }
+    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> 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)));
-    }
+    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;
 }
 
@@ -351,16 +331,14 @@ void BBP::RegenerateInds() {
 void BBP::RegenerateT() {
     // _T[i][_I]
     _T.clear();
-    _T.resize(_fg->nrVars());
+    _T.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _T[i].resize(_fg->nbV(i).size());
-        foreach(const Neighbor &I, _fg->nbV(i)) {
-            Prob prod(_fg->var(i).states(),1.0);
-            foreach(const Neighbor &J, _fg->nbV(i)) {
-                if(size_t(J)!=size_t(I)) {
-                    prod *= _bp_dual.msgM(i,J.iter);
-                }
-            }
+        _T[i].resize( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            Prob prod( _fg->var(i).states(), 1.0 );
+            foreach( const Neighbor &J, _fg->nbV(i) )
+                if( J.node != I.node )
+                    prod *= _bp_dual.msgM( i, J.iter );
             _T[i][I.iter] = prod;
         }
     }
@@ -370,20 +348,19 @@ void BBP::RegenerateT() {
 void BBP::RegenerateU() {
     // _U[I][_i]
     _U.clear();
-    _U.resize(_fg->nrFactors());
+    _U.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        _U[I].resize(_fg->nbF(I).size());
-        foreach(const Neighbor &i, _fg->nbF(I)) {
-            Prob prod(_fg->factor(I).states(), 1.0);
-            foreach(const Neighbor &j, _fg->nbF(I)) {
-                if(size_t(i) != size_t(j)) {
-                    Prob n_jI(_bp_dual.msgN(j, j.dual));
-                    const _ind_t& ind = _index(j, j.dual);
+        _U[I].resize( _fg->nbF(I).size() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            Prob prod( _fg->factor(I).states(), 1.0 );
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i.node != j.node ) {
+                    Prob n_jI( _bp_dual.msgN( j, j.dual ) );
+                    const _ind_t &ind = _index( j, j.dual );
                     // multiply prod by n_jI
-                    for(size_t x_I = 0; x_I < prod.size(); x_I++)
+                    for( size_t x_I = 0; x_I < prod.size(); x_I++ )
                         prod[x_I] *= n_jI[ind[x_I]];
                 }
-            }
             _U[I][i.iter] = prod;
         }
     }
@@ -393,29 +370,28 @@ void BBP::RegenerateU() {
 void BBP::RegenerateS() {
     // _S[i][_I][_j]
     _S.clear();
-    _S.resize(_fg->nrVars());
+    _S.resize( _fg->nrVars() );
     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        _S[i].resize(_fg->nbV(i).size());
-        foreach(const Neighbor& I, _fg->nbV(i)) {
-            _S[i][I.iter].resize(_fg->nbF(I).size());
-            foreach(const Neighbor& j, _fg->nbF(I)) {
-                if(i != j) {
-                    Factor prod(_fg->factor(I));
-                    foreach(const Neighbor& k, _fg->nbF(I)) {
-                        if(k!=i && size_t(k)!=size_t(j)) {
-                            const _ind_t& ind = _index(k,k.dual);
-                            Prob p(_bp_dual.msgN(k,k.dual));
+        _S[i].resize( _fg->nbV(i).size() );
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            _S[i][I.iter].resize( _fg->nbF(I).size() );
+            foreach( const Neighbor &j, _fg->nbF(I) )
+                if( i != j ) {
+                    Factor prod( _fg->factor(I) );
+                    foreach( const Neighbor &k, _fg->nbF(I) ) {
+                        if( k != i && k.node != j.node ) {
+                            const _ind_t &ind = _index( k, k.dual );
+                            Prob p( _bp_dual.msgN( k, k.dual ) );
                             for( size_t x_I = 0; x_I < prod.states(); x_I++ )
                                 prod.p()[x_I] *= p[ind[x_I]];
                         }
                     }
                     // "Marginalize" onto i|j (unnormalized)
-                    // XXX improve efficiency? 
+                    // XXX improve efficiency?
                     Prob marg;
-                    marg = prod.marginal(VarSet(_fg->var(i)) | VarSet(_fg->var(j)), false).p();
+                    marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p();
                     _S[i][I.iter][j.iter] = marg;
                 }
-            }
         }
     }
 }
@@ -424,20 +400,17 @@ void BBP::RegenerateS() {
 void BBP::RegenerateR() {
     // _R[I][_i][_J]
     _R.clear();
-    _R.resize(_fg->nrFactors());
+    _R.resize( _fg->nrFactors() );
     for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
-        _R[I].resize(_fg->nbF(I).size());
-        foreach(const Neighbor& i, _fg->nbF(I)) {
-            _R[I][i.iter].resize(_fg->nbV(i).size());
-            foreach(const Neighbor& J, _fg->nbV(i)) {
-                if(I != J) {
-                    Prob prod(_fg->var(i).states(), 1.0);
-                    foreach(const Neighbor& K, _fg->nbV(i)) {
-                        if(size_t(K) != size_t(I) &&
-                           size_t(K) != size_t(J)) {
+        _R[I].resize( _fg->nbF(I).size() );
+        foreach( const Neighbor &i, _fg->nbF(I) ) {
+            _R[I][i.iter].resize( _fg->nbV(i).size() );
+            foreach( const Neighbor &J, _fg->nbV(i) ) {
+                if( I != J ) {
+                    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);
-                        }
-                    }
                     _R[I][i.iter][J.iter] = prod;
                 }
             }
@@ -551,21 +524,23 @@ void BBP::RegenerateParMessageAdjoints() {
 }
 
 
-void BBP::incrSeqMsgM(size_t i, size_t _I, const Prob& p) {
-    if(props.clean_updates) {
+void BBP::incrSeqMsgM( size_t i, size_t _I, const Prob &p ) {
+    if( props.clean_updates )
         _new_adj_m[i][_I] += p;
-    else {
+    else {
         _adj_m[i][_I] += p;
         calcUnnormMsgM(i, _I);
     }
 }
 
 
+#if 0
 Real pv_thresh=1000;
+#endif
 
 
-void BBP::updateSeqMsgM(size_t i, size_t _I) {
-    if(props.clean_updates) {
+void BBP::updateSeqMsgM( size_t i, size_t _I ) {
+    if( props.clean_updates ) {
 #if 0
         if(_new_adj_m[i][_I].sumAbs() > pv_thresh ||
            _adj_m[i][_I].sumAbs() > pv_thresh) {
@@ -578,22 +553,22 @@ void BBP::updateSeqMsgM(size_t i, size_t _I) {
         }
 #endif
         _adj_m[i][_I] += _new_adj_m[i][_I];
-        calcUnnormMsgM(i, _I);
-        _new_adj_m[i][_I].fill(0.0);
+        calcUnnormMsgM( i, _I );
+        _new_adj_m[i][_I].fill( 0.0 );
     }
 }
 
 
-void BBP::setSeqMsgM(size_t i, size_t _I, const Prob& p) {
+void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) {
     _adj_m[i][_I] = p;
     calcUnnormMsgM(i, _I);
 }
 
 
-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 NeighborI = _fg->nbV(i)[_I];
-    assert(I.iter == _I);
+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);
 #if 0
     if(f_unnorm.sumAbs() > pv_thresh) {
@@ -608,8 +583,8 @@ void BBP::sendSeqMsgN(size_t i, size_t _I, const Prob &f) {
         DAI_PV(_fg->factor(I).p());
     }
 #endif
-    foreach(const Neighbor& J, _fg->nbV(i)) {
-        if(size_t(J) != size_t(I)) {
+    foreach( const Neighbor &J, _fg->nbV(i) ) {
+        if( J.node != I.node ) {
 #if 0
             if(f_unnorm.sumAbs() > pv_thresh) {
                 DAI_DMSG("in sendSeqMsgN loop");
@@ -619,16 +594,14 @@ void BBP::sendSeqMsgN(size_t i, size_t _I, const Prob &f) {
                 DAI_PV(f_unnorm * _R[J][J.dual][_I]);
             }
 #endif
-            incrSeqMsgM(i, J.iter, f_unnorm * R(J,J.dual,_I));
+            incrSeqMsgM( i, J.iter, f_unnorm * R(J, J.dual, _I) );
         }
     }
 }
 
 
-void BBP::sendSeqMsgM(size_t j, size_t _I) {
-    const Neighbor& I = _fg->nbV(j)[_I];
-    size_t _j = I.dual;
-    const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][_I];
+void BBP::sendSeqMsgM( size_t j, size_t _I ) {
+    const Neighbor &I = _fg->nbV(j)[_I];
 
 //     DAI_PV(j);
 //     DAI_PV(I);
@@ -636,26 +609,37 @@ void BBP::sendSeqMsgM(size_t j, size_t _I) {
 //     DAI_PV(_adj_m[j][_I]);
 //     DAI_PV(_bp_dual.zM(j,_I));
 
-    Prob um(U(I,_j));
-    const _ind_t& ind = _index(j, _I);
+    size_t _j = I.dual;
+    const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][_I];
+    Prob um( U(I, _j) );
+    const _ind_t &ind = _index(j, _I);
     for( size_t x_I = 0; x_I < um.size(); x_I++ )
         um[x_I] *= _adj_m_unnorm_jI[ind[x_I]];
-    um *= 1-props.damping;
+    um *= 1 - props.damping;
     _adj_psi_F[I] += um;
+    
+    /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
+    _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, _j) ) * Factor( _fg->var(j), _adj_m_unnorm[j][_I] )).p() * (1.0 - props.damping);
+    */
 
 //     DAI_DMSG("in sendSeqMsgM");
 //     DAI_PV(j);
 //     DAI_PV(I);
 //     DAI_PV(_I);
 //     DAI_PV(_fg->nbF(I).size());
-    foreach(const Neighbor& i, _fg->nbF(I)) {
-        if(size_t(i) != j) {
+    foreach( const Neighbor &i, _fg->nbF(I) ) {
+        if( i.node != j ) {
             const Prob &S = _S[i][i.dual][_j];
-            Prob msg(_fg->var(i).states(),0.0);
+            Prob msg( _fg->var(i).states(), 0.0 );
             LOOP_ij(
-                msg[xi] += S[xij]*_adj_m_unnorm_jI[xj];
+                msg[xi] += S[xij] * _adj_m_unnorm_jI[xj];
             );
-            msg *= 1-props.damping;
+            msg *= 1.0 - props.damping;
+            /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
+            Var vi = _fg->var(i);
+            Var vj = _fg->var(j);
+            msg = (Factor(VarSet(vi,vj), S) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p() * (1.0 - props.damping);
+            */
 #if 0
             if(msg.sumAbs() > pv_thresh) {
                 DAI_DMSG("in sendSeqMsgM loop");
@@ -673,12 +657,11 @@ void BBP::sendSeqMsgM(size_t j, size_t _I) {
                 DAI_PV(_fg->nbV(i).size());
             }
 #endif
-            assert(size_t(_fg->nbV(i)[i.dual]) == I);
-            sendSeqMsgN(i, i.dual, msg);
+            assert( _fg->nbV(i)[i.dual].node == I );
+            sendSeqMsgN( i, i.dual, msg );
         }
     }
-//     assert(dist(_adj_m_unnorm_jI, _adj_m_unnorm[j][_I],Prob::DISTL1)<=1.0e-5);
-    setSeqMsgM(j, _I, _adj_m[j][_I]*props.damping);
+    setSeqMsgM( j, _I, _adj_m[j][_I] * props.damping );
 }
 
 
@@ -738,115 +721,117 @@ void BBP::Regenerate() {
     RegenerateR();
     RegenerateInputs();
     RegeneratePsiAdjoints();
-    if(props.updates == Properties::UpdateType::PAR) {
+    if( props.updates == Properties::UpdateType::PAR )
         RegenerateParMessageAdjoints();
-    } else {
+    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];
+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);
+    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) {
+    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];
+                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) {
+void BBP::calcNewM( size_t i, size_t _I ) {
     const Neighbor &I = _fg->nbV(i)[_I];
-    Prob p(U(I,I.dual));
+    Prob p( U(I, I.dual) );
     const Prob &adj = _adj_m_unnorm[i][_I];
-    const _ind_tind = _index(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];
-        }
-    }
+    _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::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::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::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::upMsgM(size_t i, size_t _I) {
-    _adj_m[i][_I] = _new_adj_m[i][_I];
-    calcUnnormMsgM(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::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) ) {
+            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);
+    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;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor& I, _fg->nbV(i)) {
+    Real s = 0.0;
+    size_t e = 0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
             s += _adj_m_unnorm[i][I.iter].sumAbs();
             s += _adj_n_unnorm[i][I.iter].sumAbs();
             e++;
         }
-    }
-    return s/e;
+    return s / e;
 }
 
 
-void BBP::getMsgMags(Real &s, Real &new_s) {
-    s=0.0; new_s=0.0;
-    size_t e=0;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor& I, _fg->nbV(i)) {
+void BBP::getMsgMags( Real &s, Real &new_s ) {
+    s = 0.0; 
+    new_s = 0.0;
+    size_t e = 0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
             s += _adj_m[i][I.iter].sumAbs();
             s += _adj_n[i][I.iter].sumAbs();
             new_s += _new_adj_m[i][I.iter].sumAbs();
             new_s += _new_adj_n[i][I.iter].sumAbs();
             e++;
         }
-    }
-    s /= e; new_s /= e;
+    s /= e; 
+    new_s /= e;
 }
 
 // tuple<size_t,size_t,Real> BBP::getArgMaxPsi1Adj() {
@@ -867,64 +852,57 @@ void BBP::getMsgMags(Real &s, Real &new_s) {
 // }
 
 
+void BBP::getArgmaxMsgM( size_t &out_i, size_t &out__I, Real &mag ) {
+    bool found = false;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) ) {
+            Real thisMag = _adj_m[i][I.iter].sumAbs();
+            if( !found || mag < thisMag ) {
+                found = true;
+                mag = thisMag;
+                out_i = i;
+                out__I = I.iter;
+            }
+        }
+    assert( found );
+}
+
+
 Real BBP::getMaxMsgM() {
     size_t dummy;
     Real mag;
-    getArgmaxMsgM(dummy, dummy, mag);
+    getArgmaxMsgM( dummy, dummy, mag );
     return mag;
 }
 
 
 Real BBP::getTotalMsgM() {
-    Real mag=0.0;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor &I, _fg->nbV(i)) {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
             mag += _adj_m[i][I.iter].sumAbs();
-        }
-    }
     return mag;
 }
 
 
 Real BBP::getTotalNewMsgM() {
-    Real mag=0.0;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor &I, _fg->nbV(i)) {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
             mag += _new_adj_m[i][I.iter].sumAbs();
-        }
-    }
     return mag;
 }
 
 
 Real BBP::getTotalMsgN() {
-    Real mag=0.0;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor &I, _fg->nbV(i)) {
+    Real mag = 0.0;
+    for( size_t i = 0; i < _fg->nrVars(); i++ )
+        foreach( const Neighbor &I, _fg->nbV(i) )
             mag += _adj_n[i][I.iter].sumAbs();
-        }
-    }
     return mag;
 }
 
 
-void BBP::getArgmaxMsgM(size_t &out_i, size_t &out__I, Real &mag) {
-    bool found=false;
-    for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-        foreach(const Neighbor &I, _fg->nbV(i)) {
-            Real thisMag = _adj_m[i][I.iter].sumAbs();
-            if(!found || mag < thisMag) {
-                found = true;
-                mag = thisMag;
-                out_i = i;
-                out__I = I.iter;
-            }
-        }
-    }
-    assert(found);
-}
-
-
 void BBP::run() {
     typedef BBP::Properties::UpdateType UT;
     Real tol = props.tol;
@@ -954,16 +932,12 @@ void BBP::run() {
             mag = getTotalMsgM();
             if(mag < tol) break;
 
-            for( size_t i = 0; i < _fg->nrVars(); i++ ) {
-                foreach(const Neighbor &I, _fg->nbV(i)) {
+            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)) {
+            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) {
@@ -1214,8 +1188,8 @@ void BBP::Properties::set(const PropertySet &opts)
     bool die=false;
     for(i=keys.begin(); i!=keys.end(); i++) {
         if(*i == "verbose") continue;
-        if(*i == "tol") continue;
         if(*i == "maxiter") continue;
+        if(*i == "tol") continue;
         if(*i == "damping") continue;
         if(*i == "updates") continue;
         if(*i == "clean_updates") continue;
@@ -1229,14 +1203,14 @@ void BBP::Properties::set(const PropertySet &opts)
         cerr << "BBP: Missing property \"verbose\" for method \"BBP\"" << endl;
         die=true;
     }
-    if(!opts.hasKey("tol")) {
-        cerr << "BBP: Missing property \"tol\" for method \"BBP\"" << endl;
-        die=true;
-    }
     if(!opts.hasKey("maxiter")) {
         cerr << "BBP: Missing property \"maxiter\" for method \"BBP\"" << endl;
         die=true;
     }
+    if(!opts.hasKey("tol")) {
+        cerr << "BBP: Missing property \"tol\" for method \"BBP\"" << endl;
+        die=true;
+    }
     if(!opts.hasKey("damping")) {
         cerr << "BBP: Missing property \"damping\" for method \"BBP\"" << endl;
         die=true;
@@ -1253,8 +1227,8 @@ void BBP::Properties::set(const PropertySet &opts)
         DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
     }
     verbose = opts.getStringAs<size_t>("verbose");
-    tol = opts.getStringAs<double>("tol");
     maxiter = opts.getStringAs<size_t>("maxiter");
+    tol = opts.getStringAs<double>("tol");
     damping = opts.getStringAs<double>("damping");
     updates = opts.getStringAs<UpdateType>("updates");
     clean_updates = opts.getStringAs<bool>("clean_updates");
@@ -1262,8 +1236,8 @@ void BBP::Properties::set(const PropertySet &opts)
 PropertySet BBP::Properties::get() const {
     PropertySet opts;
     opts.Set("verbose", verbose);
-    opts.Set("tol", tol);
     opts.Set("maxiter", maxiter);
+    opts.Set("tol", tol);
     opts.Set("damping", damping);
     opts.Set("updates", updates);
     opts.Set("clean_updates", clean_updates);
@@ -1273,8 +1247,8 @@ string BBP::Properties::toString() const {
     stringstream s(stringstream::out);
     s << "[";
     s << "verbose=" << verbose << ",";
-    s << "tol=" << tol << ",";
     s << "maxiter=" << maxiter << ",";
+    s << "tol=" << tol << ",";
     s << "damping=" << damping << ",";
     s << "updates=" << updates << ",";
     s << "clean_updates=" << clean_updates;
index ef871e5..dfa4802 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=5000,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 e3ada79..7c76ea9 100755 (executable)
@@ -1,2 +1,2 @@
 #!/bin/bash
-./testdai --report-iters false --report-time false --marginals true --aliases aliases.conf --filename $1 --methods EXACT[verbose=0] JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 GBP_LOOP5 GBP_LOOP6 GBP_LOOP7 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 HAK_LOOP5 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_PAIR2CAV_SEQFIX LCBP_PAIR2CAVin_SEQFIX LCBP_PAIR2CAV_SEQRND LCBP_PAIR2CAVin_SEQRND LCBP_PAIR2CAV_NONE LCBP_PAIR2CAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND
+./testdai --report-iters false --report-time false --marginals true --aliases aliases.conf --filename $1 --methods EXACT[verbose=0] JTREE_HUGIN JTREE_SHSH BP_SEQFIX BP_SEQRND BP_SEQMAX BP_PARALL BP_SEQFIX_LOG BP_SEQRND_LOG BP_SEQMAX_LOG BP_PARALL_LOG MF_SEQRND TREEEP TREEEPWC GBP_MIN GBP_DELTA GBP_LOOP3 GBP_LOOP4 GBP_LOOP5 GBP_LOOP6 GBP_LOOP7 HAK_MIN HAK_DELTA HAK_LOOP3 HAK_LOOP4 HAK_LOOP5 MR_RESPPROP_FULL MR_CLAMPING_FULL MR_EXACT_FULL MR_RESPPROP_LINEAR MR_CLAMPING_LINEAR MR_EXACT_LINEAR LCBP_FULLCAV_SEQFIX LCBP_FULLCAVin_SEQFIX LCBP_FULLCAV_SEQRND LCBP_FULLCAVin_SEQRND LCBP_FULLCAV_NONE LCBP_FULLCAVin_NONE LCBP_PAIRCAV_SEQFIX LCBP_PAIRCAVin_SEQFIX LCBP_PAIRCAV_SEQRND LCBP_PAIRCAVin_SEQRND LCBP_PAIRCAV_NONE LCBP_PAIRCAVin_NONE LCBP_PAIR2CAV_SEQFIX LCBP_PAIR2CAVin_SEQFIX LCBP_PAIR2CAV_SEQRND LCBP_PAIR2CAVin_SEQRND LCBP_PAIR2CAV_NONE LCBP_PAIR2CAVin_NONE LCBP_UNICAV_SEQFIX LCBP_UNICAV_SEQRND BBP CBP
diff --git a/tests/testbbp.cpp b/tests/testbbp.cpp
new file mode 100644 (file)
index 0000000..502d869
--- /dev/null
@@ -0,0 +1,123 @@
+/*  Copyright (C) 2009  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
+    Max Planck Institute for Biological Cybernetics, Germany
+
+    This file is part of libDAI.
+
+    libDAI is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    libDAI is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with libDAI; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+*/
+
+
+#include <iostream>
+#include <dai/alldai.h>
+#include <dai/bbp.h>
+
+
+using namespace dai;
+using namespace std;
+
+
+int main( int argc, char *argv[] ) {
+    if ( argc != 2 ) {
+        cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
+        cout << "Reads factor graph <filename.fg> and verifies" << endl;
+        cout << "whether BBP works correctly on it." << endl << endl;
+        return 1;
+    } else {
+        // Read FactorGraph from the file specified by the first command line argument
+        FactorGraph fg;
+        fg.ReadFromFile(argv[1]);
+        
+        // Set some constants
+        size_t verbose = 0;
+        double tol = 1.0e-9;
+        size_t maxiter = 10000;
+        double damping = 0.0;
+        BBP::Properties::UpdateType updates = BBP::Properties::UpdateType::PAR;
+        bool   clean_updates = false;
+
+        // Store the constants in a PropertySet object
+        PropertySet opts;
+        opts.Set("verbose",verbose);  // Verbosity (amount of output generated)
+        opts.Set("tol",tol);          // Tolerance for convergence
+        opts.Set("maxiter",maxiter);  // Maximum number of iterations
+        opts.Set("damping",damping);  // Amount of damping applied
+
+        // Construct a BP (belief propagation) object from the FactorGraph fg
+        BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false));
+        bp.recordSentMessages = true;
+        bp.init();
+        bp.run();
+
+        vector<size_t> state( fg.nrVars(), 0 );
+
+        for( size_t t = 0; t < 90; t++ ) {
+            clean_updates = t % 2;
+            BBP::Properties::UpdateType updates;
+            switch( (t / 2) % 5 ) {
+                case BBP::Properties::UpdateType::SEQ_FIX:
+                    updates = BBP::Properties::UpdateType::SEQ_FIX;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_MAX:
+                    updates = BBP::Properties::UpdateType::SEQ_MAX;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_BP_REV:
+                    updates = BBP::Properties::UpdateType::SEQ_BP_REV;
+                    break;
+                case BBP::Properties::UpdateType::SEQ_BP_FWD:
+                    updates = BBP::Properties::UpdateType::SEQ_BP_FWD;
+                    break;
+                case BBP::Properties::UpdateType::PAR:
+                    updates = BBP::Properties::UpdateType::PAR;
+                    break;
+            }
+            bbp_cfn_t cfn;
+            switch( (t / 10) % 9 ) {
+                case 0:
+                    cfn = bbp_cfn_t::cfn_gibbs_b;
+                    break;
+                case 1:
+                    cfn = bbp_cfn_t::cfn_gibbs_b2;
+                    break;
+                case 2:
+                    cfn = bbp_cfn_t::cfn_gibbs_exp;
+                    break;
+                case 3:
+                    cfn = bbp_cfn_t::cfn_gibbs_b_factor;
+                    break;
+                case 4:
+                    cfn = bbp_cfn_t::cfn_gibbs_b2_factor;
+                    break;
+                case 5:
+                    cfn = bbp_cfn_t::cfn_gibbs_exp_factor;
+                    break;
+                case 6:
+                    cfn = bbp_cfn_t::cfn_var_ent;
+                    break;
+                case 7:
+                    cfn = bbp_cfn_t::cfn_factor_ent;
+                    break;
+                case 8:
+                    cfn = bbp_cfn_t::cfn_bethe_ent;
+                    break;
+            }
+
+            double h = 1e-4;
+            double result = numericBBPTest( bp, &state, opts("updates",updates)("clean_updates",clean_updates), cfn, h );
+            cout << "clean_updates=" << clean_updates << ", updates=" << updates << ", cfn=" << cfn << ", result: " << result << endl;
+        }
+    }
+
+    return 0;
+}
index eb632d8..61e0e54 100644 (file)
@@ -884,3 +884,37 @@ LCBP_UNICAV_SEQRND                         9.483e-02       3.078e-02       N/A             1.000e-09
 # ({x13}, (5.266e-01, 4.734e-01))
 # ({x14}, (6.033e-01, 3.967e-01))
 # ({x15}, (1.558e-01, 8.442e-01))
+BBP                                            7.049e-04       2.319e-04       +7.075e-05      1.000e-09       
+# ({x0}, (3.890e-01, 6.110e-01))
+# ({x1}, (5.555e-01, 4.445e-01))
+# ({x2}, (4.587e-01, 5.413e-01))
+# ({x3}, (5.479e-01, 4.521e-01))
+# ({x4}, (6.663e-01, 3.337e-01))
+# ({x5}, (2.105e-01, 7.895e-01))
+# ({x6}, (8.180e-01, 1.820e-01))
+# ({x7}, (2.324e-01, 7.676e-01))
+# ({x8}, (2.169e-01, 7.831e-01))
+# ({x9}, (2.050e-01, 7.950e-01))
+# ({x10}, (7.666e-01, 2.334e-01))
+# ({x11}, (1.216e-01, 8.784e-01))
+# ({x12}, (4.210e-01, 5.790e-01))
+# ({x13}, (5.351e-01, 4.649e-01))
+# ({x14}, (6.284e-01, 3.716e-01))
+# ({x15}, (1.355e-01, 8.645e-01))
+CBP                                            2.160e-05       1.082e-05       +1.720e-06      1.000e-09       
+# ({x0}, (3.888e-01, 6.112e-01))
+# ({x1}, (5.556e-01, 4.444e-01))
+# ({x2}, (4.587e-01, 5.413e-01))
+# ({x3}, (5.480e-01, 4.520e-01))
+# ({x4}, (6.660e-01, 3.340e-01))
+# ({x5}, (2.107e-01, 7.893e-01))
+# ({x6}, (8.178e-01, 1.822e-01))
+# ({x7}, (2.327e-01, 7.673e-01))
+# ({x8}, (2.171e-01, 7.829e-01))
+# ({x9}, (2.052e-01, 7.948e-01))
+# ({x10}, (7.664e-01, 2.336e-01))
+# ({x11}, (1.217e-01, 8.783e-01))
+# ({x12}, (4.214e-01, 5.786e-01))
+# ({x13}, (5.348e-01, 4.652e-01))
+# ({x14}, (6.291e-01, 3.709e-01))
+# ({x15}, (1.357e-01, 8.643e-01))