}
-JTree::JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic ) : DAIAlgRG(fg), _RTree(), _Qa(), _Qb(), _mes(), _logZ(), props() {
+JTree::JTree( const FactorGraph &fg, const PropertySet &opts, bool automatic ) : DAIAlgRG(fg), _mes(), _logZ(), RTree(), Qa(), Qb(), props() {
setProperties( opts );
if( !isConnected() )
}
// Construct maximal spanning tree using Prim's algorithm
- _RTree = MaxSpanningTreePrims( JuncGraph );
+ RTree = MaxSpanningTreePrims( JuncGraph );
// Construct corresponding region graph
RecomputeORs();
// Create inner regions and edges
- IRs.reserve( _RTree.size() );
+ IRs.reserve( RTree.size() );
vector<Edge> edges;
- edges.reserve( 2 * _RTree.size() );
- for( size_t i = 0; i < _RTree.size(); i++ ) {
- edges.push_back( Edge( _RTree[i].n1, nrIRs() ) );
- edges.push_back( Edge( _RTree[i].n2, nrIRs() ) );
+ edges.reserve( 2 * RTree.size() );
+ for( size_t i = 0; i < RTree.size(); i++ ) {
+ edges.push_back( Edge( RTree[i].n1, nrIRs() ) );
+ edges.push_back( Edge( RTree[i].n2, nrIRs() ) );
// inner clusters have counting number -1
- IRs.push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+ IRs.push_back( Region( Cliques[RTree[i].n1] & Cliques[RTree[i].n2], -1.0 ) );
}
// create bipartite graph
G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
// Create messages and beliefs
- _Qa.clear();
- _Qa.reserve( nrORs() );
+ Qa.clear();
+ Qa.reserve( nrORs() );
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- _Qa.push_back( OR(alpha) );
+ Qa.push_back( OR(alpha) );
- _Qb.clear();
- _Qb.reserve( nrIRs() );
+ Qb.clear();
+ Qb.reserve( nrIRs() );
for( size_t beta = 0; beta < nrIRs(); beta++ )
- _Qb.push_back( Factor( IR(beta), 1.0 ) );
+ Qb.push_back( Factor( IR(beta), 1.0 ) );
_mes.clear();
_mes.reserve( nrORs() );
Factor JTree::belief( const VarSet &ns ) const {
vector<Factor>::const_iterator beta;
- for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ for( beta = Qb.begin(); beta != Qb.end(); beta++ )
if( beta->vars() >> ns )
break;
- if( beta != _Qb.end() )
+ if( beta != Qb.end() )
return( beta->marginal(ns) );
else {
vector<Factor>::const_iterator alpha;
- for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ for( alpha = Qa.begin(); alpha != Qa.end(); alpha++ )
if( alpha->vars() >> ns )
break;
- assert( alpha != _Qa.end() );
+ assert( alpha != Qa.end() );
return( alpha->marginal(ns) );
}
}
vector<Factor> JTree::beliefs() const {
vector<Factor> result;
for( size_t beta = 0; beta < nrIRs(); beta++ )
- result.push_back( _Qb[beta] );
+ result.push_back( Qb[beta] );
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- result.push_back( _Qa[alpha] );
+ result.push_back( Qa[alpha] );
return result;
}
// Needs no init
void JTree::runHUGIN() {
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- _Qa[alpha] = OR(alpha);
+ Qa[alpha] = OR(alpha);
for( size_t beta = 0; beta < nrIRs(); beta++ )
- _Qb[beta].fill( 1.0 );
+ Qb[beta].fill( 1.0 );
// CollectEvidence
_logZ = 0.0;
- for( size_t i = _RTree.size(); (i--) != 0; ) {
-// Make outer region _RTree[i].n1 consistent with outer region _RTree[i].n2
-// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
- Factor new_Qb = _Qa[_RTree[i].n2].partSum( IR( i ) );
+ for( size_t i = RTree.size(); (i--) != 0; ) {
+// Make outer region RTree[i].n1 consistent with outer region RTree[i].n2
+// IR(i) = seperator OR(RTree[i].n1) && OR(RTree[i].n2)
+ Factor new_Qb = Qa[RTree[i].n2].partSum( IR( i ) );
_logZ += log(new_Qb.normalize());
- _Qa[_RTree[i].n1] *= new_Qb.divided_by( _Qb[i] );
- _Qb[i] = new_Qb;
+ Qa[RTree[i].n1] *= new_Qb.divided_by( Qb[i] );
+ Qb[i] = new_Qb;
}
- if( _RTree.empty() )
- _logZ += log(_Qa[0].normalize() );
+ if( RTree.empty() )
+ _logZ += log(Qa[0].normalize() );
else
- _logZ += log(_Qa[_RTree[0].n1].normalize());
+ _logZ += log(Qa[RTree[0].n1].normalize());
// DistributeEvidence
- for( size_t i = 0; i < _RTree.size(); i++ ) {
-// Make outer region _RTree[i].n2 consistent with outer region _RTree[i].n1
-// IR(i) = seperator OR(_RTree[i].n1) && OR(_RTree[i].n2)
- Factor new_Qb = _Qa[_RTree[i].n1].marginal( IR( i ) );
- _Qa[_RTree[i].n2] *= new_Qb.divided_by( _Qb[i] );
- _Qb[i] = new_Qb;
+ for( size_t i = 0; i < RTree.size(); i++ ) {
+// Make outer region RTree[i].n2 consistent with outer region RTree[i].n1
+// IR(i) = seperator OR(RTree[i].n1) && OR(RTree[i].n2)
+ Factor new_Qb = Qa[RTree[i].n1].marginal( IR( i ) );
+ Qa[RTree[i].n2] *= new_Qb.divided_by( Qb[i] );
+ Qb[i] = new_Qb;
}
// Normalize
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- _Qa[alpha].normalize();
+ Qa[alpha].normalize();
}
// First pass
_logZ = 0.0;
for( size_t e = nrIRs(); (e--) != 0; ) {
- // send a message from _RTree[e].n2 to _RTree[e].n1
- // or, actually, from the seperator IR(e) to _RTree[e].n1
+ // send a message from RTree[e].n2 to RTree[e].n1
+ // or, actually, from the seperator IR(e) to RTree[e].n1
- size_t i = nbIR(e)[1].node; // = _RTree[e].n2
- size_t j = nbIR(e)[0].node; // = _RTree[e].n1
+ size_t i = nbIR(e)[1].node; // = RTree[e].n2
+ size_t j = nbIR(e)[0].node; // = RTree[e].n1
size_t _e = nbIR(e)[0].dual;
Factor piet = OR(i);
// Second pass
for( size_t e = 0; e < nrIRs(); e++ ) {
- size_t i = nbIR(e)[0].node; // = _RTree[e].n1
- size_t j = nbIR(e)[1].node; // = _RTree[e].n2
+ size_t i = nbIR(e)[0].node; // = RTree[e].n1
+ size_t j = nbIR(e)[1].node; // = RTree[e].n2
size_t _e = nbIR(e)[1].dual;
Factor piet = OR(i);
piet *= message( alpha, k.iter );
if( nrIRs() == 0 ) {
_logZ += log( piet.normalize() );
- _Qa[alpha] = piet;
- } else if( alpha == nbIR(0)[0].node /*_RTree[0].n1*/ ) {
+ Qa[alpha] = piet;
+ } else if( alpha == nbIR(0)[0].node /*RTree[0].n1*/ ) {
_logZ += log( piet.normalize() );
- _Qa[alpha] = piet;
+ Qa[alpha] = piet;
} else
- _Qa[alpha] = piet.normalized();
+ Qa[alpha] = piet.normalized();
}
// Only for logZ (and for belief)...
for( size_t beta = 0; beta < nrIRs(); beta++ )
- _Qb[beta] = _Qa[nbIR(beta)[0].node].marginal( IR(beta) );
+ Qb[beta] = Qa[nbIR(beta)[0].node].marginal( IR(beta) );
}
Real JTree::logZ() const {
Real sum = 0.0;
for( size_t beta = 0; beta < nrIRs(); beta++ )
- sum += IR(beta).c() * _Qb[beta].entropy();
+ sum += IR(beta).c() * Qb[beta].entropy();
for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
- sum += OR(alpha).c() * _Qa[alpha].entropy();
- sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
+ sum += OR(alpha).c() * Qa[alpha].entropy();
+ sum += (OR(alpha).log0() * Qa[alpha]).totalSum();
}
return sum;
}
// grow new tree
Graph oldTree;
- for( DEdgeVec::const_iterator e = _RTree.begin(); e != _RTree.end(); e++ )
+ for( DEdgeVec::const_iterator e = RTree.begin(); e != RTree.end(); e++ )
oldTree.insert( UEdge(e->n1, e->n2) );
DEdgeVec newTree = GrowRootedTree( oldTree, maxalpha );
// assumes that run() has been called already
Factor JTree::calcMarginal( const VarSet& ns ) {
vector<Factor>::const_iterator beta;
- for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
+ for( beta = Qb.begin(); beta != Qb.end(); beta++ )
if( beta->vars() >> ns )
break;
- if( beta != _Qb.end() )
+ if( beta != Qb.end() )
return( beta->marginal(ns) );
else {
vector<Factor>::const_iterator alpha;
- for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
+ for( alpha = Qa.begin(); alpha != Qa.end(); alpha++ )
if( alpha->vars() >> ns )
break;
- if( alpha != _Qa.end() )
+ if( alpha != Qa.end() )
return( alpha->marginal(ns) );
else {
// Find subtree to do efficient inference
VarSet nsrem = ns / OR(T.front().n1).vars();
Factor Pns (ns, 0.0);
- // Save _Qa and _Qb on the subtree
- map<size_t,Factor> _Qa_old;
- map<size_t,Factor> _Qb_old;
+ // Save Qa and Qb on the subtree
+ map<size_t,Factor> Qa_old;
+ map<size_t,Factor> Qb_old;
vector<size_t> b(Tsize, 0);
for( size_t i = Tsize; (i--) != 0; ) {
size_t alpha1 = T[i].n1;
size_t alpha2 = T[i].n2;
size_t beta;
for( beta = 0; beta < nrIRs(); beta++ )
- if( UEdge( _RTree[beta].n1, _RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
+ if( UEdge( RTree[beta].n1, RTree[beta].n2 ) == UEdge( alpha1, alpha2 ) )
break;
assert( beta != nrIRs() );
b[i] = beta;
- if( !_Qa_old.count( alpha1 ) )
- _Qa_old[alpha1] = _Qa[alpha1];
- if( !_Qa_old.count( alpha2 ) )
- _Qa_old[alpha2] = _Qa[alpha2];
- if( !_Qb_old.count( beta ) )
- _Qb_old[beta] = _Qb[beta];
+ if( !Qa_old.count( alpha1 ) )
+ Qa_old[alpha1] = Qa[alpha1];
+ if( !Qa_old.count( alpha2 ) )
+ Qa_old[alpha2] = Qa[alpha2];
+ if( !Qb_old.count( beta ) )
+ Qb_old[beta] = Qb[beta];
}
// For all states of nsrem
// IR(i) = seperator OR(T[i].n1) && OR(T[i].n2)
for( VarSet::const_iterator n = nsrem.begin(); n != nsrem.end(); n++ )
- if( _Qa[T[i].n2].vars() >> *n ) {
+ if( Qa[T[i].n2].vars() >> *n ) {
Factor piet( *n, 0.0 );
piet[s(*n)] = 1.0;
- _Qa[T[i].n2] *= piet;
+ Qa[T[i].n2] *= piet;
}
- Factor new_Qb = _Qa[T[i].n2].partSum( IR( b[i] ) );
+ Factor new_Qb = Qa[T[i].n2].partSum( IR( b[i] ) );
logZ += log(new_Qb.normalize());
- _Qa[T[i].n1] *= new_Qb.divided_by( _Qb[b[i]] );
- _Qb[b[i]] = new_Qb;
+ Qa[T[i].n1] *= new_Qb.divided_by( Qb[b[i]] );
+ Qb[b[i]] = new_Qb;
}
- logZ += log(_Qa[T[0].n1].normalize());
+ logZ += log(Qa[T[0].n1].normalize());
Factor piet( nsrem, 0.0 );
piet[s] = exp(logZ);
- Pns += piet * _Qa[T[0].n1].partSum( ns / nsrem ); // OPTIMIZE ME
+ Pns += piet * Qa[T[0].n1].partSum( ns / nsrem ); // OPTIMIZE ME
// Restore clamped beliefs
- for( map<size_t,Factor>::const_iterator alpha = _Qa_old.begin(); alpha != _Qa_old.end(); alpha++ )
- _Qa[alpha->first] = alpha->second;
- for( map<size_t,Factor>::const_iterator beta = _Qb_old.begin(); beta != _Qb_old.end(); beta++ )
- _Qb[beta->first] = beta->second;
+ for( map<size_t,Factor>::const_iterator alpha = Qa_old.begin(); alpha != Qa_old.end(); alpha++ )
+ Qa[alpha->first] = alpha->second;
+ for( map<size_t,Factor>::const_iterator beta = Qb_old.begin(); beta != Qb_old.end(); beta++ )
+ Qb[beta->first] = beta->second;
}
return( Pns.normalized() );
}
// Construct maximal spanning tree using Prim's algorithm
- _RTree = MaxSpanningTreePrims( JuncGraph );
+ RTree = MaxSpanningTreePrims( JuncGraph );
// Construct corresponding region graph
RecomputeORs();
// Create inner regions and edges
- IRs.reserve( _RTree.size() );
+ IRs.reserve( RTree.size() );
vector<Edge> edges;
- edges.reserve( 2 * _RTree.size() );
- for( size_t i = 0; i < _RTree.size(); i++ ) {
- edges.push_back( Edge( _RTree[i].n1, IRs.size() ) );
- edges.push_back( Edge( _RTree[i].n2, IRs.size() ) );
+ edges.reserve( 2 * RTree.size() );
+ for( size_t i = 0; i < RTree.size(); i++ ) {
+ edges.push_back( Edge( RTree[i].n1, IRs.size() ) );
+ edges.push_back( Edge( RTree[i].n2, IRs.size() ) );
// inner clusters have counting number -1
- IRs.push_back( Region( Cliques[_RTree[i].n1] & Cliques[_RTree[i].n2], -1.0 ) );
+ IRs.push_back( Region( Cliques[RTree[i].n1] & Cliques[RTree[i].n2], -1.0 ) );
}
// create bipartite graph
Check_Counting_Numbers();
// Create messages and beliefs
- _Qa.clear();
- _Qa.reserve( nrORs() );
+ Qa.clear();
+ Qa.reserve( nrORs() );
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- _Qa.push_back( OR(alpha) );
+ Qa.push_back( OR(alpha) );
- _Qb.clear();
- _Qb.reserve( nrIRs() );
+ Qb.clear();
+ Qb.reserve( nrIRs() );
for( size_t beta = 0; beta < nrIRs(); beta++ )
- _Qb.push_back( Factor( IR(beta), 1.0 ) );
+ Qb.push_back( Factor( IR(beta), 1.0 ) );
// DIFF with JTree::GenerateJT: no messages
//subTree.resize( subTreeSize ); // FIXME
// cout << "subtree " << I << " has size " << subTreeSize << endl;
- TreeEPSubTree QI( subTree, _RTree, _Qa, _Qb, &factor(I) );
+ TreeEPSubTree QI( subTree, RTree, Qa, Qb, &factor(I) );
_Q[I] = QI;
}
// Previous root of first off-tree factor should be the root of the last off-tree factor
//subTree.resize( subTreeSize ); // FIXME
// cout << "subtree " << I << " has size " << subTreeSize << endl;
- TreeEPSubTree QI( subTree, _RTree, _Qa, _Qb, &factor(I) );
+ TreeEPSubTree QI( subTree, RTree, Qa, Qb, &factor(I) );
_Q[I] = QI;
break;
}
for( _iters=0; _iters < props.maxiter && diffs.maxDiff() > props.tol; _iters++ ) {
for( size_t I = 0; I < nrFactors(); I++ )
if( offtree(I) ) {
- _Q[I].InvertAndMultiply( _Qa, _Qb );
- _Q[I].HUGIN_with_I( _Qa, _Qb );
- _Q[I].InvertAndMultiply( _Qa, _Qb );
+ _Q[I].InvertAndMultiply( Qa, Qb );
+ _Q[I].HUGIN_with_I( Qa, Qb );
+ _Q[I].InvertAndMultiply( Qa, Qb );
}
// calculate new beliefs and compare with old ones
// entropy of the tree
for( size_t beta = 0; beta < nrIRs(); beta++ )
- sum -= _Qb[beta].entropy();
+ sum -= Qb[beta].entropy();
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- sum += _Qa[alpha].entropy();
+ sum += Qa[alpha].entropy();
// energy of the on-tree factors
for( size_t alpha = 0; alpha < nrORs(); alpha++ )
- sum += (OR(alpha).log0() * _Qa[alpha]).totalSum();
+ sum += (OR(alpha).log0() * Qa[alpha]).totalSum();
// energy of the off-tree factors
for( size_t I = 0; I < nrFactors(); I++ )
if( offtree(I) )
- sum += (_Q.find(I))->second.logZ( _Qa, _Qb );
+ sum += (_Q.find(I))->second.logZ( Qa, Qb );
return sum;
}