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-2010 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
15 #include <dai/exceptions.h>
24 const char *HAK::Name
= "HAK";
27 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
29 TFactor
<T
>& makePositive( TFactor
<T
> &f
, T epsilon
) {
30 for( size_t t
= 0; t
< f
.states(); t
++ )
31 if( (0 < f
[t
]) && (f
[t
] < epsilon
) )
36 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
38 TFactor
<T
>& makeZero( TFactor
<T
> &f
, T epsilon
) {
39 for( size_t t
= 0; t
< f
.states(); t
++ )
40 if( f
[t
] < epsilon
&& f
[t
] > -epsilon
)
46 void HAK::setProperties( const PropertySet
&opts
) {
47 DAI_ASSERT( opts
.hasKey("tol") );
48 DAI_ASSERT( opts
.hasKey("maxiter") );
49 DAI_ASSERT( opts
.hasKey("verbose") );
50 DAI_ASSERT( opts
.hasKey("doubleloop") );
51 DAI_ASSERT( opts
.hasKey("clusters") );
53 props
.tol
= opts
.getStringAs
<Real
>("tol");
54 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
55 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
56 props
.doubleloop
= opts
.getStringAs
<bool>("doubleloop");
57 props
.clusters
= opts
.getStringAs
<Properties::ClustersType
>("clusters");
59 if( opts
.hasKey("loopdepth") )
60 props
.loopdepth
= opts
.getStringAs
<size_t>("loopdepth");
62 DAI_ASSERT( props
.clusters
!= Properties::ClustersType::LOOP
);
63 if( opts
.hasKey("damping") )
64 props
.damping
= opts
.getStringAs
<Real
>("damping");
67 if( opts
.hasKey("init") )
68 props
.init
= opts
.getStringAs
<Properties::InitType
>("init");
70 props
.init
= Properties::InitType::UNIFORM
;
74 PropertySet
HAK::getProperties() const {
76 opts
.Set( "tol", props
.tol
);
77 opts
.Set( "maxiter", props
.maxiter
);
78 opts
.Set( "verbose", props
.verbose
);
79 opts
.Set( "doubleloop", props
.doubleloop
);
80 opts
.Set( "clusters", props
.clusters
);
81 opts
.Set( "init", props
.init
);
82 opts
.Set( "loopdepth", props
.loopdepth
);
83 opts
.Set( "damping", props
.damping
);
88 string
HAK::printProperties() const {
89 stringstream
s( stringstream::out
);
91 s
<< "tol=" << props
.tol
<< ",";
92 s
<< "maxiter=" << props
.maxiter
<< ",";
93 s
<< "verbose=" << props
.verbose
<< ",";
94 s
<< "doubleloop=" << props
.doubleloop
<< ",";
95 s
<< "clusters=" << props
.clusters
<< ",";
96 s
<< "init=" << props
.init
<< ",";
97 s
<< "loopdepth=" << props
.loopdepth
<< ",";
98 s
<< "damping=" << props
.damping
<< "]";
103 void HAK::construct() {
104 // Create outer beliefs
106 _Qa
.reserve(nrORs());
107 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
108 _Qa
.push_back( Factor( OR(alpha
) ) );
110 // Create inner beliefs
112 _Qb
.reserve(nrIRs());
113 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
114 _Qb
.push_back( Factor( IR(beta
) ) );
118 _muab
.reserve( nrORs() );
120 _muba
.reserve( nrORs() );
121 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
122 _muab
.push_back( vector
<Factor
>() );
123 _muba
.push_back( vector
<Factor
>() );
124 _muab
[alpha
].reserve( nbOR(alpha
).size() );
125 _muba
[alpha
].reserve( nbOR(alpha
).size() );
126 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
127 _muab
[alpha
].push_back( Factor( IR(beta
) ) );
128 _muba
[alpha
].push_back( Factor( IR(beta
) ) );
134 HAK::HAK( const RegionGraph
&rg
, const PropertySet
&opts
) : DAIAlgRG(rg
), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
135 setProperties( opts
);
141 void HAK::findLoopClusters( const FactorGraph
& fg
, std::set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
142 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
143 VarSet ind
= fg
.delta( fg
.findVar( *in
) );
144 if( (newcl
.size()) >= 2 && ind
.contains( root
) )
145 allcl
.insert( newcl
| *in
);
146 else if( length
> 1 )
147 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
152 HAK::HAK(const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
153 setProperties( opts
);
156 if( props
.clusters
== Properties::ClustersType::MIN
) {
158 constructCVM( fg
, cl
);
159 } else if( props
.clusters
== Properties::ClustersType::DELTA
) {
160 cl
.reserve( fg
.nrVars() );
161 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
162 cl
.push_back( fg
.Delta(i
) );
163 constructCVM( fg
, cl
);
164 } else if( props
.clusters
== Properties::ClustersType::LOOP
) {
167 for( size_t i0
= 0; i0
< fg
.nrVars(); i0
++ ) {
168 VarSet i0d
= fg
.delta(i0
);
169 if( props
.loopdepth
> 1 )
170 findLoopClusters( fg
, scl
, fg
.var(i0
), fg
.var(i0
), props
.loopdepth
- 1, fg
.delta(i0
) );
172 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
174 if( props
.verbose
>= 3 ) {
175 cerr
<< Name
<< " uses the following clusters: " << endl
;
176 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
177 cerr
<< *cli
<< endl
;
179 constructCVM( fg
, cl
);
180 } else if( props
.clusters
== Properties::ClustersType::BETHE
) {
181 // build outer regions (the cliques)
184 for( size_t c
= 0; c
< cl
.size(); c
++ )
185 nrEdges
+= cl
[c
].size();
187 // build inner regions (single variables)
189 irs
.reserve( fg
.nrVars() );
190 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
191 irs
.push_back( Region( fg
.var(i
), 1.0 ) );
193 // build edges (an outer and inner region are connected if the outer region contains the inner one)
194 // and calculate counting number for inner regions
195 vector
<std::pair
<size_t, size_t> > edges
;
196 edges
.reserve( nrEdges
);
197 for( size_t c
= 0; c
< cl
.size(); c
++ )
198 for( size_t i
= 0; i
< irs
.size(); i
++ )
199 if( cl
[c
] >> irs
[i
] ) {
200 edges
.push_back( make_pair( c
, i
) );
204 // build region graph
205 RegionGraph::construct( fg
, cl
, irs
, edges
);
207 DAI_THROW(UNKNOWN_ENUM_VALUE
);
211 if( props
.verbose
>= 3 )
212 cerr
<< Name
<< " regiongraph: " << *this << endl
;
216 string
HAK::identify() const {
217 return string(Name
) + printProperties();
221 void HAK::init( const VarSet
&ns
) {
222 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
223 if( _Qa
[alpha
].vars().intersects( ns
) ) {
224 if( props
.init
== Properties::InitType::UNIFORM
)
225 _Qa
[alpha
].setUniform();
227 _Qa
[alpha
].randomize();
228 _Qa
[alpha
] *= OR(alpha
);
229 _Qa
[alpha
].normalize();
232 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
233 if( IR(beta
).intersects( ns
) ) {
234 if( props
.init
== Properties::InitType::UNIFORM
)
235 _Qb
[beta
].fill( 1.0 );
237 _Qb
[beta
].randomize();
238 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
239 size_t _beta
= alpha
.dual
;
240 if( props
.init
== Properties::InitType::UNIFORM
) {
241 muab( alpha
, _beta
).fill( 1.0 );
242 muba( alpha
, _beta
).fill( 1.0 );
244 muab( alpha
, _beta
).randomize();
245 muba( alpha
, _beta
).randomize();
253 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
254 if( props
.init
== Properties::InitType::UNIFORM
)
255 _Qa
[alpha
].setUniform();
257 _Qa
[alpha
].randomize();
258 _Qa
[alpha
] *= OR(alpha
);
259 _Qa
[alpha
].normalize();
262 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
263 if( props
.init
== Properties::InitType::UNIFORM
)
264 _Qb
[beta
].setUniform();
266 _Qb
[beta
].randomize();
268 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
269 foreach( const Neighbor
&beta
, nbOR(alpha
) ) {
270 size_t _beta
= beta
.iter
;
271 if( props
.init
== Properties::InitType::UNIFORM
) {
272 muab( alpha
, _beta
).setUniform();
273 muba( alpha
, _beta
).setUniform();
275 muab( alpha
, _beta
).randomize();
276 muba( alpha
, _beta
).randomize();
283 if( props
.verbose
>= 1 )
284 cerr
<< "Starting " << identify() << "...";
285 if( props
.verbose
>= 3)
290 // Check whether counting numbers won't lead to problems
291 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
292 DAI_ASSERT( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
294 // Keep old beliefs to check convergence
295 vector
<Factor
> oldBeliefsV
;
296 oldBeliefsV
.reserve( nrVars() );
297 for( size_t i
= 0; i
< nrVars(); i
++ )
298 oldBeliefsV
.push_back( beliefV(i
) );
299 vector
<Factor
> oldBeliefsF
;
300 oldBeliefsF
.reserve( nrFactors() );
301 for( size_t I
= 0; I
< nrFactors(); I
++ )
302 oldBeliefsF
.push_back( beliefF(I
) );
304 // do several passes over the network until maximum number of iterations has
305 // been reached or until the maximum belief difference is smaller than tolerance
306 Real maxDiff
= INFINITY
;
307 for( _iters
= 0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
; _iters
++ ) {
308 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
309 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
310 size_t _beta
= alpha
.dual
;
311 muab( alpha
, _beta
) = _Qa
[alpha
].marginal(IR(beta
)) / muba(alpha
,_beta
);
312 /* TODO: INVESTIGATE THIS PROBLEM
314 * In some cases, the muab's can have very large entries because the muba's have very
315 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
316 * result in +inf; normalization then tries to calculate inf / inf which is NAN).
317 * A fix of this problem would consist in normalizing the messages muab.
318 * However, it is not obvious whether this is a real solution, because it has a
319 * negative performance impact and the NAN's seem to be a symptom of a fundamental
320 * numerical unstability.
322 muab(alpha
,_beta
).normalize();
326 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
327 size_t _beta
= alpha
.dual
;
328 Qb_new
*= muab(alpha
,_beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
332 if( Qb_new
.hasNaNs() ) {
333 // TODO: WHAT TO DO IN THIS CASE?
334 cerr
<< Name
<< "::doGBP: Qb_new has NaNs!" << endl
;
337 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
339 * _Qb[beta] = Qb_new.makeZero(1e-100);
342 if( props
.doubleloop
|| props
.damping
== 0.0 )
343 _Qb
[beta
] = Qb_new
; // no damping for double loop
345 _Qb
[beta
] = (Qb_new
^(1.0 - props
.damping
)) * (_Qb
[beta
]^props
.damping
);
347 foreach( const Neighbor
&alpha
, nbIR(beta
) ) {
348 size_t _beta
= alpha
.dual
;
349 muba(alpha
,_beta
) = _Qb
[beta
] / muab(alpha
,_beta
);
351 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
353 * muba(beta,*alpha).makePositive(1e-100);
357 Factor Qa_new
= OR(alpha
);
358 foreach( const Neighbor
&gamma
, nbOR(alpha
) )
359 Qa_new
*= muba(alpha
,gamma
.iter
);
360 Qa_new
^= (1.0 / OR(alpha
).c());
362 if( Qa_new
.hasNaNs() ) {
363 cerr
<< Name
<< "::doGBP: Qa_new has NaNs!" << endl
;
366 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
368 * _Qb[beta] = Qb_new.makeZero(1e-100);
371 if( props
.doubleloop
|| props
.damping
== 0.0 )
372 _Qa
[alpha
] = Qa_new
; // no damping for double loop
374 // FIXME: GEOMETRIC DAMPING IS SLOW!
375 _Qa
[alpha
] = (Qa_new
^(1.0 - props
.damping
)) * (_Qa
[alpha
]^props
.damping
);
379 // Calculate new single variable beliefs and compare with old ones
381 for( size_t i
= 0; i
< nrVars(); ++i
) {
382 Factor b
= beliefV(i
);
383 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], Prob::DISTLINF
) );
386 for( size_t I
= 0; I
< nrFactors(); ++I
) {
387 Factor b
= beliefF(I
);
388 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], Prob::DISTLINF
) );
392 if( props
.verbose
>= 3 )
393 cerr
<< Name
<< "::doGBP: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
396 if( maxDiff
> _maxdiff
)
399 if( props
.verbose
>= 1 ) {
400 if( maxDiff
> props
.tol
) {
401 if( props
.verbose
== 1 )
403 cerr
<< Name
<< "::doGBP: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
405 if( props
.verbose
>= 2 )
406 cerr
<< Name
<< "::doGBP: ";
407 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
415 Real
HAK::doDoubleLoop() {
416 if( props
.verbose
>= 1 )
417 cerr
<< "Starting " << identify() << "...";
418 if( props
.verbose
>= 3)
423 // Save original outer regions
424 vector
<FRegion
> org_ORs
= ORs
;
426 // Save original inner counting numbers and set negative counting numbers to zero
427 vector
<Real
> org_IR_cs( nrIRs(), 0.0 );
428 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
429 org_IR_cs
[beta
] = IR(beta
).c();
430 if( IR(beta
).c() < 0.0 )
434 // Keep old beliefs to check convergence
435 vector
<Factor
> oldBeliefsV
;
436 oldBeliefsV
.reserve( nrVars() );
437 for( size_t i
= 0; i
< nrVars(); i
++ )
438 oldBeliefsV
.push_back( beliefV(i
) );
439 vector
<Factor
> oldBeliefsF
;
440 oldBeliefsF
.reserve( nrFactors() );
441 for( size_t I
= 0; I
< nrFactors(); I
++ )
442 oldBeliefsF
.push_back( beliefF(I
) );
444 size_t outer_maxiter
= props
.maxiter
;
445 Real outer_tol
= props
.tol
;
446 size_t outer_verbose
= props
.verbose
;
447 Real org_maxdiff
= _maxdiff
;
449 // Set parameters for inner loop
451 props
.verbose
= outer_verbose
? outer_verbose
- 1 : 0;
453 size_t outer_iter
= 0;
454 size_t total_iter
= 0;
455 Real maxDiff
= INFINITY
;
456 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& maxDiff
> outer_tol
; outer_iter
++ ) {
457 // Calculate new outer regions
458 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
459 OR(alpha
) = org_ORs
[alpha
];
460 foreach( const Neighbor
&beta
, nbOR(alpha
) )
461 OR(alpha
) *= _Qb
[beta
] ^ ((IR(beta
).c() - org_IR_cs
[beta
]) / nbIR(beta
).size());
465 if( isnan( doGBP() ) )
468 // Calculate new single variable beliefs and compare with old ones
470 for( size_t i
= 0; i
< nrVars(); ++i
) {
471 Factor b
= beliefV(i
);
472 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], Prob::DISTLINF
) );
475 for( size_t I
= 0; I
< nrFactors(); ++I
) {
476 Factor b
= beliefF(I
);
477 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], Prob::DISTLINF
) );
481 total_iter
+= Iterations();
483 if( props
.verbose
>= 3 )
484 cerr
<< Name
<< "::doDoubleLoop: maxdiff " << maxDiff
<< " after " << total_iter
<< " passes" << endl
;
487 // restore _maxiter, _verbose and _maxdiff
488 props
.maxiter
= outer_maxiter
;
489 props
.verbose
= outer_verbose
;
490 _maxdiff
= org_maxdiff
;
493 if( maxDiff
> _maxdiff
)
496 // Restore original outer regions
499 // Restore original inner counting numbers
500 for( size_t beta
= 0; beta
< nrIRs(); ++beta
)
501 IR(beta
).c() = org_IR_cs
[beta
];
503 if( props
.verbose
>= 1 ) {
504 if( maxDiff
> props
.tol
) {
505 if( props
.verbose
== 1 )
507 cerr
<< Name
<< "::doDoubleLoop: WARNING: not converged within " << outer_maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
509 if( props
.verbose
>= 3 )
510 cerr
<< Name
<< "::doDoubleLoop: ";
511 cerr
<< "converged in " << total_iter
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
520 if( props
.doubleloop
)
521 return doDoubleLoop();
527 Factor
HAK::belief( const VarSet
&ns
) const {
528 vector
<Factor
>::const_iterator beta
;
529 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
530 if( beta
->vars() >> ns
)
532 if( beta
!= _Qb
.end() )
533 return( beta
->marginal(ns
) );
535 vector
<Factor
>::const_iterator alpha
;
536 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
537 if( alpha
->vars() >> ns
)
539 if( alpha
== _Qa
.end() )
540 DAI_THROW(BELIEF_NOT_AVAILABLE
);
541 return( alpha
->marginal(ns
) );
546 vector
<Factor
> HAK::beliefs() const {
547 vector
<Factor
> result
;
548 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
549 result
.push_back( Qb(beta
) );
550 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
551 result
.push_back( Qa(alpha
) );
556 Real
HAK::logZ() const {
558 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
559 s
+= IR(beta
).c() * Qb(beta
).entropy();
560 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
561 s
+= OR(alpha
).c() * Qa(alpha
).entropy();
562 s
+= (OR(alpha
).log(true) * Qa(alpha
)).sum();
568 } // end of namespace dai