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