Declared class Diffs as obsolete
[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
28 /// Sets factor entries that lie between 0 and \a epsilon to \a epsilon
29 template <class T>
30 TFactor<T>& makePositive( TFactor<T> &f, T epsilon ) {
31 for( size_t t = 0; t < f.states(); t++ )
32 if( (0 < f[t]) && (f[t] < epsilon) )
33 f[t] = epsilon;
34 return f;
35 }
36
37 /// Sets factor entries that are smaller (in absolute value) than \a epsilon to 0
38 template <class T>
39 TFactor<T>& makeZero( TFactor<T> &f, T epsilon ) {
40 for( size_t t = 0; t < f.states(); t++ )
41 if( f[t] < epsilon && f[t] > -epsilon )
42 f[t] = 0;
43 return f;
44 }
45
46
47 void HAK::setProperties( const PropertySet &opts ) {
48 DAI_ASSERT( opts.hasKey("tol") );
49 DAI_ASSERT( opts.hasKey("maxiter") );
50 DAI_ASSERT( opts.hasKey("verbose") );
51 DAI_ASSERT( opts.hasKey("doubleloop") );
52 DAI_ASSERT( opts.hasKey("clusters") );
53
54 props.tol = opts.getStringAs<Real>("tol");
55 props.maxiter = opts.getStringAs<size_t>("maxiter");
56 props.verbose = opts.getStringAs<size_t>("verbose");
57 props.doubleloop = opts.getStringAs<bool>("doubleloop");
58 props.clusters = opts.getStringAs<Properties::ClustersType>("clusters");
59
60 if( opts.hasKey("loopdepth") )
61 props.loopdepth = opts.getStringAs<size_t>("loopdepth");
62 else
63 DAI_ASSERT( props.clusters != Properties::ClustersType::LOOP );
64 if( opts.hasKey("damping") )
65 props.damping = opts.getStringAs<Real>("damping");
66 else
67 props.damping = 0.0;
68 if( opts.hasKey("init") )
69 props.init = opts.getStringAs<Properties::InitType>("init");
70 else
71 props.init = Properties::InitType::UNIFORM;
72 }
73
74
75 PropertySet HAK::getProperties() const {
76 PropertySet opts;
77 opts.Set( "tol", props.tol );
78 opts.Set( "maxiter", props.maxiter );
79 opts.Set( "verbose", props.verbose );
80 opts.Set( "doubleloop", props.doubleloop );
81 opts.Set( "clusters", props.clusters );
82 opts.Set( "init", props.init );
83 opts.Set( "loopdepth", props.loopdepth );
84 opts.Set( "damping", props.damping );
85 return opts;
86 }
87
88
89 string HAK::printProperties() const {
90 stringstream s( stringstream::out );
91 s << "[";
92 s << "tol=" << props.tol << ",";
93 s << "maxiter=" << props.maxiter << ",";
94 s << "verbose=" << props.verbose << ",";
95 s << "doubleloop=" << props.doubleloop << ",";
96 s << "clusters=" << props.clusters << ",";
97 s << "init=" << props.init << ",";
98 s << "loopdepth=" << props.loopdepth << ",";
99 s << "damping=" << props.damping << "]";
100 return s.str();
101 }
102
103
104 void HAK::constructMessages() {
105 // Create outer beliefs
106 _Qa.clear();
107 _Qa.reserve(nrORs());
108 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
109 _Qa.push_back( Factor( OR(alpha).vars() ) );
110
111 // Create inner beliefs
112 _Qb.clear();
113 _Qb.reserve(nrIRs());
114 for( size_t beta = 0; beta < nrIRs(); beta++ )
115 _Qb.push_back( Factor( IR(beta) ) );
116
117 // Create messages
118 _muab.clear();
119 _muab.reserve( nrORs() );
120 _muba.clear();
121 _muba.reserve( nrORs() );
122 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
123 _muab.push_back( vector<Factor>() );
124 _muba.push_back( vector<Factor>() );
125 _muab[alpha].reserve( nbOR(alpha).size() );
126 _muba[alpha].reserve( nbOR(alpha).size() );
127 foreach( const Neighbor &beta, nbOR(alpha) ) {
128 _muab[alpha].push_back( Factor( IR(beta) ) );
129 _muba[alpha].push_back( Factor( IR(beta) ) );
130 }
131 }
132 }
133
134
135 HAK::HAK( const RegionGraph &rg, const PropertySet &opts ) : DAIAlgRG(rg), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
136 setProperties( opts );
137
138 constructMessages();
139 }
140
141
142 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
143 for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
144 VarSet ind = fg.delta( fg.findVar( *in ) );
145 if( (newcl.size()) >= 2 && ind.contains( root ) ) {
146 allcl.insert( newcl | *in );
147 }
148 else if( length > 1 )
149 findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
150 }
151 }
152
153
154 HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), _Qa(), _Qb(), _muab(), _muba(), _maxdiff(0.0), _iters(0U), props() {
155 setProperties( opts );
156
157 vector<VarSet> cl;
158 if( props.clusters == Properties::ClustersType::MIN ) {
159 cl = fg.Cliques();
160 } else if( props.clusters == Properties::ClustersType::DELTA ) {
161 for( size_t i = 0; i < fg.nrVars(); i++ )
162 cl.push_back(fg.Delta(i));
163 } else if( props.clusters == Properties::ClustersType::LOOP ) {
164 cl = fg.Cliques();
165 set<VarSet> scl;
166 for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
167 VarSet i0d = fg.delta(i0);
168 if( props.loopdepth > 1 )
169 findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
170 }
171 for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
172 cl.push_back(*c);
173 if( props.verbose >= 3 ) {
174 cerr << Name << " uses the following clusters: " << endl;
175 for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
176 cerr << *cli << endl;
177 }
178 } else
179 DAI_THROW(UNKNOWN_ENUM_VALUE);
180
181 RegionGraph rg(fg,cl);
182 RegionGraph::operator=(rg);
183 constructMessages();
184
185 if( props.verbose >= 3 )
186 cerr << Name << " regiongraph: " << *this << endl;
187 }
188
189
190 string HAK::identify() const {
191 return string(Name) + printProperties();
192 }
193
194
195 void HAK::init( const VarSet &ns ) {
196 for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
197 if( alpha->vars().intersects( ns ) ) {
198 if( props.init == Properties::InitType::UNIFORM )
199 alpha->fill( 1.0 / alpha->states() );
200 else
201 alpha->randomize();
202 }
203
204 for( size_t beta = 0; beta < nrIRs(); beta++ )
205 if( IR(beta).intersects( ns ) ) {
206 if( props.init == Properties::InitType::UNIFORM )
207 _Qb[beta].fill( 1.0 );
208 else
209 _Qb[beta].randomize();
210 foreach( const Neighbor &alpha, nbIR(beta) ) {
211 size_t _beta = alpha.dual;
212 if( props.init == Properties::InitType::UNIFORM ) {
213 muab( alpha, _beta ).fill( 1.0 );
214 muba( alpha, _beta ).fill( 1.0 );
215 } else {
216 muab( alpha, _beta ).randomize();
217 muba( alpha, _beta ).randomize();
218 }
219 }
220 }
221 }
222
223
224 void HAK::init() {
225 for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
226 if( props.init == Properties::InitType::UNIFORM )
227 alpha->fill( 1.0 / alpha->states() );
228 else
229 alpha->randomize();
230
231 for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
232 if( props.init == Properties::InitType::UNIFORM )
233 beta->fill( 1.0 / beta->states() );
234 else
235 beta->randomize();
236
237 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
238 foreach( const Neighbor &beta, nbOR(alpha) ) {
239 size_t _beta = beta.iter;
240 if( props.init == Properties::InitType::UNIFORM ) {
241 muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
242 muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
243 } else {
244 muab( alpha, _beta ).randomize();
245 muba( alpha, _beta ).randomize();
246 }
247 }
248 }
249
250
251 Real HAK::doGBP() {
252 if( props.verbose >= 1 )
253 cerr << "Starting " << identify() << "...";
254 if( props.verbose >= 3)
255 cerr << endl;
256
257 double tic = toc();
258
259 // Check whether counting numbers won't lead to problems
260 for( size_t beta = 0; beta < nrIRs(); beta++ )
261 DAI_ASSERT( nbIR(beta).size() + IR(beta).c() != 0.0 );
262
263 // Keep old beliefs to check convergence
264 vector<Factor> old_beliefs;
265 old_beliefs.reserve( nrVars() );
266 for( size_t i = 0; i < nrVars(); i++ )
267 old_beliefs.push_back( belief( var(i) ) );
268
269 // Differences in single node beliefs
270 vector<Real> diffs( nrVars(), INFINITY );
271 Real maxDiff = INFINITY;
272
273 // do several passes over the network until maximum number of iterations has
274 // been reached or until the maximum belief difference is smaller than tolerance
275 for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
276 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
277 foreach( const Neighbor &alpha, nbIR(beta) ) {
278 size_t _beta = alpha.dual;
279 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)) / muba(alpha,_beta);
280 /* TODO: INVESTIGATE THIS PROBLEM
281 *
282 * In some cases, the muab's can have very large entries because the muba's have very
283 * small entries. This may cause NANs later on (e.g., multiplying large quantities may
284 * result in +inf; normalization then tries to calculate inf / inf which is NAN).
285 * A fix of this problem would consist in normalizing the messages muab.
286 * However, it is not obvious whether this is a real solution, because it has a
287 * negative performance impact and the NAN's seem to be a symptom of a fundamental
288 * numerical unstability.
289 */
290 muab(alpha,_beta).normalize();
291 }
292
293 Factor Qb_new;
294 foreach( const Neighbor &alpha, nbIR(beta) ) {
295 size_t _beta = alpha.dual;
296 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
297 }
298
299 Qb_new.normalize();
300 if( Qb_new.hasNaNs() ) {
301 // TODO: WHAT TO DO IN THIS CASE?
302 cerr << Name << "::doGBP: Qb_new has NaNs!" << endl;
303 return 1.0;
304 }
305 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
306 *
307 * _Qb[beta] = Qb_new.makeZero(1e-100);
308 */
309
310 if( props.doubleloop || props.damping == 0.0 )
311 _Qb[beta] = Qb_new; // no damping for double loop
312 else
313 _Qb[beta] = (Qb_new^(1.0 - props.damping)) * (_Qb[beta]^props.damping);
314
315 foreach( const Neighbor &alpha, nbIR(beta) ) {
316 size_t _beta = alpha.dual;
317 muba(alpha,_beta) = _Qb[beta] / muab(alpha,_beta);
318
319 /* TODO: INVESTIGATE WHETHER THIS HACK (INVENTED BY KEES) TO PREVENT NANS MAKES SENSE
320 *
321 * muba(beta,*alpha).makePositive(1e-100);
322 *
323 */
324
325 Factor Qa_new = OR(alpha);
326 foreach( const Neighbor &gamma, nbOR(alpha) )
327 Qa_new *= muba(alpha,gamma.iter);
328 Qa_new ^= (1.0 / OR(alpha).c());
329 Qa_new.normalize();
330 if( Qa_new.hasNaNs() ) {
331 cerr << Name << "::doGBP: Qa_new has NaNs!" << endl;
332 return 1.0;
333 }
334 /* TODO: WHAT IS THE PURPOSE OF THE FOLLOWING CODE?
335 *
336 * _Qb[beta] = Qb_new.makeZero(1e-100);
337 */
338
339 if( props.doubleloop || props.damping == 0.0 )
340 _Qa[alpha] = Qa_new; // no damping for double loop
341 else
342 // FIXME: GEOMETRIC DAMPING IS SLOW!
343 _Qa[alpha] = (Qa_new^(1.0 - props.damping)) * (_Qa[alpha]^props.damping);
344 }
345 }
346
347 // Calculate new single variable beliefs and compare with old ones
348 for( size_t i = 0; i < nrVars(); i++ ) {
349 Factor new_belief = belief( var( i ) );
350 diffs[i] = dist( new_belief, old_beliefs[i], Prob::DISTLINF );
351 old_beliefs[i] = new_belief;
352 }
353 maxDiff = max( diffs );
354
355 if( props.verbose >= 3 )
356 cerr << Name << "::doGBP: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
357 }
358
359 if( maxDiff > _maxdiff )
360 _maxdiff = maxDiff;
361
362 if( props.verbose >= 1 ) {
363 if( maxDiff > props.tol ) {
364 if( props.verbose == 1 )
365 cerr << endl;
366 cerr << Name << "::doGBP: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
367 } else {
368 if( props.verbose >= 2 )
369 cerr << Name << "::doGBP: ";
370 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
371 }
372 }
373
374 return maxDiff;
375 }
376
377
378 Real HAK::doDoubleLoop() {
379 if( props.verbose >= 1 )
380 cerr << "Starting " << identify() << "...";
381 if( props.verbose >= 3)
382 cerr << endl;
383
384 double tic = toc();
385
386 // Save original outer regions
387 vector<FRegion> org_ORs = ORs;
388
389 // Save original inner counting numbers and set negative counting numbers to zero
390 vector<Real> org_IR_cs( nrIRs(), 0.0 );
391 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
392 org_IR_cs[beta] = IR(beta).c();
393 if( IR(beta).c() < 0.0 )
394 IR(beta).c() = 0.0;
395 }
396
397 // Keep old beliefs to check convergence
398 vector<Factor> old_beliefs;
399 old_beliefs.reserve( nrVars() );
400 for( size_t i = 0; i < nrVars(); i++ )
401 old_beliefs.push_back( belief( var(i) ) );
402
403 // Differences in single node beliefs
404 vector<Real> diffs( nrVars(), INFINITY );
405 Real maxDiff = INFINITY;
406
407 size_t outer_maxiter = props.maxiter;
408 Real outer_tol = props.tol;
409 size_t outer_verbose = props.verbose;
410 Real org_maxdiff = _maxdiff;
411
412 // Set parameters for inner loop
413 props.maxiter = 5;
414 props.verbose = outer_verbose ? outer_verbose - 1 : 0;
415
416 size_t outer_iter = 0;
417 size_t total_iter = 0;
418 for( outer_iter = 0; outer_iter < outer_maxiter && maxDiff > outer_tol; outer_iter++ ) {
419 // Calculate new outer regions
420 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
421 OR(alpha) = org_ORs[alpha];
422 foreach( const Neighbor &beta, nbOR(alpha) )
423 OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
424 }
425
426 // Inner loop
427 if( isnan( doGBP() ) )
428 return 1.0;
429
430 // Calculate new single variable beliefs and compare with old ones
431 for( size_t i = 0; i < nrVars(); ++i ) {
432 Factor new_belief = belief( var( i ) );
433 diffs[i] = dist( new_belief, old_beliefs[i], Prob::DISTLINF );
434 old_beliefs[i] = new_belief;
435 }
436 maxDiff = max( diffs );
437
438 total_iter += Iterations();
439
440 if( props.verbose >= 3 )
441 cerr << Name << "::doDoubleLoop: maxdiff " << maxDiff << " after " << total_iter << " passes" << endl;
442 }
443
444 // restore _maxiter, _verbose and _maxdiff
445 props.maxiter = outer_maxiter;
446 props.verbose = outer_verbose;
447 _maxdiff = org_maxdiff;
448
449 _iters = total_iter;
450 if( maxDiff > _maxdiff )
451 _maxdiff = maxDiff;
452
453 // Restore original outer regions
454 ORs = org_ORs;
455
456 // Restore original inner counting numbers
457 for( size_t beta = 0; beta < nrIRs(); ++beta )
458 IR(beta).c() = org_IR_cs[beta];
459
460 if( props.verbose >= 1 ) {
461 if( maxDiff > props.tol ) {
462 if( props.verbose == 1 )
463 cerr << endl;
464 cerr << Name << "::doDoubleLoop: WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
465 } else {
466 if( props.verbose >= 3 )
467 cerr << Name << "::doDoubleLoop: ";
468 cerr << "converged in " << total_iter << " passes (" << toc() - tic << " seconds)." << endl;
469 }
470 }
471
472 return maxDiff;
473 }
474
475
476 Real HAK::run() {
477 if( props.doubleloop )
478 return doDoubleLoop();
479 else
480 return doGBP();
481 }
482
483
484 Factor HAK::belief( const VarSet &ns ) const {
485 vector<Factor>::const_iterator beta;
486 for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
487 if( beta->vars() >> ns )
488 break;
489 if( beta != _Qb.end() )
490 return( beta->marginal(ns) );
491 else {
492 vector<Factor>::const_iterator alpha;
493 for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
494 if( alpha->vars() >> ns )
495 break;
496 DAI_ASSERT( alpha != _Qa.end() );
497 return( alpha->marginal(ns) );
498 }
499 }
500
501
502 Factor HAK::belief( const Var &n ) const {
503 return belief( (VarSet)n );
504 }
505
506
507 vector<Factor> HAK::beliefs() const {
508 vector<Factor> result;
509 for( size_t beta = 0; beta < nrIRs(); beta++ )
510 result.push_back( Qb(beta) );
511 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
512 result.push_back( Qa(alpha) );
513 return result;
514 }
515
516
517 Real HAK::logZ() const {
518 Real s = 0.0;
519 for( size_t beta = 0; beta < nrIRs(); beta++ )
520 s += IR(beta).c() * Qb(beta).entropy();
521 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
522 s += OR(alpha).c() * Qa(alpha).entropy();
523 s += (OR(alpha).log(true) * Qa(alpha)).sum();
524 }
525 return s;
526 }
527
528
529 } // end of namespace dai