From ab18239d0f5f2a5f7fc7bf310b5ca71b5f74fc46 Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Tue, 3 Mar 2009 14:49:10 +0100 Subject: [PATCH] [Frederik Eaton] Added BP_Dual, BBP and CBP algorithms --- Makefile | 20 +- Makefile.CYGWIN | 2 + Makefile.LINUX | 2 + Makefile.MACOSX | 2 + Makefile.WINDOWS | 2 + include/dai/alldai.h | 9 + include/dai/bbp.h | 165 ++++++++ include/dai/bipgraph.h | 56 ++- include/dai/bp_dual.h | 147 ++++++++ include/dai/cbp.h | 155 ++++++++ include/dai/factorgraph.h | 6 + include/dai/util.h | 21 +- src/alldai.cpp | 8 + src/bbp.cpp | 770 ++++++++++++++++++++++++++++++++++++++ src/bp_dual.cpp | 444 ++++++++++++++++++++++ src/cbp.cpp | 577 ++++++++++++++++++++++++++++ tests/aliases.conf | 11 + 17 files changed, 2391 insertions(+), 6 deletions(-) create mode 100644 include/dai/bbp.h create mode 100644 include/dai/bp_dual.h create mode 100644 include/dai/cbp.h create mode 100644 src/bbp.cpp create mode 100644 src/bp_dual.cpp create mode 100644 src/cbp.cpp diff --git a/Makefile b/Makefile index bdfbdf1..f07f56b 100644 --- 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 diff --git a/Makefile.CYGWIN b/Makefile.CYGWIN index 3e2ce12..64916af 100644 --- a/Makefile.CYGWIN +++ b/Makefile.CYGWIN @@ -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? diff --git a/Makefile.LINUX b/Makefile.LINUX index bbfde5c..4c9fa51 100644 --- a/Makefile.LINUX +++ b/Makefile.LINUX @@ -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? diff --git a/Makefile.MACOSX b/Makefile.MACOSX index 9cd503e..9bcf11b 100644 --- a/Makefile.MACOSX +++ b/Makefile.MACOSX @@ -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? diff --git a/Makefile.WINDOWS b/Makefile.WINDOWS index c6ee945..9290701 100644 --- a/Makefile.WINDOWS +++ b/Makefile.WINDOWS @@ -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? diff --git a/include/dai/alldai.h b/include/dai/alldai.h index 430e8eb..91b4285 100644 --- a/include/dai/alldai.h +++ b/include/dai/alldai.h @@ -57,6 +57,12 @@ #ifdef DAI_WITH_GIBBS #include #endif +#ifdef DAI_WITH_BP_DUAL + #include +#endif +#ifdef DAI_WITH_CBP + #include +#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 index 0000000..b14b6a3 --- /dev/null +++ b/include/dai/bbp.h @@ -0,0 +1,165 @@ +#ifndef ___defined_libdai_bbp_h +#define ___defined_libdai_bbp_h + +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include + +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 getGibbsState(const BP_dual& fg, size_t iters); + +vector get_zero_adj_2(const BP_dual& bp_dual); +vector get_zero_adj_1(const BP_dual& bp_dual); + +class BBP { + // ---------------------------------------------------------------- + // inputs + const BP_dual* _bp_dual; + const FactorGraph* _fg; + vector _adj_b_1, _adj_b_2; + vector _adj_b_1_unnorm, _adj_b_2_unnorm; + // in case the objective function depends on factors as well: + vector _init_adj_psi_1; + vector _init_adj_psi_2; + // ---------------------------------------------------------------- + // ouptuts + vector _adj_psi_1, _adj_psi_2; + // the following vectors are length _fg->nr_edges() + vector _adj_n, _adj_m; + vector _adj_n_unnorm, _adj_m_unnorm; + vector _new_adj_n, _new_adj_m; + // ---------------------------------------------------------------- + // auxiliary data + typedef vector > > _trip_t; + typedef vector _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 > _trip_ind_t; + _trip_ind_t _VFV_ind; + _trip_ind_t _FVF_ind; + + // helper quantities computed from the BP messages + vector _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 &adj_b_1, const vector &adj_b_2, + const vector &adj_psi_1, const vector &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 &adj_b_1, const vector &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 &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* stateP); + +/// calculate actual value of cost function (cfn_type, stateP) +Real gibbsCostFn(const BP_dual& fg, bbp_cfn_t cfn_type, const vector *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 diff --git a/include/dai/bipgraph.h b/include/dai/bipgraph.h index f202f99..b3a85b6 100644 --- a/include/dai/bipgraph.h +++ b/include/dai/bipgraph.h @@ -101,9 +101,16 @@ class BipartiteGraph { std::vector 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 _edges; + hash_map _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 - 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& edges() const { + return _edges; + } + + size_t VV2E(size_t n1, size_t n2) const { + assert(_edge_indexed); + Edge e(n1,n2); + hash_map::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 index 0000000..aca369a --- /dev/null +++ b/include/dai/bp_dual.h @@ -0,0 +1,147 @@ +#ifndef ____defined_libdai_bp_dual_h__ +#define ____defined_libdai_bp_dual_h__ + +#include +#include +#include +#include + +namespace dai { + +using namespace std; + +struct BP_dual_messages { + // messages: + // indexed by edge index (using VV2E) + vector n; + vector Zn; + vector m; + vector Zm; +}; + +struct BP_dual_beliefs { + // beliefs: + // indexed by node + vector b1; + vector Zb1; + // indexed by factor + vector b2; + vector Zb2; +}; + +void _clamp(FactorGraph &g, const Var & n, const vector &is ); + +// clamp a factor to have one of a set of values +void _clampFactor(FactorGraph &g, size_t I, const vector &is); + +class BP_dual : public DAIAlgFG { + public: + typedef vector _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 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& 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 index 0000000..45a49f3 --- /dev/null +++ b/include/dai/cbp.h @@ -0,0 +1,155 @@ +#ifndef ____defined_libdai_cbp_h__ +#define ____defined_libdai_cbp_h__ + +#include + +#include +#include + +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 +vector concat(const vector& u, const vector& v) { + vector w; + w.reserve(u.size()+v.size()); + for(size_t i=0; i _beliefs1; + vector _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 clamped_vars_list, + set clamped_vars, + size_t &num_leaves, + double &sum_level, + Real &lz_out, vector& 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 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 &clamped_vars_list, + set &clamped_vars, + size_t &i, vector &xis, Real *maxVarOut); + + //---------------------------------------------------------------- + + void setBeliefsFrom(const BP& b); + void setBeliefs(const vector& bs, double logZ) { + size_t i=0; + _beliefs1.clear(); _beliefs1.reserve(nrVars()); + _beliefs2.clear(); _beliefs2.reserve(nrFactors()); + for(i=0; i bbpFindClampVar(BP_dual &in_bp_dual, bool clampVar, size_t numClamped, bbp_cfn_t cfn, const PropertySet &props, Real *maxVarOut); + +} + +#endif diff --git a/include/dai/factorgraph.h b/include/dai/factorgraph.h index a3f5d40..36e27e4 100644 --- a/include/dai/factorgraph.h +++ b/include/dai/factorgraph.h @@ -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& edges() const { return G.edges(); } + private: /// Part of constructors (creates edges, neighbors and adjacency matrix) void constructGraph( size_t nrEdges ); diff --git a/include/dai/util.h b/include/dai/util.h index ca3e173..45900a7 100644 --- a/include/dai/util.h +++ b/include/dai/util.h @@ -35,6 +35,7 @@ #include #include #include +#include #include @@ -50,6 +51,20 @@ /// 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 + template > class hash_map : public std::map {}; #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 - class hash_map : public std::tr1::unordered_map {}; + template > + class hash_map : public std::tr1::unordered_map {}; #endif diff --git a/src/alldai.cpp b/src/alldai.cpp index 643c917..d68d389 100644 --- a/src/alldai.cpp +++ b/src/alldai.cpp @@ -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 index 0000000..10412fd --- /dev/null +++ b/src/bbp.cpp @@ -0,0 +1,770 @@ +#include +#include +#include + +// for makeBBPGraph: { +#include +#include +#include +#include +// } + +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* stateP) { + if(cfn_type==bbp_cfn_t::bbp_cfn_t::cfn_bethe_ent) { + vector b1_adj; + vector b2_adj; + vector psi1_adj; + vector 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 b2_adj; + b2_adj.reserve(fg.nrFactors()); + for(size_t I=0; I b1_adj; + b1_adj.reserve(fg.nrVars()); + for(size_t i=0; i state; + if(stateP==NULL) { + state = getGibbsState(fg,2*fg.doneIters()); + } else { + state = *stateP; + } + assert(state.size()==fg.nrVars()); + + vector b1_adj; + b1_adj.reserve(fg.nrVars()); + for(size_t i=0; i *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 state = *stateP; + assert(state.size()==fg.nrVars()); + for(size_t i=0; i get_zero_adj_2(const BP_dual& bp_dual) { + vector adj_2; + adj_2.reserve(bp_dual.nrFactors()); + for(size_t I=0; I get_zero_adj_1(const BP_dual& bp_dual) { + vector adj_1; + adj_1.reserve(bp_dual.nrVars()); + for(size_t i=0; inr_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; xinrVars()); + _VFV_ind.clear(); + for(size_t i=0; i<_fg->nrVars(); i++) { + foreach(size_t I, _fg->nbV(i)) { + vector& 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(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& 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(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; xIvar(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 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 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 adj_est; + // for each value xi + for(size_t xi=0; xi0); + + // 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 adj_n_est; + // for each value xi + for(size_t xi=0; xi adj_m_est; + // for each value xi + for(size_t xi=0; xi adj_b_1_est; + // for each value xi + for(size_t xi=0; xi 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; i0); + 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 index 0000000..a5a9055 --- /dev/null +++ b/src/bp_dual.cpp @@ -0,0 +1,444 @@ + +#include +#include +#include +#include +#include + +#include +//#include +#include +//#include "stlutil.h" +#include + +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("tol"); + props.maxiter = opts.getStringAs("maxiter"); + props.updates = opts.getStringAs("updates"); + props.verbose = opts.getStringAs("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("tol"); +// ConvertPropertyTo("maxiter"); +// ConvertPropertyTo("verbose"); +// ConvertPropertyTo("updates"); +// } + +void BP_dual::RegenerateIndices() { + _indices.clear(); + _indices.reserve(nr_edges()); + + foreach(Edge iI, edges()) { + vector 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; ii + 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 edge_seq; + vector residuals; + + vector 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= 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 BP_dual::beliefs() const { + vector 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& 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 &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 &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 index 0000000..ca11449 --- /dev/null +++ b/src/cbp.cpp @@ -0,0 +1,577 @@ +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +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("tol"); + props.rec_tol = opts.getStringAs("rec_tol"); + props.maxiter = opts.getStringAs("maxiter"); + props.verbose = opts.getStringAs("verbose"); + props.updates = opts.getStringAs("updates"); + props.max_levels = opts.getStringAs("max_levels"); + props.min_max_adj = opts.getStringAs("min_max_adj"); + props.choose = opts.getStringAs("choose"); + props.recursion = opts.getStringAs("recursion"); + props.clamp = opts.getStringAs("clamp"); + props.bbp_cfn = opts.getStringAs("bbp_cfn"); + props.rand_seed = opts.getStringAs("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 mixBeliefs(Real p, vector b, vector c) { + vector out; + assert(b.size()==c.size()); + out.reserve(b.size()); + Real pc = 1-p; + for(size_t i=0; i("bp_alg"), *this, GetPropertyAs("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 beliefs_out; + Real lz_out; + runRecurse(bp_dual, bp_dual.logZ(), + vector(0), set(), + _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 complement(vector& xis, size_t n_states) { + vector cmp_xis(0); + size_t j=0; + for(size_t xi=0; xi=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& b1, const vector& b2, size_t nv) { + Real d=0.0; + for(size_t k=0; k clamped_vars_list, + set clamped_vars, + size_t &num_leaves, + double &sum_level, + Real &lz_out, vector& beliefs_out) { + // choose a variable/states to clamp: + size_t i; + vector 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 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 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 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 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 &clamped_vars_list, + set &clamped_vars, + size_t &i, vector &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 && t2-ary variables + DAI_IFVERB(2, endl<<"CHOOSE_RANDOM chose variable "<2-ary variables + DAI_IFVERB(2, endl<<"CHOOSE_RANDOM chose factor "<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 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 &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 bbpFindClampVar(const CBP &fg, bbp_cfn_t cfn, const Properties &props, Real *maxVarOut) { +pair 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 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("tol")); + bp_Props.Set("tol", ourTol); + bp_Props.Set("maxiter", props.GetAs("maxiter")); + bp_Props.Set("verbose", oneLess(props.GetAs("verbose"))); + // bp_Props.ConvertTo("updates"); + // DAI_PV(bp_Props.GetAs("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("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("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 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 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 diff --git a/tests/aliases.conf b/tests/aliases.conf index f031109..c7f62a1 100644 --- a/tests/aliases.conf +++ b/tests/aliases.conf @@ -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] -- 2.20.1