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