9c2ea901d5f44a9bc88593a8180203b5034f6651
1 /* This file is part of libDAI - http://www.libdai.org/
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
20 #include <dai/properties.h>
29 const char *BP::Name
= "BP";
35 void BP::setProperties( const PropertySet
&opts
) {
36 DAI_ASSERT( opts
.hasKey("tol") );
37 DAI_ASSERT( opts
.hasKey("maxiter") );
38 DAI_ASSERT( opts
.hasKey("logdomain") );
39 DAI_ASSERT( opts
.hasKey("updates") );
41 props
.tol
= opts
.getStringAs
<Real
>("tol");
42 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
43 props
.logdomain
= opts
.getStringAs
<bool>("logdomain");
44 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
46 if( opts
.hasKey("verbose") )
47 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
50 if( opts
.hasKey("damping") )
51 props
.damping
= opts
.getStringAs
<Real
>("damping");
54 if( opts
.hasKey("inference") )
55 props
.inference
= opts
.getStringAs
<Properties::InfType
>("inference");
57 props
.inference
= Properties::InfType::SUMPROD
;
61 PropertySet
BP::getProperties() const {
63 opts
.Set( "tol", props
.tol
);
64 opts
.Set( "maxiter", props
.maxiter
);
65 opts
.Set( "verbose", props
.verbose
);
66 opts
.Set( "logdomain", props
.logdomain
);
67 opts
.Set( "updates", props
.updates
);
68 opts
.Set( "damping", props
.damping
);
69 opts
.Set( "inference", props
.inference
);
74 string
BP::printProperties() const {
75 stringstream
s( stringstream::out
);
77 s
<< "tol=" << props
.tol
<< ",";
78 s
<< "maxiter=" << props
.maxiter
<< ",";
79 s
<< "verbose=" << props
.verbose
<< ",";
80 s
<< "logdomain=" << props
.logdomain
<< ",";
81 s
<< "updates=" << props
.updates
<< ",";
82 s
<< "damping=" << props
.damping
<< ",";
83 s
<< "inference=" << props
.inference
<< "]";
88 void BP::construct() {
89 // create edge properties
91 _edges
.reserve( nrVars() );
93 if( props
.updates
== Properties::UpdateType::SEQMAX
)
94 _edge2lut
.reserve( nrVars() );
95 for( size_t i
= 0; i
< nrVars(); ++i
) {
96 _edges
.push_back( vector
<EdgeProp
>() );
97 _edges
[i
].reserve( nbV(i
).size() );
98 if( props
.updates
== Properties::UpdateType::SEQMAX
) {
99 _edge2lut
.push_back( vector
<LutType::iterator
>() );
100 _edge2lut
[i
].reserve( nbV(i
).size() );
102 foreach( const Neighbor
&I
, nbV(i
) ) {
104 newEP
.message
= Prob( var(i
).states() );
105 newEP
.newMessage
= Prob( var(i
).states() );
108 newEP
.index
.reserve( factor(I
).states() );
109 for( IndexFor
k( var(i
), factor(I
).vars() ); k
.valid(); ++k
)
110 newEP
.index
.push_back( k
);
113 newEP
.residual
= 0.0;
114 _edges
[i
].push_back( newEP
);
115 if( props
.updates
== Properties::UpdateType::SEQMAX
)
116 _edge2lut
[i
].push_back( _lut
.insert( make_pair( newEP
.residual
, make_pair( i
, _edges
[i
].size() - 1 ))) );
123 Real c
= props
.logdomain
? 0.0 : 1.0;
124 for( size_t i
= 0; i
< nrVars(); ++i
) {
125 foreach( const Neighbor
&I
, nbV(i
) ) {
126 message( i
, I
.iter
).fill( c
);
127 newMessage( i
, I
.iter
).fill( c
);
128 if( props
.updates
== Properties::UpdateType::SEQMAX
)
129 updateResidual( i
, I
.iter
, 0.0 );
135 void BP::findMaxResidual( size_t &i
, size_t &_I
) {
136 DAI_ASSERT( !_lut
.empty() );
137 LutType::const_iterator largestEl
= _lut
.end();
139 i
= largestEl
->second
.first
;
140 _I
= largestEl
->second
.second
;
144 Prob
BP::calcIncomingMessageProduct( size_t I
, bool without_i
, size_t i
) const {
145 Factor
Fprod( factor(I
) );
146 Prob
&prod
= Fprod
.p();
147 if( props
.logdomain
)
150 // Calculate product of incoming messages and factor I
151 foreach( const Neighbor
&j
, nbF(I
) )
152 if( !(without_i
&& (j
== i
)) ) {
153 // prod_j will be the product of messages coming into j
154 Prob
prod_j( var(j
).states(), props
.logdomain
? 0.0 : 1.0 );
155 foreach( const Neighbor
&J
, nbV(j
) )
156 if( J
!= I
) { // for all J in nb(j) \ I
157 if( props
.logdomain
)
158 prod_j
+= message( j
, J
.iter
);
160 prod_j
*= message( j
, J
.iter
);
163 // multiply prod with prod_j
165 // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
166 if( props
.logdomain
)
167 Fprod
+= Factor( var(j
), prod_j
);
169 Fprod
*= Factor( var(j
), prod_j
);
173 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
174 const ind_t
&ind
= index(j
, _I
);
176 for( size_t r
= 0; r
< prod
.size(); ++r
)
177 if( props
.logdomain
)
178 prod
[r
] += prod_j
[ind
[r
]];
180 prod
[r
] *= prod_j
[ind
[r
]];
187 void BP::calcNewMessage( size_t i
, size_t _I
) {
188 // calculate updated message I->i
189 size_t I
= nbV(i
,_I
);
192 if( factor(I
).vars().size() == 1 ) // optimization
193 marg
= factor(I
).p();
195 Factor
Fprod( factor(I
) );
196 Prob
&prod
= Fprod
.p();
197 prod
= calcIncomingMessageProduct( I
, true, i
);
199 if( props
.logdomain
) {
204 // Marginalize onto i
206 /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
207 if( props
.inference
== Properties::InfType::SUMPROD
)
208 marg
= Fprod
.marginal( var(i
) ).p();
210 marg
= Fprod
.maxMarginal( var(i
) ).p();
212 /* OPTIMIZED VERSION */
213 marg
= Prob( var(i
).states(), 0.0 );
214 // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
215 const ind_t ind
= index(i
,_I
);
216 if( props
.inference
== Properties::InfType::SUMPROD
)
217 for( size_t r
= 0; r
< prod
.size(); ++r
)
218 marg
[ind
[r
]] += prod
[r
];
220 for( size_t r
= 0; r
< prod
.size(); ++r
)
221 if( prod
[r
] > marg
[ind
[r
]] )
222 marg
[ind
[r
]] = prod
[r
];
228 if( props
.logdomain
)
229 newMessage(i
,_I
) = marg
.log();
231 newMessage(i
,_I
) = marg
;
233 // Update the residual if necessary
234 if( props
.updates
== Properties::UpdateType::SEQMAX
)
235 updateResidual( i
, _I
, dist( newMessage( i
, _I
), message( i
, _I
), Prob::DISTLINF
) );
239 // BP::run does not check for NANs for performance reasons
240 // Somehow NaNs do not often occur in BP...
242 if( props
.verbose
>= 1 )
243 cerr
<< "Starting " << identify() << "...";
244 if( props
.verbose
>= 3)
248 Real maxDiff
= INFINITY
;
250 vector
<Factor
> oldBeliefsV
, oldBeliefsF
;
251 oldBeliefsV
.reserve( nrVars() );
252 for( size_t i
= 0; i
< nrVars(); ++i
)
253 oldBeliefsV
.push_back( beliefV(i
) );
254 oldBeliefsF
.reserve( nrFactors() );
255 for( size_t I
= 0; I
< nrFactors(); ++I
)
256 oldBeliefsF
.push_back( beliefF(I
) );
258 size_t nredges
= nrEdges();
259 vector
<Edge
> update_seq
;
260 if( props
.updates
== Properties::UpdateType::SEQMAX
) {
262 for( size_t i
= 0; i
< nrVars(); ++i
)
263 foreach( const Neighbor
&I
, nbV(i
) )
264 calcNewMessage( i
, I
.iter
);
266 update_seq
.reserve( nredges
);
267 for( size_t I
= 0; I
< nrFactors(); I
++ )
268 foreach( const Neighbor
&i
, nbF(I
) )
269 update_seq
.push_back( Edge( i
, i
.dual
) );
272 // do several passes over the network until maximum number of iterations has
273 // been reached or until the maximum belief difference is smaller than tolerance
274 for( _iters
=0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
; ++_iters
) {
275 if( props
.updates
== Properties::UpdateType::SEQMAX
) {
276 // Residuals-BP by Koller et al.
277 for( size_t t
= 0; t
< nredges
; ++t
) {
278 // update the message with the largest residual
280 findMaxResidual( i
, _I
);
281 updateMessage( i
, _I
);
283 // I->i has been updated, which means that residuals for all
284 // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
285 foreach( const Neighbor
&J
, nbV(i
) ) {
287 foreach( const Neighbor
&j
, nbF(J
) ) {
290 calcNewMessage( j
, _J
);
295 } else if( props
.updates
== Properties::UpdateType::PARALL
) {
297 for( size_t i
= 0; i
< nrVars(); ++i
)
298 foreach( const Neighbor
&I
, nbV(i
) )
299 calcNewMessage( i
, I
.iter
);
301 for( size_t i
= 0; i
< nrVars(); ++i
)
302 foreach( const Neighbor
&I
, nbV(i
) )
303 updateMessage( i
, I
.iter
);
305 // Sequential updates
306 if( props
.updates
== Properties::UpdateType::SEQRND
)
307 random_shuffle( update_seq
.begin(), update_seq
.end() );
309 foreach( const Edge
&e
, update_seq
) {
310 calcNewMessage( e
.first
, e
.second
);
311 updateMessage( e
.first
, e
.second
);
315 // calculate new beliefs and compare with old ones
317 for( size_t i
= 0; i
< nrVars(); ++i
) {
318 Factor
b( beliefV(i
) );
319 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], Prob::DISTLINF
) );
322 for( size_t I
= 0; I
< nrFactors(); ++I
) {
323 Factor
b( beliefF(I
) );
324 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], Prob::DISTLINF
) );
328 if( props
.verbose
>= 3 )
329 cerr
<< Name
<< "::run: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
332 if( maxDiff
> _maxdiff
)
335 if( props
.verbose
>= 1 ) {
336 if( maxDiff
> props
.tol
) {
337 if( props
.verbose
== 1 )
339 cerr
<< Name
<< "::run: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
341 if( props
.verbose
>= 3 )
342 cerr
<< Name
<< "::run: ";
343 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
351 void BP::calcBeliefV( size_t i
, Prob
&p
) const {
352 p
= Prob( var(i
).states(), props
.logdomain
? 0.0 : 1.0 );
353 foreach( const Neighbor
&I
, nbV(i
) )
354 if( props
.logdomain
)
355 p
+= newMessage( i
, I
.iter
);
357 p
*= newMessage( i
, I
.iter
);
361 Factor
BP::beliefV( size_t i
) const {
365 if( props
.logdomain
) {
371 return( Factor( var(i
), p
) );
375 Factor
BP::beliefF( size_t I
) const {
379 if( props
.logdomain
) {
385 return( Factor( factor(I
).vars(), p
) );
389 vector
<Factor
> BP::beliefs() const {
390 vector
<Factor
> result
;
391 for( size_t i
= 0; i
< nrVars(); ++i
)
392 result
.push_back( beliefV(i
) );
393 for( size_t I
= 0; I
< nrFactors(); ++I
)
394 result
.push_back( beliefF(I
) );
399 Factor
BP::belief( const VarSet
&ns
) const {
402 else if( ns
.size() == 1 )
403 return beliefV( findVar( *(ns
.begin() ) ) );
406 for( I
= 0; I
< nrFactors(); I
++ )
407 if( factor(I
).vars() >> ns
)
409 if( I
== nrFactors() )
410 DAI_THROW(BELIEF_NOT_AVAILABLE
);
411 return beliefF(I
).marginal(ns
);
416 Real
BP::logZ() const {
418 for( size_t i
= 0; i
< nrVars(); ++i
)
419 sum
+= (1.0 - nbV(i
).size()) * beliefV(i
).entropy();
420 for( size_t I
= 0; I
< nrFactors(); ++I
)
421 sum
-= dist( beliefF(I
), factor(I
), Prob::DISTKL
);
426 string
BP::identify() const {
427 return string(Name
) + printProperties();
431 void BP::init( const VarSet
&ns
) {
432 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); ++n
) {
433 size_t ni
= findVar( *n
);
434 foreach( const Neighbor
&I
, nbV( ni
) ) {
435 Real val
= props
.logdomain
? 0.0 : 1.0;
436 message( ni
, I
.iter
).fill( val
);
437 newMessage( ni
, I
.iter
).fill( val
);
438 if( props
.updates
== Properties::UpdateType::SEQMAX
)
439 updateResidual( ni
, I
.iter
, 0.0 );
445 void BP::updateMessage( size_t i
, size_t _I
) {
446 if( recordSentMessages
)
447 _sentMessages
.push_back(make_pair(i
,_I
));
448 if( props
.damping
== 0.0 ) {
449 message(i
,_I
) = newMessage(i
,_I
);
450 if( props
.updates
== Properties::UpdateType::SEQMAX
)
451 updateResidual( i
, _I
, 0.0 );
453 if( props
.logdomain
)
454 message(i
,_I
) = (message(i
,_I
) * props
.damping
) + (newMessage(i
,_I
) * (1.0 - props
.damping
));
456 message(i
,_I
) = (message(i
,_I
) ^ props
.damping
) * (newMessage(i
,_I
) ^ (1.0 - props
.damping
));
457 if( props
.updates
== Properties::UpdateType::SEQMAX
)
458 updateResidual( i
, _I
, dist( newMessage(i
,_I
), message(i
,_I
), Prob::DISTLINF
) );
463 void BP::updateResidual( size_t i
, size_t _I
, Real r
) {
464 EdgeProp
* pEdge
= &_edges
[i
][_I
];
467 // rearrange look-up table (delete and reinsert new key)
468 _lut
.erase( _edge2lut
[i
][_I
] );
469 _edge2lut
[i
][_I
] = _lut
.insert( make_pair( r
, make_pair(i
, _I
) ) );
473 std::vector
<size_t> BP::findMaximum() const {
474 vector
<size_t> maximum( nrVars() );
475 vector
<bool> visitedVars( nrVars(), false );
476 vector
<bool> visitedFactors( nrFactors(), false );
477 stack
<size_t> scheduledFactors
;
478 for( size_t i
= 0; i
< nrVars(); ++i
) {
481 visitedVars
[i
] = true;
483 // Maximise with respect to variable i
485 calcBeliefV( i
, prod
);
486 maximum
[i
] = prod
.argmax().first
;
488 foreach( const Neighbor
&I
, nbV(i
) )
489 if( !visitedFactors
[I
] )
490 scheduledFactors
.push(I
);
492 while( !scheduledFactors
.empty() ){
493 size_t I
= scheduledFactors
.top();
494 scheduledFactors
.pop();
495 if( visitedFactors
[I
] )
497 visitedFactors
[I
] = true;
499 // Evaluate if some neighboring variables still need to be fixed; if not, we're done
500 bool allDetermined
= true;
501 foreach( const Neighbor
&j
, nbF(I
) )
502 if( !visitedVars
[j
.node
] ) {
503 allDetermined
= false;
509 // Calculate product of incoming messages on factor I
511 calcBeliefF( I
, prod2
);
513 // The allowed configuration is restrained according to the variables assigned so far:
514 // pick the argmax amongst the allowed states
515 Real maxProb
= numeric_limits
<Real
>::min();
516 State
maxState( factor(I
).vars() );
517 for( State
s( factor(I
).vars() ); s
.valid(); ++s
){
518 // First, calculate whether this state is consistent with variables that
519 // have been assigned already
520 bool allowedState
= true;
521 foreach( const Neighbor
&j
, nbF(I
) )
522 if( visitedVars
[j
.node
] && maximum
[j
.node
] != s(var(j
.node
)) ) {
523 allowedState
= false;
526 // If it is consistent, check if its probability is larger than what we have seen so far
527 if( allowedState
&& prod2
[s
] > maxProb
) {
534 foreach( const Neighbor
&j
, nbF(I
) ) {
535 if( visitedVars
[j
.node
] ) {
536 // We have already visited j earlier - hopefully our state is consistent
537 if( maximum
[j
.node
] != maxState(var(j
.node
)) && props
.verbose
>= 1 )
538 cerr
<< "BP::findMaximum - warning: maximum not consistent due to loops." << endl
;
540 // We found a consistent state for variable j
541 visitedVars
[j
.node
] = true;
542 maximum
[j
.node
] = maxState( var(j
.node
) );
543 foreach( const Neighbor
&J
, nbV(j
) )
544 if( !visitedFactors
[J
] )
545 scheduledFactors
.push(J
);
554 } // end of namespace dai