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