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>
16 #include <dai/exceptions.h>
25 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
27 TFactor
<T
>& makePositive( TFactor
<T
> &f
, T epsilon
) {
28 for( size_t t
= 0; t
< f
.states(); t
++ )
29 if( (0 < f
[t
]) && (f
[t
] < epsilon
) )
34 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
36 TFactor
<T
>& makeZero( TFactor
<T
> &f
, T epsilon
) {
37 for( size_t t
= 0; t
< f
.states(); t
++ )
38 if( f
[t
] < epsilon
&& f
[t
] > -epsilon
)
44 void HAK::setProperties( const PropertySet
&opts
) {
45 DAI_ASSERT( opts
.hasKey("tol") );
46 DAI_ASSERT( opts
.hasKey("doubleloop") );
47 DAI_ASSERT( opts
.hasKey("clusters") );
49 props
.tol
= opts
.getStringAs
<Real
>("tol");
50 props
.doubleloop
= opts
.getStringAs
<bool>("doubleloop");
51 props
.clusters
= opts
.getStringAs
<Properties::ClustersType
>("clusters");
53 if( opts
.hasKey("maxiter") )
54 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
56 props
.maxiter
= 10000;
57 if( opts
.hasKey("maxtime") )
58 props
.maxtime
= opts
.getStringAs
<Real
>("maxtime");
60 props
.maxtime
= INFINITY
;
61 if( opts
.hasKey("verbose") )
62 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
65 if( opts
.hasKey("loopdepth") )
66 props
.loopdepth
= opts
.getStringAs
<size_t>("loopdepth");
68 DAI_ASSERT( props
.clusters
!= Properties::ClustersType::LOOP
);
69 if( opts
.hasKey("damping") )
70 props
.damping
= opts
.getStringAs
<Real
>("damping");
73 if( opts
.hasKey("init") )
74 props
.init
= opts
.getStringAs
<Properties::InitType
>("init");
76 props
.init
= Properties::InitType::UNIFORM
;
80 PropertySet
HAK::getProperties() const {
82 opts
.set( "tol", props
.tol
);
83 opts
.set( "maxiter", props
.maxiter
);
84 opts
.set( "maxtime", props
.maxtime
);
85 opts
.set( "verbose", props
.verbose
);
86 opts
.set( "doubleloop", props
.doubleloop
);
87 opts
.set( "clusters", props
.clusters
);
88 opts
.set( "init", props
.init
);
89 opts
.set( "loopdepth", props
.loopdepth
);
90 opts
.set( "damping", props
.damping
);
95 string
HAK::printProperties() const {
96 stringstream
s( stringstream::out
);
98 s
<< "tol=" << props
.tol
<< ",";
99 s
<< "maxiter=" << props
.maxiter
<< ",";
100 s
<< "maxtime=" << props
.maxtime
<< ",";
101 s
<< "verbose=" << props
.verbose
<< ",";
102 s
<< "doubleloop=" << props
.doubleloop
<< ",";
103 s
<< "clusters=" << props
.clusters
<< ",";
104 s
<< "init=" << props
.init
<< ",";
105 s
<< "loopdepth=" << props
.loopdepth
<< ",";
106 s
<< "damping=" << props
.damping
<< "]";
111 void HAK::construct() {
112 // Create outer beliefs
113 if( props
.verbose
>= 3 )
114 cerr
<< "Constructing outer beliefs" << endl
;
116 _Qa
.reserve(nrORs());
117 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
118 _Qa
.push_back( Factor( OR(alpha
) ) );
120 // Create inner beliefs
121 if( props
.verbose
>= 3 )
122 cerr
<< "Constructing inner beliefs" << endl
;
124 _Qb
.reserve(nrIRs());
125 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
126 _Qb
.push_back( Factor( IR(beta
) ) );
129 if( props
.verbose
>= 3 )
130 cerr
<< "Constructing messages" << endl
;
132 _muab
.reserve( nrORs() );
134 _muba
.reserve( nrORs() );
135 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
136 _muab
.push_back( vector
<Factor
>() );
137 _muba
.push_back( vector
<Factor
>() );
138 _muab
[alpha
].reserve( nbOR(alpha
).size() );
139 _muba
[alpha
].reserve( nbOR(alpha
).size() );
140 bforeach( const Neighbor
&beta
, nbOR(alpha
) ) {
141 _muab
[alpha
].push_back( Factor( IR(beta
) ) );
142 _muba
[alpha
].push_back( Factor( IR(beta
) ) );
148 HAK::HAK( const RegionGraph
&rg
, const PropertySet
&opts
) : DAIAlgRG(rg
), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
149 setProperties( opts
);
155 void HAK::findLoopClusters( const FactorGraph
& fg
, std::set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
156 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
157 VarSet ind
= fg
.delta( fg
.findVar( *in
) );
158 if( (newcl
.size()) >= 2 && ind
.contains( root
) )
159 allcl
.insert( newcl
| *in
);
160 else if( length
> 1 )
161 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
166 HAK::HAK(const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
167 setProperties( opts
);
169 if( props
.verbose
>= 3 )
170 cerr
<< "Constructing clusters" << endl
;
173 if( props
.clusters
== Properties::ClustersType::MIN
) {
174 cl
= fg
.maximalFactorDomains();
175 constructCVM( fg
, cl
);
176 } else if( props
.clusters
== Properties::ClustersType::DELTA
) {
177 cl
.reserve( fg
.nrVars() );
178 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
179 cl
.push_back( fg
.Delta(i
) );
180 constructCVM( fg
, cl
);
181 } else if( props
.clusters
== Properties::ClustersType::LOOP
) {
182 cl
= fg
.maximalFactorDomains();
184 if( props
.verbose
>= 2 )
185 cerr
<< "Searching loops...";
186 for( size_t i0
= 0; i0
< fg
.nrVars(); i0
++ ) {
187 VarSet i0d
= fg
.delta(i0
);
188 if( props
.loopdepth
> 1 )
189 findLoopClusters( fg
, scl
, fg
.var(i0
), fg
.var(i0
), props
.loopdepth
- 1, fg
.delta(i0
) );
191 if( props
.verbose
>= 2 )
192 cerr
<< "done" << endl
;
193 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
195 if( props
.verbose
>= 3 ) {
196 cerr
<< name() << " uses the following clusters: " << endl
;
197 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
198 cerr
<< *cli
<< endl
;
200 constructCVM( fg
, cl
);
201 } else if( props
.clusters
== Properties::ClustersType::BETHE
) {
202 // Copy factor graph structure
203 if( props
.verbose
>= 3 )
204 cerr
<< "Copying factor graph" << endl
;
205 FactorGraph::operator=( fg
);
207 // Construct inner regions (single variables)
208 if( props
.verbose
>= 3 )
209 cerr
<< "Constructing inner regions" << endl
;
210 _IRs
.reserve( fg
.nrVars() );
211 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
212 _IRs
.push_back( Region( fg
.var(i
), 1.0 ) );
215 if( props
.verbose
>= 3 )
216 cerr
<< "Constructing graph" << endl
;
217 _G
= BipartiteGraph( 0, nrIRs() );
219 // Construct outer regions:
220 // maximal factors become new outer regions
221 // non-maximal factors are assigned an outer region that contains them
222 if( props
.verbose
>= 3 )
223 cerr
<< "Construct outer regions" << endl
;
224 _fac2OR
.reserve( nrFactors() );
225 queue
<pair
<size_t, size_t> > todo
;
226 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
227 size_t J
= fg
.maximalFactor( I
);
229 // I is maximal; add it to the outer regions
230 _fac2OR
.push_back( nrORs() );
231 // Construct outer region (with counting number 1.0)
232 _ORs
.push_back( FRegion( fg
.factor(I
), 1.0 ) );
233 // Add node and edges to graph
234 SmallSet
<size_t> irs
= fg
.bipGraph().nb2Set( I
);
235 _G
.addNode1( irs
.begin(), irs
.end(), irs
.size() );
237 // J is larger and has already been assigned to an outer region
238 // so I should belong to the same outer region as J
239 _fac2OR
.push_back( _fac2OR
[J
] );
240 _ORs
[_fac2OR
[J
]] *= fg
.factor(I
);
242 // J is larger but has not yet been assigned to an outer region
243 // we handle this case later
244 _fac2OR
.push_back( -1 );
245 todo
.push( make_pair( I
, J
) );
248 // finish the construction
249 while( !todo
.empty() ) {
250 size_t I
= todo
.front().first
;
251 size_t J
= todo
.front().second
;
253 _fac2OR
[I
] = _fac2OR
[J
];
254 _ORs
[_fac2OR
[J
]] *= fg
.factor(I
);
257 // Calculate inner regions' counting numbers
258 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
259 _IRs
[beta
].c() = 1.0 - _G
.nb2(beta
).size();
261 DAI_THROW(UNKNOWN_ENUM_VALUE
);
265 if( props
.verbose
>= 3 )
266 cerr
<< name() << " regiongraph: " << *this << endl
;
270 void HAK::init( const VarSet
&ns
) {
271 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
272 if( _Qa
[alpha
].vars().intersects( ns
) ) {
273 if( props
.init
== Properties::InitType::UNIFORM
)
274 _Qa
[alpha
].setUniform();
276 _Qa
[alpha
].randomize();
277 _Qa
[alpha
] *= OR(alpha
);
278 _Qa
[alpha
].normalize();
281 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
282 if( IR(beta
).intersects( ns
) ) {
283 if( props
.init
== Properties::InitType::UNIFORM
)
284 _Qb
[beta
].fill( 1.0 );
286 _Qb
[beta
].randomize();
287 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
288 size_t _beta
= alpha
.dual
;
289 if( props
.init
== Properties::InitType::UNIFORM
) {
290 muab( alpha
, _beta
).fill( 1.0 );
291 muba( alpha
, _beta
).fill( 1.0 );
293 muab( alpha
, _beta
).randomize();
294 muba( alpha
, _beta
).randomize();
302 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
303 if( props
.init
== Properties::InitType::UNIFORM
)
304 _Qa
[alpha
].setUniform();
306 _Qa
[alpha
].randomize();
307 _Qa
[alpha
] *= OR(alpha
);
308 _Qa
[alpha
].normalize();
311 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
312 if( props
.init
== Properties::InitType::UNIFORM
)
313 _Qb
[beta
].setUniform();
315 _Qb
[beta
].randomize();
317 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
318 bforeach( const Neighbor
&beta
, nbOR(alpha
) ) {
319 size_t _beta
= beta
.iter
;
320 if( props
.init
== Properties::InitType::UNIFORM
) {
321 muab( alpha
, _beta
).setUniform();
322 muba( alpha
, _beta
).setUniform();
324 muab( alpha
, _beta
).randomize();
325 muba( alpha
, _beta
).randomize();
332 if( props
.verbose
>= 1 )
333 cerr
<< "Starting " << identify() << "...";
334 if( props
.verbose
>= 3)
339 // Check whether counting numbers won't lead to problems
340 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
341 DAI_ASSERT( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
343 // Keep old beliefs to check convergence
344 vector
<Factor
> oldBeliefsV
;
345 oldBeliefsV
.reserve( nrVars() );
346 for( size_t i
= 0; i
< nrVars(); i
++ )
347 oldBeliefsV
.push_back( beliefV(i
) );
348 vector
<Factor
> oldBeliefsF
;
349 oldBeliefsF
.reserve( nrFactors() );
350 for( size_t I
= 0; I
< nrFactors(); I
++ )
351 oldBeliefsF
.push_back( beliefF(I
) );
353 // do several passes over the network until maximum number of iterations has
354 // been reached or until the maximum belief difference is smaller than tolerance
355 Real maxDiff
= INFINITY
;
356 for( _iters
= 0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
; _iters
++ ) {
357 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
358 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
359 size_t _beta
= alpha
.dual
;
360 muab( alpha
, _beta
) = _Qa
[alpha
].marginal(IR(beta
)) / muba(alpha
,_beta
);
361 /* TODO: INVESTIGATE THIS PROBLEM
363 * In some cases, the muab's can have very large entries because the muba's have very
364 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
365 * result in +inf; normalization then tries to calculate inf / inf which is NAN).
366 * A fix of this problem would consist in normalizing the messages muab.
367 * However, it is not obvious whether this is a real solution, because it has a
368 * negative performance impact and the NAN's seem to be a symptom of a fundamental
369 * numerical unstability.
371 muab(alpha
,_beta
).normalize();
375 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
376 size_t _beta
= alpha
.dual
;
377 Qb_new
*= muab(alpha
,_beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
381 if( Qb_new
.hasNaNs() ) {
382 // TODO: WHAT TO DO IN THIS CASE?
383 cerr
<< name() << "::doGBP: Qb_new has NaNs!" << endl
;
386 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
388 * _Qb[beta] = Qb_new.makeZero(1e-100);
391 if( props
.doubleloop
|| props
.damping
== 0.0 )
392 _Qb
[beta
] = Qb_new
; // no damping for double loop
394 _Qb
[beta
] = (Qb_new
^(1.0 - props
.damping
)) * (_Qb
[beta
]^props
.damping
);
396 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
397 size_t _beta
= alpha
.dual
;
398 muba(alpha
,_beta
) = _Qb
[beta
] / muab(alpha
,_beta
);
400 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
402 * muba(beta,*alpha).makePositive(1e-100);
406 Factor Qa_new
= OR(alpha
);
407 bforeach( const Neighbor
&gamma
, nbOR(alpha
) )
408 Qa_new
*= muba(alpha
,gamma
.iter
);
409 Qa_new
^= (1.0 / OR(alpha
).c());
411 if( Qa_new
.hasNaNs() ) {
412 cerr
<< name() << "::doGBP: Qa_new has NaNs!" << endl
;
415 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
417 * _Qb[beta] = Qb_new.makeZero(1e-100);
420 if( props
.doubleloop
|| props
.damping
== 0.0 )
421 _Qa
[alpha
] = Qa_new
; // no damping for double loop
423 // FIXME: GEOMETRIC DAMPING IS SLOW!
424 _Qa
[alpha
] = (Qa_new
^(1.0 - props
.damping
)) * (_Qa
[alpha
]^props
.damping
);
428 // Calculate new single variable beliefs and compare with old ones
430 for( size_t i
= 0; i
< nrVars(); ++i
) {
431 Factor b
= beliefV(i
);
432 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], DISTLINF
) );
435 for( size_t I
= 0; I
< nrFactors(); ++I
) {
436 Factor b
= beliefF(I
);
437 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], DISTLINF
) );
441 if( props
.verbose
>= 3 )
442 cerr
<< name() << "::doGBP: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
445 if( maxDiff
> _maxdiff
)
448 if( props
.verbose
>= 1 ) {
449 if( maxDiff
> props
.tol
) {
450 if( props
.verbose
== 1 )
452 cerr
<< name() << "::doGBP: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
454 if( props
.verbose
>= 2 )
455 cerr
<< name() << "::doGBP: ";
456 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
464 Real
HAK::doDoubleLoop() {
465 if( props
.verbose
>= 1 )
466 cerr
<< "Starting " << identify() << "...";
467 if( props
.verbose
>= 3)
472 // Save original outer regions
473 vector
<FRegion
> org_ORs
= _ORs
;
475 // Save original inner counting numbers and set negative counting numbers to zero
476 vector
<Real
> org_IR_cs( nrIRs(), 0.0 );
477 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
478 org_IR_cs
[beta
] = IR(beta
).c();
479 if( IR(beta
).c() < 0.0 )
483 // Keep old beliefs to check convergence
484 vector
<Factor
> oldBeliefsV
;
485 oldBeliefsV
.reserve( nrVars() );
486 for( size_t i
= 0; i
< nrVars(); i
++ )
487 oldBeliefsV
.push_back( beliefV(i
) );
488 vector
<Factor
> oldBeliefsF
;
489 oldBeliefsF
.reserve( nrFactors() );
490 for( size_t I
= 0; I
< nrFactors(); I
++ )
491 oldBeliefsF
.push_back( beliefF(I
) );
493 size_t outer_maxiter
= props
.maxiter
;
494 Real outer_tol
= props
.tol
;
495 size_t outer_verbose
= props
.verbose
;
496 Real org_maxdiff
= _maxdiff
;
498 // Set parameters for inner loop
500 props
.verbose
= outer_verbose
? outer_verbose
- 1 : 0;
502 size_t outer_iter
= 0;
503 size_t total_iter
= 0;
504 Real maxDiff
= INFINITY
;
505 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& maxDiff
> outer_tol
&& (toc() - tic
) < props
.maxtime
; outer_iter
++ ) {
506 // Calculate new outer regions
507 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
508 OR(alpha
) = org_ORs
[alpha
];
509 bforeach( const Neighbor
&beta
, nbOR(alpha
) )
510 OR(alpha
) *= _Qb
[beta
] ^ ((IR(beta
).c() - org_IR_cs
[beta
]) / nbIR(beta
).size());
514 if( isnan( doGBP() ) )
517 // Calculate new single variable beliefs and compare with old ones
519 for( size_t i
= 0; i
< nrVars(); ++i
) {
520 Factor b
= beliefV(i
);
521 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], DISTLINF
) );
524 for( size_t I
= 0; I
< nrFactors(); ++I
) {
525 Factor b
= beliefF(I
);
526 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], DISTLINF
) );
530 total_iter
+= Iterations();
532 if( props
.verbose
>= 3 )
533 cerr
<< name() << "::doDoubleLoop: maxdiff " << maxDiff
<< " after " << total_iter
<< " passes" << endl
;
536 // restore _maxiter, _verbose and _maxdiff
537 props
.maxiter
= outer_maxiter
;
538 props
.verbose
= outer_verbose
;
539 _maxdiff
= org_maxdiff
;
542 if( maxDiff
> _maxdiff
)
545 // Restore original outer regions
548 // Restore original inner counting numbers
549 for( size_t beta
= 0; beta
< nrIRs(); ++beta
)
550 IR(beta
).c() = org_IR_cs
[beta
];
552 if( props
.verbose
>= 1 ) {
553 if( maxDiff
> props
.tol
) {
554 if( props
.verbose
== 1 )
556 cerr
<< name() << "::doDoubleLoop: WARNING: not converged after " << total_iter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
558 if( props
.verbose
>= 3 )
559 cerr
<< name() << "::doDoubleLoop: ";
560 cerr
<< "converged in " << total_iter
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
569 if( props
.doubleloop
)
570 return doDoubleLoop();
576 Factor
HAK::belief( const VarSet
&ns
) const {
577 vector
<Factor
>::const_iterator beta
;
578 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
579 if( beta
->vars() >> ns
)
581 if( beta
!= _Qb
.end() )
582 return( beta
->marginal(ns
) );
584 vector
<Factor
>::const_iterator alpha
;
585 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
586 if( alpha
->vars() >> ns
)
588 if( alpha
== _Qa
.end() )
589 DAI_THROW(BELIEF_NOT_AVAILABLE
);
590 return( alpha
->marginal(ns
) );
595 vector
<Factor
> HAK::beliefs() const {
596 vector
<Factor
> result
;
597 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
598 result
.push_back( Qb(beta
) );
599 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
600 result
.push_back( Qa(alpha
) );
605 Real
HAK::logZ() const {
607 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
608 s
+= IR(beta
).c() * Qb(beta
).entropy();
609 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
610 s
+= OR(alpha
).c() * Qa(alpha
).entropy();
611 s
+= (OR(alpha
).log(true) * Qa(alpha
)).sum();
617 } // end of namespace dai