1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
5 This file is part of libDAI.
7 libDAI is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
12 libDAI is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License
18 along with libDAI; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29 #include <dai/diffs.h>
31 #include <dai/properties.h>
40 const char *BP::Name
= "BP";
43 void BP::setProperties( const PropertySet
&opts
) {
44 assert( opts
.hasKey("tol") );
45 assert( opts
.hasKey("maxiter") );
46 assert( opts
.hasKey("logdomain") );
47 assert( opts
.hasKey("updates") );
49 props
.tol
= opts
.getStringAs
<double>("tol");
50 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
51 props
.logdomain
= opts
.getStringAs
<bool>("logdomain");
52 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
54 if( opts
.hasKey("verbose") )
55 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
58 if( opts
.hasKey("damping") )
59 props
.damping
= opts
.getStringAs
<double>("damping");
62 if( opts
.hasKey("inference") )
63 props
.inference
= opts
.getStringAs
<Properties::InfType
>("inference");
65 props
.inference
= Properties::InfType::SUMPROD
;
69 PropertySet
BP::getProperties() const {
71 opts
.Set( "tol", props
.tol
);
72 opts
.Set( "maxiter", props
.maxiter
);
73 opts
.Set( "verbose", props
.verbose
);
74 opts
.Set( "logdomain", props
.logdomain
);
75 opts
.Set( "updates", props
.updates
);
76 opts
.Set( "damping", props
.damping
);
77 opts
.Set( "inference", props
.inference
);
82 string
BP::printProperties() const {
83 stringstream
s( stringstream::out
);
85 s
<< "tol=" << props
.tol
<< ",";
86 s
<< "maxiter=" << props
.maxiter
<< ",";
87 s
<< "verbose=" << props
.verbose
<< ",";
88 s
<< "logdomain=" << props
.logdomain
<< ",";
89 s
<< "updates=" << props
.updates
<< ",";
90 s
<< "damping=" << props
.damping
<< ",";
91 s
<< "inference=" << props
.inference
<< "]";
96 void BP::construct() {
97 // create edge properties
99 _edges
.reserve( nrVars() );
100 for( size_t i
= 0; i
< nrVars(); ++i
) {
101 _edges
.push_back( vector
<EdgeProp
>() );
102 _edges
[i
].reserve( nbV(i
).size() );
103 foreach( const Neighbor
&I
, nbV(i
) ) {
105 newEP
.message
= Prob( var(i
).states() );
106 newEP
.newMessage
= Prob( var(i
).states() );
108 newEP
.index
.reserve( factor(I
).states() );
109 for( IndexFor
k( var(i
), factor(I
).vars() ); k
>= 0; ++k
)
110 newEP
.index
.push_back( k
);
112 newEP
.residual
= 0.0;
113 _edges
[i
].push_back( newEP
);
120 double c
= props
.logdomain
? 0.0 : 1.0;
121 for( size_t i
= 0; i
< nrVars(); ++i
) {
122 foreach( const Neighbor
&I
, nbV(i
) ) {
123 message( i
, I
.iter
).fill( c
);
124 newMessage( i
, I
.iter
).fill( c
);
130 void BP::findMaxResidual( size_t &i
, size_t &_I
) {
133 double maxres
= residual( i
, _I
);
134 for( size_t j
= 0; j
< nrVars(); ++j
)
135 foreach( const Neighbor
&I
, nbV(j
) )
136 if( residual( j
, I
.iter
) > maxres
) {
139 maxres
= residual( i
, _I
);
144 void BP::calcNewMessage( size_t i
, size_t _I
) {
145 // calculate updated message I->i
146 size_t I
= nbV(i
,_I
);
149 /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
150 Factor
prod( factor( I
) );
151 foreach( const Neighbor
&j
, nbF(I
) )
152 if( j
!= i
) { // for all j in I \ i
153 foreach( const Neighbor
&J
, nbV(j
) )
154 if( J
!= I
) { // for all J in nb(j) \ I
155 prod
*= Factor( var(j
), message(j
, J
.iter
) );
158 newMessage(i
,_I
) = prod
.marginal( var(i
) ).p();
160 /* OPTIMIZED VERSION */
161 Prob
prod( factor(I
).p() );
162 if( props
.logdomain
)
165 // Calculate product of incoming messages and factor I
166 foreach( const Neighbor
&j
, nbF(I
) ) {
167 if( j
!= i
) { // for all j in I \ i
169 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
170 const ind_t
&ind
= index(j
, _I
);
172 // prod_j will be the product of messages coming into j
173 Prob
prod_j( var(j
).states(), props
.logdomain
? 0.0 : 1.0 );
174 foreach( const Neighbor
&J
, nbV(j
) )
175 if( J
!= I
) { // for all J in nb(j) \ I
176 if( props
.logdomain
)
177 prod_j
+= message( j
, J
.iter
);
179 prod_j
*= message( j
, J
.iter
);
182 // multiply prod with prod_j
183 for( size_t r
= 0; r
< prod
.size(); ++r
)
184 if( props
.logdomain
)
185 prod
[r
] += prod_j
[ind
[r
]];
187 prod
[r
] *= prod_j
[ind
[r
]];
190 if( props
.logdomain
) {
191 prod
-= prod
.maxVal();
195 // Marginalize onto i
196 Prob
marg( var(i
).states(), 0.0 );
197 // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
198 const ind_t ind
= index(i
,_I
);
199 if( props
.inference
== Properties::InfType::SUMPROD
)
200 for( size_t r
= 0; r
< prod
.size(); ++r
)
201 marg
[ind
[r
]] += prod
[r
];
203 for( size_t r
= 0; r
< prod
.size(); ++r
)
204 if( prod
[r
] > marg
[ind
[r
]] )
205 marg
[ind
[r
]] = prod
[r
];
209 if( props
.logdomain
)
210 newMessage(i
,_I
) = marg
.log();
212 newMessage(i
,_I
) = marg
;
217 // BP::run does not check for NANs for performance reasons
218 // Somehow NaNs do not often occur in BP...
220 if( props
.verbose
>= 1 )
221 cout
<< "Starting " << identify() << "...";
222 if( props
.verbose
>= 3)
226 Diffs
diffs(nrVars(), 1.0);
228 vector
<Edge
> update_seq
;
230 vector
<Factor
> old_beliefs
;
231 old_beliefs
.reserve( nrVars() );
232 for( size_t i
= 0; i
< nrVars(); ++i
)
233 old_beliefs
.push_back( beliefV(i
) );
235 size_t nredges
= nrEdges();
237 if( props
.updates
== Properties::UpdateType::SEQMAX
) {
239 for( size_t i
= 0; i
< nrVars(); ++i
)
240 foreach( const Neighbor
&I
, nbV(i
) ) {
241 calcNewMessage( i
, I
.iter
);
242 // calculate initial residuals
243 residual( i
, I
.iter
) = dist( newMessage( i
, I
.iter
), message( i
, I
.iter
), Prob::DISTLINF
);
246 update_seq
.reserve( nredges
);
247 for( size_t i
= 0; i
< nrVars(); ++i
)
248 foreach( const Neighbor
&I
, nbV(i
) )
249 update_seq
.push_back( Edge( i
, I
.iter
) );
252 // do several passes over the network until maximum number of iterations has
253 // been reached or until the maximum belief difference is smaller than tolerance
254 for( _iters
=0; _iters
< props
.maxiter
&& diffs
.maxDiff() > props
.tol
; ++_iters
) {
255 if( props
.updates
== Properties::UpdateType::SEQMAX
) {
256 // Residuals-BP by Koller et al.
257 for( size_t t
= 0; t
< nredges
; ++t
) {
258 // update the message with the largest residual
260 findMaxResidual( i
, _I
);
261 updateMessage( i
, _I
);
263 // I->i has been updated, which means that residuals for all
264 // J->j with J in nb[i]\I and j in nb[J]\i have to be updated
265 foreach( const Neighbor
&J
, nbV(i
) ) {
267 foreach( const Neighbor
&j
, nbF(J
) ) {
270 calcNewMessage( j
, _J
);
271 residual( j
, _J
) = dist( newMessage( j
, _J
), message( j
, _J
), Prob::DISTLINF
);
277 } else if( props
.updates
== Properties::UpdateType::PARALL
) {
279 for( size_t i
= 0; i
< nrVars(); ++i
)
280 foreach( const Neighbor
&I
, nbV(i
) )
281 calcNewMessage( i
, I
.iter
);
283 for( size_t i
= 0; i
< nrVars(); ++i
)
284 foreach( const Neighbor
&I
, nbV(i
) )
285 updateMessage( i
, I
.iter
);
287 // Sequential updates
288 if( props
.updates
== Properties::UpdateType::SEQRND
)
289 random_shuffle( update_seq
.begin(), update_seq
.end() );
291 foreach( const Edge
&e
, update_seq
) {
292 calcNewMessage( e
.first
, e
.second
);
293 updateMessage( e
.first
, e
.second
);
297 // calculate new beliefs and compare with old ones
298 for( size_t i
= 0; i
< nrVars(); ++i
) {
299 Factor
nb( beliefV(i
) );
300 diffs
.push( dist( nb
, old_beliefs
[i
], Prob::DISTLINF
) );
304 if( props
.verbose
>= 3 )
305 cout
<< Name
<< "::run: maxdiff " << diffs
.maxDiff() << " after " << _iters
+1 << " passes" << endl
;
308 if( diffs
.maxDiff() > _maxdiff
)
309 _maxdiff
= diffs
.maxDiff();
311 if( props
.verbose
>= 1 ) {
312 if( diffs
.maxDiff() > props
.tol
) {
313 if( props
.verbose
== 1 )
315 cout
<< Name
<< "::run: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << diffs
.maxDiff() << endl
;
317 if( props
.verbose
>= 3 )
318 cout
<< Name
<< "::run: ";
319 cout
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
323 return diffs
.maxDiff();
327 Factor
BP::beliefV( size_t i
) const {
328 Prob
prod( var(i
).states(), props
.logdomain
? 0.0 : 1.0 );
329 foreach( const Neighbor
&I
, nbV(i
) )
330 if( props
.logdomain
)
331 prod
+= newMessage( i
, I
.iter
);
333 prod
*= newMessage( i
, I
.iter
);
334 if( props
.logdomain
) {
335 prod
-= prod
.maxVal();
340 return( Factor( var(i
), prod
) );
344 Factor
BP::belief (const Var
&n
) const {
345 return( beliefV( findVar( n
) ) );
349 vector
<Factor
> BP::beliefs() const {
350 vector
<Factor
> result
;
351 for( size_t i
= 0; i
< nrVars(); ++i
)
352 result
.push_back( beliefV(i
) );
353 for( size_t I
= 0; I
< nrFactors(); ++I
)
354 result
.push_back( beliefF(I
) );
359 Factor
BP::belief( const VarSet
&ns
) const {
361 return belief( *(ns
.begin()) );
364 for( I
= 0; I
< nrFactors(); I
++ )
365 if( factor(I
).vars() >> ns
)
367 assert( I
!= nrFactors() );
368 return beliefF(I
).marginal(ns
);
373 Factor
BP::beliefF (size_t I
) const {
375 /* UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
377 Factor
prod( factor(I
) );
378 foreach( const Neighbor
&j
, nbF(I
) ) {
379 foreach( const Neighbor
&J
, nbV(j
) ) {
380 if( J
!= I
) // for all J in nb(j) \ I
381 prod
*= Factor( var(j
), newMessage(j
, J
.iter
) );
384 return prod
.normalized();
386 /* OPTIMIZED VERSION */
387 Prob
prod( factor(I
).p() );
388 if( props
.logdomain
)
391 foreach( const Neighbor
&j
, nbF(I
) ) {
393 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
394 const ind_t
& ind
= index(j
, _I
);
396 // prod_j will be the product of messages coming into j
397 Prob
prod_j( var(j
).states(), props
.logdomain
? 0.0 : 1.0 );
398 foreach( const Neighbor
&J
, nbV(j
) ) {
399 if( J
!= I
) { // for all J in nb(j) \ I
400 if( props
.logdomain
)
401 prod_j
+= newMessage( j
, J
.iter
);
403 prod_j
*= newMessage( j
, J
.iter
);
407 // multiply prod with prod_j
408 for( size_t r
= 0; r
< prod
.size(); ++r
) {
409 if( props
.logdomain
)
410 prod
[r
] += prod_j
[ind
[r
]];
412 prod
[r
] *= prod_j
[ind
[r
]];
416 if( props
.logdomain
) {
417 prod
-= prod
.maxVal();
421 Factor
result( factor(I
).vars(), prod
);
429 Real
BP::logZ() const {
431 for(size_t i
= 0; i
< nrVars(); ++i
)
432 sum
+= (1.0 - nbV(i
).size()) * beliefV(i
).entropy();
433 for( size_t I
= 0; I
< nrFactors(); ++I
)
434 sum
-= KL_dist( beliefF(I
), factor(I
) );
439 string
BP::identify() const {
440 return string(Name
) + printProperties();
444 void BP::init( const VarSet
&ns
) {
445 for( VarSet::const_iterator n
= ns
.begin(); n
!= ns
.end(); ++n
) {
446 size_t ni
= findVar( *n
);
447 foreach( const Neighbor
&I
, nbV( ni
) )
448 message( ni
, I
.iter
).fill( props
.logdomain
? 0.0 : 1.0 );
453 } // end of namespace dai