[Frederik Eaton] Major cleanup of BBP and CBP code and documentation
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 1 Jun 2009 10:52:42 +0000 (12:52 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 1 Jun 2009 10:52:42 +0000 (12:52 +0200)
- BBP: change to use Neighbor interface
- tests/aliases.conf: updated with new CBP properties
- Add doxygen documentation to headers

include/dai/bbp.h
include/dai/cbp.h
include/dai/doc.h
src/bbp.cpp
src/cbp.cpp
tests/aliases.conf

index b14b6a3..4ba6498 100644 (file)
@@ -1,12 +1,37 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    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
+*/
+
+
+/// \file
+/// \brief Defines class BBP [\ref EaG09]
+/// \todo Improve documentation
+/// \todo Clean up
+
+
 #ifndef ___defined_libdai_bbp_h
 #define ___defined_libdai_bbp_h
 
+
 #include <vector>
 #include <utility>
 #include <ext/hash_map>
 
-#include <boost/tuple/tuple.hpp>
-
 #include <dai/prob.h>
 #include <dai/daialg.h>
 #include <dai/factorgraph.h>
 
 #include <dai/bp_dual.h>
 
+
 namespace dai {
 
-using namespace std;
-using namespace __gnu_cxx;
-using boost::tuple;
-using boost::get;
 
-// utility function to 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; }
+std::vector<Prob> get_zero_adj_F(const FactorGraph&);
+std::vector<Prob> get_zero_adj_V(const FactorGraph&);
 
-/// function to compute adj_w_unnorm from w, Z_w, adj_w
-Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w);
 
-vector<size_t> getGibbsState(const BP_dual& fg, size_t iters);
+/// Implements BBP (Back-belief-propagation) [\ref EaG09]
+class BBP {
+protected:
+    // ----------------------------------------------------------------
+    // inputs
+    BP_dual _bp_dual;
+    const FactorGraph* _fg;
+    const InfAlg *_ia;
 
-vector<Prob> get_zero_adj_2(const BP_dual& bp_dual);
-vector<Prob> get_zero_adj_1(const BP_dual& bp_dual);
+    // iterations
+    size_t _iters;
 
-class BBP {
-  // ----------------------------------------------------------------
-  // inputs
-  const BP_dual* _bp_dual;
-  const FactorGraph* _fg;
-  vector<Prob> _adj_b_1, _adj_b_2;
-  vector<Prob> _adj_b_1_unnorm, _adj_b_2_unnorm;
-  // in case the objective function depends on factors as well:
-  vector<Prob> _init_adj_psi_1;
-  vector<Prob> _init_adj_psi_2;
-  // ----------------------------------------------------------------
-  // ouptuts
-  vector<Prob> _adj_psi_1, _adj_psi_2;
-  // the following vectors are length _fg->nr_edges()
-  vector<Prob> _adj_n, _adj_m; 
-  vector<Prob> _adj_n_unnorm, _adj_m_unnorm;
-  vector<Prob> _new_adj_n, _new_adj_m;
-  // ----------------------------------------------------------------
-  // auxiliary data
-  typedef vector<hash_map<size_t, vector<size_t> > > _trip_t;
-  typedef vector<size_t> _trip_end_t;
-
-  // members to store indices of VFV and FVF triples
-   // XXX need to have vector going in other direction so that we can
-   // do faster looping on iIj and IiJ
-  _trip_t _VFV;
-  _trip_t _FVF;
-
-  typedef vector<tuple<size_t,size_t,size_t> > _trip_ind_t;
-  _trip_ind_t _VFV_ind;
-  _trip_ind_t _FVF_ind;
-  
-  // helper quantities computed from the BP messages
-  vector<Prob> _T,_U,_S,_R;
-  size_t _iters;
-  size_t _max_iter;
-  
- public:
-  void RegenerateInds();
-  void RegenerateT();
-  void RegenerateU();
-  void RegenerateS();
-  void RegenerateR();
-  void RegenerateInputs();
-  /// initialise members for messages and factor adjoints
-  void RegeneratePsiAdjoints();
-  void RegenerateMessageAdjoints();
-  void Regenerate();
-
-  size_t VFV(size_t i, size_t I, size_t j) const { return (*const_cast<_trip_t*>(&_VFV))[i][I][j]; }
-  size_t FVF(size_t I, size_t i, size_t J) const { return (*const_cast<_trip_t*>(&_FVF))[I][i][J]; }
-  DAI_ACCMUT(Prob & T(size_t i, size_t I), { return _T[_fg->VV2E(i,I)]; });
-  DAI_ACCMUT(Prob & U(size_t I, size_t i), { return _U[_fg->VV2E(i,I)]; });
-  DAI_ACCMUT(Prob & S(size_t i, size_t I, size_t j), { return _S[VFV(i,I,j)]; });
-  DAI_ACCMUT(Prob & R(size_t I, size_t i, size_t J), { return _R[FVF(I,i,J)]; });
-
-  DAI_ACCMUT(Prob & T(size_t iI), { return _T[iI]; });
-  DAI_ACCMUT(Prob & U(size_t iI), { return _U[iI]; });
-  DAI_ACCMUT(Prob & S(size_t iIj), { return _S[iIj]; });
-  DAI_ACCMUT(Prob & R(size_t IiJ), { return _R[IiJ]; });
+    // ----------------------------------------------------------------
+    // 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
   
-  void calcNewN(size_t iI);
-  void calcNewM(size_t iI);
-  void upMsgM(size_t iI);
-  void upMsgN(size_t iI);
-  void doParUpdate();
-  Real getUnMsgMag();
-  void getMsgMags(Real &s, Real &new_s);
-
-  void zero_adj_b_2() {
-    _adj_b_2.clear();
-    _adj_b_2.reserve(_bp_dual->nrFactors());
-    for(size_t I=0; I<_bp_dual->nrFactors(); I++) {
-      _adj_b_2.push_back(Prob(_bp_dual->factor(I).states(),Real(0.0)));
+    /// 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)));
+        }
     }
-  }
- public:
-  BBP(const BP_dual &bp_dual) :
-    _bp_dual(&bp_dual), _fg(&bp_dual), _max_iter(50) {}
-
-  DAI_ACCMUT(size_t& maxIter(), { return _max_iter; });
-
-  void init(const vector<Prob> &adj_b_1, const vector<Prob> &adj_b_2,
-            const vector<Prob> &adj_psi_1, const vector<Prob> &adj_psi_2) {
-    _adj_b_1 = adj_b_1;
-    _adj_b_2 = adj_b_2;
-    _init_adj_psi_1 = adj_psi_1;
-    _init_adj_psi_2 = adj_psi_2;
-    Regenerate(); 
-  }
-  void init(const vector<Prob> &adj_b_1, const vector<Prob> &adj_b_2) {
-    init(adj_b_1, adj_b_2, get_zero_adj_1(*_bp_dual), get_zero_adj_2(*_bp_dual));
-  }
-  void init(const vector<Prob> &adj_b_1) {
-    init(adj_b_1, get_zero_adj_2(*_bp_dual));
-  }
-
-  /// run until change is less than given tolerance
-  void run(Real tol);
-
-  size_t doneIters() { return _iters; }
-
-  DAI_ACCMUT(Prob& adj_psi_1(size_t i), { return _adj_psi_1[i]; });
-  DAI_ACCMUT(Prob& adj_psi_2(size_t I), { return _adj_psi_2[I]; });
-  DAI_ACCMUT(Prob& adj_b_1(size_t i), { return _adj_b_1[i]; });
-  DAI_ACCMUT(Prob& adj_b_2(size_t I), { return _adj_b_2[I]; });
-  DAI_ACCMUT(Prob& adj_n(size_t i, size_t I), { return _adj_n[_fg->VV2E(i,I)]; });
-  DAI_ACCMUT(Prob& adj_m(size_t I, size_t i), { return _adj_m[_fg->VV2E(i,I)]; });
+
+    //----------------------------------------------------------------
+    // 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: 
+
+/* 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;
+    }
+*/
+/* {{{ 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;
+/* }}} END OF GENERATED CODE */
 };
 
-/// Cost functions. Not used by BBP class, only used by following functions
-DAI_ENUM(bbp_cfn_t, cfn_b, cfn_b2, cfn_exp, cfn_var_ent, cfn_factor_ent, cfn_bethe_ent);
+/// 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);
+
+/// Initialise BBP using InfAlg, cost function, and stateP
+/** Calls bbp.init with adjoints calculated from ia.beliefV and
+ *  ia.beliefF. stateP is a Gibbs state and can be NULL, it will be
+ *  initialised using a Gibbs run of 2*fg.Iterations() iterations.
+ */
+void initBBPCostFnAdj(BBP& bbp, const InfAlg& ia, bbp_cfn_t cfn_type, const std::vector<size_t>* stateP);
+
+/// Answers question: Does 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);
 
-/// initialise BBP using BP_dual and (cfn_type, stateP)
-void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual& fg, bbp_cfn_t cfn_type, const 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
+
+/// 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);
 
-/// calculate actual value of cost function (cfn_type, stateP)
-Real gibbsCostFn(const BP_dual& fg, bbp_cfn_t cfn_type, const vector<size_t> *stateP);
+/// 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
-double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, double h);
 
-/// output to "./bbp-data/" a set of files which evaluate the effect
-/// on the given cost function (for a random Gibbs state) of
-/// perturbing each graph variable
-void makeBBPGraph(const BP_dual& bp_dual, bbp_cfn_t cfn);
+} // end of namespace dai
 
-}
 
 #endif
index b9a0f17..81769d3 100644 (file)
@@ -1,35 +1,65 @@
-#ifndef ____defined_libdai_cbp_h__
-#define ____defined_libdai_cbp_h__
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    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
+*/
+
+
+/// \file
+/// \brief Defines class CBP [\ref EaG09]
+/// \todo Improve documentation
+/// \todo Clean up
+
+
+#ifndef __defined_libdai_cbp_h
+#define __defined_libdai_cbp_h
+
+
+#include <fstream>
+
+#include <boost/shared_ptr.hpp>
 
 #include <dai/daialg.h>
 
 #include <dai/cbp.h>
 #include <dai/bbp.h>
 
-using namespace std;
 
-/* Class for "clamped belief propagation".
+namespace dai {
 
-'chooseNextClamp' can be overridden (e.g. in BPC); this version
-randomly chooses variables to clamp.
-*/
 
-namespace dai {
+/// Find a variable to clamp using BBP (goes with maximum adjoint) 
+/// \see BBP
+std::pair<size_t, size_t> bbpFindClampVar(const InfAlg &in_bp, bool clampingVar,
+            const PropertySet &bbp_props, bbp_cfn_t cfn, 
+            Real *maxVarOut);
 
-template<class T>
-vector<T> concat(const vector<T>& u, const vector<T>& v) {
-    vector<T> w;
-    w.reserve(u.size()+v.size());
-    for(size_t i=0; i<u.size(); i++)
-        w.push_back(u[i]);
-    for(size_t i=0; i<v.size(); i++)
-        w.push_back(v[i]);
-    return w;
-}
 
+/// Class for CBP (Clamped Belief Propagation)
+/** This algorithm uses configurable heuristics to choose a variable
+ *  x_i and a state x_i*. Inference is done with x_i "clamped" to x_i*
+ *  (e.g. conditional on x_i==x_i*), and also with the negation of this
+ *  condition. Clamping is done recursively up to a fixed number of
+ *  levels (other stopping criteria are also implemented, see
+ *  "recursion" property). The resulting approximate marginals are
+ *  combined using logZ estimates.
+ */
 class CBP : public DAIAlgFG {
-    vector<Factor> _beliefs1;
-    vector<Factor> _beliefs2;
+    std::vector<Factor> _beliefsV;
+    std::vector<Factor> _beliefsF;
     double _logZ;
     double _est_logZ;
     double _old_est_logZ;
@@ -37,10 +67,7 @@ class CBP : public DAIAlgFG {
     double _sum_level;
     size_t _num_leaves;
 
-    size_t maxClampLevel() { return props.max_levels; }
-    double minMaxAdj() { return props.min_max_adj; }
-    double recTol() { return props.rec_tol; }
-
+    boost::shared_ptr<std::ofstream> _clamp_ofstream;
 
     bbp_cfn_t BBP_cost_function() {
       return props.bbp_cfn;
@@ -48,67 +75,70 @@ class CBP : public DAIAlgFG {
 
     void printDebugInfo();
 
-    void runRecurse(BP_dual &bp_dual,
+    /// 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,
-                    vector<size_t> clamped_vars_list,
-                    set<size_t> clamped_vars,
+                    std::vector<size_t> clamped_vars_list,
                     size_t &num_leaves,
+                    size_t &choose_count,
                     double &sum_level,
-                    Real &lz_out, vector<Factor>& beliefs_out);
-
-    BP getBP();
-    BP_dual getBP_dual();
+                    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();
 
-  public:
     size_t _iters;
     double _maxdiff;
 
-    struct Properties {
-      typedef BP::Properties::UpdateType UpdateType;
-      DAI_ENUM(RecurseType,REC_FIXED,REC_LOGZ,REC_BDIFF);
-      DAI_ENUM(ChooseMethodType,CHOOSE_RANDOM,CHOOSE_BBP,CHOOSE_BP_ALL);
-      DAI_ENUM(ClampType,CLAMP_VAR,CLAMP_FACTOR);
-
-      double tol;
-      double rec_tol;
-      size_t maxiter;
-      size_t verbose;
-      UpdateType updates;
-      size_t max_levels;
-      double min_max_adj;
-
-      ChooseMethodType choose;
-      RecurseType recursion;
-      ClampType clamp;
-      bbp_cfn_t bbp_cfn;
-      size_t rand_seed;
-    } props;
+    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;
+    }
 
-    Properties::ChooseMethodType ChooseMethod() { return props.choose; }
-    Properties::RecurseType Recursion() { return props.recursion; }
-    Properties::ClampType Clamping() { return props.clamp; }
+    void construct();
+
+  public:
 
     //----------------------------------------------------------------
 
-    // construct CBP object from FactorGraph
+    /// construct CBP object from FactorGraph
     CBP(const FactorGraph &fg, const PropertySet &opts) : DAIAlgFG(fg) {
-        setProperties(opts);
+        props.set(opts);
         construct();
     }
 
     static const char *Name;
-    /// List of property names
-    static const char *PropertyList[];
 
     /// @name General InfAlg interface
     //@{
     virtual CBP* clone() const { return new CBP(*this); }
-//     virtual CBP* create() const { return new CBP(); }
     virtual CBP* create() const { DAI_THROW(NOT_IMPLEMENTED); }
-    virtual std::string identify() const { return string(Name) + printProperties(); }
-    virtual Factor belief (const Var &n) const { return _beliefs1[findVar(n)]; }
+    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 { assert(0); }
-    virtual vector<Factor> beliefs() const { return concat(_beliefs1, _beliefs2); }
+    virtual std::vector<Factor> beliefs() const { return concat(_beliefsV, _beliefsF); }
     virtual Real logZ() const { return _logZ; }
     virtual void init() {};
     virtual void init( const VarSet & ) {};
@@ -117,39 +147,119 @@ class CBP : public DAIAlgFG {
     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]; }
 
     //----------------------------------------------------------------
 
-    virtual bool chooseNextClampVar(BP_dual& bp,
-                                    vector<size_t> &clamped_vars_list,
-                                    set<size_t> &clamped_vars,
-                                    size_t &i, vector<size_t> &xis, Real *maxVarOut);
+ public:
+/* 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 = "";
+   }
+*/
+/* {{{ 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;
+/* }}} END OF GENERATED CODE */
 
-    //----------------------------------------------------------------
+    Properties::ChooseMethodType ChooseMethod() { return props.choose; }
+    Properties::RecurseType Recursion() { return props.recursion; }
+    Properties::ClampType Clamping() { return props.clamp; }
+    size_t maxClampLevel() { return props.max_levels; }
+    double minMaxAdj() { return props.min_max_adj; }
+    double recTol() { return props.rec_tol; }
+};
 
-    void setBeliefsFrom(const BP& b);
-    void setBeliefs(const vector<Factor>& bs, double logZ) {
-        size_t i=0;
-        _beliefs1.clear(); _beliefs1.reserve(nrVars());
-        _beliefs2.clear(); _beliefs2.reserve(nrFactors());
-        for(i=0; i<nrVars(); i++) _beliefs1.push_back(bs[i]);
-        for(; i<nrVars()+nrFactors(); i++) _beliefs2.push_back(bs[i]);
-        _logZ = logZ;
-    }
-    Factor belief1 (size_t i) const { return _beliefs1[i]; }
-    Factor belief2 (size_t I) const { return _beliefs2[I]; }
+/// Given a sorted vector of states \a xis and total state count \a n_states, return a vector of states not in \a xis
+std::vector<size_t> complement( std::vector<size_t>& xis, size_t n_states );
 
-    //----------------------------------------------------------------
+/// Computes \f$\frac{\exp(a)}{\exp(a)+\exp(b)}\f$
+Real unSoftMax( Real a, Real b );
 
-    void construct();
+/// Computes log of sum of exponents, i.e., \f$\log\left(\exp(a) + \exp(b)\right)\f$
+Real logSumExp( Real a, Real b );
+
+/// Compute sum of pairwise L-infinity distances of the first \a nv factors in each vector
+Real dist( const std::vector<Factor>& b1, const std::vector<Factor>& b2, size_t nv );
 
-    void setProperties( const PropertySet &opts );
-    PropertySet getProperties() const;
-    std::string printProperties() const;
-};
 
-pair<size_t, size_t> bbpFindClampVar(BP_dual &in_bp_dual, bool clampVar, size_t numClamped, bbp_cfn_t cfn, const PropertySet &props, Real *maxVarOut);
+} // end of namespace dai
 
-}
 
 #endif
index d201fe5..670b662 100644 (file)
  *  "Sufficient Conditions for Convergence of the Sum-Product Algorithm",
  *  <em>IEEE Transactions on Information Theory</em> 53(12):4422-4437.
  *  http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=4385778
+ *
+ *  \anchor EaG09 \ref EaG09
+ *  F. Eaton and Z. Ghahramani (2009):
+ *  "Choosing a Variable to Clamp",
+ *  <em>Proceedings of the Twelfth International Conference on Artificial Intelligence and Statistics (AISTATS 2009)</em> 5:145-152
+ *  http://jmlr.csail.mit.edu/proceedings/papers/v5/eaton09a/eaton09a.pdf
  */
 
 
index 35d502a..73dfac7 100644 (file)
@@ -1,21 +1,46 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    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 <dai/bp.h>
 #include <dai/bbp.h>
 #include <dai/gibbs.h>
+#include <dai/util.h>
+#include <dai/bipgraph.h>
 
-// for makeBBPGraph: {
+// for makeBBPPlot: {
 #include <sys/stat.h>
 #include <sys/types.h>
 #include <iostream>
 #include <fstream>
 // }
 
+
 namespace dai {
 
+
 using namespace std;
 
-#define rnd_multi(x) rnd_int(0,(x)-1)
 
-/// function to compute \sld\~w from w, Z_w, \sld w
+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);
@@ -29,10 +54,25 @@ Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w) {
     return adj_w_unnorm;
 }
 
-/// Function to turn Gibbs state into b1_adj
-/// calls bbp.init with calculated adjoints
-void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, const vector<size_t>* stateP) {
-    if(cfn_type==bbp_cfn_t::bbp_cfn_t::cfn_bethe_ent) {
+
+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
+        // of variables as iterating over factor(I).vars()
+        size_t j = fg.nbF(I)[_j];
+        f_entry *= fg.var(j).states();
+        f_entry += state[j];
+    }
+    return f_entry;
+}
+
+
+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;
@@ -46,12 +86,12 @@ void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, cons
             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(fg.belief1(i)[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] = -fg.belief1(i)[xi];
+                p[xi] = -ia.beliefV(i)[xi];
             }
             psi1_adj.push_back(p);
         }
@@ -59,24 +99,26 @@ void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, cons
             size_t dim = fg.factor(I).states();
             Prob p(dim,0.0);
             for(size_t xI=0; xI<dim; xI++) {
-                p[xI] = 1 + log(fg.belief2(I)[xI]/fg.factor(I).p()[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] = -fg.belief2(I)[xI]/fg.factor(I).p()[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);
-    } else if(cfn_type==bbp_cfn_t::cfn_factor_ent) {
+        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 = fg.belief2(I)[xI];
+                double bIxI = ia.beliefF(I)[xI];
                 if(bIxI<1.0e-15) {
                     p[xI] = -1.0e10;
                 } else {
@@ -85,15 +127,17 @@ void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, cons
             }
             b2_adj.push_back(p);
         }
-        bbp.init(get_zero_adj_1(fg), b2_adj);
-    } else if(cfn_type==bbp_cfn_t::cfn_var_ent) {
+        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 = fg.belief1(i)[xi];
+                double bixi = ia.beliefV(i)[xi];
                 if(bixi<1.0e-15) {
                     p[xi] = -1.0e10;
                 } else {
@@ -103,10 +147,16 @@ void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, cons
             b1_adj.push_back(p);
         }
         bbp.init(b1_adj);
-    } else { // cost functions that use Gibbs sample: cfn_b, cfn_b2, cfn_exp
+        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(fg,2*fg.doneIters());
+            state = getGibbsState(ia,2*ia.Iterations());
         } else {
             state = *stateP;
         }
@@ -118,82 +168,143 @@ void gibbsInitBBPCostFnAdj(BBP& bbp, const BP_dual &fg, bbp_cfn_t cfn_type, cons
             size_t n = fg.var(i).states();
             Prob delta(n,0.0);
             assert(/*0<=state[i] &&*/ state[i] < n);
-            double b = fg.belief1(i)[state[i]];
-            if(cfn_type==bbp_cfn_t::cfn_b) {
-                delta[state[i]] = 1.0;
-            } else if(cfn_type==bbp_cfn_t::cfn_b2) {
-                delta[state[i]] = b;
-            } else if(cfn_type==bbp_cfn_t::cfn_exp) {
-                delta[state[i]] = exp(b);
-            } else { abort(); }
+            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: abort();
+            }
             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: abort();
+            }
+            b2_adj.push_back(delta);
+        }
+        bbp.init(get_zero_adj_V(fg), b2_adj);
+        break;
+    }
+    default: abort();
     }
 }
 
-/// 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 gibbsCostFn(const BP_dual &fg, bbp_cfn_t cfn_type, const vector<size_t> *stateP) {
+
+Real getCostFn(const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP) {
     double cf=0.0;
-    if(cfn_type==bbp_cfn_t::cfn_bethe_ent) { // ignores state
-        cf = -fg.logZ();
-    } else if(cfn_type==bbp_cfn_t::cfn_var_ent) { // ignores state
+    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 += -fg.belief1(i).entropy();
+            cf += -ia.beliefV(i).entropy();
         }
-    } else {
+        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 = fg.belief1(i)[state[i]];
-            if(cfn_type==bbp_cfn_t::cfn_b) {
-                cf += b;
-            } else if(cfn_type==bbp_cfn_t::cfn_b2) {
-                cf += b*b/2;
-            } else if(cfn_type==bbp_cfn_t::cfn_exp) {
-                cf += exp(b);
-            } else { abort(); }
+            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: abort();
+            }
         }
+        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: abort();
+            }
+        }
+        break;
+    }
+    default: 
+        cerr << "Unknown cost function " << cfn_type << endl;
+        abort();
     }
     return cf;
 }
 
-vector<Prob> get_zero_adj_2(const BP_dual& bp_dual) {
+
+vector<Prob> get_zero_adj_F(const FactorGraph& fg) {
     vector<Prob> adj_2;
-    adj_2.reserve(bp_dual.nrFactors());
-    for(size_t I=0; I<bp_dual.nrFactors(); I++) {
-        adj_2.push_back(Prob(bp_dual.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_1(const BP_dual& bp_dual) {
+
+vector<Prob> get_zero_adj_V(const FactorGraph& fg) {
     vector<Prob> adj_1;
-    adj_1.reserve(bp_dual.nrVars());
-    for(size_t i=0; i<bp_dual.nrVars(); i++) {
-        adj_1.push_back(Prob(bp_dual.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;
 }
 
-#define foreach_iI(_fg)                                                 \
-    for(size_t i=0, I=0, iI=0; (iI<_fg->nr_edges()) && ((i=_fg->edge(iI).first) || 1) && ((I=_fg->edge(iI).second) || 1); iI++)
-    
-#define foreach_iIj(_fg)                                        \
-    for(size_t i=0, I=0, j=0, iIj=0; (iIj<_VFV_ind.size()) &&   \
-            ((i=get<0>(_VFV_ind[iIj])) || 1) &&                 \
-            ((I=get<1>(_VFV_ind[iIj])) || 1) &&                 \
-            ((j=get<2>(_VFV_ind[iIj])) || 1); iIj++)
-
-#define foreach_IiJ(_fg)                                        \
-    for(size_t I=0, i=0, J=0, IiJ=0; (IiJ<_FVF_ind.size()) &&   \
-            ((I=get<0>(_FVF_ind[IiJ])) || 1) &&                 \
-            ((i=get<1>(_FVF_ind[IiJ])) || 1) &&                 \
-            ((J=get<2>(_FVF_ind[IiJ])) || 1); IiJ++)
-
 
 #define LOOP_ij(body) {                                 \
         size_t i_states = _fg->var(i).states();         \
@@ -217,205 +328,408 @@ vector<Prob> get_zero_adj_1(const BP_dual& bp_dual) {
         }                                               \
     }
 
+
 void BBP::RegenerateInds() {
-    // initialise _VFV and _VFV_ind
-    {
-        size_t ind=0;
-        _VFV.resize(_fg->nrVars());
-        _VFV_ind.clear();
-        for(size_t i=0; i<_fg->nrVars(); i++) {
-            foreach(size_t I, _fg->nbV(i)) {
-                vector<size_t>& v = _VFV[i][I];
-                v.resize(_fg->nrVars());
-                foreach(size_t j, _fg->nbF(I)) {
-                    if(j!=i) {
-                        v[j] = ind++;
-                        _VFV_ind.push_back(tuple<size_t,size_t,size_t>(i,I,j));
-                    }
-                }
-            }
-        }
-    }
-    // initialise _FVF
-    {
-        size_t ind=0;
-        _FVF.resize(_fg->nrFactors());
-        _FVF_ind.clear();
-        for(size_t I=0; I<_fg->nrFactors(); I++) {
-            foreach(size_t i, _fg->nbF(I)) {
-                vector<size_t>& v = _FVF[I][i];
-                v.resize(_fg->nrFactors());
-                foreach(size_t J, _fg->nbV(i)) {
-                    if(J!=I) {
-                        v[J] = ind++;
-                        _FVF_ind.push_back(tuple<size_t,size_t,size_t>(I,i,J));
-                    }
-                }
+    // initialise _indices
+    //     typedef std::vector<size_t>  _ind_t;
+    //     std::vector<std::vector<_ind_t> >  _indices; 
+    _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)) {
+            _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);
         }
     }
 }
 
+
 void BBP::RegenerateT() {
+    // _T[i][_I]
     _T.clear();
-    foreach_iI(_fg) {
-        Prob prod(_fg->var(i).states(),1.0);
-        foreach(size_t J, _fg->nbV(i)) {
-            if(J!=I) {
-                prod *= _bp_dual->msgM(J,i);
+    _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][I.iter] = prod;
         }
-        _T.push_back(prod);
     }
 }
 
+
 void BBP::RegenerateU() {
+    // _U[I][_i]
     _U.clear();
-    foreach_iI(_fg) {
-        //     Prob prod(_fg->var(i).states(), 1.0);
-        Prob prod(_fg->factor(I).states(), 1.0);
-        foreach(size_t j, _fg->nbF(I)) {
-            if(i!=j) {
-                Prob n_jI(_bp_dual->msgN(j,I));
-                const BP_dual::_ind_t* ind = &(_bp_dual->index(j,I));
-                // 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]];
-                //         prod *= _bp_dual->msgN(j,I);
+    _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);
+                    // multiply prod by n_jI
+                    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;
         }
-        _U.push_back(prod);
     }
 }
 
+
 void BBP::RegenerateS() {
+    // _S[i][_I][_j]
     _S.clear();
-    foreach_iIj(_fg) {
-        Factor prod(_fg->factor(I));
-        foreach(size_t k, _fg->nbF(I)) {
-            if(k!=i && k!=j) {
-                if(0) {
-                    prod *= Factor(_fg->var(k), _bp_dual->msgN(k,I));
-                } else {
-                    const BP_dual::_ind_t& ind = _bp_dual->index(k,I);
-                    Prob p(_bp_dual->msgN(k,I));
-                    for(size_t xI=0; xI<prod.states(); xI++) {
-                        prod.p()[xI] *= p[ind[xI]];
+    _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));
+                            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? 
+                    Prob marg;
+                    marg = prod.marginal(VarSet(_fg->var(i)) | VarSet(_fg->var(j)), false).p();
+                    _S[i][I.iter][j.iter] = marg;
                 }
             }
         }
-        // XXX improve efficiency? 
-        // "Marginalize" onto i|j (unnormalized)
-        Prob marg;
-        marg = prod.marginal(VarSet(_fg->var(i)) | VarSet(_fg->var(j)), false).p();
-        _S.push_back(marg);
     }
 }
 
+
 void BBP::RegenerateR() {
+    // _R[I][_i][_J]
     _R.clear();
-    foreach_IiJ(_fg) {
-        Prob prod(_fg->var(i).states(), 1.0);
-        foreach(size_t K, _fg->nbV(i)) {
-            if(K!=I && K!=J) {
-                prod *= _bp_dual->msgM(K,i);
+    _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)) {
+                            prod *= _bp_dual.msgM(i,K.iter);
+                        }
+                    }
+                    _R[I][i.iter][J.iter] = prod;
+                }
             }
         }
-        _R.push_back(prod);
     }
 }
 
+
 void BBP::RegenerateInputs() {
-    _adj_b_1_unnorm.clear();
-    for(size_t i=0; i<_fg->nrVars(); i++) {
-        _adj_b_1_unnorm.push_back(unnormAdjoint(_bp_dual->belief1(i).p(), _bp_dual->belief1Z(i), _adj_b_1[i]));
+    _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_2_unnorm.clear();
-    for(size_t I=0; I<_fg->nrFactors(); I++) {
-        _adj_b_2_unnorm.push_back(unnormAdjoint(_bp_dual->belief2(I).p(), _bp_dual->belief2Z(I), _adj_b_2[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]));
     }
 }
 
-/// initialise members for factor adjoints (call after RegenerateInputs)
+
 void BBP::RegeneratePsiAdjoints() {
-    _adj_psi_1.clear();
-    _adj_psi_1.reserve(_fg->nrVars());
-    for(size_t i=0; i<_fg->nrVars(); i++) {
-        Prob p(_adj_b_1_unnorm[i]);
+    _adj_psi_V.clear();
+    _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(size_t I, _fg->nbV(i)) {
-            p *= _bp_dual->msgM(I,i);
+        foreach(const Neighbor& I, _fg->nbV(i)) {
+            p *= _bp_dual.msgM(i,I.iter);
         }
-        p += _init_adj_psi_1[i];
-        _adj_psi_1.push_back(p);
+        p += _init_adj_psi_V[i];
+        _adj_psi_V.push_back(p);
     }
-    _adj_psi_2.clear();
-    _adj_psi_2.reserve(_fg->nrFactors());
-    for(size_t I=0; I<_fg->nrFactors(); I++) {
-        Prob p(_adj_b_2_unnorm[I]);
+    _adj_psi_F.clear();
+    _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(size_t i, _fg->nbF(I)) {
-            Prob n_iI(_bp_dual->msgN(i,I));
-            const BP_dual::_ind_t* ind = &(_bp_dual->index(i,I));
+        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++)
-                p[x_I] *= n_iI[(*ind)[x_I]];
+                p[x_I] *= n_iI[ind[x_I]];
         }
-        p += _init_adj_psi_2[I];
-        _adj_psi_2.push_back(p);
+        p += _init_adj_psi_F[I];
+        _adj_psi_F.push_back(p);
     }
 }
 
-/// initialise members for messages adjoints (call after RegenerateInputs)
-void BBP::RegenerateMessageAdjoints() {
-    _adj_n.resize(_fg->nr_edges());
-    _adj_m.resize(_fg->nr_edges());
-    _adj_n_unnorm.resize(_fg->nr_edges());
-    _adj_m_unnorm.resize(_fg->nr_edges());
+
+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.reserve(_fg->nr_edges());
-    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
-        size_t i = _fg->edge(iI).first;
-        size_t I = _fg->edge(iI).second;
-        Prob prod(_fg->factor(I).p());
-        prod *= _adj_b_2_unnorm[I];
-        foreach(size_t j, _fg->nbF(I)) {
-            if(j!=i) { // for all j in I\i
-                Prob n_jI(_bp_dual->msgN(j,I));; 
-                const BP_dual::_ind_t* ind = &(_bp_dual->index(j,I));
-                // 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]];
+    _new_adj_n.resize(nv);
+    _new_adj_m.clear();
+    _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());
+                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];
+                _new_adj_n[i][I.iter] = marg;
+                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)) {
+                        prod *= _bp_dual.msgM(i,J.iter);
+                    }
+                }
+                _new_adj_m[i][I.iter] = prod;
+                upMsgM(i,I.iter);
             }
         }
-        Prob marg( _fg->var(i).states(), 0.0 );
-        // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
-        const BP_dual::_ind_t* ind = &(_bp_dual->index(i,I));
-        for( size_t r = 0; r < prod.size(); r++ )
-            marg[(*ind)[r]] += prod[r];
+    }
+}
 
-        _new_adj_n.push_back(marg);
-        upMsgN(iI); // propagate new _new_adj_n to _adj_n and _adj_n_unnorm
+
+void BBP::incrSeqMsgM(size_t i, size_t _I, const Prob& p) {
+    if(props.clean_updates) {
+        _new_adj_m[i][_I] += p;
+    } else {
+        _adj_m[i][_I] += p;
+        calcUnnormMsgM(i, _I);
     }
-    _new_adj_m.clear();
-    _new_adj_m.reserve(_fg->nr_edges());
-    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
-        size_t i = _fg->edge(iI).first;
-        size_t I = _fg->edge(iI).second;
-        //     Prob prod(_fg->var(i).states(),1.0);
-        Prob prod(_adj_b_1_unnorm[i]);
-        assert(prod.size()==_fg->var(i).states());
-        foreach(size_t J, _fg->nbV(i)) {
-            if(J!=I) { // for all J in i\I
-                prod *= _bp_dual->msgM(J,i);
+}
+
+
+Real pv_thresh=1000;
+
+
+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) {
+
+            DAI_DMSG("in updateSeqMsgM");
+            DAI_PV(i);
+            DAI_PV(_I);
+            DAI_PV(_adj_m[i][_I]);
+            DAI_PV(_new_adj_m[i][_I]);
+        }
+#endif
+        _adj_m[i][_I] += _new_adj_m[i][_I];
+        calcUnnormMsgM(i, _I);
+        _new_adj_m[i][_I].fill(0.0);
+    }
+}
+
+
+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 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) {
+        DAI_DMSG("in sendSeqMsgN");
+        DAI_PV(i);
+        DAI_PV(I);
+        DAI_PV(_I);
+        DAI_PV(_bp_dual.msgN(i,_I));
+        DAI_PV(_bp_dual.zN(i,_I));
+        DAI_PV(_bp_dual.msgM(i,_I));
+        DAI_PV(_bp_dual.zM(i,_I));
+        DAI_PV(_fg->factor(I).p());
+    }
+#endif
+    foreach(const Neighbor& J, _fg->nbV(i)) {
+        if(size_t(J) != size_t(I)) {
+#if 0
+            if(f_unnorm.sumAbs() > pv_thresh) {
+                DAI_DMSG("in sendSeqMsgN loop");
+                DAI_PV(J);
+                DAI_PV(f_unnorm);
+                DAI_PV(_R[J][J.dual][_I]);
+                DAI_PV(f_unnorm * _R[J][J.dual][_I]);
+            }
+#endif
+            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];
+
+//     DAI_PV(j);
+//     DAI_PV(I);
+//     DAI_PV(_adj_m_unnorm_jI);
+//     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);
+    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;
+    _adj_psi_F[I] += um;
+
+//     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) {
+            const Prob &S = _S[i][i.dual][_j];
+            Prob msg(_fg->var(i).states(),0.0);
+            LOOP_ij(
+                msg[xi] += S[xij]*_adj_m_unnorm_jI[xj];
+            );
+            msg *= 1-props.damping;
+#if 0
+            if(msg.sumAbs() > pv_thresh) {
+                DAI_DMSG("in sendSeqMsgM loop");
+
+                DAI_PV(j);
+                DAI_PV(I);
+                DAI_PV(_I);
+                DAI_PV(_fg->nbF(I).size());
+                DAI_PV(_fg->factor(I).p());
+                DAI_PV(_S[i][i.dual][_j]);
+
+                DAI_PV(i);
+                DAI_PV(i.dual);
+                DAI_PV(msg);
+                DAI_PV(_fg->nbV(i).size());
+            }
+#endif
+            assert(size_t(_fg->nbV(i)[i.dual]) == 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);
+}
+
+
+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);
         }
-        _new_adj_m.push_back(prod);
-        upMsgM(iI); // propagate new _new_adj_m to _adj_m and _adj_m_unnorm
     }
 }
 
+
 void BBP::Regenerate() {
     RegenerateInds();
     RegenerateT();
@@ -424,184 +738,377 @@ void BBP::Regenerate() {
     RegenerateR();
     RegenerateInputs();
     RegeneratePsiAdjoints();
-    RegenerateMessageAdjoints();
+    if(props.updates == Properties::UpdateType::PAR) {
+        RegenerateParMessageAdjoints();
+    } else {
+        RegenerateSeqMessageAdjoints();
+    }
     _iters = 0;
 }
-  
-void BBP::calcNewN(size_t iI) {
-    size_t i=_fg->edge(iI).first;
-    size_t I=_fg->edge(iI).second;
-    _adj_psi_1[i] += T(i,I)*_adj_n_unnorm[iI];
-    _trip_end_t vfv_iI = _VFV[i][I];
-    Prob &new_adj_n_iI = _new_adj_n[iI];
+
+
+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);
-    foreach(size_t j, _fg->nbF(I)) {
+    size_t I = _fg->nbV(i)[_I];
+    foreach(const Neighbor& j, _fg->nbF(I)) {
         if(j!=i) {
-            size_t iIj = vfv_iI[j];
-            Prob &p = _S[iIj];
-            size_t jI = _fg->VV2E(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];
             );
         }
     }
 }
 
-void BBP::calcNewM(size_t iI) {
-    size_t i=_fg->edge(iI).first;
-    size_t I=_fg->edge(iI).second;
 
-    Prob p(U(I,i));
-    Prob adj(_adj_m_unnorm[iI]);
-    const BP_dual::_ind_t* ind = &(_bp_dual->index(i,I));
-    for(size_t x_I = 0; x_I < p.size(); x_I++)
-        p[x_I] *= adj[(*ind)[x_I]];
-    _adj_psi_2[I] += 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;
 
-    _new_adj_m[iI] = Prob(_fg->var(i).states(),0.0);
-    _trip_end_t fvf_Ii = _FVF[I][i];
-    foreach(size_t J, _fg->nbV(i)) {
+    _new_adj_m[i][_I] = Prob(_fg->var(i).states(),0.0);
+    foreach(const Neighbor& J, _fg->nbV(i)) {
         if(J!=I) {
-            size_t iJ = _fg->VV2E(i,J);
-            _new_adj_m[iI] += _R[fvf_Ii[J]]*_adj_n_unnorm[iJ];
+            _new_adj_m[i][_I] += _R[I][I.dual][J.iter]*_adj_n_unnorm[i][J.iter];
         }
     }
 }
 
-void BBP::upMsgM(size_t iI) {
-    _adj_m[iI] = _new_adj_m[iI];
-    _adj_m_unnorm[iI] = unnormAdjoint(_bp_dual->msgM(iI), _bp_dual->zM(iI), _adj_m[iI]);
+
+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::upMsgM(size_t i, size_t _I) {
+    _adj_m[i][_I] = _new_adj_m[i][_I];
+    calcUnnormMsgM(i,_I);
 }
 
-void BBP::upMsgN(size_t iI) {
-    _adj_n[iI] = _new_adj_n[iI];
-    _adj_n_unnorm[iI] = unnormAdjoint(_bp_dual->msgN(iI), _bp_dual->zN(iI), _adj_n[iI]);
+
+void BBP::upMsgN(size_t i, size_t _I) {
+    _adj_n[i][_I] = _new_adj_n[i][_I];
+    calcUnnormMsgN(i,_I);
 }
 
+
 void BBP::doParUpdate() {
-    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
-        calcNewM(iI);
-        calcNewN(iI);
+    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 iI=0; iI<_fg->nr_edges(); iI++) {
-        upMsgM(iI);
-        upMsgN(iI);
+    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;
-    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
-        s += _adj_m_unnorm[iI].sumAbs();
-        s += _adj_n_unnorm[iI].sumAbs();
+    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/_fg->nr_edges();
+    return s/e;
 }
 
+
 void BBP::getMsgMags(Real &s, Real &new_s) {
     s=0.0; new_s=0.0;
-    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
-        s += _adj_m[iI].sumAbs();
-        s += _adj_n[iI].sumAbs();
-        new_s += _new_adj_m[iI].sumAbs();
-        new_s += _new_adj_n[iI].sumAbs();
-    }
-    Real t = _fg->nr_edges();
-    s /= t; new_s /= t;
-}
-
-/// run until change is less than given tolerance
-void BBP::run(Real tol) {
-    size_t minIters = 2;
-    do {
-        _iters++;
-        doParUpdate();
-    } while((_iters < minIters || getUnMsgMag() > tol) && _iters < maxIter());
-    if(_iters==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;
-    }
-}
-
-vector<size_t> getGibbsState(const BP_dual& fg, size_t iters) {
+    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;
+}
+
+// tuple<size_t,size_t,Real> BBP::getArgMaxPsi1Adj() {
+//     size_t argmax_var=0;
+//     size_t argmax_var_state=0;
+//     Real max_var=0;
+//     for( size_t i = 0; i < _fg->nrVars(); i++ ) {
+//         pair<size_t,Real> argmax_state = adj_psi_V(i).argmax();
+//         if(i==0 || argmax_state.second>max_var) {
+//             argmax_var = i;
+//             max_var = argmax_state.second;
+//             argmax_var_state = argmax_state.first;
+//         }
+//     }
+//     assert(/*0 <= argmax_var_state &&*/
+//            argmax_var_state < _fg->var(argmax_var).states());
+//     return tuple<size_t,size_t,Real>(argmax_var,argmax_var_state,max_var);
+// }
+
+
+Real BBP::getMaxMsgM() {
+    size_t dummy;
+    Real 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)) {
+            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)) {
+            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)) {
+            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;
+    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;
+    }
+}
+
+
+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;
+    }
+}
+
+
+vector<size_t> getGibbsState(const InfAlg& ia, size_t iters) {
     PropertySet gibbsProps;
     gibbsProps.Set("iters", iters);
-    gibbsProps.Set("verbose", oneLess(fg.Verbose()));
-    Gibbs gibbs(fg, gibbsProps);
+    gibbsProps.Set("verbose", size_t(0));
+    Gibbs gibbs(ia.fg(), gibbsProps);
     gibbs.run();
     return gibbs.state();
 }
 
-void testUnnormAdjoint(int n);
 
-/// 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_1.
-/// 'h' controls size of perturbation.
-/// 'bbpTol' controls tolerance of BBP run.
-double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, double h) {
-    //   testUnnormAdjoint(4);
-
-    vector<size_t> state = getGibbsState(bp_dual,2*bp_dual.doneIters());
+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 = gibbsCostFn(bp_dual, cfn, &state);
+    Real cf0 = getCostFn(bp, cfn, state);
 
     // run BBP to estimate adjoints
-    BBP bbp(bp_dual);
-    gibbsInitBBPCostFnAdj(bbp, bp_dual, cfn, &state);
-    bbp.run(bbpTol);
+    BBP bbp(&bp,bbp_props);
+    initBBPCostFnAdj(bbp, bp, cfn, state);
+    bbp.run();
 
     Real d=0;
+    const FactorGraph& fg = bp.fg();
 
     if(1) {
-        // verify bbp.adj_psi_1
+        // verify bbp.adj_psi_V
 
         // for each variable i
-        for(size_t i=0; i<bp_dual.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<bp_dual.var(i).states(); xi++) {
-                // create a copy of bp_dual which multiplies this into some factor containing the variable
-                BP_dual bp_dual_prb(bp_dual);
-                assert(bp_dual_prb.nbV(i).size()>0);
+            for(size_t xi=0; xi<fg.var(i).states(); xi++) {
+                // Clone 'bp' (which may be any InfAlg)
+                InfAlg *bp_prb = bp.clone();
 
-                // create a perturbed psi_1[i] (there is no single-variable
-                // psi in libDAI so we initialize to all 1 and then multiply
-                // into nearest multivariate factor)
-                size_t n = bp_dual.var(i).states();
-                Prob psi_1_prb(n, 1.0);
+                // perturb it
+                size_t n = bp_prb->fg().var(i).states();
+                Prob psi_1_prb(n,1.0);
                 psi_1_prb[xi] += h;
-                psi_1_prb.normalize();
+//                 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);
+                
+                // call 'init' on the perturbed variables
+                bp_prb->init(bp_prb->fg().var(i));
+                
+                // run copy to convergence
+                bp_prb->run();
 
-                // update bp_dual_prb
-                size_t I = bp_dual_prb.nbV(i)[0]; // use the first factor in list of neighbours of i
-                bp_dual_prb.factor(I) *= Factor(bp_dual.var(i), psi_1_prb);
-                bp_dual_prb.init(bp_dual.var(i)); // reset messages involving i
-                //       bp_dual_prb.init();
-    
-                // run the copy to convergence
-                bp_dual_prb.run();
-    
-                // calculate the new value of the cost function
-                Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                // calculate new value of cost function
+                Real cf_prb = getCostFn(*bp_prb, cfn, state);
 
-                // use this to estimate the adjoint for i
-                // XXX why is it off by a factor of 2?
+                // use to estimate adjoint for i
                 adj_est.push_back((cf_prb-cf0)/h);
+                
+                // free cloned InfAlg
+                delete bp_prb;
             }
             Prob p_adj_est(adj_est.begin(), adj_est.end());
             // compare this numerical estimate to the BBP estimate; sum the distances
             cerr << "i: " << i
                  << ", p_adj_est: " << p_adj_est
-                 << ", bbp.adj_psi_1(i): " << bbp.adj_psi_1(i) << endl;
-            d += dist(p_adj_est, bbp.adj_psi_1(i), Prob::DISTL1);
+                 << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
+            d += dist(p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1);
         }
     }
-    if(1) {
+    /*    if(1) {
         // verify bbp.adj_n and bbp.adj_m
 
         // We actually want to check the responsiveness of objective
@@ -624,7 +1131,7 @@ double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, doub
                     // recalculate beliefs
                     bp_dual_prb.CalcBeliefs();
                     // get cost function value
-                    Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                    Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
                     // add it to list of adjoints
                     adj_n_est.push_back((cf_prb-cf0)/h);
                 }
@@ -638,7 +1145,7 @@ double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, doub
                     // recalculate beliefs
                     bp_dual_prb.CalcBeliefs();
                     // get cost function value
-                    Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                    Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
                     // add it to list of adjoints
                     adj_m_est.push_back((cf_prb-cf0)/h);
                 }
@@ -659,10 +1166,11 @@ double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, doub
             }
         }
     }
-    if(0) {
-        // verify bbp.adj_b_1
+    */
+    /*    if(0) {
+        // verify bbp.adj_b_V
         for(size_t i=0; i<bp_dual.nrVars(); i++) {
-            vector<double> adj_b_1_est;
+            vector<double> adj_b_V_est;
             // for each value xi
             for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
                 BP_dual bp_dual_prb(bp_dual);
@@ -671,100 +1179,107 @@ double numericBBPTest(const BP_dual& bp_dual, bbp_cfn_t cfn, double bbpTol, doub
                 bp_dual_prb._beliefs.b1[i][xi] += h;
 
                 // get cost function value
-                Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
 
                 // add it to list of adjoints
-                adj_b_1_est.push_back((cf_prb-cf0)/h);
+                adj_b_V_est.push_back((cf_prb-cf0)/h);
             }
-            Prob p_adj_b_1_est(adj_b_1_est.begin(), adj_b_1_est.end());
+            Prob p_adj_b_V_est(adj_b_V_est.begin(), adj_b_V_est.end());
             // compare this numerical estimate to the BBP estimate; sum the distances
             cerr << "i: " << i
-                 << ", adj_b_1_est: " << p_adj_b_1_est
-                 << ", bbp.adj_b_1(i): " << bbp.adj_b_1(i) << endl;
-            d += dist(p_adj_b_1_est, bbp.adj_b_1(i), Prob::DISTL1);
+                 << ", adj_b_V_est: " << p_adj_b_V_est
+                 << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
+            d += dist(p_adj_b_V_est, bbp.adj_b_V(i), Prob::DISTL1);
         }
     }
+    */
 
     // return total of distances
     return d;
 }
 
-void testUnnormAdjoint(int n) {
-    double h = 1.0e-5;
-    Prob q_unnorm(n);
-    // create a random q_unnorm
-    for(int i=0; i<n; i++) {
-        q_unnorm[i] = exp(4*rnd_stdnormal());
-    }
-    // normalize it to get q and Z_q
-    Prob q(q_unnorm);
-    double Z_q = q.normalize();
-    // cost function V just selects one (random) element of q
-    int vk = rnd_multi(n);
-    double V = q[vk];
-    // calculate adj_q for this V
-    Prob adj_q(n, 0.0);
-    adj_q[vk] = 1;
-    // use unnormAdjoint to get adj_q_unnorm
-    Prob adj_q_unnorm = unnormAdjoint(q, Z_q, adj_q);
-    // with perturbations of q_unnorm, test that adj_q_unnorm is correct
-    Prob adj_q_unnorm_numeric(n, 0.0);
-    for(int i=0; i<n; i++) {
-        Prob q_unnorm_prb = q_unnorm;
-        q_unnorm_prb[i] += h;
-        Prob q_prb(q_unnorm_prb);
-        q_prb.normalize();
-        DAI_PV(q_unnorm_prb);
-        DAI_PV(q_prb);
-        double V_prb = q_prb[vk];
-        adj_q_unnorm_numeric[i] = (V_prb-V)/h;
-    }
-    DAI_PV(q_unnorm);
-    DAI_PV(q);
-    DAI_PV(Z_q);
-    DAI_PV(vk);
-    DAI_PV(V);
-    DAI_PV(adj_q_unnorm);
-    DAI_PV(adj_q_unnorm_numeric);
-}
-
-void makeBBPGraph(const BP_dual& bp_dual, bbp_cfn_t cfn) {
-    double tiny=1.0e-7;
-    vector<size_t> state = getGibbsState(bp_dual,2*bp_dual.doneIters());
-    const char *graphs_dir = "./bbp-data";
-    mkdir(graphs_dir, -1);
-    size_t num_samples=500;
-    for(size_t i=0; i<bp_dual.nrVars(); i++) {
-        for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
-            if(bp_dual.belief1(i)[xi]<tiny || abs(bp_dual.belief1(i)[xi]-1)<tiny)
-                continue;
-            char *fn;
-            asprintf(&fn,"%s/i:%d-xi:%d.out",graphs_dir,(int)i,(int)xi);
-            ofstream os(fn, ios_base::out|ios_base::trunc);
-      
-            for(size_t n=0; n<num_samples; n++) {
-                double c=((double)n)/(num_samples-1);
-                Prob psi_1_prb(bp_dual.var(i).states(),1.0-c);
-                psi_1_prb[xi] = 1.0;
 
-                BP_dual bp_dual_prb(bp_dual);
-                assert(bp_dual_prb.nbV(i).size()>0);
-                size_t I = bp_dual_prb.nbV(i)[0]; // use the first factor in list of neighbours of i
-                bp_dual_prb.factor(I) *= Factor(bp_dual.var(i), psi_1_prb);
-                bp_dual_prb.init(bp_dual.var(i)); // reset messages involving i
-    
-                // run the copy to convergence
-                bp_dual_prb.run();
-        
-                Real cf = gibbsCostFn(bp_dual_prb, cfn, &state);
+} // end of namespace dai
 
-                os << n << "\t" << cf << endl;
-            }
-      
-            os.close();
-            free(fn);
-        }
+
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp 
+*/
+namespace dai {
+
+void BBP::Properties::set(const PropertySet &opts)
+{
+    const std::set<PropertyKey> &keys = opts.allKeys();
+    std::set<PropertyKey>::const_iterator i;
+    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 == "damping") continue;
+        if(*i == "updates") continue;
+        if(*i == "clean_updates") continue;
+        cerr << "BBP: Unknown property " << *i << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
+    }
+    if(!opts.hasKey("verbose")) {
+        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("damping")) {
+        cerr << "BBP: Missing property \"damping\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("updates")) {
+        cerr << "BBP: Missing property \"updates\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("clean_updates")) {
+        cerr << "BBP: Missing property \"clean_updates\" for method \"BBP\"" << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
+    }
+    verbose = opts.getStringAs<size_t>("verbose");
+    tol = opts.getStringAs<double>("tol");
+    maxiter = opts.getStringAs<size_t>("maxiter");
+    damping = opts.getStringAs<double>("damping");
+    updates = opts.getStringAs<UpdateType>("updates");
+    clean_updates = opts.getStringAs<bool>("clean_updates");
+}
+PropertySet BBP::Properties::get() const {
+    PropertySet opts;
+    opts.Set("verbose", verbose);
+    opts.Set("tol", tol);
+    opts.Set("maxiter", maxiter);
+    opts.Set("damping", damping);
+    opts.Set("updates", updates);
+    opts.Set("clean_updates", clean_updates);
+    return opts;
+}
+string BBP::Properties::toString() const {
+    stringstream s(stringstream::out);
+    s << "[";
+    s << "verbose=" << verbose << ",";
+    s << "tol=" << tol << ",";
+    s << "maxiter=" << maxiter << ",";
+    s << "damping=" << damping << ",";
+    s << "updates=" << updates << ",";
+    s << "clean_updates=" << clean_updates;
+    s << "]";
+    return s.str();
 }
-
 } // end of namespace dai
+/* }}} END OF GENERATED CODE */
index ca11449..aab36e6 100644 (file)
@@ -1,3 +1,23 @@
+/*  Copyright (C) 2009  Frederik Eaton [frederik at ofb dot net]
+
+    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 <sstream>
 #include <map>
 #include <dai/cbp.h>
 #include <dai/bbp.h>
 
-using namespace std;
 
 namespace dai {
 
-const char *CBP::Name = "CBP";
-
-const char *CBP::PropertyList[] = {"updates","tol","rec_tol","maxiter","verbose","max_levels","min_max_adj","choose","clamp","recursion","bbp_cfn","rand_seed"};
 
-#define rnd_multi(x) rnd_int(0,(x)-1)
+using namespace std;
+using boost::shared_ptr;
 
-void CBP::setProperties(const PropertySet &opts) {
-//     DAI_DMSG("in CBP::setProperties");
-//     DAI_PV(opts);
-    foreach(const char* p, PropertyList) {
-        if(!opts.hasKey(p)) {
-            // XXX probably leaks pointer?
-            throw (string("CBP: Missing property ")+p).c_str();
-        }
-    }
-    
-    props.tol = opts.getStringAs<double>("tol");
-    props.rec_tol = opts.getStringAs<double>("rec_tol");
-    props.maxiter = opts.getStringAs<size_t>("maxiter");
-    props.verbose = opts.getStringAs<size_t>("verbose");
-    props.updates = opts.getStringAs<Properties::UpdateType>("updates");
-    props.max_levels = opts.getStringAs<size_t>("max_levels");
-    props.min_max_adj = opts.getStringAs<double>("min_max_adj");
-    props.choose = opts.getStringAs<Properties::ChooseMethodType>("choose");
-    props.recursion = opts.getStringAs<Properties::RecurseType>("recursion");
-    props.clamp = opts.getStringAs<Properties::ClampType>("clamp");
-    props.bbp_cfn = opts.getStringAs<bbp_cfn_t>("bbp_cfn");
-    props.rand_seed = opts.getStringAs<size_t>("rand_seed");
-}
 
-PropertySet CBP::getProperties() const {
-    PropertySet opts;
-    opts.Set( "tol", props.tol );
-    opts.Set( "rec_tol", props.rec_tol );
-    opts.Set( "maxiter", props.maxiter );
-    opts.Set( "verbose", props.verbose );
-    opts.Set( "updates", props.updates );
-    opts.Set( "max_levels", props.max_levels );
-    opts.Set( "min_max_adj", props.min_max_adj );
-    opts.Set( "choose", props.choose );
-    opts.Set( "recursion", props.recursion );
-    opts.Set( "clamp", props.clamp );
-    opts.Set( "bbp_cfn", props.bbp_cfn );
-    opts.Set( "rand_seed", props.rand_seed );
-    return opts;
-}
+const char *CBP::Name = "CBP";
 
-std::string CBP::printProperties() const {
-    stringstream s( stringstream::out );
-    s << "[";
-    s << "tol=" << props.tol << ",";
-    s << "rec_tol=" << props.rec_tol << ",";
-    s << "maxiter=" << props.maxiter << ",";
-    s << "verbose=" << props.verbose << ",";
-    s << "updates=" << props.updates << ",";
-    s << "max_levels=" << props.max_levels << ",";
-    s << "min_max_adj=" << props.min_max_adj << ",";
-    s << "choose=" << props.choose << ",";
-    s << "recursion=" << props.recursion << ",";
-    s << "clamp=" << props.clamp << ",";
-    s << "bbp_cfn=" << props.bbp_cfn << ",";
-    s << "rand_seed=" << props.rand_seed << ",";
-    s << "]";
-    return s.str();
-}
 
 void CBP::construct() {
-//     DAIAlgFG::Regenerate();
     indexEdges();
 
-    _beliefs1.clear(); _beliefs1.reserve(nrVars());
+    _beliefsV.clear(); _beliefsV.reserve(nrVars());
     for( size_t i = 0; i < nrVars(); i++ )
-        _beliefs1.push_back( Factor(var(i)).normalized() );
+        _beliefsV.push_back( Factor(var(i)).normalized() );
 
-    _beliefs2.clear(); _beliefs2.reserve(nrFactors());
+    _beliefsF.clear(); _beliefsF.reserve(nrFactors());
     for( size_t I = 0; I < nrFactors(); I++ ) {
         Factor f = factor(I);
         f.fill(1); f.normalize();
-        _beliefs2.push_back(f);
+        _beliefsF.push_back(f);
     }
 
     // to compute average level
@@ -102,8 +62,14 @@ void CBP::construct() {
 
     _maxdiff = 0;
     _iters = 0;
+
+    if(props.clamp_outfile.length()>0) {
+        _clamp_ofstream = shared_ptr<ofstream>(new ofstream(props.clamp_outfile.c_str(), ios_base::out|ios_base::trunc));
+        *_clamp_ofstream << "# COUNT LEVEL VAR STATE" << endl;
+    }
 }
 
+
 static
 vector<Factor> mixBeliefs(Real p, vector<Factor> b, vector<Factor> c) {
   vector<Factor> out;
@@ -111,29 +77,29 @@ vector<Factor> mixBeliefs(Real p, vector<Factor> b, vector<Factor> c) {
   out.reserve(b.size());
   Real pc = 1-p;
   for(size_t i=0; i<b.size(); i++) {
-      // XXX probably already normalized
+    // probably already normalized, but do it again just in case
     out.push_back(b[i].normalized()*p+
                   c[i].normalized()*pc);
   }
   return out;
 }
 
+
 double CBP::run() {
-//     BP bp(getBP());
-//     InfAlg *bp = newInfAlg( GetPropertyAs<string>("bp_alg"), *this, GetPropertyAs<Properties>("bp_opts") );
     size_t seed = props.rand_seed;
     if(seed>0) rnd_seed(seed);
 
-    BP_dual bp_dual(getBP_dual());
-    bp_dual.init();
-    bp_dual.run();
-    _iters += bp_dual.Iterations();
+    InfAlg *bp = getInfAlg();
+    bp->init();
+    bp->run();
+    _iters += bp->Iterations();
 
     vector<Factor> beliefs_out;
     Real lz_out;
-    runRecurse(bp_dual, bp_dual.logZ(), 
-               vector<size_t>(0), set<size_t>(),
-               _num_leaves, _sum_level,
+    size_t choose_count=0;
+    runRecurse(bp, bp->logZ(), 
+               vector<size_t>(0),
+               _num_leaves, choose_count, _sum_level,
                lz_out, beliefs_out);
     if(props.verbose>=1) {
       cerr << "CBP average levels = " << (_sum_level/_num_leaves) << ", leaves = " << _num_leaves << endl;
@@ -142,69 +108,28 @@ double CBP::run() {
     return 0;
 }
 
-BP CBP::getBP() {
-    PropertySet bpProps;
-    bpProps.Set("updates", string("PARALL"));
-    bpProps.Set("tol", props.tol);
-    bpProps.Set("maxiter", props.maxiter);
-    bpProps.Set("verbose", oneLess(props.verbose));
-    BP bp(*this,bpProps);
-    bp.init();
-    return bp;
-}
 
-BP_dual CBP::getBP_dual() {
+/// \todo Eventually this allow other inference algorithms that BP to be selected, via a property
+InfAlg* CBP::getInfAlg() {
     PropertySet bpProps;
-    bpProps.Set("updates", string("PARALL"));
+    bpProps.Set("updates", props.updates);
     bpProps.Set("tol", props.tol);
     bpProps.Set("maxiter", props.maxiter);
-    bpProps.Set("verbose", oneLess(props.verbose));
-//     cerr << "In getBP_dual" << endl;
-//     DAI_PV(bpProps);
-    BP_dual bp_dual(*this,bpProps);
-    return bp_dual;
-}
-
-vector<size_t> complement(vector<size_t>& xis, size_t n_states) {
-    vector<size_t> cmp_xis(0);
-    size_t j=0;
-    for(size_t xi=0; xi<n_states; xi++) {
-        while(j<xis.size() && xis[j]<xi) j++;
-        if(j>=xis.size() || xis[j]>xi) cmp_xis.push_back(xi);
-    }
-    assert( xis.size()+cmp_xis.size() == n_states );
-    return cmp_xis;
-}
-
-Real max(Real x,Real y) { return x>y?x:y; }
-
-Real unSoftMax(Real lz, Real cmp_lz) {
-    double m = max(lz, cmp_lz);
-    lz -= m; cmp_lz -= m;
-    double p = exp(lz)/(exp(lz)+exp(cmp_lz));
-    return p;
-}
-
-Real logSumExp(Real lz, Real cmp_lz) {
-    double m = max(lz, cmp_lz);
-    lz -= m; cmp_lz -= m;
-    return m+log(exp(lz)+exp(cmp_lz));
-}
-
-Real dist(const vector<Factor>& b1, const vector<Factor>& b2, size_t nv) {
-    Real d=0.0;
-    for(size_t k=0; k<nv; k++) {
-        d += dist( b1[k], b2[k], Prob::DISTLINF );
-    }
-    return d;
+    bpProps.Set("verbose", props.verbose);
+    bpProps.Set("logdomain", false);
+    bpProps.Set("damping", 0.0);
+    BP *bp = new BP(*this,bpProps);
+    bp->recordSentMessages = true;
+    bp->init();
+    return bp;
 }
 
 
-void CBP::runRecurse(BP_dual &bp_dual,
+void CBP::runRecurse(InfAlg *bp,
                      double orig_logZ,
                      vector<size_t> clamped_vars_list,
-                     set<size_t> clamped_vars,
                      size_t &num_leaves,
+                     size_t &choose_count,
                      double &sum_level,
                      Real &lz_out, vector<Factor>& beliefs_out) {
     // choose a variable/states to clamp:
@@ -212,29 +137,31 @@ void CBP::runRecurse(BP_dual &bp_dual,
     vector<size_t> xis;
     Real maxVar=0.0;
     bool found;
-    bool clampVar = (Clamping()==Properties::ClampType::CLAMP_VAR);
-
-    // XXX fix to just pass orig_logZ
+    bool clampingVar = (Clamping()==Properties::ClampType::CLAMP_VAR);
 
     if(Recursion()==Properties::RecurseType::REC_LOGZ && recTol()>0 &&
-       exp(bp_dual.logZ()-orig_logZ) < recTol()) {
+       exp(bp->logZ()-orig_logZ) < recTol()) {
         found = false;
     } else {
-        found = chooseNextClampVar(bp_dual,
+        found = chooseNextClampVar(bp,
                                    clamped_vars_list,
-                                   clamped_vars,
                                    i, xis, &maxVar);
     }
 
     if(!found) {
         num_leaves++;
         sum_level += clamped_vars_list.size();
-        beliefs_out = bp_dual.beliefs();
-        lz_out = bp_dual.logZ();
+        beliefs_out = bp->beliefs();
+        lz_out = bp->logZ();
         return;
     }
+
+    choose_count++;
+    if(props.clamp_outfile.length()>0) {
+        *_clamp_ofstream << choose_count << "\t" << clamped_vars_list.size() << "\t" << i << "\t" << xis[0] << endl;
+    }
     
-    if(clampVar) {
+    if(clampingVar) {
         foreach(size_t xi, xis) { assert(/*0<=xi &&*/ xi<var(i).states()); }
     } else {
         foreach(size_t xI, xis) { assert(/*0<=xI &&*/ xI<factor(i).states()); }
@@ -243,39 +170,39 @@ void CBP::runRecurse(BP_dual &bp_dual,
     // clamp setting. afterwards, combine estimates.
 
     // compute complement of 'xis'
-    vector<size_t> cmp_xis=complement(xis, clampVar?var(i).states():factor(i).states());
+    vector<size_t> cmp_xis=complement(xis, clampingVar?var(i).states():factor(i).states());
 
-    // XXX could do this more efficiently with a nesting version of
-    // saveProbs/undoProbs
+    /// \todo could do this more efficiently with a nesting version of
+    /// saveProbs/undoProbs
     Real lz; vector<Factor> b;
-    BP_dual bp_dual_c(bp_dual);
-    if(clampVar) {
-        _clamp((FactorGraph&)bp_dual_c, var(i), xis);
-        bp_dual_c.init(var(i));
+    InfAlg *bp_c = bp->clone();
+    if(clampingVar) {
+        bp_c->fg().clampVar(i, xis);
+        bp_c->init(var(i));
     } else {
-        _clampFactor((FactorGraph&)bp_dual_c, i, xis);
-        bp_dual_c.init(factor(i).vars());
+        bp_c->fg().clampFactor(i, xis);
+        bp_c->init(factor(i).vars());
     }
-    bp_dual_c.run();
-    _iters += bp_dual_c.Iterations();
+    bp_c->run();
+    _iters += bp_c->Iterations();
 
-    lz = bp_dual_c.logZ();
-    b = bp_dual_c.beliefs();
+    lz = bp_c->logZ();
+    b = bp_c->beliefs();
 
     Real cmp_lz; vector<Factor> cmp_b;
-    BP_dual cmp_bp_dual_c(bp_dual);
-    if(clampVar) {
-        _clamp(cmp_bp_dual_c,var(i),cmp_xis);
-        cmp_bp_dual_c.init(var(i));
+    InfAlg *cmp_bp_c = bp->clone();
+    if(clampingVar) {
+        cmp_bp_c->fg().clampVar(i,cmp_xis);
+        cmp_bp_c->init(var(i));
     } else {
-        _clampFactor(cmp_bp_dual_c,i,cmp_xis);
-        cmp_bp_dual_c.init(factor(i).vars());
+        cmp_bp_c->fg().clampFactor(i,cmp_xis);
+        cmp_bp_c->init(factor(i).vars());
     }
-    cmp_bp_dual_c.run();
-    _iters += cmp_bp_dual_c.Iterations();
+    cmp_bp_c->run();
+    _iters += cmp_bp_c->Iterations();
 
-    cmp_lz = cmp_bp_dual_c.logZ();
-    cmp_b = cmp_bp_dual_c.beliefs();
+    cmp_lz = cmp_bp_c->logZ();
+    cmp_b = cmp_bp_c->beliefs();
 
     double p = unSoftMax(lz, cmp_lz);
     Real bp__d=0.0;
@@ -283,7 +210,7 @@ void CBP::runRecurse(BP_dual &bp_dual,
     if(Recursion()==Properties::RecurseType::REC_BDIFF && recTol() > 0) {
         vector<Factor> combined_b(mixBeliefs(p,b,cmp_b));
         Real new_lz = logSumExp(lz,cmp_lz);
-        bp__d = dist(bp_dual.beliefs(),combined_b,nrVars());
+        bp__d = dist(bp->beliefs(),combined_b,nrVars());
         if(exp(new_lz-orig_logZ)*bp__d < recTol()) {
             num_leaves++;
             sum_level += clamped_vars_list.size();
@@ -295,15 +222,12 @@ void CBP::runRecurse(BP_dual &bp_dual,
 
     // either we are not doing REC_BDIFF or the distance was large
     // enough to recurse:
-
-    runRecurse(bp_dual_c, orig_logZ, 
+    runRecurse(bp_c, orig_logZ, 
                clamped_vars_list,
-               clamped_vars,
-               num_leaves, sum_level, lz, b);
-    runRecurse(cmp_bp_dual_c, orig_logZ, 
+               num_leaves, choose_count, sum_level, lz, b);
+    runRecurse(cmp_bp_c, orig_logZ, 
                clamped_vars_list,
-               clamped_vars,
-               num_leaves, sum_level, cmp_lz, cmp_b);
+               num_leaves, choose_count, sum_level, cmp_lz, cmp_b);
 
     p = unSoftMax(lz,cmp_lz);
 
@@ -311,18 +235,22 @@ void CBP::runRecurse(BP_dual &bp_dual,
     lz_out = logSumExp(lz,cmp_lz);
 
     if(props.verbose>=2) {
-        Real d = dist( bp_dual.beliefs(), beliefs_out, nrVars() );
+        Real d = dist( bp->beliefs(), beliefs_out, nrVars() );
         cerr << "Distance (clamping " << i << "): " << d;
         if(Recursion()==Properties::RecurseType::REC_BDIFF)
             cerr << "; bp_dual predicted " << bp__d;
-        cerr << "; max adjoint = " << maxVar << "; level = " << clamped_vars_list.size() << endl;
+        fprintf(stderr, "; max_adjoint = %f; logZ = %f (in %f) (orig %f); p = %f; level = %zd\n",
+                maxVar, lz_out, bp->logZ(), orig_logZ, p,
+                clamped_vars_list.size());
     }
+
+    delete bp_c; delete cmp_bp_c;
 }
 
+
 // 'xis' must be sorted
-bool CBP::chooseNextClampVar(BP_dual& bp,
+bool CBP::chooseNextClampVar(InfAlg* bp,
                              vector<size_t> &clamped_vars_list,
-                             set<size_t> &clamped_vars,
                              size_t &i, vector<size_t> &xis, Real *maxVarOut) {
     Real tiny=1.0e-14;
     if(props.verbose>=3) {
@@ -335,9 +263,9 @@ bool CBP::chooseNextClampVar(BP_dual& bp,
         if(Clamping()==Properties::ClampType::CLAMP_VAR) {
             int t=0, t1=100;
             do {
-                i = rnd_multi(nrVars());
+                i = rnd(nrVars());
                 t++;
-            } while(abs(bp.belief1(i).p().max()-1) < tiny &&
+            } while(abs(bp->beliefV(i).p().max()-1) < tiny &&
                     t < t1);
             if(t==t1) {
                 return false;
@@ -346,19 +274,20 @@ bool CBP::chooseNextClampVar(BP_dual& bp,
             // only pick probable values for variable
             size_t xi;
             do {
-                xi = rnd_multi(var(i).states());
+                xi = rnd(var(i).states());
                 t++;
-            } while(bp.belief1(i).p()[xi] < tiny && t<t1);
+            } while(bp->beliefV(i).p()[xi] < tiny && t<t1);
             assert(t<t1);
             xis.resize(1, xi);
             //         assert(!_clamped_vars.count(i)); // not true for >2-ary variables
-            DAI_IFVERB(2, endl<<"CHOOSE_RANDOM chose variable "<<i<<" state "<<xis[0]<<endl);
+            DAI_IFVERB(2, "CHOOSE_RANDOM at level " << clamped_vars_list.size()
+                       << " chose variable " << i << " state " << xis[0] << endl);
         } else {
             int t=0, t1=100;
             do {
-                i = rnd_multi(nrFactors());
+                i = rnd(nrFactors());
                 t++;
-            } while(abs(bp.belief2(i).p().max()-1) < tiny &&
+            } while(abs(bp->beliefF(i).p().max()-1) < tiny &&
                     t < t1);
             if(t==t1) {
                 return false;
@@ -367,179 +296,178 @@ bool CBP::chooseNextClampVar(BP_dual& bp,
             // only pick probable values for variable
             size_t xi;
             do {
-                xi = rnd_multi(factor(i).states());
+                xi = rnd(factor(i).states());
                 t++;
-            } while(bp.belief2(i).p()[xi] < tiny && t<t1);
+            } while(bp->beliefF(i).p()[xi] < tiny && t<t1);
             assert(t<t1);
             xis.resize(1, xi);
             //         assert(!_clamped_vars.count(i)); // not true for >2-ary variables
             DAI_IFVERB(2, endl<<"CHOOSE_RANDOM chose factor "<<i<<" state "<<xis[0]<<endl);
         }
-    } else if(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_ALL) {
-      // try clamping each variable manually
-      assert(Clamping()==Properties::ClampType::CLAMP_VAR);
-      Real max_diff=0.0;
-      int win_k=-1, win_xk=-1;
-      for(size_t k=0; k<nrVars(); k++) {
-        for(size_t xk=0; xk<var(k).states(); xk++) {
-          if(bp.belief1(k)[xk]<tiny) continue;
-          BP_dual bp1(bp);
-          bp1.clamp(var(k), xk);
-          bp1.init(var(k));
-          bp1.run();
-          Real diff=0;
-          for(size_t j=0; j<nrVars(); j++) {
-            diff += dist(bp.belief1(j), bp1.belief1(j), Prob::DISTL1);
-          }
-          if(diff>max_diff) {
-            max_diff=diff; win_k=k; win_xk=xk;
-          }
+    } else if(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_MAXENT) {
+        if(Clamping()==Properties::ClampType::CLAMP_VAR) {
+            Real max_ent=-1.0;
+            int win_k=-1, win_xk=-1;
+            for(size_t k=0; k<nrVars(); k++) {
+                Real ent=bp->beliefV(k).entropy();
+                if(max_ent < ent) {
+                    max_ent = ent;
+                    win_k = k;
+                    win_xk = bp->beliefV(k).p().argmax().first;
+                }
+            }
+            assert(win_k>=0); assert(win_xk>=0);
+            i = win_k; xis.resize(1, win_xk);
+            DAI_IFVERB(2, endl<<"CHOOSE_MAXENT chose variable "<<i<<" state "<<xis[0]<<endl);
+            if(bp->beliefV(i).p()[xis[0]] < tiny) {
+                DAI_IFVERB(2, "Warning: CHOOSE_MAXENT found unlikely state, not recursing");
+                return false;
+            }
+        } else {
+            Real max_ent=-1.0;
+            int win_k=-1, win_xk=-1;
+            for(size_t k=0; k<nrFactors(); k++) {
+                Real ent=bp->beliefF(k).entropy();
+                if(max_ent < ent) {
+                    max_ent = ent;
+                    win_k = k;
+                    win_xk = bp->beliefF(k).p().argmax().first;
+                }
+            }
+            assert(win_k>=0); assert(win_xk>=0);
+            i = win_k; xis.resize(1, win_xk);
+            DAI_IFVERB(2, endl<<"CHOOSE_MAXENT chose factor "<<i<<" state "<<xis[0]<<endl);
+            if(bp->beliefF(i).p()[xis[0]] < tiny) {
+                DAI_IFVERB(2, "Warning: CHOOSE_MAXENT found unlikely state, not recursing");
+                return false;
+            }
+        }
+    } else if(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_L1 ||
+              ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_CFN) {
+        bool doL1=(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_L1);
+        vector<size_t> state;
+        if(!doL1 && needGibbsState(BBP_cost_function())) {
+            state = getGibbsState(*bp,2*bp->Iterations());
         }
-      }
-      assert(win_k>=0); assert(win_xk>=0);
-      i = win_k; xis.resize(1, win_xk);
+        // try clamping each variable manually
+        assert(Clamping()==Properties::ClampType::CLAMP_VAR);
+        Real max_cost=0.0;
+        int win_k=-1, win_xk=-1;
+        for(size_t k=0; k<nrVars(); k++) {
+            for(size_t xk=0; xk<var(k).states(); xk++) {
+                if(bp->beliefV(k)[xk]<tiny) continue;
+                InfAlg *bp1 = bp->clone();
+                bp1->clamp(var(k), xk);
+                bp1->init(var(k));
+                bp1->run();
+                Real cost=0;
+                if(doL1) {
+                    for(size_t j=0; j<nrVars(); j++) {
+                        cost += dist(bp->beliefV(j), bp1->beliefV(j), Prob::DISTL1);
+                    }
+                } else {
+                    cost = getCostFn(*bp1, BBP_cost_function(), &state);
+                }
+                if(cost>max_cost || win_k==-1) {
+                    max_cost=cost; win_k=k; win_xk=xk;
+                }
+                delete bp1;
+            }
+        }
+        assert(win_k>=0); assert(win_xk>=0);
+        i = win_k; xis.resize(1, win_xk);
     } else if(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BBP) {
         Real mvo; if(!maxVarOut) maxVarOut = &mvo;
-        bool clampVar = (Clamping()==Properties::ClampType::CLAMP_VAR);
+        bool clampingVar = (Clamping()==Properties::ClampType::CLAMP_VAR);
         pair<size_t, size_t> cv =
-            bbpFindClampVar(bp,
-                            clampVar,
-                            clamped_vars_list.size(),
-                            BBP_cost_function(),getProperties(),maxVarOut);
+            bbpFindClampVar(*bp,
+                            clampingVar,
+                            props.bbp_props,
+                            BBP_cost_function(),&mvo);
 
         // if slope isn't big enough then don't clamp
-        if(*maxVarOut < minMaxAdj()) return false;
+        if(mvo < minMaxAdj()) return false;
 
         size_t xi=cv.second;
         i = cv.first;
-#define VAR_INFO (clampVar?"variable ":"factor ")                       \
+#define VAR_INFO (clampingVar?"variable ":"factor ")                       \
             << i << " state " << xi                                     \
-            << " (p=" << (clampVar?bp.belief1(i)[xi]:bp.belief2(i)[xi]) \
-            << ", entropy = " << (clampVar?bp.belief1(i):bp.belief2(i)).entropy() \
-            << ", maxVar = "<< *maxVarOut << ")" 
-        Prob b = (clampVar?bp.belief1(i).p():bp.belief2(i).p());
+            << " (p=" << (clampingVar?bp->beliefV(i)[xi]:bp->beliefF(i)[xi]) \
+            << ", entropy = " << (clampingVar?bp->beliefV(i):bp->beliefF(i)).entropy() \
+            << ", maxVar = "<< mvo << ")" 
+        Prob b = (clampingVar?bp->beliefV(i).p():bp->beliefF(i).p());
         if(b[xi] < tiny) {
-            cerr << "Warning, bbpFindClampVar found unlikely "
+            cerr << "Warning, at level "
+                 << clamped_vars_list.size()
+                 << ", bbpFindClampVar found unlikely "
                  << VAR_INFO << endl;
             return false;
         }
         if(abs(b[xi]-1) < tiny) {
-            cerr << "Warning, bbpFindClampVar found overly likely "
+            cerr << "Warning at level " 
+                 << clamped_vars_list.size()
+                 << ", bbpFindClampVar found overly likely "
                  << VAR_INFO << endl;
             return false;
         }
 
         xis.resize(1,xi);
-        if(clampVar) {
+        if(clampingVar) {
             assert(/*0<=xi &&*/ xi<var(i).states());
         } else {
             assert(/*0<=xi &&*/ xi<factor(i).states());
         }
         DAI_IFVERB(2, "CHOOSE_BBP (num clamped = " << clamped_vars_list.size()
-               << ") chose " << i << " state " << xi << endl);
+                   << ") chose " << i << " state " << xi << endl);
     } else {
         abort();
     }
     clamped_vars_list.push_back(i);
-    clamped_vars.insert(i);
     return true;
 }
 
-// void CBP::clamp(const Var & n, size_t i) {
-//   FactorGraph::clamp(n,i);
-//   _clamped_vars.insert(findVar(n));
-//   _clamped_vars_list.push_back(findVar(n));
-// }
-
-// void CBP::clamp(const Var & n, const vector<size_t> &is) {
-//   FactorGraph::clamp(n,is);
-//   _clamped_vars.insert(findVar(n));
-//   _clamped_vars_list.push_back(findVar(n));
-// }
 
 void CBP::printDebugInfo() {
-    DAI_PV(_beliefs1);
-    DAI_PV(_beliefs2);
+    DAI_PV(_beliefsV);
+    DAI_PV(_beliefsF);
     DAI_PV(_logZ);
 }
 
+
 //----------------------------------------------------------------
 
-bool doBBPTest=false;
-bool doBBPGraph=false;
-size_t bbpGraphLevel=3;
-
-#define BPP_INIT_GIBBS 1
-
-/// function which takes a factor graph as input, runs Gibbs and BP_dual,
-/// creates and runs a BBP object, finds best variable, returns
-/// (variable,state) pair for clamping
-// pair<size_t, size_t> bbpFindClampVar(const CBP &fg, bbp_cfn_t cfn, const Properties &props, Real *maxVarOut) {
-pair<size_t, size_t> bbpFindClampVar(BP_dual &in_bp_dual, bool clampVar,
-    size_t numClamped, bbp_cfn_t cfn, const PropertySet &props, Real *maxVarOut) {
-#if BPP_INIT_GIBBS
-    vector<size_t> state = getGibbsState(in_bp_dual, 100);
-    in_bp_dual.init(state);
-    in_bp_dual.run();
-#endif
-  
-    Real ourTol = doBBPTest ? 1.0e-11 : 1.0e-3;
-    if(0) {
-        PropertySet bp_Props;
-        bp_Props.Set("updates", string("PARALL"));
-        //   bp_Props.Set("tol", props.GetAs<double>("tol"));
-        bp_Props.Set("tol", ourTol);
-        bp_Props.Set("maxiter", props.GetAs<size_t>("maxiter"));
-        bp_Props.Set("verbose", oneLess(props.GetAs<size_t>("verbose")));
-        //   bp_Props.ConvertTo<BP_dual::UpdateType>("updates");
-        //   DAI_PV(bp_Props.GetAs<BP_dual::UpdateType>("updates"));
-        BP_dual bp_dual(in_bp_dual, bp_Props);
-#if BPP_INIT_GIBBS
-        bp_dual.init(state);
-#endif
-        bp_dual.run();
-    }
-
-    if(doBBPGraph && numClamped == bbpGraphLevel) {
-        cerr << "Writing BBP graph data" << endl;
-        makeBBPGraph(in_bp_dual,cfn);
-        doBBPGraph=false; // only do it once
-        cerr << "Done writing BBP graph data" << endl;
-    }
-    if(doBBPTest) {
-        double err = numericBBPTest(in_bp_dual, cfn, /*bbp tol*/ ourTol, /*h*/ 1.0e-5);
-        cerr << "Error from numericBBPTest: " << err << endl;
-    }
-    Real tic1=toc();
-    BBP bbp(in_bp_dual);
-    bbp.maxIter() = props.GetAs<size_t>("maxiter");
-#if BPP_INIT_GIBBS
-    gibbsInitBBPCostFnAdj(bbp, in_bp_dual, cfn, &state);
-#else
-    gibbsInitBBPCostFnAdj(bbp, in_bp_dual, cfn, NULL);
-#endif
-    Real tic2=toc();
-    bbp.run(ourTol);
-    if(props.GetAs<size_t>("verbose") >= 3) {
-        cerr << "BBP took " << toc()-tic1 << " seconds (BBP.run = " << toc()-tic2 << " seconds), "
-             << bbp.doneIters() << " iterations" << endl;
-    }
+
+/** function which takes a factor graph as input, runs Gibbs and BP_dual,
+ *  creates and runs a BBP object, finds best variable, returns
+ *  (variable,state) pair for clamping
+ */
+pair<size_t, size_t> bbpFindClampVar(const InfAlg &in_bp, bool clampingVar,
+            const PropertySet &bbp_props, bbp_cfn_t cfn, 
+            Real *maxVarOut)
+{
+    BBP bbp(&in_bp,bbp_props);
+    initBBPCostFnAdj(bbp, in_bp, cfn, NULL);
+    bbp.run();
   
-    // find and return the (variable,state) with the largest adj_psi_1
+    // find and return the (variable,state) with the largest adj_psi_V
     size_t argmax_var=0;
     size_t argmax_var_state=0;
     Real max_var=0;
-    if(clampVar) {
-        for(size_t i=0; i<in_bp_dual.nrVars(); i++) {
-            Prob adj_psi_1 = bbp.adj_psi_1(i);
+    if(clampingVar) {
+        for(size_t i=0; i<in_bp.fg().nrVars(); i++) {
+            Prob adj_psi_V = bbp.adj_psi_V(i);
             if(0) {
                 // helps to account for amount of movement possible in variable
                 // i's beliefs? seems not..
-                adj_psi_1 *= in_bp_dual.belief1(i).entropy();
+                adj_psi_V *= in_bp.beliefV(i).entropy();
+            }
+            if(0) {
+//                 adj_psi_V *= Prob(in_bp.fg().var(i).states(),1.0)-in_bp.beliefV(i).p();
+                adj_psi_V *= in_bp.beliefV(i).p();
             }
             // try to compensate for effect on same variable (doesn't work)
-            //     adj_psi_1[gibbs.state()[i]] -= bp_dual.belief1(i)[gibbs.state()[i]]/10;
-            pair<size_t,Real> argmax_state = adj_psi_1.argmax();
+            //     adj_psi_V[gibbs.state()[i]] -= bp_dual.beliefV(i)[gibbs.state()[i]]/10;
+            pair<size_t,Real> argmax_state = adj_psi_V.argmax();
 
             if(i==0 || argmax_state.second>max_var) {
                 argmax_var = i;
@@ -548,18 +476,18 @@ pair<size_t, size_t> bbpFindClampVar(BP_dual &in_bp_dual, bool clampVar,
             }
         }
         assert(/*0 <= argmax_var_state &&*/
-               argmax_var_state < in_bp_dual.var(argmax_var).states());
+               argmax_var_state < in_bp.fg().var(argmax_var).states());
     } else {
-        for(size_t I=0; I<in_bp_dual.nrFactors(); I++) {
-            Prob adj_psi_2 = bbp.adj_psi_2(I);
+        for(size_t I=0; I<in_bp.fg().nrFactors(); I++) {
+            Prob adj_psi_F = bbp.adj_psi_F(I);
             if(0) {
                 // helps to account for amount of movement possible in variable
                 // i's beliefs? seems not..
-                adj_psi_2 *= in_bp_dual.belief2(I).entropy();
+                adj_psi_F *= in_bp.beliefF(I).entropy();
             }
             // try to compensate for effect on same variable (doesn't work)
-            //     adj_psi_1[gibbs.state()[i]] -= bp_dual.belief1(i)[gibbs.state()[i]]/10;
-            pair<size_t,Real> argmax_state = adj_psi_2.argmax();
+            //     adj_psi_V[gibbs.state()[i]] -= bp_dual.beliefV(i)[gibbs.state()[i]]/10;
+            pair<size_t,Real> argmax_state = adj_psi_F.argmax();
 
             if(I==0 || argmax_state.second>max_var) {
                 argmax_var = I;
@@ -568,10 +496,196 @@ pair<size_t, size_t> bbpFindClampVar(BP_dual &in_bp_dual, bool clampVar,
             }
         }
         assert(/*0 <= argmax_var_state &&*/
-               argmax_var_state < in_bp_dual.factor(argmax_var).states());
+               argmax_var_state < in_bp.fg().factor(argmax_var).states());
     }
     if(maxVarOut) *maxVarOut = max_var;
     return make_pair(argmax_var,argmax_var_state);
 }
 
+
+vector<size_t> complement( vector<size_t>& xis, size_t n_states ) {
+    vector<size_t> cmp_xis(0);
+    size_t j = 0;
+    for( size_t xi = 0; xi < n_states; xi++ ) {
+        while( j < xis.size() && xis[j] < xi )
+            j++;
+        if( j >= xis.size() || xis[j] > xi )
+            cmp_xis.push_back(xi);
+    }
+    assert( xis.size()+cmp_xis.size() == n_states );
+    return cmp_xis;
+}
+
+
+Real unSoftMax( Real a, Real b ) {
+    if( a > b )
+        return 1.0 / (1.0 + exp(b-a));
+    else
+        return exp(a-b) / (exp(a-b) + 1.0);
+}
+
+
+Real logSumExp( Real a, Real b ) {
+    if( a > b )
+        return a + log1p(exp(b-a));
+    else
+        return b + log1p(exp(a-b));
+}
+
+
+Real dist( const vector<Factor>& b1, const vector<Factor>& b2, size_t nv ) {
+    Real d = 0.0;
+    for( size_t k = 0; k < nv; k++ )
+        d += dist( b1[k], b2[k], Prob::DISTLINF );
+    return d;
+}
+
+
+} // end of namespace dai
+
+
+/* {{{ GENERATED CODE: DO NOT EDIT. Created by 
+    ./scripts/regenerate-properties include/dai/cbp.h src/cbp.cpp 
+*/
+namespace dai {
+
+void CBP::Properties::set(const PropertySet &opts)
+{
+    const std::set<PropertyKey> &keys = opts.allKeys();
+    std::set<PropertyKey>::const_iterator i;
+    bool die=false;
+    for(i=keys.begin(); i!=keys.end(); i++) {
+        if(*i == "verbose") continue;
+        if(*i == "tol") continue;
+        if(*i == "updates") continue;
+        if(*i == "maxiter") continue;
+        if(*i == "rec_tol") continue;
+        if(*i == "max_levels") continue;
+        if(*i == "min_max_adj") continue;
+        if(*i == "choose") continue;
+        if(*i == "recursion") continue;
+        if(*i == "clamp") continue;
+        if(*i == "bbp_props") continue;
+        if(*i == "bbp_cfn") continue;
+        if(*i == "rand_seed") continue;
+        if(*i == "clamp_outfile") continue;
+        cerr << "CBP: Unknown property " << *i << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(UNKNOWN_PROPERTY_TYPE);
+    }
+    if(!opts.hasKey("tol")) {
+        cerr << "CBP: Missing property \"tol\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("updates")) {
+        cerr << "CBP: Missing property \"updates\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("maxiter")) {
+        cerr << "CBP: Missing property \"maxiter\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("rec_tol")) {
+        cerr << "CBP: Missing property \"rec_tol\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("min_max_adj")) {
+        cerr << "CBP: Missing property \"min_max_adj\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("choose")) {
+        cerr << "CBP: Missing property \"choose\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("recursion")) {
+        cerr << "CBP: Missing property \"recursion\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("clamp")) {
+        cerr << "CBP: Missing property \"clamp\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("bbp_props")) {
+        cerr << "CBP: Missing property \"bbp_props\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(!opts.hasKey("bbp_cfn")) {
+        cerr << "CBP: Missing property \"bbp_cfn\" for method \"CBP\"" << endl;
+        die=true;
+    }
+    if(die) {
+        DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
+    }
+    if(opts.hasKey("verbose")) {
+        verbose = opts.getStringAs<size_t>("verbose");
+    } else {
+        verbose = 0;
+    }
+    tol = opts.getStringAs<double>("tol");
+    updates = opts.getStringAs<UpdateType>("updates");
+    maxiter = opts.getStringAs<size_t>("maxiter");
+    rec_tol = opts.getStringAs<double>("rec_tol");
+    if(opts.hasKey("max_levels")) {
+        max_levels = opts.getStringAs<size_t>("max_levels");
+    } else {
+        max_levels = 10;
+    }
+    min_max_adj = opts.getStringAs<double>("min_max_adj");
+    choose = opts.getStringAs<ChooseMethodType>("choose");
+    recursion = opts.getStringAs<RecurseType>("recursion");
+    clamp = opts.getStringAs<ClampType>("clamp");
+    bbp_props = opts.getStringAs<PropertySet>("bbp_props");
+    bbp_cfn = opts.getStringAs<bbp_cfn_t>("bbp_cfn");
+    if(opts.hasKey("rand_seed")) {
+        rand_seed = opts.getStringAs<size_t>("rand_seed");
+    } else {
+        rand_seed = 0;
+    }
+    if(opts.hasKey("clamp_outfile")) {
+        clamp_outfile = opts.getStringAs<std::string>("clamp_outfile");
+    } else {
+        clamp_outfile = "";
+    }
+}
+PropertySet CBP::Properties::get() const {
+    PropertySet opts;
+    opts.Set("verbose", verbose);
+    opts.Set("tol", tol);
+    opts.Set("updates", updates);
+    opts.Set("maxiter", maxiter);
+    opts.Set("rec_tol", rec_tol);
+    opts.Set("max_levels", max_levels);
+    opts.Set("min_max_adj", min_max_adj);
+    opts.Set("choose", choose);
+    opts.Set("recursion", recursion);
+    opts.Set("clamp", clamp);
+    opts.Set("bbp_props", bbp_props);
+    opts.Set("bbp_cfn", bbp_cfn);
+    opts.Set("rand_seed", rand_seed);
+    opts.Set("clamp_outfile", clamp_outfile);
+    return opts;
+}
+string CBP::Properties::toString() const {
+    stringstream s(stringstream::out);
+    s << "[";
+    s << "verbose=" << verbose << ",";
+    s << "tol=" << tol << ",";
+    s << "updates=" << updates << ",";
+    s << "maxiter=" << maxiter << ",";
+    s << "rec_tol=" << rec_tol << ",";
+    s << "max_levels=" << max_levels << ",";
+    s << "min_max_adj=" << min_max_adj << ",";
+    s << "choose=" << choose << ",";
+    s << "recursion=" << recursion << ",";
+    s << "clamp=" << clamp << ",";
+    s << "bbp_props=" << bbp_props << ",";
+    s << "bbp_cfn=" << bbp_cfn << ",";
+    s << "rand_seed=" << rand_seed << ",";
+    s << "clamp_outfile=" << clamp_outfile;
+    s << "]";
+    return s.str();
+}
 } // end of namespace dai
+/* }}} END OF GENERATED CODE */
index c7f62a1..e40ffbd 100644 (file)
@@ -109,5 +109,5 @@ BP_dual_SEQFIX:                 BP_dual[updates=SEQFIX,tol=1e-9,maxiter=10000,ve
 
 # --- CBP ---------------------
 
-CBP:                            CBP[max_levels=12,updates=SEQMAX,tol=1e-9,rec_tol=1e-9,maxiter=100,choose=CHOOSE_RANDOM,recursion=REC_FIXED,clamp=CLAMP_VAR,min_max_adj=1.0e-9,bbp_cfn=cfn_b,verbose=0,rand_seed=0]
+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=]
 BBP:                            CBP[choose=CHOOSE_BBP]