Improve documentation:
- bbp.h
cbp.h
evidence.h
emalg.h
# The macro definition that is found in the sources will be used.
# Use the PREDEFINED tag if you want to use a different macro definition.
-EXPAND_AS_DEFINED =
+EXPAND_AS_DEFINED = DAI_ENUM
# If the SKIP_FUNCTION_MACROS tag is set to YES (the default) then
# doxygen's preprocessor will remove all function-like macros that are alone
/// \file
-/// \brief Defines class BBP [\ref EaG09]
-/// \todo Fit more closely into libDAI framework
-/// \todo Improve documentation
-/// \author Frederik Eaton
+/// \brief Defines class BBP, which implements Back-Belief-Propagation
#ifndef ___defined_libdai_bbp_h
namespace dai {
-/// 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 );
+/// Enumeration of several cost functions that can be used with BBP
+/** \note This class is meant as a base class for BBPCostFunction, which provides additional functionality.
+ */
+DAI_ENUM(BBPCostFunctionBase,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
-/// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state
-std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters );
+
+/// Predefined cost functions that can be used with BBP
+class BBPCostFunction : public BBPCostFunctionBase {
+ public:
+ /// Returns whether this cost function depends on having a Gibbs state
+ bool needGibbsState() const;
+
+ /// Evaluates cost function in state \a stateP using the information in inference algorithm \a ia
+ Real evaluate( const InfAlg &ia, const std::vector<size_t> *stateP ) const;
+
+ /// Assignment operator
+ BBPCostFunction& operator=( const BBPCostFunctionBase &x ) {
+ if( this != &x ) {
+ (BBPCostFunctionBase)*this = x;
+ }
+ return *this;
+ }
+};
/// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
/** \author Frederik Eaton
*/
class BBP {
- protected:
- /// \name Inputs
+ private:
+ /// \name Input variables
//@{
+ /// Stores a BP_dual helper object
BP_dual _bp_dual;
+ /// Pointer to the factor graph
const FactorGraph *_fg;
+ /// Pointer to the approximate inference algorithm
const InfAlg *_ia;
//@}
- /// Number of iterations done
- size_t _iters;
-
- /// \name Outputs
+ /// \name Output variables
//@{
/// Variable factor adjoints
std::vector<Prob> _adj_psi_V;
std::vector<Prob> _adj_b_F;
//@}
- /// \name Helper quantities computed from the BP messages
+ /// \name Internal state variables
//@{
- /// _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<std::vector<Prob> > _new_adj_n;
/// Updated normalized factor->variable message adjoint (indexed [i][_I])
std::vector<std::vector<Prob> > _new_adj_m;
+ /// Unnormalized variable belief adjoints
+ std::vector<Prob> _adj_b_V_unnorm;
+ /// Unnormalized factor belief adjoints
+ std::vector<Prob> _adj_b_F_unnorm;
- /// \name Optimized indexing (for performance)
- //@{
- /// Calculates _indices, which is a cache of IndexFor @see bp.cpp
- void RegenerateInds();
+ /// _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;
+
+ /// Number of iterations done
+ size_t _iters;
+ //@}
+ /// \name Index cache management (for performance)
+ //@{
/// Index type
typedef std::vector<size_t> _ind_t;
/// Cached indices (indexed [i][_I])
std::vector<std::vector<_ind_t> > _indices;
+ /// Prepares index cache _indices
+ /** \see bp.cpp
+ */
+ void RegenerateInds();
/// Returns an index from the cache
const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
//@}
- /// \name Initialization
+ /// \name Initialization helper functions
//@{
/// Calculate T values; see eqn. (41) in [\ref EaG09]
void RegenerateT();
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)
+ /// Initialise members for factor adjoints
+ /** \pre RegenerateInputs() should be called first
+ */
void RegeneratePsiAdjoints();
- /// Initialise members for message adjoints (call after RegenerateInputs) for parallel algorithm
+ /// Initialise members for message adjoints for parallel algorithm
+ /** \pre RegenerateInputs() should be called first
+ */
void RegenerateParMessageAdjoints();
- /// Initialise members for message adjoints (call after RegenerateInputs) for sequential algorithm
+ /// Initialise members for message adjoints for sequential algorithm
/** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
* than updating _adj_n (and friends) which are unused in the sequential algorithm.
+ * \pre RegenerateInputs() should be called first
*/
void RegenerateSeqMessageAdjoints();
+ /// Called by \a init, recalculates intermediate values
+ void Regenerate();
//@}
+ /// \name Accessors/mutators
+ //@{
/// Returns reference to T value; see eqn. (41) in [\ref EaG09]
Prob & T(size_t i, size_t _I) { return _T[i][_I]; }
/// Returns constant reference to T value; see eqn. (41) in [\ref EaG09]
/// Returns constant reference to R value; see eqn. (44) in [\ref EaG09]
const Prob & R(size_t I, size_t _i, size_t _J) const { return _R[I][_i][_J]; }
+ /// Returns reference to variable->factor message adjoint
+ Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
+ /// Returns constant reference to variable->factor message adjoint
+ const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
+ /// Returns reference to factor->variable message adjoint
+ Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
+ /// Returns constant reference to factor->variable message adjoint
+ const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
+ //@}
+
/// \name Parallel algorithm
//@{
/// Calculates new variable->factor message adjoint
/// \name Sequential algorithm
//@{
- /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
+ /// Helper function for sendSeqMsgM(): increases factor->variable message adjoint by \a p and calculates the corresponding unnormalized adjoint
void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
// DISABLED BECAUSE IT IS BUGGY:
// void updateSeqMsgM( size_t i, size_t _I );
void sendSeqMsgM( size_t i, size_t _I );
//@}
- /// Calculates averaged L-1 norm of unnormalized message adjoints
+ /// 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 );
+
+ /// Calculates averaged L1 norm of unnormalized message adjoints
Real getUnMsgMag();
- /// Calculates averaged L-1 norms of current and new normalized message adjoints
+ /// Calculates averaged L1 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 ) ) );
- }
-
/// 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
/// Calculates sum of L1 norms of all normalized variable->factor message adjoints
Real getTotalMsgN();
- public:
- /// Called by \a init, recalculates intermediate values
- void Regenerate();
+ /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph \a fg
+ std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
+ /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph \a fg
+ std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
- /// Constructor
+ public:
+ /// \name Constructors/destructors
+ //@{
+ /// Construct from a InfAlg \a ia and a PropertySet \a opts
BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
props.set(opts);
}
+ //@}
- /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg
- std::vector<Prob> getZeroAdjF( 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> getZeroAdjV( const FactorGraph &fg );
-
- /// Initializes belief adjoints and initial factor adjoints and regenerates
+ /// \name Initialization
+ //@{
+ /// Initializes from given belief adjoints \a adj_b_V, \a adj_b_F and initial factor adjoints \a adj_psi_V, \a adj_psi_F
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;
Regenerate();
}
- /// Initializes belief adjoints and with zero initial factor adjoints and regenerates
+ /// Initializes from given belief adjoints \a adj_b_V and \a adj_b_F (setting initial factor adjoints to zero)
void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
init( adj_b_V, adj_b_F, getZeroAdjV(*_fg), getZeroAdjF(*_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, getZeroAdjF(*_fg));
+ /// Initializes variable belief adjoints \a adj_b_V (and sets factor belief adjoints and initial factor adjoints to zero)
+ void init
_V( const std::vector<Prob> &adj_b_V ) {
+ init( adj_b_V, getZeroAdjF(*_fg) );
}
- /// Run until change is less than given tolerance
- void run();
+ /// Initializes factor belief adjoints \a adj_b_F (and sets variable belief adjoints and initial factor adjoints to zero)
+ void init_F( const std::vector<Prob> &adj_b_F ) {
+ init( getZeroAdjV(*_fg), adj_b_F );
+ }
- /// Return number of iterations done so far
- size_t doneIters() { return _iters; }
+ /// Initializes with adjoints calculated from cost function \a cfn, and state \a stateP
+ /** Uses the internal pointer to an inference algorithm in combination with the cost function and state for initialization.
+ * \param cfn Cost function used for initialization;
+ * \param stateP is a Gibbs state and can be NULL; it will be initialised using a Gibbs run.
+ */
+ void initCostFnAdj( const BBPCostFunction &cfn, const std::vector<size_t> *stateP );
+ //@}
+ /// \name BBP Algorithm
+ //@{
+ /// Perform iterative updates until change is less than given tolerance
+ void run();
+ //@}
+
+ /// \name Query results
+ //@{
/// Returns reference to variable factor adjoint
Prob& adj_psi_V(size_t i) { return _adj_psi_V[i]; }
/// Returns constant reference to variable factor adjoint
Prob& adj_b_F(size_t I) { return _adj_b_F[I]; }
/// Returns constant reference to factor belief adjoint
const Prob& adj_b_F(size_t I) const { return _adj_b_F[I]; }
+ /// Return number of iterations done so far
+ size_t Iterations() { return _iters; }
+ //@}
- protected:
- /// Returns reference to variable->factor message adjoint
- Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
- /// Returns constant reference to variable->factor message adjoint
- const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
- /// Returns reference to factor->variable message adjoint
- Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
- /// Returns constant reference to factor->variable message adjoint
- const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
-
- public:
+ public:
/// Parameters of this algorithm
/* PROPERTIES(props,BBP) {
/// Enumeration of possible update schedules
};
-/// Enumeration of several cost functions that can be used with BBP.
-DAI_ENUM(bbp_cfn_t,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
-
-/// Initialise BBP using InfAlg, cost function, and stateP
-/** Calls bbp.init with adjoints calculated from ia.beliefV and
- * ia.beliefF. stateP is a Gibbs state and can be NULL, it will be
- * initialised using a Gibbs run of 2*fg.Iterations() iterations.
- */
-void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
-
-/// Answers question: does the given cost function depend on having a Gibbs state?
-bool needGibbsState( bbp_cfn_t cfn );
-
-/// Calculate actual value of cost function (cfn_type, stateP)
-/** This function returns the actual value of the cost function whose
- * gradient with respect to singleton beliefs is given by
- * gibbsToB1Adj on the same arguments
- */
-Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
-
-/// Function to test the validity of adjoints computed by BBP given a state for each variable using numerical derivatives.
-/** Factors containing a variable are multiplied by psi_1 adjustments to verify accuracy of _adj_psi_V.
- * \param bp BP object.
- * \param state Global state of all variables.
- * \param bbp_props BBP Properties.
- * \param cfn Cost function to be used.
- * \param h controls size of perturbation.
+/// Function to verify the validity of adjoints computed by BBP using numerical differentiation.
+/** Factors containing a variable are multiplied by small adjustments to verify accuracy of calculated variable factor adjoints.
+ * \param bp BP object;
+ * \param state Global state of all variables;
+ * \param bbp_props BBP parameters;
+ * \param cfn Cost function to be used;
+ * \param h Size of perturbation.
+ * \relates BBP
*/
-Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h );
+Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h );
} // end of namespace dai
/// Find a variable to clamp using BBP (goes with maximum adjoint)
/// \see BBP
-std::pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut );
+std::pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real *maxVarOut );
/// Class for CBP (Clamped Belief Propagation)
boost::shared_ptr<std::ofstream> _clamp_ofstream;
/// Returns BBP cost function used
- bbp_cfn_t BBP_cost_function() { return props.bbp_cfn; }
+ BBPCostFunction BBP_cost_function() { return props.bbp_cfn; }
/// Prints beliefs, variables and partition sum, in case of a debugging build
void printDebugInfo();
/// Properties to pass to BBP
PropertySet bbp_props;
/// Cost function to use for BBP
- bbp_cfn_t bbp_cfn;
+ BBPCostFunction bbp_cfn;
/// Random seed
size_t rand_seed = 0;
/// Properties to pass to BBP
PropertySet bbp_props;
/// Cost function to use for BBP
- bbp_cfn_t bbp_cfn;
+ BBPCostFunction bbp_cfn;
/// Random seed
size_t rand_seed;
/// If non-empty, write clamping choices to this file
return os;\
}\
\
- private:\
+ protected:\
value v;\
};
};
+/// Runs Gibbs sampling for \a iters iterations (of which \a burnin for burn-in) on FactorGraph \a fg, and returns the resulting state
+/** \relates Gibbs
+ */
+std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t iters, size_t burnin=0 );
+
+
} // end of namespace dai
/// \file
/// \brief Defines class TreeEP, which implements Tree Expectation Propagation
+/// \todo Clean up the TreeEP code
#ifndef __defined_libdai_treeep_h
typedef BipartiteGraph::Neighbor Neighbor;
-Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
- DAI_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;
-}
-
-
-std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters ) {
- PropertySet gibbsProps;
- gibbsProps.Set("iters", iters);
- gibbsProps.Set("verbose", size_t(0));
- Gibbs gibbs( ia.fg(), gibbsProps );
- gibbs.run();
- return gibbs.state();
-}
-
-
/// Returns the entry of the I'th factor corresponding to a global state
size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
size_t f_entry = 0;
}
+bool BBPCostFunction::needGibbsState() const {
+ switch( (size_t)(*this) ) {
+ case CFN_GIBBS_B:
+ case CFN_GIBBS_B2:
+ case CFN_GIBBS_EXP:
+ case CFN_GIBBS_B_FACTOR:
+ case CFN_GIBBS_B2_FACTOR:
+ case CFN_GIBBS_EXP_FACTOR:
+ return true;
+ default:
+ return false;
+ }
+}
+
+
+Real BBPCostFunction::evaluate( const InfAlg &ia, const vector<size_t> *stateP ) const {
+ Real cf = 0.0;
+ const FactorGraph &fg = ia.fg();
+
+ switch( (size_t)(*this) ) {
+ case CFN_BETHE_ENT: // ignores state
+ cf = -ia.logZ();
+ break;
+ case CFN_VAR_ENT: // ignores state
+ for( size_t i = 0; i < fg.nrVars(); i++ )
+ cf += -ia.beliefV(i).entropy();
+ break;
+ case CFN_FACTOR_ENT: // ignores state
+ for( size_t I = 0; I < fg.nrFactors(); I++ )
+ cf += -ia.beliefF(I).entropy();
+ break;
+ case CFN_GIBBS_B:
+ case CFN_GIBBS_B2:
+ case CFN_GIBBS_EXP: {
+ DAI_ASSERT( stateP != NULL );
+ vector<size_t> state = *stateP;
+ DAI_ASSERT( state.size() == fg.nrVars() );
+ for( size_t i = 0; i < fg.nrVars(); i++ ) {
+ Real b = ia.beliefV(i)[state[i]];
+ switch( (size_t)(*this) ) {
+ case CFN_GIBBS_B:
+ cf += b;
+ break;
+ case CFN_GIBBS_B2:
+ cf += b * b / 2.0;
+ break;
+ case CFN_GIBBS_EXP:
+ cf += exp( b );
+ break;
+ default:
+ DAI_THROW(UNKNOWN_ENUM_VALUE);
+ }
+ }
+ break;
+ } case CFN_GIBBS_B_FACTOR:
+ case CFN_GIBBS_B2_FACTOR:
+ case CFN_GIBBS_EXP_FACTOR: {
+ DAI_ASSERT( stateP != NULL );
+ vector<size_t> state = *stateP;
+ DAI_ASSERT( state.size() == fg.nrVars() );
+ for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+ size_t x_I = getFactorEntryForState( fg, I, state );
+ Real b = ia.beliefF(I)[x_I];
+ switch( (size_t)(*this) ) {
+ case CFN_GIBBS_B_FACTOR:
+ cf += b;
+ break;
+ case CFN_GIBBS_B2_FACTOR:
+ cf += b * b / 2.0;
+ break;
+ case CFN_GIBBS_EXP_FACTOR:
+ cf += exp( b );
+ break;
+ default:
+ DAI_THROW(UNKNOWN_ENUM_VALUE);
+ }
+ }
+ break;
+ } default:
+ DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(*this));
+ }
+ return cf;
+}
+
+
#define LOOP_ij(body) { \
size_t i_states = _fg->var(i).states(); \
size_t j_states = _fg->var(j).states(); \
}
+void BBP::Regenerate() {
+ RegenerateInds();
+ RegenerateT();
+ RegenerateU();
+ RegenerateS();
+ RegenerateR();
+ RegenerateInputs();
+ RegeneratePsiAdjoints();
+ if( props.updates == Properties::UpdateType::PAR )
+ RegenerateParMessageAdjoints();
+ 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];
Prob &new_adj_n_iI = _new_adj_n[i][_I];
}
+Prob BBP::unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
+ DAI_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;
+}
+
+
Real BBP::getUnMsgMag() {
Real s = 0.0;
size_t e = 0;
}
-void BBP::Regenerate() {
- RegenerateInds();
- RegenerateT();
- RegenerateU();
- RegenerateS();
- RegenerateR();
- RegenerateInputs();
- RegeneratePsiAdjoints();
- if( props.updates == Properties::UpdateType::PAR )
- RegenerateParMessageAdjoints();
- else
- RegenerateSeqMessageAdjoints();
- _iters = 0;
-}
-
-
std::vector<Prob> BBP::getZeroAdjF( const FactorGraph &fg ) {
vector<Prob> adj_2;
adj_2.reserve( fg.nrFactors() );
}
+void BBP::initCostFnAdj( const BBPCostFunction &cfn, const vector<size_t> *stateP ) {
+ const FactorGraph &fg = _ia->fg();
+
+ switch( (size_t)cfn ) {
+ case BBPCostFunction::CFN_BETHE_ENT: {
+ vector<Prob> b1_adj;
+ vector<Prob> b2_adj;
+ vector<Prob> psi1_adj;
+ vector<Prob> psi2_adj;
+ b1_adj.reserve( fg.nrVars() );
+ psi1_adj.reserve( fg.nrVars() );
+ b2_adj.reserve( fg.nrFactors() );
+ psi2_adj.reserve( fg.nrFactors() );
+ for( size_t i = 0; i < fg.nrVars(); i++ ) {
+ size_t dim = fg.var(i).states();
+ int c = fg.nbV(i).size();
+ Prob p(dim,0.0);
+ for( size_t xi = 0; xi < dim; xi++ )
+ p[xi] = (1 - c) * (1 + log( _ia->beliefV(i)[xi] ));
+ b1_adj.push_back( p );
+
+ 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++ ) {
+ size_t dim = fg.factor(I).states();
+ Prob p( dim, 0.0 );
+ 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++ )
+ p[xI] = -_ia->beliefF(I)[xI] / fg.factor(I).p()[xI];
+ psi2_adj.push_back( p );
+ }
+ init( b1_adj, b2_adj, psi1_adj, psi2_adj );
+ break;
+ } case BBPCostFunction::CFN_FACTOR_ENT: {
+ vector<Prob> b2_adj;
+ b2_adj.reserve( fg.nrFactors() );
+ for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+ size_t dim = fg.factor(I).states();
+ Prob p( dim, 0.0 );
+ for( size_t xI = 0; xI < dim; xI++ ) {
+ Real bIxI = _ia->beliefF(I)[xI];
+ if( bIxI < 1.0e-15 )
+ p[xI] = -1.0e10;
+ else
+ p[xI] = 1 + log( bIxI );
+ }
+ b2_adj.push_back(p);
+ }
+ init_F( b2_adj );
+ break;
+ } case BBPCostFunction::CFN_VAR_ENT: {
+ vector<Prob> b1_adj;
+ b1_adj.reserve( fg.nrVars() );
+ for( size_t i = 0; i < fg.nrVars(); i++ ) {
+ size_t dim = fg.var(i).states();
+ Prob p( dim, 0.0 );
+ for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
+ Real bixi = _ia->beliefV(i)[xi];
+ if( bixi < 1.0e-15 )
+ p[xi] = -1.0e10;
+ else
+ p[xi] = 1 + log( bixi );
+ }
+ b1_adj.push_back( p );
+ }
+ init_V( b1_adj );
+ break;
+ } case BBPCostFunction::CFN_GIBBS_B:
+ case BBPCostFunction::CFN_GIBBS_B2:
+ case BBPCostFunction::CFN_GIBBS_EXP: {
+ // cost functions that use Gibbs sample, summing over variable marginals
+ vector<size_t> state;
+ if( stateP == NULL )
+ state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+ else
+ state = *stateP;
+ DAI_ASSERT( state.size() == fg.nrVars() );
+
+ vector<Prob> b1_adj;
+ b1_adj.reserve(fg.nrVars());
+ for( size_t i = 0; i < state.size(); i++ ) {
+ size_t n = fg.var(i).states();
+ Prob delta( n, 0.0 );
+ DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
+ Real b = _ia->beliefV(i)[state[i]];
+ switch( (size_t)cfn ) {
+ case BBPCostFunction::CFN_GIBBS_B:
+ delta[state[i]] = 1.0;
+ break;
+ case BBPCostFunction::CFN_GIBBS_B2:
+ delta[state[i]] = b;
+ break;
+ case BBPCostFunction::CFN_GIBBS_EXP:
+ delta[state[i]] = exp(b);
+ break;
+ default:
+ DAI_THROW(UNKNOWN_ENUM_VALUE);
+ }
+ b1_adj.push_back( delta );
+ }
+ init_V( b1_adj );
+ break;
+ } case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+ case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+ case BBPCostFunction::CFN_GIBBS_EXP_FACTOR: {
+ // cost functions that use Gibbs sample, summing over factor marginals
+ vector<size_t> state;
+ if( stateP == NULL )
+ state = getGibbsState( _ia->fg(), 2*_ia->Iterations() );
+ else
+ state = *stateP;
+ DAI_ASSERT( state.size() == fg.nrVars() );
+
+ vector<Prob> b2_adj;
+ b2_adj.reserve( fg.nrVars() );
+ for( size_t I = 0; I < fg.nrFactors(); I++ ) {
+ size_t n = fg.factor(I).states();
+ Prob delta( n, 0.0 );
+
+ size_t x_I = getFactorEntryForState( fg, I, state );
+ DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
+
+ Real b = _ia->beliefF(I)[x_I];
+ switch( (size_t)cfn ) {
+ case BBPCostFunction::CFN_GIBBS_B_FACTOR:
+ delta[x_I] = 1.0;
+ break;
+ case BBPCostFunction::CFN_GIBBS_B2_FACTOR:
+ delta[x_I] = b;
+ break;
+ case BBPCostFunction::CFN_GIBBS_EXP_FACTOR:
+ delta[x_I] = exp( b );
+ break;
+ default:
+ DAI_THROW(UNKNOWN_ENUM_VALUE);
+ }
+ b2_adj.push_back( delta );
+ }
+ init_F( b2_adj );
+ break;
+ } default:
+ DAI_THROW(UNKNOWN_ENUM_VALUE);
+ }
+}
+
+
void BBP::run() {
typedef BBP::Properties::UpdateType UT;
Real tol = props.tol;
}
}
if( props.verbose >= 3 )
- cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl;
+ cerr << "BBP::run() took " << toc()-tic << " seconds " << Iterations() << " iterations" << endl;
}
-Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h ) {
+Real numericBBPTest( const InfAlg &bp, const vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h ) {
+ BBP bbp( &bp, bbp_props );
// calculate the value of the unperturbed cost function
- Real cf0 = getCostFn( bp, cfn, state );
-
+ Real cf0 = cfn.evaluate( bp, state );
// run BBP to estimate adjoints
- BBP bbp( &bp, bbp_props );
- initBBPCostFnAdj( bbp, bp, cfn, state );
+ bbp.initCostFnAdj( cfn, state );
bbp.run();
Real d = 0;
bp_prb->run();
// calculate new value of cost function
- Real cf_prb = getCostFn( *bp_prb, cfn, state );
+ Real cf_prb = cfn.evaluate( *bp_prb, state );
// use to estimate adjoint for i
adj_est.push_back( (cf_prb - cf0) / h );
}
-bool needGibbsState( bbp_cfn_t cfn ) {
- switch( (size_t)cfn ) {
- case bbp_cfn_t::CFN_GIBBS_B:
- case bbp_cfn_t::CFN_GIBBS_B2:
- case bbp_cfn_t::CFN_GIBBS_EXP:
- case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
- return true;
- default:
- return false;
- }
-}
-
-
-void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
- const FactorGraph &fg = ia.fg();
-
- switch( (size_t)cfn_type ) {
- case bbp_cfn_t::CFN_BETHE_ENT: {
- vector<Prob> b1_adj;
- vector<Prob> b2_adj;
- vector<Prob> psi1_adj;
- vector<Prob> psi2_adj;
- b1_adj.reserve( fg.nrVars() );
- psi1_adj.reserve( fg.nrVars() );
- b2_adj.reserve( fg.nrFactors() );
- psi2_adj.reserve( fg.nrFactors() );
- for( size_t i = 0; i < fg.nrVars(); i++ ) {
- size_t dim = fg.var(i).states();
- int c = fg.nbV(i).size();
- Prob p(dim,0.0);
- for( size_t xi = 0; xi < dim; xi++ )
- p[xi] = (1 - c) * (1 + log( ia.beliefV(i)[xi] ));
- b1_adj.push_back( p );
-
- 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++ ) {
- size_t dim = fg.factor(I).states();
- Prob p( dim, 0.0 );
- 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++ )
- 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 );
- break;
- } case bbp_cfn_t::CFN_FACTOR_ENT: {
- vector<Prob> b2_adj;
- b2_adj.reserve( fg.nrFactors() );
- for( size_t I = 0; I < fg.nrFactors(); I++ ) {
- size_t dim = fg.factor(I).states();
- Prob p( dim, 0.0 );
- for( size_t xI = 0; xI < dim; xI++ ) {
- Real bIxI = ia.beliefF(I)[xI];
- if( bIxI < 1.0e-15 )
- p[xI] = -1.0e10;
- else
- p[xI] = 1 + log( bIxI );
- }
- b2_adj.push_back(p);
- }
- bbp.init( bbp.getZeroAdjV(fg), b2_adj );
- break;
- } case bbp_cfn_t::CFN_VAR_ENT: {
- vector<Prob> b1_adj;
- b1_adj.reserve( fg.nrVars() );
- for( size_t i = 0; i < fg.nrVars(); i++ ) {
- size_t dim = fg.var(i).states();
- Prob p( dim, 0.0 );
- for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
- Real bixi = ia.beliefV(i)[xi];
- if( bixi < 1.0e-15 )
- p[xi] = -1.0e10;
- else
- p[xi] = 1 + log( bixi );
- }
- b1_adj.push_back( p );
- }
- bbp.init( b1_adj );
- break;
- } case bbp_cfn_t::CFN_GIBBS_B:
- case bbp_cfn_t::CFN_GIBBS_B2:
- case bbp_cfn_t::CFN_GIBBS_EXP: {
- // cost functions that use Gibbs sample, summing over variable marginals
- vector<size_t> state;
- if( stateP == NULL )
- state = getGibbsState( ia, 2*ia.Iterations() );
- else
- state = *stateP;
- DAI_ASSERT( state.size() == fg.nrVars() );
-
- vector<Prob> b1_adj;
- b1_adj.reserve(fg.nrVars());
- for( size_t i = 0; i < state.size(); i++ ) {
- size_t n = fg.var(i).states();
- Prob delta( n, 0.0 );
- DAI_ASSERT(/*0<=state[i] &&*/ state[i] < n);
- Real b = ia.beliefV(i)[state[i]];
- switch( (size_t)cfn_type ) {
- case bbp_cfn_t::CFN_GIBBS_B:
- delta[state[i]] = 1.0;
- break;
- case bbp_cfn_t::CFN_GIBBS_B2:
- delta[state[i]] = b;
- break;
- case bbp_cfn_t::CFN_GIBBS_EXP:
- delta[state[i]] = exp(b);
- break;
- default:
- DAI_THROW(UNKNOWN_ENUM_VALUE);
- }
- b1_adj.push_back( delta );
- }
- bbp.init( b1_adj );
- break;
- } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
- // cost functions that use Gibbs sample, summing over factor marginals
- vector<size_t> state;
- if( stateP == NULL )
- state = getGibbsState( ia, 2*ia.Iterations() );
- else
- state = *stateP;
- DAI_ASSERT( state.size() == fg.nrVars() );
-
- vector<Prob> b2_adj;
- b2_adj.reserve( fg.nrVars() );
- for( size_t I = 0; I < fg.nrFactors(); I++ ) {
- size_t n = fg.factor(I).states();
- Prob delta( n, 0.0 );
-
- size_t x_I = getFactorEntryForState( fg, I, state );
- DAI_ASSERT(/*0<=x_I &&*/ x_I < n);
-
- Real b = ia.beliefF(I)[x_I];
- switch( (size_t)cfn_type ) {
- case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
- delta[x_I] = 1.0;
- break;
- case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
- delta[x_I] = b;
- break;
- case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
- delta[x_I] = exp( b );
- break;
- default:
- DAI_THROW(UNKNOWN_ENUM_VALUE);
- }
- b2_adj.push_back( delta );
- }
- bbp.init( bbp.getZeroAdjV(fg), b2_adj );
- break;
- } default:
- DAI_THROW(UNKNOWN_ENUM_VALUE);
- }
-}
-
-
-Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
- Real cf = 0.0;
- const FactorGraph &fg = ia.fg();
-
- switch( (size_t)cfn_type ) {
- case bbp_cfn_t::CFN_BETHE_ENT: // ignores state
- cf = -ia.logZ();
- break;
- case bbp_cfn_t::CFN_VAR_ENT: // ignores state
- for( size_t i = 0; i < fg.nrVars(); i++ )
- cf += -ia.beliefV(i).entropy();
- break;
- case bbp_cfn_t::CFN_FACTOR_ENT: // ignores state
- for( size_t I = 0; I < fg.nrFactors(); I++ )
- cf += -ia.beliefF(I).entropy();
- break;
- case bbp_cfn_t::CFN_GIBBS_B:
- case bbp_cfn_t::CFN_GIBBS_B2:
- case bbp_cfn_t::CFN_GIBBS_EXP: {
- DAI_ASSERT( stateP != NULL );
- vector<size_t> state = *stateP;
- DAI_ASSERT( state.size() == fg.nrVars() );
- for( size_t i = 0; i < fg.nrVars(); i++ ) {
- Real b = ia.beliefV(i)[state[i]];
- switch( (size_t)cfn_type ) {
- case bbp_cfn_t::CFN_GIBBS_B:
- cf += b;
- break;
- case bbp_cfn_t::CFN_GIBBS_B2:
- cf += b * b / 2.0;
- break;
- case bbp_cfn_t::CFN_GIBBS_EXP:
- cf += exp( b );
- break;
- default:
- DAI_THROW(UNKNOWN_ENUM_VALUE);
- }
- }
- break;
- } case bbp_cfn_t::CFN_GIBBS_B_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_B2_FACTOR:
- case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: {
- DAI_ASSERT( stateP != NULL );
- vector<size_t> state = *stateP;
- DAI_ASSERT( state.size() == fg.nrVars() );
- for( size_t I = 0; I < fg.nrFactors(); I++ ) {
- size_t x_I = getFactorEntryForState( fg, I, state );
- Real 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.0;
- break;
- case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR:
- cf += exp( b );
- break;
- default:
- DAI_THROW(UNKNOWN_ENUM_VALUE);
- }
- }
- break;
- } default:
- DAI_THROWE(UNKNOWN_ENUM_VALUE, "Unknown cost function " + std::string(cfn_type));
- }
- return cf;
-}
-
-
} // end of namespace dai
#include <dai/util.h>
#include <dai/properties.h>
-
+#include <dai/gibbs.h>
#include <dai/bp.h>
#include <dai/cbp.h>
#include <dai/bbp.h>
ChooseMethod()==Properties::ChooseMethodType::CHOOSE_BP_CFN ) {
bool doL1 = (ChooseMethod() == Properties::ChooseMethodType::CHOOSE_BP_L1);
vector<size_t> state;
- if( !doL1 && needGibbsState( BBP_cost_function() ) )
- state = getGibbsState( *bp, 2*bp->Iterations() );
+ if( !doL1 && BBP_cost_function().needGibbsState() )
+ state = getGibbsState( bp->fg(), 2*bp->Iterations() );
// try clamping each variable manually
DAI_ASSERT( Clamping() == Properties::ClampType::CLAMP_VAR );
Real max_cost = 0.0;
for( size_t j = 0; j < nrVars(); j++ )
cost += dist( bp->beliefV(j), bp1->beliefV(j), Prob::DISTL1 );
else
- cost = getCostFn( *bp1, BBP_cost_function(), &state );
+ cost = BBP_cost_function().evaluate( *bp1, &state );
if( cost > max_cost || win_k == -1 ) {
max_cost = cost;
win_k = k;
* creates and runs a BBP object, finds best variable, returns
* (variable,state) pair for clamping
*/
-pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, bbp_cfn_t cfn, Real *maxVarOut ) {
+pair<size_t, size_t> bbpFindClampVar( const InfAlg &in_bp, bool clampingVar, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real *maxVarOut ) {
BBP bbp( &in_bp, bbp_props );
- initBBPCostFnAdj( bbp, in_bp, cfn, NULL );
+ bbp.initCostFnAdj( cfn, NULL );
bbp.run();
// find and return the (variable,state) with the largest adj_psi_V
recursion = opts.getStringAs<RecurseType>("recursion");
clamp = opts.getStringAs<ClampType>("clamp");
bbp_props = opts.getStringAs<PropertySet>("bbp_props");
- bbp_cfn = opts.getStringAs<bbp_cfn_t>("bbp_cfn");
+ bbp_cfn = opts.getStringAs<BBPCostFunction>("bbp_cfn");
if(opts.hasKey("rand_seed")) {
rand_seed = opts.getStringAs<size_t>("rand_seed");
} else {
}
+std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t iters, size_t burnin ) {
+ PropertySet gibbsProps;
+ gibbsProps.Set("iters", iters);
+ gibbsProps.Set("burnin", burnin);
+ gibbsProps.Set("verbose", size_t(0));
+ Gibbs gibbs( fg, gibbsProps );
+ gibbs.run();
+ return gibbs.state();
+}
+
+
} // end of namespace dai
updates = BBP::Properties::UpdateType::PAR;
break;
}
- bbp_cfn_t cfn;
+ BBPCostFunction cfn;
switch( (t / 5) % 9 ) {
case 0:
- cfn = bbp_cfn_t::CFN_GIBBS_B;
+ cfn = BBPCostFunction::CFN_GIBBS_B;
break;
case 1:
- cfn = bbp_cfn_t::CFN_GIBBS_B2;
+ cfn = BBPCostFunction::CFN_GIBBS_B2;
break;
case 2:
- cfn = bbp_cfn_t::CFN_GIBBS_EXP;
+ cfn = BBPCostFunction::CFN_GIBBS_EXP;
break;
case 3:
- cfn = bbp_cfn_t::CFN_GIBBS_B_FACTOR;
+ cfn = BBPCostFunction::CFN_GIBBS_B_FACTOR;
break;
case 4:
- cfn = bbp_cfn_t::CFN_GIBBS_B2_FACTOR;
+ cfn = BBPCostFunction::CFN_GIBBS_B2_FACTOR;
break;
case 5:
- cfn = bbp_cfn_t::CFN_GIBBS_EXP_FACTOR;
+ cfn = BBPCostFunction::CFN_GIBBS_EXP_FACTOR;
break;
case 6:
- cfn = bbp_cfn_t::CFN_VAR_ENT;
+ cfn = BBPCostFunction::CFN_VAR_ENT;
break;
case 7:
- cfn = bbp_cfn_t::CFN_FACTOR_ENT;
+ cfn = BBPCostFunction::CFN_FACTOR_ENT;
break;
case 8:
- cfn = bbp_cfn_t::CFN_BETHE_ENT;
+ cfn = BBPCostFunction::CFN_BETHE_ENT;
break;
}
projects / libdai.git / commitdiff