4b79981471794340b2b28d500e1dac4d60816df2
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.
12 #include <dai/exceptions.h>
21 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
23 TFactor
<T
>& makePositive( TFactor
<T
> &f
, T epsilon
) {
24 for( size_t t
= 0; t
< f
.states(); t
++ )
25 if( (0 < f
[t
]) && (f
[t
] < epsilon
) )
30 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
32 TFactor
<T
>& makeZero( TFactor
<T
> &f
, T epsilon
) {
33 for( size_t t
= 0; t
< f
.states(); t
++ )
34 if( f
[t
] < epsilon
&& f
[t
] > -epsilon
)
40 void HAK::setProperties( const PropertySet
&opts
) {
41 DAI_ASSERT( opts
.hasKey("tol") );
42 DAI_ASSERT( opts
.hasKey("doubleloop") );
43 DAI_ASSERT( opts
.hasKey("clusters") );
45 props
.tol
= opts
.getStringAs
<Real
>("tol");
46 props
.doubleloop
= opts
.getStringAs
<bool>("doubleloop");
47 props
.clusters
= opts
.getStringAs
<Properties::ClustersType
>("clusters");
49 if( opts
.hasKey("maxiter") )
50 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
52 props
.maxiter
= 10000;
53 if( opts
.hasKey("maxtime") )
54 props
.maxtime
= opts
.getStringAs
<Real
>("maxtime");
56 props
.maxtime
= INFINITY
;
57 if( opts
.hasKey("verbose") )
58 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
61 if( opts
.hasKey("loopdepth") )
62 props
.loopdepth
= opts
.getStringAs
<size_t>("loopdepth");
64 DAI_ASSERT( props
.clusters
!= Properties::ClustersType::LOOP
);
65 if( opts
.hasKey("damping") )
66 props
.damping
= opts
.getStringAs
<Real
>("damping");
69 if( opts
.hasKey("init") )
70 props
.init
= opts
.getStringAs
<Properties::InitType
>("init");
72 props
.init
= Properties::InitType::UNIFORM
;
76 PropertySet
HAK::getProperties() const {
78 opts
.set( "tol", props
.tol
);
79 opts
.set( "maxiter", props
.maxiter
);
80 opts
.set( "maxtime", props
.maxtime
);
81 opts
.set( "verbose", props
.verbose
);
82 opts
.set( "doubleloop", props
.doubleloop
);
83 opts
.set( "clusters", props
.clusters
);
84 opts
.set( "init", props
.init
);
85 opts
.set( "loopdepth", props
.loopdepth
);
86 opts
.set( "damping", props
.damping
);
91 string
HAK::printProperties() const {
92 stringstream
s( stringstream::out
);
94 s
<< "tol=" << props
.tol
<< ",";
95 s
<< "maxiter=" << props
.maxiter
<< ",";
96 s
<< "maxtime=" << props
.maxtime
<< ",";
97 s
<< "verbose=" << props
.verbose
<< ",";
98 s
<< "doubleloop=" << props
.doubleloop
<< ",";
99 s
<< "clusters=" << props
.clusters
<< ",";
100 s
<< "init=" << props
.init
<< ",";
101 s
<< "loopdepth=" << props
.loopdepth
<< ",";
102 s
<< "damping=" << props
.damping
<< "]";
107 void HAK::construct() {
108 // Create outer beliefs
109 if( props
.verbose
>= 3 )
110 cerr
<< "Constructing outer beliefs" << endl
;
112 _Qa
.reserve(nrORs());
113 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
114 _Qa
.push_back( Factor( OR(alpha
) ) );
116 // Create inner beliefs
117 if( props
.verbose
>= 3 )
118 cerr
<< "Constructing inner beliefs" << endl
;
120 _Qb
.reserve(nrIRs());
121 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
122 _Qb
.push_back( Factor( IR(beta
) ) );
125 if( props
.verbose
>= 3 )
126 cerr
<< "Constructing messages" << endl
;
128 _muab
.reserve( nrORs() );
130 _muba
.reserve( nrORs() );
131 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
132 _muab
.push_back( vector
<Factor
>() );
133 _muba
.push_back( vector
<Factor
>() );
134 _muab
[alpha
].reserve( nbOR(alpha
).size() );
135 _muba
[alpha
].reserve( nbOR(alpha
).size() );
136 bforeach( const Neighbor
&beta
, nbOR(alpha
) ) {
137 _muab
[alpha
].push_back( Factor( IR(beta
) ) );
138 _muba
[alpha
].push_back( Factor( IR(beta
) ) );
144 HAK::HAK( const RegionGraph
&rg
, const PropertySet
&opts
) : DAIAlgRG(rg
), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
145 setProperties( opts
);
151 void HAK::findLoopClusters( const FactorGraph
& fg
, std::set
<VarSet
> &allcl
, VarSet newcl
, const Var
& root
, size_t length
, VarSet vars
) {
152 for( VarSet::const_iterator in
= vars
.begin(); in
!= vars
.end(); in
++ ) {
153 VarSet ind
= fg
.delta( fg
.findVar( *in
) );
154 if( (newcl
.size()) >= 2 && ind
.contains( root
) )
155 allcl
.insert( newcl
| *in
);
156 else if( length
> 1 )
157 findLoopClusters( fg
, allcl
, newcl
| *in
, root
, length
- 1, ind
/ newcl
);
162 HAK::HAK(const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
163 setProperties( opts
);
165 if( props
.verbose
>= 3 )
166 cerr
<< "Constructing clusters" << endl
;
169 if( props
.clusters
== Properties::ClustersType::MIN
) {
170 cl
= fg
.maximalFactorDomains();
171 constructCVM( fg
, cl
);
172 } else if( props
.clusters
== Properties::ClustersType::DELTA
) {
173 cl
.reserve( fg
.nrVars() );
174 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
175 cl
.push_back( fg
.Delta(i
) );
176 constructCVM( fg
, cl
);
177 } else if( props
.clusters
== Properties::ClustersType::LOOP
) {
178 cl
= fg
.maximalFactorDomains();
180 if( props
.verbose
>= 2 )
181 cerr
<< "Searching loops...";
182 for( size_t i0
= 0; i0
< fg
.nrVars(); i0
++ ) {
183 VarSet i0d
= fg
.delta(i0
);
184 if( props
.loopdepth
> 1 )
185 findLoopClusters( fg
, scl
, fg
.var(i0
), fg
.var(i0
), props
.loopdepth
- 1, fg
.delta(i0
) );
187 if( props
.verbose
>= 2 )
188 cerr
<< "done" << endl
;
189 for( set
<VarSet
>::const_iterator c
= scl
.begin(); c
!= scl
.end(); c
++ )
191 if( props
.verbose
>= 3 ) {
192 cerr
<< name() << " uses the following clusters: " << endl
;
193 for( vector
<VarSet
>::const_iterator cli
= cl
.begin(); cli
!= cl
.end(); cli
++ )
194 cerr
<< *cli
<< endl
;
196 constructCVM( fg
, cl
);
197 } else if( props
.clusters
== Properties::ClustersType::BETHE
) {
198 // Copy factor graph structure
199 if( props
.verbose
>= 3 )
200 cerr
<< "Copying factor graph" << endl
;
201 FactorGraph::operator=( fg
);
203 // Construct inner regions (single variables)
204 if( props
.verbose
>= 3 )
205 cerr
<< "Constructing inner regions" << endl
;
206 _IRs
.reserve( fg
.nrVars() );
207 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
208 _IRs
.push_back( Region( fg
.var(i
), 1.0 ) );
211 if( props
.verbose
>= 3 )
212 cerr
<< "Constructing graph" << endl
;
213 _G
= BipartiteGraph( 0, nrIRs() );
215 // Construct outer regions:
216 // maximal factors become new outer regions
217 // non-maximal factors are assigned an outer region that contains them
218 if( props
.verbose
>= 3 )
219 cerr
<< "Construct outer regions" << endl
;
220 _fac2OR
.reserve( nrFactors() );
221 queue
<pair
<size_t, size_t> > todo
;
222 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
223 size_t J
= fg
.maximalFactor( I
);
225 // I is maximal; add it to the outer regions
226 _fac2OR
.push_back( nrORs() );
227 // Construct outer region (with counting number 1.0)
228 _ORs
.push_back( FRegion( fg
.factor(I
), 1.0 ) );
229 // Add node and edges to graph
230 SmallSet
<size_t> irs
= fg
.bipGraph().nb2Set( I
);
231 _G
.addNode1( irs
.begin(), irs
.end(), irs
.size() );
233 // J is larger and has already been assigned to an outer region
234 // so I should belong to the same outer region as J
235 _fac2OR
.push_back( _fac2OR
[J
] );
236 _ORs
[_fac2OR
[J
]] *= fg
.factor(I
);
238 // J is larger but has not yet been assigned to an outer region
239 // we handle this case later
240 _fac2OR
.push_back( -1 );
241 todo
.push( make_pair( I
, J
) );
244 // finish the construction
245 while( !todo
.empty() ) {
246 size_t I
= todo
.front().first
;
247 size_t J
= todo
.front().second
;
249 _fac2OR
[I
] = _fac2OR
[J
];
250 _ORs
[_fac2OR
[J
]] *= fg
.factor(I
);
253 // Calculate inner regions' counting numbers
254 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
255 _IRs
[beta
].c() = 1.0 - _G
.nb2(beta
).size();
257 DAI_THROW(UNKNOWN_ENUM_VALUE
);
261 if( props
.verbose
>= 3 )
262 cerr
<< name() << " regiongraph: " << *this << endl
;
266 void HAK::init( const VarSet
&ns
) {
267 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
268 if( _Qa
[alpha
].vars().intersects( ns
) ) {
269 if( props
.init
== Properties::InitType::UNIFORM
)
270 _Qa
[alpha
].setUniform();
272 _Qa
[alpha
].randomize();
273 _Qa
[alpha
] *= OR(alpha
);
274 _Qa
[alpha
].normalize();
277 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
278 if( IR(beta
).intersects( ns
) ) {
279 if( props
.init
== Properties::InitType::UNIFORM
)
280 _Qb
[beta
].fill( 1.0 );
282 _Qb
[beta
].randomize();
283 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
284 size_t _beta
= alpha
.dual
;
285 if( props
.init
== Properties::InitType::UNIFORM
) {
286 muab( alpha
, _beta
).fill( 1.0 );
287 muba( alpha
, _beta
).fill( 1.0 );
289 muab( alpha
, _beta
).randomize();
290 muba( alpha
, _beta
).randomize();
298 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
299 if( props
.init
== Properties::InitType::UNIFORM
)
300 _Qa
[alpha
].setUniform();
302 _Qa
[alpha
].randomize();
303 _Qa
[alpha
] *= OR(alpha
);
304 _Qa
[alpha
].normalize();
307 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
308 if( props
.init
== Properties::InitType::UNIFORM
)
309 _Qb
[beta
].setUniform();
311 _Qb
[beta
].randomize();
313 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
314 bforeach( const Neighbor
&beta
, nbOR(alpha
) ) {
315 size_t _beta
= beta
.iter
;
316 if( props
.init
== Properties::InitType::UNIFORM
) {
317 muab( alpha
, _beta
).setUniform();
318 muba( alpha
, _beta
).setUniform();
320 muab( alpha
, _beta
).randomize();
321 muba( alpha
, _beta
).randomize();
328 if( props
.verbose
>= 1 )
329 cerr
<< "Starting " << identify() << "...";
330 if( props
.verbose
>= 3)
335 // Check whether counting numbers won't lead to problems
336 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
337 DAI_ASSERT( nbIR(beta
).size() + IR(beta
).c() != 0.0 );
339 // Keep old beliefs to check convergence
340 vector
<Factor
> oldBeliefsV
;
341 oldBeliefsV
.reserve( nrVars() );
342 for( size_t i
= 0; i
< nrVars(); i
++ )
343 oldBeliefsV
.push_back( beliefV(i
) );
344 vector
<Factor
> oldBeliefsF
;
345 oldBeliefsF
.reserve( nrFactors() );
346 for( size_t I
= 0; I
< nrFactors(); I
++ )
347 oldBeliefsF
.push_back( beliefF(I
) );
349 // do several passes over the network until maximum number of iterations has
350 // been reached or until the maximum belief difference is smaller than tolerance
351 Real maxDiff
= INFINITY
;
352 for( _iters
= 0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
; _iters
++ ) {
353 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
354 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
355 size_t _beta
= alpha
.dual
;
356 muab( alpha
, _beta
) = _Qa
[alpha
].marginal(IR(beta
)) / muba(alpha
,_beta
);
357 /* TODO: INVESTIGATE THIS PROBLEM
359 * In some cases, the muab's can have very large entries because the muba's have very
360 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
361 * result in +inf; normalization then tries to calculate inf / inf which is NAN).
362 * A fix of this problem would consist in normalizing the messages muab.
363 * However, it is not obvious whether this is a real solution, because it has a
364 * negative performance impact and the NAN's seem to be a symptom of a fundamental
365 * numerical unstability.
367 muab(alpha
,_beta
).normalize();
371 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
372 size_t _beta
= alpha
.dual
;
373 Qb_new
*= muab(alpha
,_beta
) ^ (1 / (nbIR(beta
).size() + IR(beta
).c()));
377 if( Qb_new
.hasNaNs() ) {
378 // TODO: WHAT TO DO IN THIS CASE?
379 cerr
<< name() << "::doGBP: Qb_new has NaNs!" << endl
;
382 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
384 * _Qb[beta] = Qb_new.makeZero(1e-100);
387 if( props
.doubleloop
|| props
.damping
== 0.0 )
388 _Qb
[beta
] = Qb_new
; // no damping for double loop
390 _Qb
[beta
] = (Qb_new
^(1.0 - props
.damping
)) * (_Qb
[beta
]^props
.damping
);
392 bforeach( const Neighbor
&alpha
, nbIR(beta
) ) {
393 size_t _beta
= alpha
.dual
;
394 muba(alpha
,_beta
) = _Qb
[beta
] / muab(alpha
,_beta
);
396 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
398 * muba(beta,*alpha).makePositive(1e-100);
402 Factor Qa_new
= OR(alpha
);
403 bforeach( const Neighbor
&gamma
, nbOR(alpha
) )
404 Qa_new
*= muba(alpha
,gamma
.iter
);
405 Qa_new
^= (1.0 / OR(alpha
).c());
407 if( Qa_new
.hasNaNs() ) {
408 cerr
<< name() << "::doGBP: Qa_new has NaNs!" << endl
;
411 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
413 * _Qb[beta] = Qb_new.makeZero(1e-100);
416 if( props
.doubleloop
|| props
.damping
== 0.0 )
417 _Qa
[alpha
] = Qa_new
; // no damping for double loop
419 // FIXME: GEOMETRIC DAMPING IS SLOW!
420 _Qa
[alpha
] = (Qa_new
^(1.0 - props
.damping
)) * (_Qa
[alpha
]^props
.damping
);
424 // Calculate new single variable beliefs and compare with old ones
426 for( size_t i
= 0; i
< nrVars(); ++i
) {
427 Factor b
= beliefV(i
);
428 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], DISTLINF
) );
431 for( size_t I
= 0; I
< nrFactors(); ++I
) {
432 Factor b
= beliefF(I
);
433 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], DISTLINF
) );
437 if( props
.verbose
>= 3 )
438 cerr
<< name() << "::doGBP: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
441 if( maxDiff
> _maxdiff
)
444 if( props
.verbose
>= 1 ) {
445 if( maxDiff
> props
.tol
) {
446 if( props
.verbose
== 1 )
448 cerr
<< name() << "::doGBP: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
450 if( props
.verbose
>= 2 )
451 cerr
<< name() << "::doGBP: ";
452 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
460 Real
HAK::doDoubleLoop() {
461 if( props
.verbose
>= 1 )
462 cerr
<< "Starting " << identify() << "...";
463 if( props
.verbose
>= 3)
468 // Save original outer regions
469 vector
<FRegion
> org_ORs
= _ORs
;
471 // Save original inner counting numbers and set negative counting numbers to zero
472 vector
<Real
> org_IR_cs( nrIRs(), 0.0 );
473 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
474 org_IR_cs
[beta
] = IR(beta
).c();
475 if( IR(beta
).c() < 0.0 )
479 // Keep old beliefs to check convergence
480 vector
<Factor
> oldBeliefsV
;
481 oldBeliefsV
.reserve( nrVars() );
482 for( size_t i
= 0; i
< nrVars(); i
++ )
483 oldBeliefsV
.push_back( beliefV(i
) );
484 vector
<Factor
> oldBeliefsF
;
485 oldBeliefsF
.reserve( nrFactors() );
486 for( size_t I
= 0; I
< nrFactors(); I
++ )
487 oldBeliefsF
.push_back( beliefF(I
) );
489 size_t outer_maxiter
= props
.maxiter
;
490 Real outer_tol
= props
.tol
;
491 size_t outer_verbose
= props
.verbose
;
492 Real org_maxdiff
= _maxdiff
;
494 // Set parameters for inner loop
496 props
.verbose
= outer_verbose
? outer_verbose
- 1 : 0;
498 size_t outer_iter
= 0;
499 size_t total_iter
= 0;
500 Real maxDiff
= INFINITY
;
501 for( outer_iter
= 0; outer_iter
< outer_maxiter
&& maxDiff
> outer_tol
&& (toc() - tic
) < props
.maxtime
; outer_iter
++ ) {
502 // Calculate new outer regions
503 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
504 OR(alpha
) = org_ORs
[alpha
];
505 bforeach( const Neighbor
&beta
, nbOR(alpha
) )
506 OR(alpha
) *= _Qb
[beta
] ^ ((IR(beta
).c() - org_IR_cs
[beta
]) / nbIR(beta
).size());
510 if( isnan( doGBP() ) )
513 // Calculate new single variable beliefs and compare with old ones
515 for( size_t i
= 0; i
< nrVars(); ++i
) {
516 Factor b
= beliefV(i
);
517 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsV
[i
], DISTLINF
) );
520 for( size_t I
= 0; I
< nrFactors(); ++I
) {
521 Factor b
= beliefF(I
);
522 maxDiff
= std::max( maxDiff
, dist( b
, oldBeliefsF
[I
], DISTLINF
) );
526 total_iter
+= Iterations();
528 if( props
.verbose
>= 3 )
529 cerr
<< name() << "::doDoubleLoop: maxdiff " << maxDiff
<< " after " << total_iter
<< " passes" << endl
;
532 // restore _maxiter, _verbose and _maxdiff
533 props
.maxiter
= outer_maxiter
;
534 props
.verbose
= outer_verbose
;
535 _maxdiff
= org_maxdiff
;
538 if( maxDiff
> _maxdiff
)
541 // Restore original outer regions
544 // Restore original inner counting numbers
545 for( size_t beta
= 0; beta
< nrIRs(); ++beta
)
546 IR(beta
).c() = org_IR_cs
[beta
];
548 if( props
.verbose
>= 1 ) {
549 if( maxDiff
> props
.tol
) {
550 if( props
.verbose
== 1 )
552 cerr
<< name() << "::doDoubleLoop: WARNING: not converged after " << total_iter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
554 if( props
.verbose
>= 3 )
555 cerr
<< name() << "::doDoubleLoop: ";
556 cerr
<< "converged in " << total_iter
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
565 if( props
.doubleloop
)
566 return doDoubleLoop();
572 Factor
HAK::belief( const VarSet
&ns
) const {
573 vector
<Factor
>::const_iterator beta
;
574 for( beta
= _Qb
.begin(); beta
!= _Qb
.end(); beta
++ )
575 if( beta
->vars() >> ns
)
577 if( beta
!= _Qb
.end() )
578 return( beta
->marginal(ns
) );
580 vector
<Factor
>::const_iterator alpha
;
581 for( alpha
= _Qa
.begin(); alpha
!= _Qa
.end(); alpha
++ )
582 if( alpha
->vars() >> ns
)
584 if( alpha
== _Qa
.end() )
585 DAI_THROW(BELIEF_NOT_AVAILABLE
);
586 return( alpha
->marginal(ns
) );
591 vector
<Factor
> HAK::beliefs() const {
592 vector
<Factor
> result
;
593 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
594 result
.push_back( Qb(beta
) );
595 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ )
596 result
.push_back( Qa(alpha
) );
601 Real
HAK::logZ() const {
603 for( size_t beta
= 0; beta
< nrIRs(); beta
++ )
604 s
+= IR(beta
).c() * Qb(beta
).entropy();
605 for( size_t alpha
= 0; alpha
< nrORs(); alpha
++ ) {
606 s
+= OR(alpha
).c() * Qa(alpha
).entropy();
607 s
+= (OR(alpha
).log(true) * Qa(alpha
)).sum();
613 } // end of namespace dai