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