[Frederik Eaton] Added BP_Dual, BBP and CBP algorithms
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Mar 2009 13:49:10 +0000 (14:49 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Mar 2009 13:53:58 +0000 (14:53 +0100)
17 files changed:
Makefile
Makefile.CYGWIN
Makefile.LINUX
Makefile.MACOSX
Makefile.WINDOWS
include/dai/alldai.h
include/dai/bbp.h [new file with mode: 0644]
include/dai/bipgraph.h
include/dai/bp_dual.h [new file with mode: 0644]
include/dai/cbp.h [new file with mode: 0644]
include/dai/factorgraph.h
include/dai/util.h
src/alldai.cpp
src/bbp.cpp [new file with mode: 0644]
src/bp_dual.cpp [new file with mode: 0644]
src/cbp.cpp [new file with mode: 0644]
tests/aliases.conf

index bdfbdf1..f07f56b 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -78,9 +78,18 @@ ifdef WITH_GIBBS
   CCFLAGS:=$(CCFLAGS) -DDAI_WITH_GIBBS
   OBJECTS:=$(OBJECTS) gibbs$(OE)
 endif
+ifdef WITH_BP_DUAL
+  CCFLAGS:=$(CCFLAGS) -DDAI_WITH_BP_DUAL
+  OBJECTS:=$(OBJECTS) bp_dual$(OE)
+endif
+ifdef WITH_CBP
+  CCFLAGS:=$(CCFLAGS) -DDAI_WITH_CBP
+  OBJECTS:=$(OBJECTS) bbp$(OE) cbp$(OE)
+endif
+
 
 # Define standard libDAI header dependencies
-HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h
+HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h $(INC)/util.h
 
 # Setup final command for C++ compiler and MEX
 ifdef DEBUG
@@ -126,6 +135,15 @@ exactinf$(OE) : $(SRC)/exactinf.cpp $(INC)/exactinf.h $(HEADERS)
 bp$(OE) : $(SRC)/bp.cpp $(INC)/bp.h $(HEADERS)
        $(CC) -c $(SRC)/bp.cpp
 
+bp_dual$(OE) : $(SRC)/bp_dual.cpp $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bp_dual.cpp
+
+bbp$(OE) : $(SRC)/bbp.cpp $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/bbp.cpp
+
+cbp$(OE) : $(SRC)/cbp.cpp $(INC)/cbp.h $(INC)/bbp.h $(INC)/bp_dual.h $(HEADERS)
+       $(CC) -c $(SRC)/cbp.cpp
+
 lc$(OE) : $(SRC)/lc.cpp $(INC)/lc.h $(HEADERS)
        $(CC) -c $(SRC)/lc.cpp
 
index 3e2ce12..64916af 100644 (file)
@@ -15,6 +15,8 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_BP_DUAL=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index bbfde5c..4c9fa51 100644 (file)
@@ -16,6 +16,8 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_BP_DUAL=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index 9cd503e..9bcf11b 100644 (file)
@@ -15,6 +15,8 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_BP_DUAL=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=true
 # Build with debug info?
index c6ee945..9290701 100644 (file)
@@ -16,6 +16,8 @@ WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
 WITH_GIBBS=true
+WITH_BP_DUAL=true
+WITH_CBP=true
 # Build doxygen documentation?
 WITH_DOC=
 # Build with debug info?
index 430e8eb..91b4285 100644 (file)
 #ifdef DAI_WITH_GIBBS
     #include <dai/gibbs.h>
 #endif
+#ifdef DAI_WITH_BP_DUAL
+    #include <dai/bp_dual.h>
+#endif
+#ifdef DAI_WITH_CBP
+    #include <dai/cbp.h>
+#endif
 
 
 /// Namespace for libDAI
@@ -98,6 +104,9 @@ static const char* DAINames[] = {
 #endif
 #ifdef DAI_WITH_GIBBS
     Gibbs::Name,
+#endif
+#ifdef DAI_WITH_CBP
+    CBP::Name,
 #endif
     ""
 };
diff --git a/include/dai/bbp.h b/include/dai/bbp.h
new file mode 100644 (file)
index 0000000..b14b6a3
--- /dev/null
@@ -0,0 +1,165 @@
+#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/enum.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; }
+
+/// 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);
+
+vector<Prob> get_zero_adj_2(const BP_dual& bp_dual);
+vector<Prob> get_zero_adj_1(const BP_dual& bp_dual);
+
+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]; });
+  
+  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)));
+    }
+  }
+ 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)]; });
+};
+
+/// 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);
+
+/// 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);
+
+/// 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);
+
+/// 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);
+
+}
+
+#endif
index f202f99..b3a85b6 100644 (file)
@@ -101,9 +101,16 @@ class BipartiteGraph {
             std::vector<size_t> ind2;       // indices of nodes of type 2
         };
 
+        /// Support for some backwards compatibility with old interface
+        /** Call indexEdges() first to initialize these members
+         */
+        bool _edge_indexed;
+        std::vector<Edge> _edges;
+        hash_map<Edge,size_t> _vv2e;
+
     public:
         /// Default constructor (creates an empty bipartite graph)
-        BipartiteGraph() : _nb1(), _nb2() {}
+        BipartiteGraph() : _nb1(), _nb2(), _edge_indexed(false) {}
 
         /// Constructs BipartiteGraph from a range of edges.
         /** \tparam EdgeInputIterator Iterator that iterates over instances of BipartiteGraph::Edge.
@@ -113,7 +120,7 @@ class BipartiteGraph {
          *  \param end Points just beyond the last edge.
          */
         template<typename EdgeInputIterator>
-        BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
+        BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ), _edge_indexed(false) {
             construct( nr1, nr2, begin, end );
         }
 
@@ -302,6 +309,51 @@ class BipartiteGraph {
         /// Writes this BipartiteGraph to an output stream in GraphViz .dot syntax
         void printDot( std::ostream& os ) const;
 
+        // ----------------------------------------------------------------
+        // backwards compatibility layer
+        void indexEdges() {
+            _edges.clear();
+            _vv2e.clear();
+            size_t i=0;
+            foreach(const Neighbors &nb1s, _nb1) {
+                foreach(const Neighbor &n2, nb1s) {
+                    Edge e(i, n2.node);
+                    _edges.push_back(e);
+                }
+                i++;
+            }
+            sort(_edges.begin(), _edges.end()); // unnecessary?
+
+            i=0;
+            foreach(const Edge& e, _edges) {
+                _vv2e[e] = i++;
+            }
+
+            _edge_indexed = true;
+        }
+        
+        const Edge& edge(size_t e) const {
+            assert(_edge_indexed);
+            return _edges[e];
+        }
+        
+        const std::vector<Edge>& edges() const {
+            return _edges;
+        }
+        
+        size_t VV2E(size_t n1, size_t n2) const {
+            assert(_edge_indexed);
+            Edge e(n1,n2);
+            hash_map<Edge,size_t>::const_iterator i = _vv2e.find(e);
+            assert(i != _vv2e.end());
+            return i->second;
+        }
+
+        size_t nr_edges() const {
+            assert(_edge_indexed);
+            return _edges.size();
+        }
+
     private:
         /// Checks internal consistency
         void check() const;
diff --git a/include/dai/bp_dual.h b/include/dai/bp_dual.h
new file mode 100644 (file)
index 0000000..aca369a
--- /dev/null
@@ -0,0 +1,147 @@
+#ifndef ____defined_libdai_bp_dual_h__
+#define ____defined_libdai_bp_dual_h__
+
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+#include <dai/enum.h>
+#include <dai/bp.h>
+
+namespace dai {
+
+using namespace std;
+
+struct BP_dual_messages {
+  // messages:
+  // indexed by edge index (using VV2E)
+  vector<Prob> n;
+  vector<Real> Zn;
+  vector<Prob> m;
+  vector<Real> Zm;
+};
+
+struct BP_dual_beliefs {
+  // beliefs:
+  // indexed by node
+  vector<Prob> b1;
+  vector<Real> Zb1;
+  // indexed by factor
+  vector<Prob> b2;
+  vector<Real> Zb2;
+};
+
+void _clamp(FactorGraph &g, const Var & n, const vector<size_t> &is );
+
+// clamp a factor to have one of a set of values
+void _clampFactor(FactorGraph &g, size_t I, const vector<size_t> &is);
+
+class BP_dual : public DAIAlgFG {
+ public:
+        typedef vector<size_t>  _ind_t;
+
+ protected:
+
+        // indexed by edge index. for each edge i->I, contains a
+        // vector whose entries correspond to those of I, and the
+        // value of each entry is the corresponding entry of i
+        vector<_ind_t>          _indices; 
+
+        BP_dual_messages _msgs;
+        BP_dual_messages _new_msgs;
+ public:
+        BP_dual_beliefs _beliefs;
+
+        size_t _iters;
+        double _maxdiff;
+
+        struct Properties {
+          typedef BP::Properties::UpdateType UpdateType;
+          UpdateType updates;
+          double tol;
+          size_t maxiter;
+          size_t verbose;
+        } props;
+
+        /// List of property names
+        static const char *PropertyList[];
+        /// Name of this inference algorithm
+        static const char *Name;
+
+    public:
+        void Regenerate(); // used by constructor
+        void RegenerateIndices();
+        void RegenerateMessages();
+        void RegenerateBeliefs();
+
+        void CalcBelief1(size_t i);
+        void CalcBelief2(size_t I);
+        void CalcBeliefs(); // called after run()
+
+        void calcNewM(size_t iI);
+        void calcNewN(size_t iI);
+        void upMsgM(size_t iI);
+        void upMsgN(size_t iI);
+
+/*             DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL) */
+        typedef BP::Properties::UpdateType UpdateType;
+        UpdateType Updates() const { return props.updates; }
+        size_t Verbose() const { return props.verbose; }
+
+        /// construct BP_dual object from FactorGraph
+        BP_dual(const FactorGraph & fg, const PropertySet &opts) : DAIAlgFG(fg) {
+            setProperties(opts);
+            Regenerate();
+        }
+        
+        DAI_ACCMUT(Prob & msgM(size_t I, size_t i), { return _msgs.m[VV2E(i,I)]; });
+        DAI_ACCMUT(Prob & msgN(size_t i, size_t I), { return _msgs.n[VV2E(i,I)]; });
+        DAI_ACCMUT(Prob & msgM(size_t iI), { return _msgs.m[iI]; });
+        DAI_ACCMUT(Prob & msgN(size_t iI), { return _msgs.n[iI]; });
+        DAI_ACCMUT(Real & zM(size_t I, size_t i), { return _msgs.Zm[VV2E(i,I)]; });
+        DAI_ACCMUT(Real & zN(size_t i, size_t I), { return _msgs.Zn[VV2E(i,I)]; });
+        DAI_ACCMUT(Real & zM(size_t iI), { return _msgs.Zm[iI]; });
+        DAI_ACCMUT(Real & zN(size_t iI), { return _msgs.Zn[iI]; });
+        DAI_ACCMUT(Prob & newMsgM(size_t I, size_t i), { return _new_msgs.m[VV2E(i,I)]; });
+        DAI_ACCMUT(Prob & newMsgN(size_t i, size_t I), { return _new_msgs.n[VV2E(i,I)]; });
+        DAI_ACCMUT(Real & newZM(size_t I, size_t i), { return _new_msgs.Zm[VV2E(i,I)]; });
+        DAI_ACCMUT(Real & newZN(size_t i, size_t I), { return _new_msgs.Zn[VV2E(i,I)]; });
+
+        DAI_ACCMUT(_ind_t & index(size_t i, size_t I), { return( _indices[VV2E(i,I)] ); });
+
+        Real belief1Z(size_t i) const { return _beliefs.Zb1[i]; }
+        Real belief2Z(size_t I) const { return _beliefs.Zb2[I]; }
+
+        size_t doneIters() const { return _iters; }
+
+
+        /// @name General InfAlg interface
+        //@{
+        virtual BP_dual* clone() const { return new BP_dual(*this); }
+/*         virtual BP_dual* create() const { return new BP_dual(); } */
+        virtual BP_dual* create() const { return NULL; }
+        virtual std::string identify() const;
+        virtual Factor belief (const Var &n) const { return( belief1( findVar( n ) ) ); }
+        virtual Factor belief (const VarSet &n) const;
+        virtual std::vector<Factor> beliefs() const;
+        virtual Real logZ() const;
+        virtual void init();
+        virtual void init( const VarSet &ns );
+        virtual double run();
+        virtual double maxDiff() const { return _maxdiff; }
+        virtual size_t Iterations() const { return _iters; }
+        //@}
+
+        void init(const vector<size_t>& state);
+        Factor belief1 (size_t i) const { return Factor(var(i), _beliefs.b1[i]); }
+        Factor belief2 (size_t I) const { return Factor(factor(I).vars(), _beliefs.b2[I]); }
+
+        void updateMaxDiff( double maxdiff ) { if( maxdiff > _maxdiff ) _maxdiff = maxdiff; }
+
+        /// Set Props according to the PropertySet opts, where the values can be stored as std::strings or as the type of the corresponding Props member
+        void setProperties( const PropertySet &opts );
+        PropertySet getProperties() const;
+        std::string printProperties() const;
+};
+
+}
+
+#endif
diff --git a/include/dai/cbp.h b/include/dai/cbp.h
new file mode 100644 (file)
index 0000000..45a49f3
--- /dev/null
@@ -0,0 +1,155 @@
+#ifndef ____defined_libdai_cbp_h__
+#define ____defined_libdai_cbp_h__
+
+#include <dai/daialg.h>
+
+#include <dai/cbp.h>
+#include <dai/bbp.h>
+
+using namespace std;
+
+/* Class for "clamped belief propagation".
+
+'chooseNextClamp' can be overridden (e.g. in BPC); this version
+randomly chooses variables to clamp.
+*/
+
+namespace dai {
+
+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 CBP : public DAIAlgFG {
+    vector<Factor> _beliefs1;
+    vector<Factor> _beliefs2;
+    double _logZ;
+    double _est_logZ;
+    double _old_est_logZ;
+
+    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; }
+
+
+    bbp_cfn_t BBP_cost_function() {
+      return props.bbp_cfn;
+    }
+
+    void printDebugInfo();
+
+    void runRecurse(BP_dual &bp_dual,
+                    double orig_logZ,
+                    vector<size_t> clamped_vars_list,
+                    set<size_t> clamped_vars,
+                    size_t &num_leaves,
+                    double &sum_level,
+                    Real &lz_out, vector<Factor>& beliefs_out);
+
+    BP getBP();
+    BP_dual getBP_dual();
+
+  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;
+
+    Properties::ChooseMethodType ChooseMethod() { return props.choose; }
+    Properties::RecurseType Recursion() { return props.recursion; }
+    Properties::ClampType Clamping() { return props.clamp; }
+
+    //----------------------------------------------------------------
+
+    // construct CBP object from FactorGraph
+    CBP(const FactorGraph &fg, const PropertySet &opts) : DAIAlgFG(fg) {
+        setProperties(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 { throw "Unimplemented"; }
+    virtual std::string identify() const { return string(Name) + printProperties(); }
+    virtual Factor belief (const Var &n) const { return _beliefs1[findVar(n)]; }
+    virtual Factor belief (const VarSet &) const { assert(0); }
+    virtual vector<Factor> beliefs() const { return concat(_beliefs1, _beliefs2); }
+    virtual Real logZ() const { return _logZ; }
+    virtual void init() {};
+    virtual void init( const VarSet & ) {};
+    virtual double run();
+    virtual double maxDiff() const { return _maxdiff; }
+    virtual size_t Iterations() const { return _iters; }
+    //@}
+
+
+    //----------------------------------------------------------------
+
+    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);
+
+    //----------------------------------------------------------------
+
+    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]; }
+
+    //----------------------------------------------------------------
+
+    void construct();
+
+    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);
+
+}
+
+#endif
index a3f5d40..36e27e4 100644 (file)
@@ -276,6 +276,12 @@ class FactorGraph {
         friend std::ostream& operator << (std::ostream& os, const FactorGraph& fg);
         friend std::istream& operator >> (std::istream& is, FactorGraph& fg);
 
+        size_t VV2E(size_t n1, size_t n2) const { return G.VV2E(n1,n2); }
+        const Edge& edge(size_t e) const { return G.edge(e); }
+        void indexEdges() { G.indexEdges(); }
+        size_t nr_edges() const { return G.nr_edges(); }
+        const std::vector<Edge>& edges() const { return G.edges(); }
+
     private:
         /// Part of constructors (creates edges, neighbors and adjacency matrix)
         void constructGraph( size_t nrEdges );
index ca3e173..45900a7 100644 (file)
@@ -35,6 +35,7 @@
 #include <iostream>
 #include <cstdio>
 #include <boost/foreach.hpp>
+#include <boost/functional/hash.hpp>
 #include <algorithm>
 
 
 /// An alias to the BOOST_FOREACH macro from the boost::foreach library
 #define foreach BOOST_FOREACH
 
+#ifdef DAI_DEBUG
+#define DAI_PV(x) do {std::cerr << #x "= " << x << std::endl;} while(0)
+#define DAI_DMSG(str) do {std::cerr << str << std::endl;} while(0)
+#else
+#define DAI_PV(x) do {} while(0)
+#define DAI_DMSG(str) do {} while(0)
+#endif
+
+#define DAI_ACCMUT(x,y)                     \
+      x y;                                  \
+      const x const y;
+
+#define DAI_IFVERB(n, stmt) if(props.verbose>=n) { cerr << stmt; }
+
 
 /// Real number (alias for double, which could be changed to long double if necessary)
 typedef double Real;
@@ -77,14 +92,14 @@ namespace dai {
     /// hash_map is an alias for std::map.
     /** Since there is no TR1 unordered_map implementation available yet, we fall back on std::map.
      */
-    template <typename T, typename U>
+    template <typename T, typename U, typename H = boost::hash<T> >
         class hash_map : public std::map<T,U> {};
 #else
     /// hash_map is an alias for std::tr1::unordered_map.
     /** We use the (experimental) TR1 unordered_map implementation included in modern GCC distributions.
      */
-    template <typename T, typename U>
-        class hash_map : public std::tr1::unordered_map<T,U> {};
+    template <typename T, typename U, typename H = boost::hash<T> >
+        class hash_map : public std::tr1::unordered_map<T,U,H> {};
 #endif
 
 
index 643c917..d68d389 100644 (file)
@@ -66,6 +66,14 @@ InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const PropertySet
 #ifdef DAI_WITH_GIBBS
     if( name == Gibbs::Name )
         return new Gibbs (fg, opts);
+#endif
+#ifdef DAI_WITH_BP_DUAL
+    if( name == BP_dual::Name )
+        return new BP_dual (fg, opts);
+#endif
+#ifdef DAI_WITH_CBP
+    if( name == CBP::Name )
+        return new CBP (fg, opts);
 #endif
     DAI_THROW(UNKNOWN_DAI_ALGORITHM);
 }
diff --git a/src/bbp.cpp b/src/bbp.cpp
new file mode 100644 (file)
index 0000000..10412fd
--- /dev/null
@@ -0,0 +1,770 @@
+#include <dai/bp.h>
+#include <dai/bbp.h>
+#include <dai/gibbs.h>
+
+// for makeBBPGraph: {
+#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
+Prob unnormAdjoint(const Prob &w, Real Z_w, const Prob &adj_w) {
+    assert(w.size()==adj_w.size());
+    Prob adj_w_unnorm(w.size(),0.0);
+    Real s=0.0;
+    for(size_t i=0; i<w.size(); i++) {
+        s += w[i]*adj_w[i];
+    }
+    for(size_t i=0; i<w.size(); i++) {
+        adj_w_unnorm[i] = (adj_w[i]-s)/Z_w;
+    }
+    return adj_w_unnorm;
+}
+
+/// 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) {
+        vector<Prob> b1_adj;
+        vector<Prob> b2_adj;
+        vector<Prob> psi1_adj;
+        vector<Prob> psi2_adj;
+        b1_adj.reserve(fg.nrVars());
+        psi1_adj.reserve(fg.nrVars());
+        b2_adj.reserve(fg.nrFactors());
+        psi2_adj.reserve(fg.nrFactors());
+        for(size_t i=0; i<fg.nrVars(); i++) {
+            size_t dim = fg.var(i).states();
+            int c = fg.nbV(i).size();
+            Prob p(dim,0.0);
+            for(size_t xi=0; xi<dim; xi++) {
+                p[xi] = (1-c)*(1+log(fg.belief1(i)[xi]));
+            }
+            b1_adj.push_back(p);
+
+            for(size_t xi=0; xi<dim; xi++) {
+                p[xi] = -fg.belief1(i)[xi];
+            }
+            psi1_adj.push_back(p);
+        }
+        for(size_t I=0; I<fg.nrFactors(); I++) {
+            size_t dim = fg.factor(I).states();
+            Prob p(dim,0.0);
+            for(size_t xI=0; xI<dim; xI++) {
+                p[xI] = 1 + log(fg.belief2(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];
+            }
+            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) {
+        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];
+                if(bIxI<1.0e-15) {
+                    p[xI] = -1.0e10;
+                } else {
+                    p[xI] = 1+log(bIxI);
+                }
+            }
+            b2_adj.push_back(p);
+        }
+        bbp.init(get_zero_adj_1(fg), b2_adj);
+    } else if(cfn_type==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];
+                if(bixi<1.0e-15) {
+                    p[xi] = -1.0e10;
+                } else {
+                    p[xi] = 1+log(bixi);
+                }
+            }
+            b1_adj.push_back(p);
+        }
+        bbp.init(b1_adj);
+    } else { // cost functions that use Gibbs sample: cfn_b, cfn_b2, cfn_exp
+        vector<size_t> state;
+        if(stateP==NULL) {
+            state = getGibbsState(fg,2*fg.doneIters());
+        } else {
+            state = *stateP;
+        }
+        assert(state.size()==fg.nrVars());
+
+        vector<Prob> b1_adj;
+        b1_adj.reserve(fg.nrVars());
+        for(size_t i=0; i<state.size(); i++) {
+            size_t n = fg.var(i).states();
+            Prob delta(n,0.0);
+            assert(/*0<=state[i] &&*/ state[i] < n);
+            double b = 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(); }
+            b1_adj.push_back(delta);
+        }
+        bbp.init(b1_adj);
+    }
+}
+
+/// 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) {
+    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
+        for(size_t i=0; i<fg.nrVars(); i++) {
+            cf += -fg.belief1(i).entropy();
+        }
+    } else {
+        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(); }
+        }
+    }
+    return cf;
+}
+
+vector<Prob> get_zero_adj_2(const BP_dual& bp_dual) {
+    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)));
+    }
+    return adj_2;
+}
+
+vector<Prob> get_zero_adj_1(const BP_dual& bp_dual) {
+    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)));
+    }
+    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();         \
+        size_t j_states = _fg->var(j).states();         \
+        if(_fg->var(i) > _fg->var(j)) {                 \
+            size_t xij=0;                               \
+            for(size_t xi=0; xi<i_states; xi++) {       \
+                for(size_t xj=0; xj<j_states; xj++) {   \
+                    body;                               \
+                    xij++;                              \
+                }                                       \
+            }                                           \
+        } else {                                        \
+            size_t xij=0;                               \
+            for(size_t xj=0; xj<j_states; xj++) {       \
+                for(size_t xi=0; xi<i_states; xi++) {   \
+                    body;                               \
+                    xij++;                              \
+                }                                       \
+            }                                           \
+        }                                               \
+    }
+
+void BBP::RegenerateInds() {
+    // initialise _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));
+                    }
+                }
+            }
+        }
+    }
+}
+
+void BBP::RegenerateT() {
+    _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.push_back(prod);
+    }
+}
+
+void BBP::RegenerateU() {
+    _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.push_back(prod);
+    }
+}
+
+void BBP::RegenerateS() {
+    _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]];
+                    }
+                }
+            }
+        }
+        // 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.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.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_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]));
+    }
+}
+
+/// 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]);
+        assert(p.size()==_fg->var(i).states());
+        foreach(size_t I, _fg->nbV(i)) {
+            p *= _bp_dual->msgM(I,i);
+        }
+        p += _init_adj_psi_1[i];
+        _adj_psi_1.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]);
+        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));
+            // 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 += _init_adj_psi_2[I];
+        _adj_psi_2.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());
+    _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]];
+            }
+        }
+        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
+    }
+    _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);
+            }
+        }
+        _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();
+    RegenerateU();
+    RegenerateS();
+    RegenerateR();
+    RegenerateInputs();
+    RegeneratePsiAdjoints();
+    RegenerateMessageAdjoints();
+    _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];
+    new_adj_n_iI = Prob(_fg->var(i).states(),0.0);
+    foreach(size_t j, _fg->nbF(I)) {
+        if(j!=i) {
+            size_t iIj = vfv_iI[j];
+            Prob &p = _S[iIj];
+            size_t jI = _fg->VV2E(j,I);
+            LOOP_ij(
+                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;
+
+    _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)) {
+        if(J!=I) {
+            size_t iJ = _fg->VV2E(i,J);
+            _new_adj_m[iI] += _R[fvf_Ii[J]]*_adj_n_unnorm[iJ];
+        }
+    }
+}
+
+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::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::doParUpdate() {
+    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
+        calcNewM(iI);
+        calcNewN(iI);
+    }
+    for(size_t iI=0; iI<_fg->nr_edges(); iI++) {
+        upMsgM(iI);
+        upMsgN(iI);
+    }
+}
+
+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();
+    }
+    return s/_fg->nr_edges();
+}
+
+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) {
+    PropertySet gibbsProps;
+    gibbsProps.Set("iters", iters);
+    gibbsProps.Set("verbose", oneLess(fg.Verbose()));
+    Gibbs gibbs(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());
+    // calculate the value of the unperturbed cost function
+    Real cf0 = gibbsCostFn(bp_dual, cfn, &state);
+
+    // run BBP to estimate adjoints
+    BBP bbp(bp_dual);
+    gibbsInitBBPCostFnAdj(bbp, bp_dual, cfn, &state);
+    bbp.run(bbpTol);
+
+    Real d=0;
+
+    if(1) {
+        // verify bbp.adj_psi_1
+
+        // for each variable i
+        for(size_t i=0; i<bp_dual.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);
+
+                // 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);
+                psi_1_prb[xi] += h;
+                psi_1_prb.normalize();
+
+                // 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);
+
+                // use this to estimate the adjoint for i
+                // XXX why is it off by a factor of 2?
+                adj_est.push_back((cf_prb-cf0)/h);
+            }
+            Prob p_adj_est(adj_est);
+            // 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);
+        }
+    }
+    if(1) {
+        // verify bbp.adj_n and bbp.adj_m
+
+        // We actually want to check the responsiveness of objective
+        // function to changes in the final messages. But at the end of a
+        // BBP run, the message adjoints are for the initial messages (and
+        // they should be close to zero, see paper). So this resets the
+        // BBP adjoints to the refer to the desired final messages
+        bbp.RegenerateMessageAdjoints();
+
+        // for each variable i
+        for(size_t i=0; i<bp_dual.nrVars(); i++) {
+            // for each factor I ~ i
+            foreach(size_t I, bp_dual.nbV(i)) {
+                vector<double> adj_n_est;
+                // for each value xi
+                for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                    BP_dual bp_dual_prb(bp_dual);
+                    // make h-sized change to newMsgN
+                    bp_dual_prb.newMsgN(i,I)[xi] += h;
+                    // recalculate beliefs
+                    bp_dual_prb.CalcBeliefs();
+                    // get cost function value
+                    Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                    // add it to list of adjoints
+                    adj_n_est.push_back((cf_prb-cf0)/h);
+                }
+        
+                vector<double> adj_m_est;
+                // for each value xi
+                for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                    BP_dual bp_dual_prb(bp_dual);
+                    // make h-sized change to newMsgM
+                    bp_dual_prb.newMsgM(I,i)[xi] += h;
+                    // recalculate beliefs
+                    bp_dual_prb.CalcBeliefs();
+                    // get cost function value
+                    Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+                    // add it to list of adjoints
+                    adj_m_est.push_back((cf_prb-cf0)/h);
+                }
+
+                Prob p_adj_n_est(adj_n_est);
+                // compare this numerical estimate to the BBP estimate; sum the distances
+                cerr << "i: " << i << ", I: " << I
+                     << ", adj_n_est: " << p_adj_n_est
+                     << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
+                d += dist(p_adj_n_est, bbp.adj_n(i,I), Prob::DISTL1);
+
+                Prob p_adj_m_est(adj_m_est);
+                // compare this numerical estimate to the BBP estimate; sum the distances
+                cerr << "i: " << i << ", I: " << I
+                     << ", adj_m_est: " << p_adj_m_est
+                     << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
+                d += dist(p_adj_m_est, bbp.adj_m(I,i), Prob::DISTL1);
+            }
+        }
+    }
+    if(0) {
+        // verify bbp.adj_b_1
+        for(size_t i=0; i<bp_dual.nrVars(); i++) {
+            vector<double> adj_b_1_est;
+            // for each value xi
+            for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
+                BP_dual bp_dual_prb(bp_dual);
+
+                // make h-sized change to b_1(i)[x_i]
+                bp_dual_prb._beliefs.b1[i][xi] += h;
+
+                // get cost function value
+                Real cf_prb = gibbsCostFn(bp_dual_prb, cfn, &state);
+
+                // add it to list of adjoints
+                adj_b_1_est.push_back((cf_prb-cf0)/h);
+            }
+            Prob p_adj_b_1_est(adj_b_1_est);
+            // 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);
+        }
+    }
+
+    // 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);
+
+                os << n << "\t" << cf << endl;
+            }
+      
+            os.close();
+            free(fn);
+        }
+    }
+}
+
+} // end of namespace dai
diff --git a/src/bp_dual.cpp b/src/bp_dual.cpp
new file mode 100644 (file)
index 0000000..a5a9055
--- /dev/null
@@ -0,0 +1,444 @@
+
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include <algorithm>
+
+#include <dai/bbp.h>
+//#include <dai/diffs.h>
+#include <dai/util.h>
+//#include "stlutil.h"
+#include <dai/properties.h>
+
+namespace dai {
+
+using namespace std;
+
+const char *BP_dual::Name = "BP_dual";
+
+const char *BP_dual::PropertyList[] = {"tol","maxiter","updates","verbose"};
+
+void BP_dual::setProperties( const PropertySet &opts ) {
+//     DAI_DMSG("in BP_dual::setProperties");
+//     DAI_PV(opts);
+
+    bool die=false;
+    foreach(const char *p, PropertyList) {
+        if( !opts.hasKey(p) ) {
+            cerr << "BP_dual: missing property " << p << endl;
+            die=true;
+        }
+    }
+    if(die) throw "BP_dual: Couldn't set properties";
+    
+    props.tol = opts.getStringAs<double>("tol");
+    props.maxiter = opts.getStringAs<size_t>("maxiter");
+    props.updates = opts.getStringAs<Properties::UpdateType>("updates");
+    props.verbose = opts.getStringAs<size_t>("verbose");
+
+//     DAI_PV(printProperties());
+}
+
+PropertySet BP_dual::getProperties() const {
+    PropertySet opts;
+    opts.Set( "tol", props.tol );
+    opts.Set( "maxiter", props.maxiter );
+    opts.Set( "updates", props.updates );
+    opts.Set( "verbose", props.verbose );
+    return opts;
+}
+
+std::string BP_dual::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "tol=" << props.tol << ",";
+    s << "maxiter=" << props.maxiter << ",";
+    s << "updates=" << props.updates << ",";
+    s << "verbose=" << props.verbose;
+    s << "]";
+    return s.str();
+}
+
+// void BP_dual::checkProperties() {
+//     const char *props[] = {"updates","tol","maxiter","verbose"};
+//     for(size_t i=0; i<sizeof(props)/sizeof(*props); i++) {
+//         if(!HasProperty(props[i]))
+//             die("BP_dual: Missing property \"%s\"", props[i]);
+//     }
+    
+//     ConvertPropertyTo<double>("tol");
+//     ConvertPropertyTo<size_t>("maxiter");
+//     ConvertPropertyTo<size_t>("verbose");
+//     ConvertPropertyTo<UpdateType>("updates");
+// }
+
+void BP_dual::RegenerateIndices() {
+    _indices.clear();
+    _indices.reserve(nr_edges());
+
+    foreach(Edge iI, edges()) {
+        vector<size_t> ind( factor(iI.second).states(), 0 );
+        IndexFor i (var(iI.first), factor(iI.second).vars() );
+        for( size_t j = 0; i >= 0; ++i,++j )
+            ind[j] = i; 
+        _indices.push_back( ind );
+    }
+}
+
+void BP_dual::RegenerateMessages() {
+    _msgs.Zn.resize(nr_edges(),1.0);
+    _msgs.Zm.resize(nr_edges(),1.0);
+
+    // clear messages
+    _msgs.m.clear();
+    _msgs.m.reserve(nr_edges());
+    _msgs.n.clear();
+    _msgs.n.reserve(nr_edges());
+
+    // create messages and indices
+    foreach(Edge iI, edges()) {
+        // initialize to uniform distributions
+        _msgs.m.push_back( Prob( var(iI.first).states() ) );
+        _msgs.n.push_back( Prob( var(iI.first).states() ) );
+    }
+
+    // create new_messages
+    _new_msgs = _msgs;
+}
+
+void BP_dual::RegenerateBeliefs() {
+    _beliefs.b1.clear();
+    _beliefs.b1.reserve(nrVars());
+    _beliefs.Zb1.resize(nrVars(), 1.0);
+    _beliefs.b2.clear();
+    _beliefs.b2.reserve(nrFactors());
+    _beliefs.Zb2.resize(nrFactors(), 1.0);
+
+    for(size_t i=0; i<nrVars(); i++) {
+        _beliefs.b1.push_back( Prob( var(i).states() ).setUniform() );
+    }
+    for(size_t I=0; I<nrFactors(); I++) {
+        _beliefs.b2.push_back( Prob( factor(I).states() ).setUniform() );
+    }
+}
+
+// called by constructor, called before 'init'
+void BP_dual::Regenerate() {
+
+    indexEdges(); // so we can use compatibility interface
+
+//     DAIAlgFG::Regenerate(); // located in BipartiteGraph
+  
+    RegenerateIndices();
+    RegenerateMessages();
+    RegenerateBeliefs();
+
+    _maxdiff = 0;
+    _iters = 0;
+}
+
+void BP_dual::CalcBelief1(size_t i) {
+    Prob prod( var(i).states(), 1.0 );
+    foreach(size_t I, nbV(i)) {
+        prod *= newMsgM(I,i);
+    }
+    _beliefs.Zb1[i] = prod.normalize();
+    _beliefs.b1[i] = prod;
+}
+
+void BP_dual::CalcBelief2(size_t I) {
+    Prob prod( factor(I).p() );
+    foreach(size_t j, nbF(I)) {
+        const _ind_t *ind = &(index(j, I));
+        for(size_t r=0; r<prod.size(); r++) {
+            Prob n(newMsgN(j,I));
+            prod[r] *= n[(*ind)[r]];
+        }
+    }
+    _beliefs.Zb2[I] = prod.normalize();
+    _beliefs.b2[I] = prod;
+}
+
+// called after run()
+void BP_dual::CalcBeliefs() {
+    for(size_t i=0; i<nrVars(); i++) {
+        // calculate b_i
+        CalcBelief1(i);
+    }
+    for(size_t I=0; I<nrFactors(); I++) {
+        // calculate b_I
+        CalcBelief2(I);
+    }
+}
+
+void BP_dual::calcNewM(size_t iI) {
+    // calculate updated message I->i
+    size_t i = edge(iI).first;
+    size_t I = edge(iI).second;
+
+    Prob prod( factor(I).p() );
+
+    foreach(size_t j, nbF(I)) {
+        if( j != i ) {     // for all j in I \ i
+            _ind_t* ind = &(index(j,I));
+            Prob n(msgN(j,I));
+            for( size_t r = 0; r < prod.size(); r++ )
+                prod[r] *= n[(*ind)[r]];
+        }
+    }
+
+    // Marginalize onto i
+    Prob marg( var(i).states(), 0.0 );
+    // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
+    _ind_t* ind = &(index(i,I));
+    for( size_t r = 0; r < prod.size(); r++ )
+        marg[(*ind)[r]] += prod[r];
+    
+    _new_msgs.Zm[iI] = marg.normalize();
+    _new_msgs.m[iI] = marg;
+}
+
+void BP_dual::calcNewN(size_t iI) {
+    // XXX optimize
+    // calculate updated message i->I
+    size_t i = edge(iI).first;
+    size_t I = edge(iI).second;
+
+    Prob prod(var(i).states(), 1.0);
+    foreach(size_t J, nbV(i)) {
+        if(J != I) { // for all J in i \ I
+            prod *= msgM(J,i);
+        }
+    }
+    _new_msgs.Zn[iI] = prod.normalize();
+    _new_msgs.n[iI] = prod;
+}
+
+void BP_dual::upMsgM(size_t iI) {
+    _msgs.m[iI] = _new_msgs.m[iI];
+    _msgs.Zm[iI] = _new_msgs.Zm[iI];
+}
+
+void BP_dual::upMsgN(size_t iI) {
+    _msgs.n[iI] = _new_msgs.n[iI];
+    _msgs.Zn[iI] = _new_msgs.Zn[iI];
+}
+
+double BP_dual::run() {
+    DAI_IFVERB(1, "Starting " << identify() << "..." << endl);
+
+    double tic = toc();
+    // for some reason we need 2* here, where orig BP doesn't
+    Diffs diffs(2*nrVars(), 1.0);
+    
+    vector<size_t> edge_seq;
+    vector<double> residuals;
+
+    vector<Factor> old_beliefs;
+    old_beliefs.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ ) {
+        CalcBelief1(i);
+        old_beliefs.push_back(belief1(i));
+    }
+
+    size_t iter = 0;
+
+    if( Updates() == UpdateType::SEQMAX ) {
+        // do the first pass
+        for(size_t iI = 0; iI < nr_edges(); iI++ ) {
+            calcNewM(iI);
+            calcNewN(iI);
+        }
+
+        // calculate initial residuals
+        residuals.reserve(nr_edges());
+        for( size_t iI = 0; iI < nr_edges(); iI++ )
+            residuals.push_back( dist( _new_msgs.m[iI], _msgs.m[iI], Prob::DISTLINF ) );
+    } else {
+        edge_seq.reserve( nr_edges() );
+        for( size_t i = 0; i < nr_edges(); i++ )
+            edge_seq.push_back( i );
+    }
+
+    // do several passes over the network until maximum number of iterations has
+    // been reached or until the maximum belief difference is smaller than tolerance
+    for( iter=0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
+        if( Updates() == UpdateType::SEQMAX ) {
+            // Residuals-BP by Koller et al.
+            for( size_t t = 0; t < nr_edges(); t++ ) {
+                // update the message with the largest residual
+                size_t iI = max_element(residuals.begin(), residuals.end()) - residuals.begin();
+                upMsgM(iI);
+                residuals[iI] = 0;
+
+                // I->i has been updated, which means that residuals for all
+                // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
+                size_t i = edge(iI).first;
+                size_t I = edge(iI).second;
+                foreach(size_t J, nbV(i)) {
+                    if(J != I) {
+                        size_t iJ = VV2E(i,J);
+                        calcNewN(iJ);
+                        upMsgN(iJ);
+                        foreach(size_t j, nbF(J)) {
+                            if(j != i) {
+                                size_t jJ = VV2E(j,J);
+                                calcNewM(jJ);
+                                residuals[jJ] = dist( _new_msgs.m[jJ], _msgs.m[jJ], Prob::DISTLINF );
+                            }
+                        }
+                    }
+                }
+            }
+        } else if( Updates() == UpdateType::PARALL ) {
+            // Parallel updates 
+            for( size_t t = 0; t < nr_edges(); t++ ) {
+                calcNewM(t);
+                calcNewN(t);
+            }
+            if(0) {
+                for(size_t t=0; t<nr_edges(); t++) {
+                    upMsgM(t); upMsgN(t);
+                }
+            } else {
+                _msgs = _new_msgs;
+            }
+        } else {
+            // Sequential updates
+            if( Updates() == UpdateType::SEQRND )
+                random_shuffle( edge_seq.begin(), edge_seq.end() );
+
+            foreach(size_t k, edge_seq) {
+                calcNewM(k);
+                calcNewN(k);
+                upMsgM(k);
+                upMsgN(k);
+            }
+        }
+
+        // calculate new beliefs and compare with old ones
+        for( size_t i = 0; i < nrVars(); i++ ) {
+            CalcBelief1(i);
+            Factor nb( belief1(i) );
+            diffs.push( dist( nb, old_beliefs[i], Prob::DISTLINF ) );
+            old_beliefs[i] = nb;
+        }
+
+        DAI_IFVERB(3,"BP_dual::run:  maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl);
+
+        _iters++;
+    }
+
+    updateMaxDiff( diffs.maxDiff() );
+
+    if( props.verbose >= 1 ) {
+        if( diffs.maxDiff() > props.tol ) {
+            DAI_IFVERB(1, endl << "BP_dual::run:  WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl);
+        } else {
+            DAI_IFVERB(3, "BP_dual::run:  converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl);
+        }
+    }
+
+    CalcBeliefs();
+    
+    return diffs.maxDiff();
+}
+
+string BP_dual::identify() const { 
+    return string(Name) + printProperties();
+}
+
+vector<Factor> BP_dual::beliefs() const {
+    vector<Factor> result;
+    for( size_t i = 0; i < nrVars(); i++ )
+        result.push_back( belief1(i) );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        result.push_back( belief2(I) );
+    return result;
+}
+
+Factor BP_dual::belief( const VarSet &ns ) const {
+    if( ns.size() == 1 )
+        return belief( *(ns.begin()) );
+    else {
+        size_t I;
+        for( I = 0; I < nrFactors(); I++ )
+            if( factor(I).vars() >> ns )
+                break;
+        assert( I != nrFactors() );
+        return belief2(I).marginal(ns);
+    }
+}
+
+Real BP_dual::logZ() const {
+    Real sum = 0.0;
+    for(size_t i = 0; i < nrVars(); i++ )
+        sum += Real(1.0 - nbV(i).size()) * belief1(i).entropy();
+    for( size_t I = 0; I < nrFactors(); I++ )
+        sum -= dist( belief2(I), factor(I), Prob::DISTKL );
+    return sum;
+}
+
+// reset only messages to/from certain variables
+void BP_dual::init(const VarSet &ns) {
+    _iters=0;
+    foreach(Var n, ns) {
+        size_t ni = findVar(n);
+        size_t st = n.states();
+        foreach(Neighbor I, nbV(ni)) {
+            msgM(I.node,ni).fill(1.0/st);
+            zM(I.node,ni) = 1.0;
+            msgN(ni,I.node).fill(1.0/st);
+            zN(ni,I.node) = 1.0;
+        }
+    }
+}
+
+void BP_dual::init() {
+    _iters=0;
+    for(size_t iI = 0; iI < nr_edges(); iI++ ) {
+        _msgs.m[iI].setUniform();
+        _msgs.Zm[iI] = 1;
+        _msgs.n[iI].setUniform();
+        _msgs.Zn[iI] = 1;
+    }
+    _new_msgs = _msgs;
+}
+
+
+void BP_dual::init(const vector<size_t>& state) {
+    _iters=0;
+    for(size_t iI = 0; iI < nr_edges(); iI++ ) {
+        size_t i = edge(iI).first;
+        _msgs.m[iI].fill(0.1);
+        _msgs.m[iI][state[i]]=1;
+        _msgs.Zm[iI] = _msgs.m[iI].normalize();
+        _msgs.n[iI].fill(0.1);
+        _msgs.n[iI][state[i]]=1;
+        _msgs.Zn[iI] = _msgs.n[iI].normalize();
+    }
+    _new_msgs = _msgs;
+}
+
+void _clamp(FactorGraph &g, const Var & n, const vector<size_t> &is ) {
+    Factor mask_n(n,0.0);
+
+    foreach(size_t i, is) { assert( i <= n.states() ); mask_n[i] = 1.0; }
+
+    for( size_t I = 0; I < g.nrFactors(); I++ ) 
+        if( g.factor(I).vars() && n )
+          g.factor(I) *= mask_n;
+}
+
+// clamp a factor to have one of a set of values
+void _clampFactor(FactorGraph &g, size_t I, const vector<size_t> &is) {
+    size_t st = g.factor(I).states();
+    Prob mask_n(st,0.0);
+
+    foreach(size_t i, is) { assert( i <= st ); mask_n[i] = 1.0; }
+
+    g.factor(I).p() *= mask_n;
+}
+
+} // end of namespace dai
diff --git a/src/cbp.cpp b/src/cbp.cpp
new file mode 100644 (file)
index 0000000..ca11449
--- /dev/null
@@ -0,0 +1,577 @@
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include <algorithm>
+
+#include <dai/util.h>
+#include <dai/properties.h>
+
+#include <dai/bp.h>
+#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)
+
+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;
+}
+
+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());
+    for( size_t i = 0; i < nrVars(); i++ )
+        _beliefs1.push_back( Factor(var(i)).normalized() );
+
+    _beliefs2.clear(); _beliefs2.reserve(nrFactors());
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        Factor f = factor(I);
+        f.fill(1); f.normalize();
+        _beliefs2.push_back(f);
+    }
+
+    // to compute average level
+    _sum_level = 0;
+    _num_leaves = 0;
+
+    _maxdiff = 0;
+    _iters = 0;
+}
+
+static
+vector<Factor> mixBeliefs(Real p, vector<Factor> b, vector<Factor> c) {
+  vector<Factor> out;
+  assert(b.size()==c.size());
+  out.reserve(b.size());
+  Real pc = 1-p;
+  for(size_t i=0; i<b.size(); i++) {
+      // XXX probably already normalized
+    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();
+
+    vector<Factor> beliefs_out;
+    Real lz_out;
+    runRecurse(bp_dual, bp_dual.logZ(), 
+               vector<size_t>(0), set<size_t>(),
+               _num_leaves, _sum_level,
+               lz_out, beliefs_out);
+    if(props.verbose>=1) {
+      cerr << "CBP average levels = " << (_sum_level/_num_leaves) << ", leaves = " << _num_leaves << endl;
+    }
+    setBeliefs(beliefs_out, lz_out);
+    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() {
+    PropertySet bpProps;
+    bpProps.Set("updates", string("PARALL"));
+    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;
+}
+
+
+void CBP::runRecurse(BP_dual &bp_dual,
+                     double orig_logZ,
+                     vector<size_t> clamped_vars_list,
+                     set<size_t> clamped_vars,
+                     size_t &num_leaves,
+                     double &sum_level,
+                     Real &lz_out, vector<Factor>& beliefs_out) {
+    // choose a variable/states to clamp:
+    size_t i;
+    vector<size_t> xis;
+    Real maxVar=0.0;
+    bool found;
+    bool clampVar = (Clamping()==Properties::ClampType::CLAMP_VAR);
+
+    // XXX fix to just pass orig_logZ
+
+    if(Recursion()==Properties::RecurseType::REC_LOGZ && recTol()>0 &&
+       exp(bp_dual.logZ()-orig_logZ) < recTol()) {
+        found = false;
+    } else {
+        found = chooseNextClampVar(bp_dual,
+                                   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();
+        return;
+    }
+    
+    if(clampVar) {
+        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()); }
+    }
+    // - otherwise, clamp and recurse, saving margin estimates for each
+    // clamp setting. afterwards, combine estimates.
+
+    // compute complement of 'xis'
+    vector<size_t> cmp_xis=complement(xis, clampVar?var(i).states():factor(i).states());
+
+    // XXX 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));
+    } else {
+        _clampFactor((FactorGraph&)bp_dual_c, i, xis);
+        bp_dual_c.init(factor(i).vars());
+    }
+    bp_dual_c.run();
+    _iters += bp_dual_c.Iterations();
+
+    lz = bp_dual_c.logZ();
+    b = bp_dual_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));
+    } else {
+        _clampFactor(cmp_bp_dual_c,i,cmp_xis);
+        cmp_bp_dual_c.init(factor(i).vars());
+    }
+    cmp_bp_dual_c.run();
+    _iters += cmp_bp_dual_c.Iterations();
+
+    cmp_lz = cmp_bp_dual_c.logZ();
+    cmp_b = cmp_bp_dual_c.beliefs();
+
+    double p = unSoftMax(lz, cmp_lz);
+    Real bp__d=0.0;
+    
+    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());
+        if(exp(new_lz-orig_logZ)*bp__d < recTol()) {
+            num_leaves++;
+            sum_level += clamped_vars_list.size();
+            beliefs_out = combined_b;
+            lz_out = new_lz;
+            return;
+        }
+    }
+
+    // either we are not doing REC_BDIFF or the distance was large
+    // enough to recurse:
+
+    runRecurse(bp_dual_c, orig_logZ, 
+               clamped_vars_list,
+               clamped_vars,
+               num_leaves, sum_level, lz, b);
+    runRecurse(cmp_bp_dual_c, orig_logZ, 
+               clamped_vars_list,
+               clamped_vars,
+               num_leaves, sum_level, cmp_lz, cmp_b);
+
+    p = unSoftMax(lz,cmp_lz);
+
+    beliefs_out = mixBeliefs(p, b, cmp_b);
+    lz_out = logSumExp(lz,cmp_lz);
+
+    if(props.verbose>=2) {
+        Real d = dist( bp_dual.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;
+    }
+}
+
+// 'xis' must be sorted
+bool CBP::chooseNextClampVar(BP_dual& 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) {
+        cerr << "clamped_vars_list" << clamped_vars_list << endl;
+    }
+    if(clamped_vars_list.size() >= maxClampLevel()) {
+        return false;
+    }
+    if(ChooseMethod()==Properties::ChooseMethodType::CHOOSE_RANDOM) {
+        if(Clamping()==Properties::ClampType::CLAMP_VAR) {
+            int t=0, t1=100;
+            do {
+                i = rnd_multi(nrVars());
+                t++;
+            } while(abs(bp.belief1(i).p().max()-1) < tiny &&
+                    t < t1);
+            if(t==t1) {
+                return false;
+                //             die("Too many levels requested in CBP");
+            }
+            // only pick probable values for variable
+            size_t xi;
+            do {
+                xi = rnd_multi(var(i).states());
+                t++;
+            } while(bp.belief1(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);
+        } else {
+            int t=0, t1=100;
+            do {
+                i = rnd_multi(nrFactors());
+                t++;
+            } while(abs(bp.belief2(i).p().max()-1) < tiny &&
+                    t < t1);
+            if(t==t1) {
+                return false;
+                //             die("Too many levels requested in CBP");
+            }
+            // only pick probable values for variable
+            size_t xi;
+            do {
+                xi = rnd_multi(factor(i).states());
+                t++;
+            } while(bp.belief2(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;
+          }
+        }
+      }
+      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);
+        pair<size_t, size_t> cv =
+            bbpFindClampVar(bp,
+                            clampVar,
+                            clamped_vars_list.size(),
+                            BBP_cost_function(),getProperties(),maxVarOut);
+
+        // if slope isn't big enough then don't clamp
+        if(*maxVarOut < minMaxAdj()) return false;
+
+        size_t xi=cv.second;
+        i = cv.first;
+#define VAR_INFO (clampVar?"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());
+        if(b[xi] < tiny) {
+            cerr << "Warning, bbpFindClampVar found unlikely "
+                 << VAR_INFO << endl;
+            return false;
+        }
+        if(abs(b[xi]-1) < tiny) {
+            cerr << "Warning, bbpFindClampVar found overly likely "
+                 << VAR_INFO << endl;
+            return false;
+        }
+
+        xis.resize(1,xi);
+        if(clampVar) {
+            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);
+    } 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(_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;
+    }
+  
+    // find and return the (variable,state) with the largest adj_psi_1
+    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(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();
+            }
+            // 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();
+
+            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 < in_bp_dual.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);
+            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();
+            }
+            // 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();
+
+            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 < in_bp_dual.factor(argmax_var).states());
+    }
+    if(maxVarOut) *maxVarOut = max_var;
+    return make_pair(argmax_var,argmax_var_state);
+}
+
+} // end of namespace dai
index f031109..c7f62a1 100644 (file)
@@ -1,5 +1,6 @@
 # --- BP ----------------------
 
+BP:                             BP[updates=SEQFIX,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQFIX:                      BP[updates=SEQFIX,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQRND:                      BP[updates=SEQRND,tol=1e-9,maxiter=10000,logdomain=0]
 BP_SEQMAX:                      BP[updates=SEQMAX,tol=1e-9,maxiter=10000,logdomain=0]
@@ -100,3 +101,13 @@ GIBBS_1e6:                      GIBBS[iters=1000000]
 GIBBS_1e7:                      GIBBS[iters=10000000]
 GIBBS_1e8:                      GIBBS[iters=100000000]
 GIBBS_1e9:                      GIBBS[iters=1000000000]
+
+# --- BP_dual ---------------------
+
+BP_dual:                        BP_dual[updates=SEQFIX,tol=1e-9,maxiter=10000,verbose=0]
+BP_dual_SEQFIX:                 BP_dual[updates=SEQFIX,tol=1e-9,maxiter=10000,verbose=0]
+
+# --- 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]
+BBP:                            CBP[choose=CHOOSE_BBP]