1 /* This file is part of libDAI - http://www.libdai.org/
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
9 #include <dai/dai_config.h>
10 #ifdef DAI_WITH_TREEEP
16 #include <dai/jtree.h>
17 #include <dai/treeep.h>
27 void TreeEP::setProperties( const PropertySet
&opts
) {
28 DAI_ASSERT( opts
.hasKey("tol") );
29 DAI_ASSERT( opts
.hasKey("type") );
31 props
.tol
= opts
.getStringAs
<Real
>("tol");
32 props
.type
= opts
.getStringAs
<Properties::TypeType
>("type");
33 if( opts
.hasKey("maxiter") )
34 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
36 props
.maxiter
= 10000;
37 if( opts
.hasKey("maxtime") )
38 props
.maxtime
= opts
.getStringAs
<Real
>("maxtime");
40 props
.maxtime
= INFINITY
;
41 if( opts
.hasKey("verbose") )
42 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
48 PropertySet
TreeEP::getProperties() const {
50 opts
.set( "tol", props
.tol
);
51 opts
.set( "maxiter", props
.maxiter
);
52 opts
.set( "maxtime", props
.maxtime
);
53 opts
.set( "verbose", props
.verbose
);
54 opts
.set( "type", props
.type
);
59 string
TreeEP::printProperties() const {
60 stringstream
s( stringstream::out
);
62 s
<< "tol=" << props
.tol
<< ",";
63 s
<< "maxiter=" << props
.maxiter
<< ",";
64 s
<< "maxtime=" << props
.maxtime
<< ",";
65 s
<< "verbose=" << props
.verbose
<< ",";
66 s
<< "type=" << props
.type
<< "]";
71 TreeEP::TreeEP( const FactorGraph
&fg
, const PropertySet
&opts
) : JTree(fg
, opts("updates",string("HUGIN")), false), _maxdiff(0.0), _iters(0), props(), _Q() {
72 setProperties( opts
);
74 if( opts
.hasKey("tree") ) {
75 construct( fg
, opts
.getAs
<RootedTree
>("tree") );
77 if( props
.type
== Properties::TypeType::ORG
|| props
.type
== Properties::TypeType::ALT
) {
78 // ORG: construct weighted graph with as weights a crude estimate of the
79 // mutual information between the nodes
80 // ALT: construct weighted graph with as weights an upper bound on the
81 // effective interaction strength between pairs of nodes
83 WeightedGraph
<Real
> wg
;
84 // in order to get a connected weighted graph, we start
85 // by connecting every variable to the zero'th variable with weight 0
86 for( size_t i
= 1; i
< fg
.nrVars(); i
++ )
88 for( size_t i
= 0; i
< fg
.nrVars(); i
++ ) {
89 SmallSet
<size_t> delta_i
= fg
.bipGraph().delta1( i
, false );
90 const Var
& v_i
= fg
.var(i
);
91 bforeach( size_t j
, delta_i
)
93 const Var
& v_j
= fg
.var(j
);
94 VarSet
v_ij( v_i
, v_j
);
95 SmallSet
<size_t> nb_ij
= fg
.bipGraph().nb1Set( i
) | fg
.bipGraph().nb1Set( j
);
97 bforeach( size_t I
, nb_ij
) {
98 const VarSet
& Ivars
= fg
.factor(I
).vars();
99 if( props
.type
== Properties::TypeType::ORG
) {
100 if( (Ivars
== v_i
) || (Ivars
== v_j
) )
101 piet
*= fg
.factor(I
);
102 else if( Ivars
>> v_ij
)
103 piet
*= fg
.factor(I
).marginal( v_ij
);
106 piet
*= fg
.factor(I
);
109 if( props
.type
== Properties::TypeType::ORG
) {
110 if( piet
.vars() >> v_ij
) {
111 piet
= piet
.marginal( v_ij
);
112 Factor pietf
= piet
.marginal(v_i
) * piet
.marginal(v_j
);
113 wg
[UEdge(i
,j
)] = dist( piet
, pietf
, DISTKL
);
115 // this should never happen...
116 DAI_ASSERT( 0 == 1 );
120 wg
[UEdge(i
,j
)] = piet
.strength(v_i
, v_j
);
124 // find maximal spanning tree
125 if( props
.verbose
>= 3 )
126 cerr
<< "WeightedGraph: " << wg
<< endl
;
127 RootedTree t
= MaxSpanningTree( wg
, /*true*/false ); // WORKAROUND FOR BUG IN BOOST GRAPH LIBRARY VERSION 1.54
128 if( props
.verbose
>= 3 )
129 cerr
<< "Spanningtree: " << t
<< endl
;
132 DAI_THROW(UNKNOWN_ENUM_VALUE
);
137 void TreeEP::construct( const FactorGraph
& fg
, const RootedTree
& tree
) {
138 // Copy the factor graph
139 FactorGraph::operator=( fg
);
142 for( size_t i
= 0; i
< tree
.size(); i
++ )
143 cl
.push_back( VarSet( var(tree
[i
].first
), var(tree
[i
].second
) ) );
145 // If no outer region can be found subsuming that factor, label the
146 // factor as off-tree.
147 JTree::construct( *this, cl
, false );
149 if( props
.verbose
>= 1 )
150 cerr
<< "TreeEP::construct: The tree has size " << JTree::RTree
.size() << endl
;
151 if( props
.verbose
>= 3 )
152 cerr
<< " it is " << JTree::RTree
<< " with cliques " << cl
<< endl
;
154 // Create factor approximations
156 size_t PreviousRoot
= (size_t)-1;
157 // Second repetition: previous root of first off-tree factor should be the root of the last off-tree factor
158 for( size_t repeats
= 0; repeats
< 2; repeats
++ )
159 for( size_t I
= 0; I
< nrFactors(); I
++ )
161 // find efficient subtree
163 size_t subTreeSize
= findEfficientTree( factor(I
).vars(), subTree
, PreviousRoot
);
164 PreviousRoot
= subTree
[0].first
;
165 subTree
.resize( subTreeSize
);
166 if( props
.verbose
>= 1 )
167 cerr
<< "Subtree " << I
<< " has size " << subTreeSize
<< endl
;
168 if( props
.verbose
>= 3 )
169 cerr
<< " it is " << subTree
<< endl
;
170 _Q
[I
] = TreeEPSubTree( subTree
, RTree
, Qa
, Qb
, &factor(I
) );
175 if( props
.verbose
>= 3 )
176 cerr
<< "Resulting regiongraph: " << *this << endl
;
180 void TreeEP::init() {
183 // Init factor approximations
184 for( size_t I
= 0; I
< nrFactors(); I
++ )
191 if( props
.verbose
>= 1 )
192 cerr
<< "Starting " << identify() << "...";
193 if( props
.verbose
>= 3 )
198 vector
<Factor
> oldBeliefs
= beliefs();
200 // do several passes over the network until maximum number of iterations has
201 // been reached or until the maximum belief difference is smaller than tolerance
202 Real maxDiff
= INFINITY
;
203 for( _iters
= 0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
&& (toc() - tic
) < props
.maxtime
; _iters
++ ) {
204 for( size_t I
= 0; I
< nrFactors(); I
++ )
206 _Q
[I
].InvertAndMultiply( Qa
, Qb
);
207 _Q
[I
].HUGIN_with_I( Qa
, Qb
);
208 _Q
[I
].InvertAndMultiply( Qa
, Qb
);
211 // calculate new beliefs and compare with old ones
212 vector
<Factor
> newBeliefs
= beliefs();
214 for( size_t t
= 0; t
< oldBeliefs
.size(); t
++ )
215 maxDiff
= std::max( maxDiff
, dist( newBeliefs
[t
], oldBeliefs
[t
], DISTLINF
) );
216 swap( newBeliefs
, oldBeliefs
);
218 if( props
.verbose
>= 3 )
219 cerr
<< name() << "::run: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
222 if( maxDiff
> _maxdiff
)
225 if( props
.verbose
>= 1 ) {
226 if( maxDiff
> props
.tol
) {
227 if( props
.verbose
== 1 )
229 cerr
<< name() << "::run: WARNING: not converged after " << _iters
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
231 if( props
.verbose
>= 3 )
232 cerr
<< name() << "::run: ";
233 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
241 Real
TreeEP::logZ() const {
244 // entropy of the tree
245 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
246 s
-= Qb
[beta
].entropy();
247 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
248 s
+= Qa
[alpha
].entropy();
250 // energy of the on-tree factors
251 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
252 s
+= (OR(alpha
).log(true) * Qa
[alpha
]).sum();
254 // energy of the off-tree factors
255 for( size_t I
= 0; I
< nrFactors(); I
++ )
257 s
+= (_Q
.find(I
))->second
.logZ( Qa
, Qb
);
263 TreeEP::TreeEPSubTree::TreeEPSubTree( const RootedTree
&subRTree
, const RootedTree
&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) {
266 // Make _Qa, _Qb, _a and _b corresponding to the subtree
267 _b
.reserve( subRTree
.size() );
268 _Qb
.reserve( subRTree
.size() );
269 _RTree
.reserve( subRTree
.size() );
270 for( size_t i
= 0; i
< subRTree
.size(); i
++ ) {
271 size_t alpha1
= subRTree
[i
].first
; // old index 1
272 size_t alpha2
= subRTree
[i
].second
; // old index 2
273 size_t beta
; // old sep index
274 for( beta
= 0; beta
< jt_RTree
.size(); beta
++ )
275 if( UEdge( jt_RTree
[beta
].first
, jt_RTree
[beta
].second
) == UEdge( alpha1
, alpha2
) )
277 DAI_ASSERT( beta
!= jt_RTree
.size() );
279 size_t newalpha1
= find(_a
.begin(), _a
.end(), alpha1
) - _a
.begin();
280 if( newalpha1
== _a
.size() ) {
281 _Qa
.push_back( Factor( jt_Qa
[alpha1
].vars(), 1.0 ) );
282 _a
.push_back( alpha1
); // save old index in index conversion table
285 size_t newalpha2
= find(_a
.begin(), _a
.end(), alpha2
) - _a
.begin();
286 if( newalpha2
== _a
.size() ) {
287 _Qa
.push_back( Factor( jt_Qa
[alpha2
].vars(), 1.0 ) );
288 _a
.push_back( alpha2
); // save old index in index conversion table
291 _RTree
.push_back( DEdge( newalpha1
, newalpha2
) );
292 _Qb
.push_back( Factor( jt_Qb
[beta
].vars(), 1.0 ) );
293 _b
.push_back( beta
);
296 // Find remaining variables (which are not in the new root)
297 _nsrem
= _ns
/ _Qa
[0].vars();
301 void TreeEP::TreeEPSubTree::init() {
302 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
303 _Qa
[alpha
].fill( 1.0 );
304 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
305 _Qb
[beta
].fill( 1.0 );
309 void TreeEP::TreeEPSubTree::InvertAndMultiply( const std::vector
<Factor
> &Qa
, const std::vector
<Factor
> &Qb
) {
310 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
311 _Qa
[alpha
] = Qa
[_a
[alpha
]] / _Qa
[alpha
];
313 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
314 _Qb
[beta
] = Qb
[_b
[beta
]] / _Qb
[beta
];
318 void TreeEP::TreeEPSubTree::HUGIN_with_I( std::vector
<Factor
> &Qa
, std::vector
<Factor
> &Qb
) {
319 // Backup _Qa and _Qb
320 vector
<Factor
> _Qa_old(_Qa
);
321 vector
<Factor
> _Qb_old(_Qb
);
324 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
325 Qa
[_a
[alpha
]].fill( 0.0 );
326 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
327 Qb
[_b
[beta
]].fill( 0.0 );
329 // For all states of _nsrem
330 for( State
s(_nsrem
); s
.valid(); s
++ ) {
331 // Multiply root with slice of I
332 _Qa
[0] *= _I
->slice( _nsrem
, s
);
335 for( size_t i
= _RTree
.size(); (i
--) != 0; ) {
336 // clamp variables in nsrem
337 for( VarSet::const_iterator n
= _nsrem
.begin(); n
!= _nsrem
.end(); n
++ )
338 if( _Qa
[_RTree
[i
].second
].vars() >> *n
)
339 _Qa
[_RTree
[i
].second
] *= createFactorDelta( *n
, s(*n
) );
340 Factor new_Qb
= _Qa
[_RTree
[i
].second
].marginal( _Qb
[i
].vars(), false );
341 _Qa
[_RTree
[i
].first
] *= new_Qb
/ _Qb
[i
];
345 // DistributeEvidence
346 for( size_t i
= 0; i
< _RTree
.size(); i
++ ) {
347 Factor new_Qb
= _Qa
[_RTree
[i
].first
].marginal( _Qb
[i
].vars(), false );
348 _Qa
[_RTree
[i
].second
] *= new_Qb
/ _Qb
[i
];
352 // Store Qa's and Qb's
353 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
354 Qa
[_a
[alpha
]].p() += _Qa
[alpha
].p();
355 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
356 Qb
[_b
[beta
]].p() += _Qb
[beta
].p();
358 // Restore _Qa and _Qb
363 // Normalize Qa and Qb
365 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ ) {
366 _logZ
+= log(Qa
[_a
[alpha
]].sum());
367 Qa
[_a
[alpha
]].normalize();
369 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ ) {
370 _logZ
-= log(Qb
[_b
[beta
]].sum());
371 Qb
[_b
[beta
]].normalize();
376 Real
TreeEP::TreeEPSubTree::logZ( const std::vector
<Factor
> &Qa
, const std::vector
<Factor
> &Qb
) const {
378 for( size_t alpha
= 0; alpha
< _Qa
.size(); alpha
++ )
379 s
+= (Qa
[_a
[alpha
]] * _Qa
[alpha
].log(true)).sum();
380 for( size_t beta
= 0; beta
< _Qb
.size(); beta
++ )
381 s
-= (Qb
[_b
[beta
]] * _Qb
[beta
].log(true)).sum();
386 } // end of namespace dai