1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
4 This file is part of libDAI.
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 #include <dai/jtree.h>
26 #include <dai/treeep.h>
28 #include <dai/diffs.h>
37 const char *TreeEP::Name
= "TREEEP";
40 bool TreeEP::checkProperties() {
41 if( !HasProperty("type") )
43 if( !HasProperty("tol") )
45 if (!HasProperty("maxiter") )
47 if (!HasProperty("verbose") )
50 ConvertPropertyTo
<TypeType
>("type");
51 ConvertPropertyTo
<double>("tol");
52 ConvertPropertyTo
<size_t>("maxiter");
53 ConvertPropertyTo
<size_t>("verbose");
59 TreeEPSubTree::TreeEPSubTree( const DEdgeVec
&subRTree
, const DEdgeVec
&jt_RTree
, const std::vector
<Factor
> &jt_Qa
, const std::vector
<Factor
> &jt_Qb
, const Factor
*I
) : _Qa(), _Qb(), _RTree(), _a(), _b(), _I(I
), _ns(), _nsrem(), _logZ(0.0) {
62 // Make _Qa, _Qb, _a and _b corresponding to the subtree
63 _b
.reserve( subRTree
.size() );
64 _Qb
.reserve( subRTree
.size() );
65 _RTree
.reserve( subRTree
.size() );
66 for( size_t i
= 0; i
< subRTree
.size(); i
++ ) {
67 size_t alpha1
= subRTree
[i
].n1
; // old index 1
68 size_t alpha2
= subRTree
[i
].n2
; // old index 2
69 size_t beta
; // old sep index
70 for( beta
= 0; beta
< jt_RTree
.size(); beta
++ )
71 if( UEdge( jt_RTree
[beta
].n1
, jt_RTree
[beta
].n2
) == UEdge( alpha1
, alpha2
) )
73 assert( beta
!= jt_RTree
.size() );
75 size_t newalpha1
= find(_a
.begin(), _a
.end(), alpha1
) - _a
.begin();
76 if( newalpha1
== _a
.size() ) {
77 _Qa
.push_back( Factor( jt_Qa
[alpha1
].vars(), 1.0 ) );
78 _a
.push_back( alpha1
); // save old index in index conversion table
81 size_t newalpha2
= find(_a
.begin(), _a
.end(), alpha2
) - _a
.begin();
82 if( newalpha2
== _a
.size() ) {
83 _Qa
.push_back( Factor( jt_Qa
[alpha2
].vars(), 1.0 ) );
84 _a
.push_back( alpha2
); // save old index in index conversion table
87 _RTree
.push_back( DEdge( newalpha1
, newalpha2
) );
88 _Qb
.push_back( Factor( jt_Qb
[beta
].vars(), 1.0 ) );
92 // Find remaining variables (which are not in the new root)
93 _nsrem
= _ns
/ _Qa
[0].vars();
97 void TreeEPSubTree::init() {
98 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
99 _Qa
[alpha
].fill( 1.0 );
100 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
101 _Qb
[beta
].fill( 1.0 );
105 void TreeEPSubTree::InvertAndMultiply( const std::vector
<Factor
> &Qa
, const std::vector
<Factor
> &Qb
) {
106 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
107 _Qa
[alpha
] = Qa
[_a
[alpha
]].divided_by( _Qa
[alpha
] );
109 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
110 _Qb
[beta
] = Qb
[_b
[beta
]].divided_by( _Qb
[beta
] );
114 void TreeEPSubTree::HUGIN_with_I( std::vector
<Factor
> &Qa
, std::vector
<Factor
> &Qb
) {
115 // Backup _Qa and _Qb
116 vector
<Factor
> _Qa_old(_Qa
);
117 vector
<Factor
> _Qb_old(_Qb
);
120 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
121 Qa
[_a
[alpha
]].fill( 0.0 );
122 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
123 Qb
[_b
[beta
]].fill( 0.0 );
125 // For all states of _nsrem
126 for( State
s(_nsrem
); s
.valid(); s
++ ) {
127 // Multiply root with slice of I
128 _Qa
[0] *= _I
->slice( _nsrem
, s
);
131 for( size_t i
= _RTree
.size(); (i
--) != 0; ) {
132 // clamp variables in nsrem
133 for( VarSet::const_iterator n
= _nsrem
.begin(); n
!= _nsrem
.end(); n
++ )
134 if( _Qa
[_RTree
[i
].n2
].vars() >> *n
) {
135 Factor
delta( *n
, 0.0 );
137 _Qa
[_RTree
[i
].n2
] *= delta
;
139 Factor new_Qb
= _Qa
[_RTree
[i
].n2
].part_sum( _Qb
[i
].vars() );
140 _Qa
[_RTree
[i
].n1
] *= new_Qb
.divided_by( _Qb
[i
] );
144 // DistributeEvidence
145 for( size_t i
= 0; i
< _RTree
.size(); i
++ ) {
146 Factor new_Qb
= _Qa
[_RTree
[i
].n1
].part_sum( _Qb
[i
].vars() );
147 _Qa
[_RTree
[i
].n2
] *= new_Qb
.divided_by( _Qb
[i
] );
151 // Store Qa's and Qb's
152 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
153 Qa
[_a
[alpha
]].p() += _Qa
[alpha
].p();
154 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
155 Qb
[_b
[beta
]].p() += _Qb
[beta
].p();
157 // Restore _Qa and _Qb
162 // Normalize Qa and Qb
164 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ ) {
165 _logZ
+= log(Qa
[_a
[alpha
]].totalSum());
166 Qa
[_a
[alpha
]].normalize( Prob::NORMPROB
);
168 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ ) {
169 _logZ
-= log(Qb
[_b
[beta
]].totalSum());
170 Qb
[_b
[beta
]].normalize( Prob::NORMPROB
);
175 double TreeEPSubTree::logZ( const std::vector
<Factor
> &Qa
, const std::vector
<Factor
> &Qb
) const {
177 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
178 sum
+= (Qa
[_a
[alpha
]] * _Qa
[alpha
].log0()).totalSum();
179 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
180 sum
-= (Qb
[_b
[beta
]] * _Qb
[beta
].log0()).totalSum();
185 TreeEP::TreeEP( const FactorGraph
&fg
, const Properties
&opts
) : JTree(fg
, opts("updates",string("HUGIN")), false) {
186 assert( checkProperties() );
188 assert( fg
.G
.isConnected() );
190 if( opts
.hasKey("tree") ) {
191 ConstructRG( opts
.GetAs
<DEdgeVec
>("tree") );
193 if( Type() == TypeType::ORG
) {
194 // construct weighted graph with as weights a crude estimate of the
195 // mutual information between the nodes
196 WeightedGraph
<double> wg
;
197 for( size_t i
= 0; i
< nrVars(); ++i
) {
199 VarSet di
= delta(i
);
200 for( VarSet::const_iterator j
= di
.begin(); j
!= di
.end(); j
++ )
203 for( size_t I
= 0; I
< nrFactors(); I
++ ) {
204 VarSet Ivars
= factor(I
).vars();
205 if( (Ivars
== v_i
) || (Ivars
== *j
) )
207 else if( Ivars
>> (v_i
| *j
) )
208 piet
*= factor(I
).marginal( v_i
| *j
);
210 if( piet
.vars() >> (v_i
| *j
) ) {
211 piet
= piet
.marginal( v_i
| *j
);
212 Factor pietf
= piet
.marginal(v_i
) * piet
.marginal(*j
);
213 wg
[UEdge(i
,findVar(*j
))] = real( KL_dist( piet
, pietf
) );
215 wg
[UEdge(i
,findVar(*j
))] = 0;
219 // find maximal spanning tree
220 ConstructRG( MaxSpanningTreePrims( wg
) );
222 // cout << "Constructing maximum spanning tree..." << endl;
223 // DEdgeVec MST = MaxSpanningTreePrims( wg );
224 // cout << "Maximum spanning tree:" << endl;
225 // for( DEdgeVec::const_iterator e = MST.begin(); e != MST.end(); e++ )
226 // cout << *e << endl;
227 // ConstructRG( MST );
228 } else if( Type() == TypeType::ALT
) {
229 // construct weighted graph with as weights an upper bound on the
230 // effective interaction strength between pairs of nodes
231 WeightedGraph
<double> wg
;
232 for( size_t i
= 0; i
< nrVars(); ++i
) {
234 VarSet di
= delta(i
);
235 for( VarSet::const_iterator j
= di
.begin(); j
!= di
.end(); j
++ )
238 for( size_t I
= 0; I
< nrFactors(); I
++ ) {
239 VarSet Ivars
= factor(I
).vars();
240 if( Ivars
>> (v_i
| *j
) )
243 wg
[UEdge(i
,findVar(*j
))] = piet
.strength(v_i
, *j
);
247 // find maximal spanning tree
248 ConstructRG( MaxSpanningTreePrims( wg
) );
256 void TreeEP::ConstructRG( const DEdgeVec
&tree
) {
257 vector
<VarSet
> Cliques
;
258 for( size_t i
= 0; i
< tree
.size(); i
++ )
259 Cliques
.push_back( var(tree
[i
].n1
) | var(tree
[i
].n2
) );
261 // Construct a weighted graph (each edge is weighted with the cardinality
262 // of the intersection of the nodes, where the nodes are the elements of
264 WeightedGraph
<int> JuncGraph
;
265 for( size_t i
= 0; i
< Cliques
.size(); i
++ )
266 for( size_t j
= i
+1; j
< Cliques
.size(); j
++ ) {
267 size_t w
= (Cliques
[i
] & Cliques
[j
]).size();
268 JuncGraph
[UEdge(i
,j
)] = w
;
271 // Construct maximal spanning tree using Prim's algorithm
272 _RTree
= MaxSpanningTreePrims( JuncGraph
);
274 // Construct corresponding region graph
276 // Create outer regions
277 ORs
.reserve( Cliques
.size() );
278 for( size_t i
= 0; i
< Cliques
.size(); i
++ )
279 ORs
.push_back( FRegion( Factor(Cliques
[i
], 1.0), 1.0 ) );
281 // For each factor, find an outer region that subsumes that factor.
282 // Then, multiply the outer region with that factor.
283 // If no outer region can be found subsuming that factor, label the
284 // factor as off-tree.
286 fac2OR
.resize( nrFactors(), -1U );
287 for( size_t I
= 0; I
< nrFactors(); I
++ ) {
289 for( alpha
= 0; alpha
< nrORs(); alpha
++ )
290 if( OR(alpha
).vars() >> factor(I
).vars() ) {
294 // DIFF WITH JTree::GenerateJT: assert
298 // Create inner regions and edges
299 IRs
.reserve( _RTree
.size() );
301 edges
.reserve( 2 * _RTree
.size() );
302 for( size_t i
= 0; i
< _RTree
.size(); i
++ ) {
303 edges
.push_back( Edge( _RTree
[i
].n1
, IRs
.size() ) );
304 edges
.push_back( Edge( _RTree
[i
].n2
, IRs
.size() ) );
305 // inner clusters have counting number -1
306 IRs
.push_back( Region( Cliques
[_RTree
[i
].n1
] & Cliques
[_RTree
[i
].n2
], -1.0 ) );
309 // create bipartite graph
310 G
.create( nrORs(), nrIRs(), edges
.begin(), edges
.end() );
312 // Check counting numbers
313 Check_Counting_Numbers();
315 // Create messages and beliefs
317 _Qa
.reserve( nrORs() );
318 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
319 _Qa
.push_back( OR(alpha
) );
322 _Qb
.reserve( nrIRs() );
323 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
324 _Qb
.push_back( Factor( IR(beta
), 1.0 ) );
326 // DIFF with JTree::GenerateJT: no messages
328 // DIFF with JTree::GenerateJT:
329 // Create factor approximations
331 size_t PreviousRoot
= (size_t)-1;
332 for( size_t I
= 0; I
< nrFactors(); I
++ )
334 // find efficient subtree
336 /*size_t subTreeSize =*/ findEfficientTree( factor(I
).vars(), subTree
, PreviousRoot
);
337 PreviousRoot
= subTree
[0].n1
;
338 //subTree.resize( subTreeSize ); // FIXME
339 // cout << "subtree " << I << " has size " << subTreeSize << endl;
343 sprintf( fn, "/tmp/subtree_%d.dot", I );
344 std::ofstream dots(fn);
345 dots << "graph G {" << endl;
346 dots << "graph[size=\"9,9\"];" << endl;
347 dots << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
348 for( size_t i = 0; i < nrVars(); i++ )
349 dots << "\tx" << var(i).label() << ((factor(I).vars() >> var(i)) ? "[color=blue];" : ";") << endl;
350 dots << "node[shape=box,style=filled,color=lightgrey,width=0.3,height=0.3,fixedsize=true];" << endl;
351 for( size_t J = 0; J < nrFactors(); J++ )
352 dots << "\tp" << J << ";" << endl;
353 for( size_t iI = 0; iI < FactorGraph::nr_edges(); iI++ )
354 dots << "\tx" << var(FactorGraph::edge(iI).first).label() << " -- p" << FactorGraph::edge(iI).second << ";" << endl;
355 for( size_t a = 0; a < tree.size(); a++ )
356 dots << "\tx" << var(tree[a].n1).label() << " -- x" << var(tree[a].n2).label() << " [color=red];" << endl;
361 TreeEPSubTree
QI( subTree
, _RTree
, _Qa
, _Qb
, &factor(I
) );
364 // Previous root of first off-tree factor should be the root of the last off-tree factor
365 for( size_t I
= 0; I
< nrFactors(); I
++ )
368 /*size_t subTreeSize =*/ findEfficientTree( factor(I
).vars(), subTree
, PreviousRoot
);
369 PreviousRoot
= subTree
[0].n1
;
370 //subTree.resize( subTreeSize ); // FIXME
371 // cout << "subtree " << I << " has size " << subTreeSize << endl;
373 TreeEPSubTree
QI( subTree
, _RTree
, _Qa
, _Qb
, &factor(I
) );
378 if( Verbose() >= 3 ) {
379 cout
<< "Resulting regiongraph: " << *this << endl
;
384 string
TreeEP::identify() const {
385 stringstream
result (stringstream::out
);
386 result
<< Name
<< GetProperties();
391 void TreeEP::init() {
392 assert( checkProperties() );
396 // Init factor approximations
397 for( size_t I
= 0; I
< nrFactors(); I
++ )
403 double TreeEP::run() {
405 cout
<< "Starting " << identify() << "...";
410 Diffs
diffs(nrVars(), 1.0);
412 vector
<Factor
> old_beliefs
;
413 old_beliefs
.reserve( nrVars() );
414 for( size_t i
= 0; i
< nrVars(); i
++ )
415 old_beliefs
.push_back(belief(var(i
)));
419 // do several passes over the network until maximum number of iterations has
420 // been reached or until the maximum belief difference is smaller than tolerance
421 for( iter
=0; iter
< MaxIter() && diffs
.maxDiff() > Tol(); iter
++ ) {
422 for( size_t I
= 0; I
< nrFactors(); I
++ )
424 _Q
[I
].InvertAndMultiply( _Qa
, _Qb
);
425 _Q
[I
].HUGIN_with_I( _Qa
, _Qb
);
426 _Q
[I
].InvertAndMultiply( _Qa
, _Qb
);
429 // calculate new beliefs and compare with old ones
430 for( size_t i
= 0; i
< nrVars(); i
++ ) {
431 Factor
nb( belief(var(i
)) );
432 diffs
.push( dist( nb
, old_beliefs
[i
], Prob::DISTLINF
) );
437 cout
<< "TreeEP::run: maxdiff " << diffs
.maxDiff() << " after " << iter
+1 << " passes" << endl
;
440 updateMaxDiff( diffs
.maxDiff() );
442 if( Verbose() >= 1 ) {
443 if( diffs
.maxDiff() > Tol() ) {
446 cout
<< "TreeEP::run: WARNING: not converged within " << MaxIter() << " passes (" << toc() - tic
<< " clocks)...final maxdiff:" << diffs
.maxDiff() << endl
;
449 cout
<< "TreeEP::run: ";
450 cout
<< "converged in " << iter
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
454 return diffs
.maxDiff();
458 Complex
TreeEP::logZ() const {
461 // entropy of the tree
462 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
463 sum
-= real(_Qb
[beta
].entropy());
464 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
465 sum
+= real(_Qa
[alpha
].entropy());
467 // energy of the on-tree factors
468 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
469 sum
+= (OR(alpha
).log0() * _Qa
[alpha
]).totalSum();
471 // energy of the off-tree factors
472 for( size_t I
= 0; I
< nrFactors(); I
++ )
474 sum
+= (_Q
.find(I
))->second
.logZ( _Qa
, _Qb
);
480 } // end of namespace dai