From b2d044defebc5a9104765d53c81a4bc291898fcf Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Tue, 4 Aug 2009 15:36:37 +0200 Subject: [PATCH] Even more cleanup of BBP code --- include/dai/bbp.h | 80 ++- src/bbp.cpp | 1172 ++++++++++++++++++++++---------------------- tests/aliases.conf | 2 +- tests/testbbp.cpp | 18 +- 4 files changed, 621 insertions(+), 651 deletions(-) diff --git a/include/dai/bbp.h b/include/dai/bbp.h index bf68121..2b41568 100644 --- a/include/dai/bbp.h +++ b/include/dai/bbp.h @@ -21,7 +21,6 @@ /// \file /// \brief Defines class BBP [\ref EaG09] /// \todo Improve documentation -/// \todo Clean up /// \todo Debug clean_updates @@ -42,11 +41,11 @@ namespace dai { -/// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg -std::vector get_zero_adj_F( const FactorGraph& fg ); +/// 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 ); -/// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph fg -std::vector get_zero_adj_V( const FactorGraph& fg ); +/// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state +std::vector getGibbsState( const InfAlg &ia, size_t iters ); /// Implements BBP (Back-Belief-Propagation) [\ref EaG09] @@ -178,6 +177,19 @@ class BBP { void doParUpdate(); //@} + /// @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 ); + /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint + void setSeqMsgM( size_t i, size_t _I, const Prob &p ); + /// 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 ); + //@} + /// Calculates averaged L-1 norm of unnormalized message adjoints Real getUnMsgMag(); /// Calculates averaged L-1 norms of current and new normalized message adjoints @@ -191,22 +203,8 @@ class BBP { _adj_b_F.push_back( Prob( _fg->factor(I).states(), Real( 0.0 ) ) ); } - /// @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 @@ -225,6 +223,11 @@ class BBP { 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 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 getZeroAdjV( const FactorGraph &fg ); + /// Initializes belief adjoints and initial factor adjoints and regenerates void init( const std::vector &adj_b_V, const std::vector &adj_b_F, const std::vector &adj_psi_V, const std::vector &adj_psi_F ) { _adj_b_V = adj_b_V; @@ -236,12 +239,12 @@ class BBP { /// Initializes belief adjoints and with zero initial factor adjoints and regenerates void init( const std::vector &adj_b_V, const std::vector &adj_b_F ) { - init( adj_b_V, adj_b_F, get_zero_adj_V(*_fg), get_zero_adj_F(*_fg) ); + 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 &adj_b_V ) { - init(adj_b_V, get_zero_adj_F(*_fg)); + init(adj_b_V, getZeroAdjF(*_fg)); } /// Run until change is less than given tolerance @@ -317,43 +320,32 @@ class BBP { /* }}} END OF GENERATED CODE */ }; -/// Cost functions. Not used by BBP class, only used by following functions. -DAI_ENUM(bbp_cfn_t,cfn_gibbs_b,cfn_gibbs_b2,cfn_gibbs_exp,cfn_gibbs_b_factor,cfn_gibbs_b2_factor,cfn_gibbs_exp_factor,cfn_var_ent,cfn_factor_ent,cfn_bethe_ent); + +/// 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* stateP); +void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector *stateP ); -/// Answers question: Does given cost function depend on having a Gibbs state? -bool needGibbsState(bbp_cfn_t cfn); +/// 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 *stateP); - -/// Function to test the validity of adjoints computed by BBP -/** given a state for each variable, use numerical derivatives - * (multiplying a factor containing a variable by psi_1 adjustments) - * to verify accuracy of _adj_psi_V. - * 'h' controls size of perturbation. - * 'bbpTol' controls tolerance of BBP run. - */ -double numericBBPTest(const InfAlg& bp, const std::vector *state, const PropertySet& bbp_props, bbp_cfn_t cfn, double h); - -// ---------------------------------------------------------------- -// Utility functions, some of which are used elsewhere +Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector *stateP ); -/// 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 getGibbsState(const InfAlg& ia, size_t iters); +/// 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 h controls size of perturbation. + */ +double numericBBPTest( const InfAlg &bp, const std::vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ); } // end of namespace dai diff --git a/src/bbp.cpp b/src/bbp.cpp index 06750db..1e55c2e 100644 --- a/src/bbp.cpp +++ b/src/bbp.cpp @@ -47,6 +47,16 @@ Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) { } +std::vector 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(); +} + + size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector &state ) { size_t f_entry = 0; for( int _j = fg.nbF(I).size() - 1; _j >= 0; _j-- ) { @@ -60,269 +70,43 @@ size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector *stateP ) { - const FactorGraph &fg = ia.fg(); - - switch( (size_t)cfn_type ) { - case bbp_cfn_t::cfn_bethe_ent: { - vector b1_adj; - vector b2_adj; - vector psi1_adj; - vector psi2_adj; - b1_adj.reserve(fg.nrVars()); - psi1_adj.reserve(fg.nrVars()); - b2_adj.reserve(fg.nrFactors()); - psi2_adj.reserve(fg.nrFactors()); - for( size_t i = 0; i < 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 b2_adj; - b2_adj.reserve(fg.nrFactors()); - for( size_t I = 0; I < fg.nrFactors(); I++ ) { - size_t dim = fg.factor(I).states(); - Prob p(dim,0.0); - for( size_t xI = 0; xI < dim; xI++ ) { - double bIxI = 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(get_zero_adj_V(fg), b2_adj); - break; - } - case bbp_cfn_t::cfn_var_ent: { - vector b1_adj; - b1_adj.reserve(fg.nrVars()); - for(size_t i=0; i state; - if(stateP==NULL) - state = getGibbsState(ia,2*ia.Iterations()); - else - state = *stateP; - assert(state.size()==fg.nrVars()); - - vector b1_adj; - b1_adj.reserve(fg.nrVars()); - for( size_t i = 0; i < state.size(); i++ ) { - size_t n = fg.var(i).states(); - Prob delta(n,0.0); - assert(/*0<=state[i] &&*/ state[i] < n); - double b = 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(INTERNAL_ERROR); - } - 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 state; - if(stateP==NULL) - state = getGibbsState(ia,2*ia.Iterations()); - else - state = *stateP; - assert(state.size()==fg.nrVars()); - - vector b2_adj; - b2_adj.reserve(fg.nrVars()); - for( size_t I = 0; I < fg.nrFactors(); I++ ) { - size_t n = fg.factor(I).states(); - Prob delta(n,0.0); - - size_t x_I = getFactorEntryForState(fg, I, state); - assert(/*0<=x_I &&*/ x_I < n); - - double b = ia.beliefF(I)[x_I]; - switch( (size_t)cfn_type ) { - case bbp_cfn_t::cfn_gibbs_b_factor: - delta[x_I] = 1.0; break; - case bbp_cfn_t::cfn_gibbs_b2_factor: - delta[x_I] = b; break; - case bbp_cfn_t::cfn_gibbs_exp_factor: - delta[x_I] = exp(b); break; - default: DAI_THROW(INTERNAL_ERROR); - } - b2_adj.push_back(delta); - } - bbp.init(get_zero_adj_V(fg), b2_adj); - break; - } - default: DAI_THROW(UNKNOWN_ENUM_VALUE); - } -} - - -Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector *stateP ) { - double 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 state = *stateP; - assert(state.size()==fg.nrVars()); - for( size_t i = 0; i < fg.nrVars(); i++ ) { - double b = ia.beliefV(i)[state[i]]; - switch( (size_t)cfn_type ) { - case bbp_cfn_t::cfn_gibbs_b: cf += b; break; - case bbp_cfn_t::cfn_gibbs_b2: cf += b*b/2; break; - case bbp_cfn_t::cfn_gibbs_exp: cf += exp(b); break; - default: DAI_THROW(INTERNAL_ERROR); - } - } - break; - } - case bbp_cfn_t::cfn_gibbs_b_factor: - case bbp_cfn_t::cfn_gibbs_b2_factor: - case bbp_cfn_t::cfn_gibbs_exp_factor: { - assert(stateP != NULL); - vector state = *stateP; - assert(state.size()==fg.nrVars()); - for( size_t I = 0; I < fg.nrFactors(); I++ ) { - size_t x_I = getFactorEntryForState(fg, I, state); - double b = ia.beliefF(I)[x_I]; - switch((size_t)cfn_type) { - case bbp_cfn_t::cfn_gibbs_b_factor: cf += b; break; - case bbp_cfn_t::cfn_gibbs_b2_factor: cf += b*b/2; break; - case bbp_cfn_t::cfn_gibbs_exp_factor: cf += exp(b); break; - default: DAI_THROW(INTERNAL_ERROR); - } - } - break; - } - default: - cerr << "Unknown cost function " << cfn_type << endl; - DAI_THROW(UNKNOWN_ENUM_VALUE); - } - return cf; -} - - -vector get_zero_adj_F( const FactorGraph &fg ) { - vector 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) ) ); - return adj_2; -} - - -vector get_zero_adj_V( const FactorGraph &fg ) { - vector 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) ) ); - return adj_1; +#define LOOP_ij(body) { \ + size_t i_states = _fg->var(i).states(); \ + size_t j_states = _fg->var(j).states(); \ + if(_fg->var(i) > _fg->var(j)) { \ + size_t xij=0; \ + for(size_t xi=0; xivar(i).states(); \ - size_t j_states = _fg->var(j).states(); \ - if(_fg->var(i) > _fg->var(j)) { \ - size_t xij=0; \ - for(size_t xi=0; xi _ind_t; + // typedef std::vector _ind_t; // std::vector > _indices; - _indices.resize(_fg->nrVars()); + _indices.resize( _fg->nrVars() ); for( size_t i = 0; i < _fg->nrVars(); i++ ) { - _indices[i].reserve(_fg->nbV(i).size()); - foreach(const Neighbor &I, _fg->nbV(i)) { + _indices[i].clear(); + _indices[i].reserve( _fg->nbV(i).size() ); + foreach( const Neighbor &I, _fg->nbV(i) ) { _ind_t index; - index.reserve(_fg->factor(I).states()); - for(IndexFor k( _fg->var(i), _fg->factor(I).vars() ); k >= 0; ++k) { - index.push_back(k); - } - _indices[i].push_back(index); + index.reserve( _fg->factor(I).states() ); + for( IndexFor k(_fg->var(i), _fg->factor(I).vars()); k >= 0; ++k ) + index.push_back( k ); + _indices[i].push_back( index ); } } } @@ -330,7 +114,6 @@ void BBP::RegenerateInds() { void BBP::RegenerateT() { // _T[i][_I] - _T.clear(); _T.resize( _fg->nrVars() ); for( size_t i = 0; i < _fg->nrVars(); i++ ) { _T[i].resize( _fg->nbV(i).size() ); @@ -347,7 +130,6 @@ void BBP::RegenerateT() { void BBP::RegenerateU() { // _U[I][_i] - _U.clear(); _U.resize( _fg->nrFactors() ); for( size_t I = 0; I < _fg->nrFactors(); I++ ) { _U[I].resize( _fg->nbF(I).size() ); @@ -369,7 +151,6 @@ void BBP::RegenerateU() { void BBP::RegenerateS() { // _S[i][_I][_j] - _S.clear(); _S.resize( _fg->nrVars() ); for( size_t i = 0; i < _fg->nrVars(); i++ ) { _S[i].resize( _fg->nbV(i).size() ); @@ -387,7 +168,6 @@ void BBP::RegenerateS() { } } // "Marginalize" onto i|j (unnormalized) - // XXX improve efficiency? Prob marg; marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p(); _S[i][I.iter][j.iter] = marg; @@ -399,7 +179,6 @@ void BBP::RegenerateS() { void BBP::RegenerateR() { // _R[I][_i][_J] - _R.clear(); _R.resize( _fg->nrFactors() ); for( size_t I = 0; I < _fg->nrFactors(); I++ ) { _R[I].resize( _fg->nbF(I).size() ); @@ -410,7 +189,7 @@ void BBP::RegenerateR() { 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); + prod *= _bp_dual.msgM( i, K.iter ); _R[I][i.iter][J.iter] = prod; } } @@ -421,109 +200,216 @@ void BBP::RegenerateR() { void BBP::RegenerateInputs() { _adj_b_V_unnorm.clear(); - _adj_b_V_unnorm.reserve(_fg->nrVars()); - for( size_t i = 0; i < _fg->nrVars(); i++ ) { - _adj_b_V_unnorm.push_back( - unnormAdjoint(_bp_dual.beliefV(i).p(), - _bp_dual.beliefVZ(i), _adj_b_V[i])); - } + _adj_b_V_unnorm.reserve( _fg->nrVars() ); + for( size_t i = 0; i < _fg->nrVars(); i++ ) + _adj_b_V_unnorm.push_back( unnormAdjoint( _bp_dual.beliefV(i).p(), _bp_dual.beliefVZ(i), _adj_b_V[i] ) ); _adj_b_F_unnorm.clear(); - _adj_b_F_unnorm.reserve(_fg->nrFactors()); - for( size_t I = 0; I < _fg->nrFactors(); I++ ) { - _adj_b_F_unnorm.push_back( - unnormAdjoint(_bp_dual.beliefF(I).p(), - _bp_dual.beliefFZ(I), _adj_b_F[I])); - } + _adj_b_F_unnorm.reserve( _fg->nrFactors() ); + for( size_t I = 0; I < _fg->nrFactors(); I++ ) + _adj_b_F_unnorm.push_back( unnormAdjoint( _bp_dual.beliefF(I).p(), _bp_dual.beliefFZ(I), _adj_b_F[I] ) ); } void BBP::RegeneratePsiAdjoints() { _adj_psi_V.clear(); - _adj_psi_V.reserve(_fg->nrVars()); + _adj_psi_V.reserve( _fg->nrVars() ); for( size_t i = 0; i < _fg->nrVars(); i++ ) { - Prob p(_adj_b_V_unnorm[i]); - assert(p.size()==_fg->var(i).states()); - foreach(const Neighbor& I, _fg->nbV(i)) { - p *= _bp_dual.msgM(i,I.iter); - } + Prob p( _adj_b_V_unnorm[i] ); + assert( p.size() == _fg->var(i).states() ); + foreach( const Neighbor &I, _fg->nbV(i) ) + p *= _bp_dual.msgM( i, I.iter ); p += _init_adj_psi_V[i]; - _adj_psi_V.push_back(p); + _adj_psi_V.push_back( p ); } _adj_psi_F.clear(); - _adj_psi_F.reserve(_fg->nrFactors()); + _adj_psi_F.reserve( _fg->nrFactors() ); for( size_t I = 0; I < _fg->nrFactors(); I++ ) { - Prob p(_adj_b_F_unnorm[I]); - assert(p.size()==_fg->factor(I).states()); - foreach(const Neighbor& i, _fg->nbF(I)) { - Prob n_iI(_bp_dual.msgN(i,i.dual)); - const _ind_t& ind = _index(i,i.dual); + Prob p( _adj_b_F_unnorm[I] ); + assert( p.size() == _fg->factor(I).states() ); + foreach( const Neighbor &i, _fg->nbF(I) ) { + Prob n_iI( _bp_dual.msgN( i, i.dual ) ); + const _ind_t& ind = _index( i, i.dual ); // multiply prod with n_jI - for(size_t x_I = 0; x_I < p.size(); x_I++) + for( size_t x_I = 0; x_I < p.size(); x_I++ ) p[x_I] *= n_iI[ind[x_I]]; } p += _init_adj_psi_F[I]; - _adj_psi_F.push_back(p); + _adj_psi_F.push_back( p ); } } void BBP::RegenerateParMessageAdjoints() { size_t nv = _fg->nrVars(); - _adj_n.resize(nv); - _adj_m.resize(nv); - _adj_n_unnorm.resize(nv); - _adj_m_unnorm.resize(nv); - _new_adj_n.clear(); - _new_adj_n.resize(nv); - _new_adj_m.clear(); - _new_adj_m.resize(nv); + _adj_n.resize( nv ); + _adj_m.resize( nv ); + _adj_n_unnorm.resize( nv ); + _adj_m_unnorm.resize( nv ); + _new_adj_n.resize( nv ); + _new_adj_m.resize( nv ); for( size_t i = 0; i < _fg->nrVars(); i++ ) { size_t n_i = _fg->nbV(i).size(); - _adj_n[i].resize(n_i); - _adj_m[i].resize(n_i); - _adj_n_unnorm[i].resize(n_i); - _adj_m_unnorm[i].resize(n_i); - _new_adj_n[i].resize(n_i); - _new_adj_m[i].resize(n_i); - foreach(const Neighbor& I, _fg->nbV(i)) { - // calculate adj_n - { - Prob prod(_fg->factor(I).p()); + _adj_n[i].resize( n_i ); + _adj_m[i].resize( n_i ); + _adj_n_unnorm[i].resize( n_i ); + _adj_m_unnorm[i].resize( n_i ); + _new_adj_n[i].resize( n_i ); + _new_adj_m[i].resize( n_i ); + foreach( const Neighbor &I, _fg->nbV(i) ) { + { // calculate adj_n + Prob prod( _fg->factor(I).p() ); prod *= _adj_b_F_unnorm[I]; - foreach(const Neighbor& j, _fg->nbF(I)) { - if(i != j) { - Prob n_jI(_bp_dual.msgN(j,j.dual)); - const _ind_t& ind = _index(j,j.dual); + foreach( const Neighbor &j, _fg->nbF(I) ) + if( i != j ) { + Prob n_jI( _bp_dual.msgN( j, j.dual ) ); + const _ind_t &ind = _index( j, j.dual ); // multiply prod with n_jI for( size_t x_I = 0; x_I < prod.size(); x_I++ ) prod[x_I] *= n_jI[ind[x_I]]; } - } - Prob marg(_fg->var(i).states(), 0.0); - const _ind_t &ind = _index(i,I.iter); + Prob marg( _fg->var(i).states(), 0.0 ); + const _ind_t &ind = _index( i, I.iter ); for( size_t r = 0; r < prod.size(); r++ ) marg[ind[r]] += prod[r]; _new_adj_n[i][I.iter] = marg; - upMsgN(i,I.iter); + upMsgN( i, I.iter ); } - // calculate adj_m - { - Prob prod(_adj_b_V_unnorm[i]); - assert(prod.size()==_fg->var(i).states()); - foreach(const Neighbor& J, _fg->nbV(i)) { - if(size_t(J) != size_t(I)) { + { // calculate adj_m + Prob prod( _adj_b_V_unnorm[i] ); + assert( prod.size() == _fg->var(i).states() ); + foreach( const Neighbor &J, _fg->nbV(i) ) + if( J.node != I.node ) prod *= _bp_dual.msgM(i,J.iter); - } - } _new_adj_m[i][I.iter] = prod; - upMsgM(i,I.iter); + upMsgM( i, I.iter ); } } } } +void BBP::RegenerateSeqMessageAdjoints() { + size_t nv = _fg->nrVars(); + _adj_m.resize( nv ); + _adj_m_unnorm.resize( nv ); + _new_adj_m.resize( nv ); + for( size_t i = 0; i < _fg->nrVars(); i++ ) { + size_t n_i = _fg->nbV(i).size(); + _adj_m[i].resize( n_i ); + _adj_m_unnorm[i].resize( n_i ); + _new_adj_m[i].resize( n_i ); + foreach( const Neighbor &I, _fg->nbV(i) ) { + // calculate adj_m + Prob prod( _adj_b_V_unnorm[i] ); + assert( prod.size() == _fg->var(i).states() ); + foreach( const Neighbor &J, _fg->nbV(i) ) + if( J.node != I.node ) + prod *= _bp_dual.msgM( i, J.iter ); + _adj_m[i][I.iter] = prod; + calcUnnormMsgM( i, I.iter ); + _new_adj_m[i][I.iter] = Prob( _fg->var(i).states(), 0.0 ); + } + } + for( size_t i = 0; i < _fg->nrVars(); i++ ) { + foreach( const Neighbor &I, _fg->nbV(i) ) { + // calculate adj_n + Prob prod( _fg->factor(I).p() ); + prod *= _adj_b_F_unnorm[I]; + foreach( const Neighbor &j, _fg->nbF(I) ) + if( i != j ) { + Prob n_jI( _bp_dual.msgN( j, j.dual) ); + const _ind_t& ind = _index( j, j.dual ); + // multiply prod with n_jI + for( size_t x_I = 0; x_I < prod.size(); x_I++ ) + prod[x_I] *= n_jI[ind[x_I]]; + } + Prob marg( _fg->var(i).states(), 0.0 ); + const _ind_t &ind = _index( i, I.iter ); + for( size_t r = 0; r < prod.size(); r++ ) + marg[ind[r]] += prod[r]; + sendSeqMsgN( i, I.iter,marg ); + } + } +} + + +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 ); + size_t I = _fg->nbV(i)[_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]; + ); + /* 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 ) { + const Neighbor &I = _fg->nbV(i)[_I]; + Prob p( U(I, I.dual) ); + const Prob &adj = _adj_m_unnorm[i][_I]; + const _ind_t &ind = _index(i,_I); + for( size_t x_I = 0; x_I < p.size(); x_I++ ) + p[x_I] *= adj[ind[x_I]]; + _adj_psi_F[I] += p; + /* 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]; +} + + +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::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) ) { + upMsgM( i, I.iter ); + upMsgN( i, I.iter ); + } +} + + void BBP::incrSeqMsgM( size_t i, size_t _I, const Prob &p ) { if( props.clean_updates ) _new_adj_m[i][_I] += p; @@ -561,7 +447,7 @@ void BBP::updateSeqMsgM( size_t i, size_t _I ) { void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) { _adj_m[i][_I] = p; - calcUnnormMsgM(i, _I); + calcUnnormMsgM( i, _I ); } @@ -569,7 +455,7 @@ 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); + _adj_psi_V[i] += f_unnorm * T( i, _I ); #if 0 if(f_unnorm.sumAbs() > pv_thresh) { DAI_DMSG("in sendSeqMsgN"); @@ -665,146 +551,6 @@ void BBP::sendSeqMsgM( size_t j, size_t _I ) { } -void BBP::RegenerateSeqMessageAdjoints() { - size_t nv = _fg->nrVars(); - _adj_m.resize(nv); - _adj_m_unnorm.resize(nv); - _new_adj_m.resize(nv); - for( size_t i = 0; i < _fg->nrVars(); i++ ) { - size_t n_i = _fg->nbV(i).size(); - _adj_m[i].resize(n_i); - _adj_m_unnorm[i].resize(n_i); - _new_adj_m[i].resize(n_i); - foreach(const Neighbor& I, _fg->nbV(i)) { - // calculate adj_m - Prob prod(_adj_b_V_unnorm[i]); - assert(prod.size()==_fg->var(i).states()); - foreach(const Neighbor& J, _fg->nbV(i)) { - if(size_t(J) != size_t(I)) { - prod *= _bp_dual.msgM(i,J.iter); - } - } - _adj_m[i][I.iter] = prod; - calcUnnormMsgM(i, I.iter); - _new_adj_m[i][I.iter] = Prob(_fg->var(i).states(),0.0); - } - } - for( size_t i = 0; i < _fg->nrVars(); i++ ) { - foreach(const Neighbor& I, _fg->nbV(i)) { - // calculate adj_n - Prob prod(_fg->factor(I).p()); - prod *= _adj_b_F_unnorm[I]; - foreach(const Neighbor& j, _fg->nbF(I)) { - if(i != j) { - Prob n_jI(_bp_dual.msgN(j,j.dual)); - const _ind_t& ind = _index(j,j.dual); - // multiply prod with n_jI - for( size_t x_I = 0; x_I < prod.size(); x_I++ ) - prod[x_I] *= n_jI[ind[x_I]]; - } - } - Prob marg(_fg->var(i).states(), 0.0); - const _ind_t &ind = _index(i,I.iter); - for( size_t r = 0; r < prod.size(); r++ ) - marg[ind[r]] += prod[r]; - sendSeqMsgN(i,I.iter,marg); - } - } -} - - -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]; - 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 ) { - 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]; - ); - /* 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 ) { - const Neighbor &I = _fg->nbV(i)[_I]; - Prob p( U(I, I.dual) ); - const Prob &adj = _adj_m_unnorm[i][_I]; - const _ind_t &ind = _index(i,_I); - for( size_t x_I = 0; x_I < p.size(); x_I++ ) - p[x_I] *= adj[ind[x_I]]; - _adj_psi_F[I] += p; - /* 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]; -} - - -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::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) ) { - upMsgM( i, I.iter ); - upMsgN( i, I.iter ); - } -} - - Real BBP::getUnMsgMag() { Real s = 0.0; size_t e = 0; @@ -903,183 +649,178 @@ Real BBP::getTotalMsgN() { } -void BBP::run() { - typedef BBP::Properties::UpdateType UT; - Real tol = props.tol; - UT updates = props.updates; - Real tic=toc(); - switch((size_t)updates) { - case UT::SEQ_MAX: { - size_t i, _I; - Real mag; - do { - _iters++; - getArgmaxMsgM(i,_I,mag); - sendSeqMsgM(i,_I); - } while(mag > tol && _iters < props.maxiter); - - if(_iters >= props.maxiter) { - cerr << "Warning: BBP didn't converge in " << _iters - << " iterations (greatest message magnitude = " << mag << ")" - << endl; - } - break; - } - case UT::SEQ_FIX: { - Real mag; - do { - _iters++; - mag = getTotalMsgM(); - if(mag < tol) break; - - for( size_t i = 0; i < _fg->nrVars(); i++ ) - foreach(const Neighbor &I, _fg->nbV(i)) - sendSeqMsgM(i, I.iter); - for( size_t i = 0; i < _fg->nrVars(); i++ ) - foreach(const Neighbor &I, _fg->nbV(i)) - updateSeqMsgM(i, I.iter); - } while(mag > tol && _iters < props.maxiter); - - if(_iters >= props.maxiter) { - cerr << "Warning: BBP didn't converge in " << _iters - << " iterations (greatest message magnitude = " << mag << ")" - << endl; - break; - } - break; - } - case UT::SEQ_BP_REV: - case UT::SEQ_BP_FWD: { - const BP *bp = static_cast(_ia); - vector > sentMessages = bp->getSentMessages(); - size_t totalMessages = sentMessages.size(); - if(totalMessages==0) { - cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl; - DAI_THROW(INTERNAL_ERROR); - } - if(updates==UT::SEQ_BP_FWD) { - reverse(sentMessages.begin(), sentMessages.end()); - } -// DAI_PV(sentMessages.size()); -// DAI_PV(_iters); -// DAI_PV(props.maxiter); - while(sentMessages.size()>0 && _iters < props.maxiter) { -// DAI_PV(sentMessages.size()); -// DAI_PV(_iters); - _iters++; - pair e = sentMessages.back(); - sentMessages.pop_back(); - size_t i = e.first, _I = e.second; - sendSeqMsgM(i,_I); - } - if(_iters >= props.maxiter) { - cerr << "Warning: BBP updates limited to " << props.maxiter - << " iterations, but using UpdateType " << updates - << " with " << totalMessages << " messages" - << endl; - } - break; - } - case UT::PAR: { - do { - _iters++; - doParUpdate(); - } while((_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter); - if(_iters==props.maxiter) { - Real s, new_s; - getMsgMags(s,new_s); - cerr << "Warning: BBP didn't converge in " << _iters - << " iterations (unnorm message magnitude = " << getUnMsgMag() - << ", norm message mags = " << s << " -> " << new_s - << ")" << endl; - } - break; - } - } - if(props.verbose >= 3) { - cerr << "BBP::run() took " << toc()-tic << " seconds " - << doneIters() << " iterations" << endl; - } +void BBP::Regenerate() { + RegenerateInds(); + RegenerateT(); + RegenerateU(); + RegenerateS(); + RegenerateR(); + RegenerateInputs(); + RegeneratePsiAdjoints(); + if( props.updates == Properties::UpdateType::PAR ) + RegenerateParMessageAdjoints(); + else + RegenerateSeqMessageAdjoints(); + _iters = 0; } -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; - } +std::vector BBP::getZeroAdjF( const FactorGraph &fg ) { + vector 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(), 0.0 ) ); + return adj_2; } -vector 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(); +std::vector BBP::getZeroAdjV( const FactorGraph &fg ) { + vector 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(), 0.0 ) ); + return adj_1; } -double numericBBPTest(const InfAlg& bp, const vector *state, const PropertySet& bbp_props, bbp_cfn_t cfn, double h) { +void BBP::run() { + typedef BBP::Properties::UpdateType UT; + Real &tol = props.tol; + UT &updates = props.updates; + + Real tic = toc(); + switch( (size_t)updates ) { + case UT::SEQ_MAX: { + size_t i, _I; + Real mag; + do { + _iters++; + getArgmaxMsgM( i, _I, mag ); + sendSeqMsgM( i, _I ); + } while( mag > tol && _iters < props.maxiter ); + + if( _iters >= props.maxiter ) + if( props.verbose >= 1 ) + cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl; + break; + } case UT::SEQ_FIX: { + Real mag; + do { + _iters++; + mag = getTotalMsgM(); + if( mag < tol ) + break; + + for( size_t i = 0; i < _fg->nrVars(); i++ ) + foreach( const Neighbor &I, _fg->nbV(i) ) + sendSeqMsgM( i, I.iter ); + for( size_t i = 0; i < _fg->nrVars(); i++ ) + foreach( const Neighbor &I, _fg->nbV(i) ) + updateSeqMsgM( i, I.iter ); + } while( mag > tol && _iters < props.maxiter ); + + if( _iters >= props.maxiter ) + if( props.verbose >= 1 ) + cerr << "Warning: BBP didn't converge in " << _iters << " iterations (greatest message magnitude = " << mag << ")" << endl; + break; + } case UT::SEQ_BP_REV: + case UT::SEQ_BP_FWD: { + const BP *bp = static_cast(_ia); + vector > sentMessages = bp->getSentMessages(); + size_t totalMessages = sentMessages.size(); + if( totalMessages == 0 ) { + cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl; + DAI_THROW(INTERNAL_ERROR); + } + if( updates==UT::SEQ_BP_FWD ) + reverse( sentMessages.begin(), sentMessages.end() ); +// DAI_PV(sentMessages.size()); +// DAI_PV(_iters); +// DAI_PV(props.maxiter); + while( sentMessages.size() > 0 && _iters < props.maxiter ) { +// DAI_PV(sentMessages.size()); +// DAI_PV(_iters); + _iters++; + pair e = sentMessages.back(); + sentMessages.pop_back(); + size_t i = e.first, _I = e.second; + sendSeqMsgM( i, _I ); + } + if( _iters >= props.maxiter ) + if( props.verbose >= 1 ) + cerr << "Warning: BBP updates limited to " << props.maxiter << " iterations, but using UpdateType " << updates << " with " << totalMessages << " messages" << endl; + break; + } case UT::PAR: { + do { + _iters++; + doParUpdate(); + } while( (_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter ); + if( _iters == props.maxiter ) { + Real s, new_s; + getMsgMags( s, new_s ); + if( props.verbose >= 1 ) + cerr << "Warning: BBP didn't converge in " << _iters << " iterations (unnorm message magnitude = " << getUnMsgMag() << ", norm message mags = " << s << " -> " << new_s << ")" << endl; + } + break; + } + } + if( props.verbose >= 3 ) + cerr << "BBP::run() took " << toc()-tic << " seconds " << doneIters() << " iterations" << endl; +} + + +double numericBBPTest( const InfAlg &bp, const vector *state, const PropertySet &bbp_props, bbp_cfn_t cfn, double h ) { // calculate the value of the unperturbed cost function - Real cf0 = getCostFn(bp, cfn, state); + Real cf0 = getCostFn( bp, cfn, state ); // run BBP to estimate adjoints - BBP bbp(&bp,bbp_props); - initBBPCostFnAdj(bbp, bp, cfn, state); + BBP bbp( &bp, bbp_props ); + initBBPCostFnAdj( bbp, bp, cfn, state ); bbp.run(); - Real d=0; + Real d = 0; const FactorGraph& fg = bp.fg(); - if(1) { + if( 1 ) { // verify bbp.adj_psi_V // for each variable i - for(size_t i=0; i adj_est; // for each value xi - for(size_t xi=0; xifg().var(i).states(); - Prob psi_1_prb(n,1.0); + Prob psi_1_prb( n, 1.0 ); psi_1_prb[xi] += h; // psi_1_prb.normalize(); size_t I = bp_prb->fg().nbV(i)[0]; // use first factor in list of neighbors of i - bp_prb->fg().factor(I) *= Factor(bp_prb->fg().var(i), psi_1_prb); + bp_prb->fg().factor(I) *= Factor( bp_prb->fg().var(i), psi_1_prb ); // call 'init' on the perturbed variables - bp_prb->init(bp_prb->fg().var(i)); + bp_prb->init( bp_prb->fg().var(i) ); // run copy to convergence bp_prb->run(); // calculate new value of cost function - Real cf_prb = getCostFn(*bp_prb, cfn, state); + Real cf_prb = getCostFn( *bp_prb, cfn, state ); // use to estimate adjoint for i - adj_est.push_back((cf_prb-cf0)/h); + adj_est.push_back( (cf_prb - cf0) / h ); // free cloned InfAlg delete bp_prb; } - Prob p_adj_est(adj_est.begin(), adj_est.end()); + Prob p_adj_est( adj_est.begin(), adj_est.end() ); // compare this numerical estimate to the BBP estimate; sum the distances - cerr << "i: " << i + cout << "i: " << i << ", p_adj_est: " << p_adj_est << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl; - d += dist(p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1); + d += dist( p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1 ); } } /* if(1) { @@ -1173,6 +914,243 @@ double numericBBPTest(const InfAlg& bp, const vector *state, const Prope } +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 *stateP ) { + const FactorGraph &fg = ia.fg(); + + switch( (size_t)cfn_type ) { + case bbp_cfn_t::CFN_BETHE_ENT: { + vector b1_adj; + vector b2_adj; + vector psi1_adj; + vector psi2_adj; + b1_adj.reserve( fg.nrVars() ); + psi1_adj.reserve( fg.nrVars() ); + b2_adj.reserve( fg.nrFactors() ); + psi2_adj.reserve( fg.nrFactors() ); + for( size_t i = 0; i < 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 b2_adj; + b2_adj.reserve( fg.nrFactors() ); + for( size_t I = 0; I < fg.nrFactors(); I++ ) { + size_t dim = fg.factor(I).states(); + Prob p( dim, 0.0 ); + for( size_t xI = 0; xI < dim; xI++ ) { + double bIxI = 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 b1_adj; + b1_adj.reserve( fg.nrVars() ); + for( size_t i = 0; i < fg.nrVars(); i++ ) { + size_t dim = fg.var(i).states(); + Prob p( dim, 0.0 ); + for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) { + double bixi = 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 state; + if( stateP == NULL ) + state = getGibbsState( ia, 2*ia.Iterations() ); + else + state = *stateP; + assert( state.size() == fg.nrVars() ); + + vector b1_adj; + b1_adj.reserve(fg.nrVars()); + for( size_t i = 0; i < state.size(); i++ ) { + size_t n = fg.var(i).states(); + Prob delta( n, 0.0 ); + assert(/*0<=state[i] &&*/ state[i] < n); + double b = 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(INTERNAL_ERROR); + } + 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 state; + if( stateP == NULL ) + state = getGibbsState( ia, 2*ia.Iterations() ); + else + state = *stateP; + assert( state.size() == fg.nrVars() ); + + vector b2_adj; + b2_adj.reserve( fg.nrVars() ); + for( size_t I = 0; I < fg.nrFactors(); I++ ) { + size_t n = fg.factor(I).states(); + Prob delta( n, 0.0 ); + + size_t x_I = getFactorEntryForState( fg, I, state ); + assert(/*0<=x_I &&*/ x_I < n); + + double b = ia.beliefF(I)[x_I]; + switch( (size_t)cfn_type ) { + case bbp_cfn_t::CFN_GIBBS_B_FACTOR: + delta[x_I] = 1.0; + break; + case bbp_cfn_t::CFN_GIBBS_B2_FACTOR: + delta[x_I] = b; + break; + case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: + delta[x_I] = exp( b ); + break; + default: + DAI_THROW(INTERNAL_ERROR); + } + 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 *stateP ) { + double 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: { + assert( stateP != NULL ); + vector state = *stateP; + assert( state.size() == fg.nrVars() ); + for( size_t i = 0; i < fg.nrVars(); i++ ) { + double b = ia.beliefV(i)[state[i]]; + switch( (size_t)cfn_type ) { + case bbp_cfn_t::CFN_GIBBS_B: + cf += b; + break; + case bbp_cfn_t::CFN_GIBBS_B2: + cf += b * b / 2.0; + break; + case bbp_cfn_t::CFN_GIBBS_EXP: + cf += exp( b ); + break; + default: + DAI_THROW(INTERNAL_ERROR); + } + } + break; + } case bbp_cfn_t::CFN_GIBBS_B_FACTOR: + case bbp_cfn_t::CFN_GIBBS_B2_FACTOR: + case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: { + assert( stateP != NULL ); + vector state = *stateP; + assert( state.size() == fg.nrVars() ); + for( size_t I = 0; I < fg.nrFactors(); I++ ) { + size_t x_I = getFactorEntryForState( fg, I, state ); + double b = ia.beliefF(I)[x_I]; + switch( (size_t)cfn_type ) { + case bbp_cfn_t::CFN_GIBBS_B_FACTOR: + cf += b; + break; + case bbp_cfn_t::CFN_GIBBS_B2_FACTOR: + cf += b * b / 2.0; + break; + case bbp_cfn_t::CFN_GIBBS_EXP_FACTOR: + cf += exp( b ); + break; + default: + DAI_THROW(INTERNAL_ERROR); + } + } + break; + } default: + cerr << "Unknown cost function " << cfn_type << endl; + DAI_THROW(UNKNOWN_ENUM_VALUE); + } + return cf; +} + + } // end of namespace dai diff --git a/tests/aliases.conf b/tests/aliases.conf index dfa4802..dc9a71b 100644 --- a/tests/aliases.conf +++ b/tests/aliases.conf @@ -104,5 +104,5 @@ GIBBS_1e9: GIBBS[iters=1000000000] # --- 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=10000,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] diff --git a/tests/testbbp.cpp b/tests/testbbp.cpp index 502d869..8c061cb 100644 --- a/tests/testbbp.cpp +++ b/tests/testbbp.cpp @@ -85,31 +85,31 @@ int main( int argc, char *argv[] ) { bbp_cfn_t cfn; switch( (t / 10) % 9 ) { case 0: - cfn = bbp_cfn_t::cfn_gibbs_b; + cfn = bbp_cfn_t::CFN_GIBBS_B; break; case 1: - cfn = bbp_cfn_t::cfn_gibbs_b2; + cfn = bbp_cfn_t::CFN_GIBBS_B2; break; case 2: - cfn = bbp_cfn_t::cfn_gibbs_exp; + cfn = bbp_cfn_t::CFN_GIBBS_EXP; break; case 3: - cfn = bbp_cfn_t::cfn_gibbs_b_factor; + cfn = bbp_cfn_t::CFN_GIBBS_B_FACTOR; break; case 4: - cfn = bbp_cfn_t::cfn_gibbs_b2_factor; + cfn = bbp_cfn_t::CFN_GIBBS_B2_FACTOR; break; case 5: - cfn = bbp_cfn_t::cfn_gibbs_exp_factor; + cfn = bbp_cfn_t::CFN_GIBBS_EXP_FACTOR; break; case 6: - cfn = bbp_cfn_t::cfn_var_ent; + cfn = bbp_cfn_t::CFN_VAR_ENT; break; case 7: - cfn = bbp_cfn_t::cfn_factor_ent; + cfn = bbp_cfn_t::CFN_FACTOR_ENT; break; case 8: - cfn = bbp_cfn_t::cfn_bethe_ent; + cfn = bbp_cfn_t::CFN_BETHE_ENT; break; } -- 2.20.1