Added BETHE method for GBP/HAK cluster choice
[libdai.git] / src / hak.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
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.
6 *
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
9 */
10
11
12 #include <map>
13 #include <dai/hak.h>
14 #include <dai/util.h>
15 #include <dai/exceptions.h>
16
17
18 namespace dai {
19
20
21 using namespace std;
22
23
24 const char *HAK::Name = "HAK";
25
26
27 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
28 template <class T>
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) )
32 f[t] = epsilon;
33 return f;
34 }
35
36 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
37 template <class T>
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 )
41 f[t] = 0;
42 return f;
43 }
44
45
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") );
52
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");
58
59 if( opts.hasKey("loopdepth") )
60 props.loopdepth = opts.getStringAs<size_t>("loopdepth");
61 else
62 DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP );
63 if( opts.hasKey("damping") )
64 props.damping = opts.getStringAs<Real>("damping");
65 else
66 props.damping = 0.0;
67 if( opts.hasKey("init") )
68 props.init = opts.getStringAs<Properties::InitType>("init");
69 else
70 props.init = Properties::InitType::UNIFORM;
71 }
72
73
74 PropertySet HAK::getProperties() const {
75 PropertySet opts;
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 );
84 return opts;
85 }
86
87
88 string HAK::printProperties() const {
89 stringstream s( stringstream::out );
90 s << "[";
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 << "]";
99 return s.str();
100 }
101
102
103 void HAK::construct() {
104 // Create outer beliefs
105 _Qa.clear();
106 _Qa.reserve(nrORs());
107 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
108 _Qa.push_back( Factor( OR(alpha) ) );
109
110 // Create inner beliefs
111 _Qb.clear();
112 _Qb.reserve(nrIRs());
113 for( size_t beta = 0; beta < nrIRs(); beta++ )
114 _Qb.push_back( Factor( IR(beta) ) );
115
116 // Create messages
117 _muab.clear();
118 _muab.reserve( nrORs() );
119 _muba.clear();
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) ) );
129 }
130 }
131 }
132
133
134 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
135 setProperties( opts );
136
137 construct();
138 }
139
140
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 );
148 }
149 }
150
151
152 HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
153 setProperties( opts );
154
155 vector<VarSet> cl;
156 if( props.clusters == Properties::ClustersType::MIN ) {
157 cl = fg.Cliques();
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 ) {
165 cl = fg.Cliques();
166 set<VarSet> scl;
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) );
171 }
172 for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
173 cl.push_back(*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;
178 }
179 constructCVM( fg, cl );
180 } else if( props.clusters == Properties::ClustersType::BETHE ) {
181 // build outer regions (the cliques)
182 cl = fg.Cliques();
183 size_t nrEdges = 0;
184 for( size_t c = 0; c < cl.size(); c++ )
185 nrEdges += cl[c].size();
186
187 // build inner regions (single variables)
188 vector<Region> irs;
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 ) );
192
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 ) );
201 irs[i].c() -= 1.0;
202 }
203
204 // build region graph
205 RegionGraph::construct( fg, cl, irs, edges );
206 } else
207 DAI_THROW(UNKNOWN_ENUM_VALUE);
208
209 construct();
210
211 if( props.verbose >= 3 )
212 cerr << Name << " regiongraph: " << *this << endl;
213 }
214
215
216 string HAK::identify() const {
217 return string(Name) + printProperties();
218 }
219
220
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();
226 else
227 _Qa[alpha].randomize();
228 _Qa[alpha] *= OR(alpha);
229 _Qa[alpha].normalize();
230 }
231
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 );
236 else
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 );
243 } else {
244 muab( alpha, _beta ).randomize();
245 muba( alpha, _beta ).randomize();
246 }
247 }
248 }
249 }
250
251
252 void HAK::init() {
253 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
254 if( props.init == Properties::InitType::UNIFORM )
255 _Qa[alpha].setUniform();
256 else
257 _Qa[alpha].randomize();
258 _Qa[alpha] *= OR(alpha);
259 _Qa[alpha].normalize();
260 }
261
262 for( size_t beta = 0; beta < nrIRs(); beta++ )
263 if( props.init == Properties::InitType::UNIFORM )
264 _Qb[beta].setUniform();
265 else
266 _Qb[beta].randomize();
267
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();
274 } else {
275 muab( alpha, _beta ).randomize();
276 muba( alpha, _beta ).randomize();
277 }
278 }
279 }
280
281
282 Real HAK::doGBP() {
283 if( props.verbose >= 1 )
284 cerr << "Starting " << identify() << "...";
285 if( props.verbose >= 3)
286 cerr << endl;
287
288 double tic = toc();
289
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 );
293
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) );
303
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
313 *
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.
321 */
322 muab(alpha,_beta).normalize();
323 }
324
325 Factor Qb_new;
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()));
329 }
330
331 Qb_new.normalize();
332 if( Qb_new.hasNaNs() ) {
333 // TODO: WHAT TO DO IN THIS CASE?
334 cerr << Name << "::doGBP: Qb_new has NaNs!" << endl;
335 return 1.0;
336 }
337 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
338 *
339 * _Qb[beta] = Qb_new.makeZero(1e-100);
340 */
341
342 if( props.doubleloop || props.damping == 0.0 )
343 _Qb[beta] = Qb_new; // no damping for double loop
344 else
345 _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
346
347 foreach( const Neighbor &alpha, nbIR(beta) ) {
348 size_t _beta = alpha.dual;
349 muba(alpha,_beta) = _Qb[beta] / muab(alpha,_beta);
350
351 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
352 *
353 * muba(beta,*alpha).makePositive(1e-100);
354 *
355 */
356
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());
361 Qa_new.normalize();
362 if( Qa_new.hasNaNs() ) {
363 cerr << Name << "::doGBP: Qa_new has NaNs!" << endl;
364 return 1.0;
365 }
366 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
367 *
368 * _Qb[beta] = Qb_new.makeZero(1e-100);
369 */
370
371 if( props.doubleloop || props.damping == 0.0 )
372 _Qa[alpha] = Qa_new; // no damping for double loop
373 else
374 // FIXME: GEOMETRIC DAMPING IS SLOW!
375 _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
376 }
377 }
378
379 // Calculate new single variable beliefs and compare with old ones
380 maxDiff = -INFINITY;
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 ) );
384 oldBeliefsV[i] = b;
385 }
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 ) );
389 oldBeliefsF[I] = b;
390 }
391
392 if( props.verbose >= 3 )
393 cerr << Name << "::doGBP: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
394 }
395
396 if( maxDiff > _maxdiff )
397 _maxdiff = maxDiff;
398
399 if( props.verbose >= 1 ) {
400 if( maxDiff > props.tol ) {
401 if( props.verbose == 1 )
402 cerr << endl;
403 cerr << Name << "::doGBP: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
404 } else {
405 if( props.verbose >= 2 )
406 cerr << Name << "::doGBP: ";
407 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
408 }
409 }
410
411 return maxDiff;
412 }
413
414
415 Real HAK::doDoubleLoop() {
416 if( props.verbose >= 1 )
417 cerr << "Starting " << identify() << "...";
418 if( props.verbose >= 3)
419 cerr << endl;
420
421 double tic = toc();
422
423 // Save original outer regions
424 vector<FRegion> org_ORs = ORs;
425
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 )
431 IR(beta).c() = 0.0;
432 }
433
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) );
443
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;
448
449 // Set parameters for inner loop
450 props.maxiter = 5;
451 props.verbose = outer_verbose ? outer_verbose - 1 : 0;
452
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());
462 }
463
464 // Inner loop
465 if( isnan( doGBP() ) )
466 return 1.0;
467
468 // Calculate new single variable beliefs and compare with old ones
469 maxDiff = -INFINITY;
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 ) );
473 oldBeliefsV[i] = b;
474 }
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 ) );
478 oldBeliefsF[I] = b;
479 }
480
481 total_iter += Iterations();
482
483 if( props.verbose >= 3 )
484 cerr << Name << "::doDoubleLoop: maxdiff " << maxDiff << " after " << total_iter << " passes" << endl;
485 }
486
487 // restore _maxiter, _verbose and _maxdiff
488 props.maxiter = outer_maxiter;
489 props.verbose = outer_verbose;
490 _maxdiff = org_maxdiff;
491
492 _iters = total_iter;
493 if( maxDiff > _maxdiff )
494 _maxdiff = maxDiff;
495
496 // Restore original outer regions
497 ORs = org_ORs;
498
499 // Restore original inner counting numbers
500 for( size_t beta = 0; beta < nrIRs(); ++beta )
501 IR(beta).c() = org_IR_cs[beta];
502
503 if( props.verbose >= 1 ) {
504 if( maxDiff > props.tol ) {
505 if( props.verbose == 1 )
506 cerr << endl;
507 cerr << Name << "::doDoubleLoop: WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
508 } else {
509 if( props.verbose >= 3 )
510 cerr << Name << "::doDoubleLoop: ";
511 cerr << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
512 }
513 }
514
515 return maxDiff;
516 }
517
518
519 Real HAK::run() {
520 if( props.doubleloop )
521 return doDoubleLoop();
522 else
523 return doGBP();
524 }
525
526
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 )
531 break;
532 if( beta != _Qb.end() )
533 return( beta->marginal(ns) );
534 else {
535 vector<Factor>::const_iterator alpha;
536 for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
537 if( alpha->vars() >> ns )
538 break;
539 if( alpha == _Qa.end() )
540 DAI_THROW(BELIEF_NOT_AVAILABLE);
541 return( alpha->marginal(ns) );
542 }
543 }
544
545
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) );
552 return result;
553 }
554
555
556 Real HAK::logZ() const {
557 Real s = 0.0;
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();
563 }
564 return s;
565 }
566
567
568 } // end of namespace dai