Merge branch 'master' of git@git.tuebingen.mpg.de:libdai
[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 void HAK::constructMessages() {
71 // Create outer beliefs
72 _Qa.clear();
73 _Qa.reserve(nrORs());
74 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
75 _Qa.push_back( Factor( OR(alpha).vars() ) );
76
77 // Create inner beliefs
78 _Qb.clear();
79 _Qb.reserve(nrIRs());
80 for( size_t beta = 0; beta < nrIRs(); beta++ )
81 _Qb.push_back( Factor( IR(beta) ) );
82
83 // Create messages
84 _muab.clear();
85 _muab.reserve( nrORs() );
86 _muba.clear();
87 _muba.reserve( nrORs() );
88 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
89 _muab.push_back( vector<Factor>() );
90 _muba.push_back( vector<Factor>() );
91 _muab[alpha].reserve( nbOR(alpha).size() );
92 _muba[alpha].reserve( nbOR(alpha).size() );
93 foreach( const Neighbor &beta, nbOR(alpha) ) {
94 _muab[alpha].push_back( Factor( IR(beta) ) );
95 _muba[alpha].push_back( Factor( IR(beta) ) );
96 }
97 }
98 }
99
100
101 HAK::HAK(const RegionGraph & rg, const PropertySet &opts ) : DAIAlgRG(rg) {
102 setProperties( opts );
103
104 constructMessages();
105 }
106
107
108 void HAK::findLoopClusters( const FactorGraph & fg, std::set<VarSet> &allcl, VarSet newcl, const Var & root, size_t length, VarSet vars ) {
109 for( VarSet::const_iterator in = vars.begin(); in != vars.end(); in++ ) {
110 VarSet ind = fg.delta( fg.findVar( *in ) );
111 if( (newcl.size()) >= 2 && (ind >> root) ) {
112 allcl.insert( newcl | *in );
113 }
114 else if( length > 1 )
115 findLoopClusters( fg, allcl, newcl | *in, root, length - 1, ind / newcl );
116 }
117 }
118
119
120 HAK::HAK(const FactorGraph & fg, const PropertySet &opts) : DAIAlgRG(), props(), maxdiff(0.0) {
121 setProperties( opts );
122
123 vector<VarSet> cl;
124 if( props.clusters == Properties::ClustersType::MIN ) {
125 cl = fg.Cliques();
126 } else if( props.clusters == Properties::ClustersType::DELTA ) {
127 for( size_t i = 0; i < fg.nrVars(); i++ )
128 cl.push_back(fg.Delta(i));
129 } else if( props.clusters == Properties::ClustersType::LOOP ) {
130 cl = fg.Cliques();
131 set<VarSet> scl;
132 for( size_t i0 = 0; i0 < fg.nrVars(); i0++ ) {
133 VarSet i0d = fg.delta(i0);
134 if( props.loopdepth > 1 )
135 findLoopClusters( fg, scl, fg.var(i0), fg.var(i0), props.loopdepth - 1, fg.delta(i0) );
136 }
137 for( set<VarSet>::const_iterator c = scl.begin(); c != scl.end(); c++ )
138 cl.push_back(*c);
139 if( props.verbose >= 3 ) {
140 cout << "HAK uses the following clusters: " << endl;
141 for( vector<VarSet>::const_iterator cli = cl.begin(); cli != cl.end(); cli++ )
142 cout << *cli << endl;
143 }
144 } else
145 DAI_THROW(INTERNAL_ERROR);
146
147 RegionGraph rg(fg,cl);
148 RegionGraph::operator=(rg);
149 constructMessages();
150
151 if( props.verbose >= 3 )
152 cout << "HAK regiongraph: " << *this << endl;
153 }
154
155
156 string HAK::identify() const {
157 stringstream result (stringstream::out);
158 result << Name << getProperties();
159 return result.str();
160 }
161
162
163 void HAK::init( const VarSet &ns ) {
164 for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
165 if( alpha->vars().intersects( ns ) )
166 alpha->fill( 1.0 / alpha->states() );
167
168 for( size_t beta = 0; beta < nrIRs(); beta++ )
169 if( IR(beta).intersects( ns ) ) {
170 _Qb[beta].fill( 1.0 );
171 foreach( const Neighbor &alpha, nbIR(beta) ) {
172 size_t _beta = alpha.dual;
173 muab( alpha, _beta ).fill( 1.0 / IR(beta).states() );
174 muba( alpha, _beta ).fill( 1.0 / IR(beta).states() );
175 }
176 }
177 }
178
179
180 void HAK::init() {
181 for( vector<Factor>::iterator alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
182 alpha->fill( 1.0 / alpha->states() );
183
184 for( vector<Factor>::iterator beta = _Qb.begin(); beta != _Qb.end(); beta++ )
185 beta->fill( 1.0 / beta->states() );
186
187 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
188 foreach( const Neighbor &beta, nbOR(alpha) ) {
189 size_t _beta = beta.iter;
190 muab( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
191 muba( alpha, _beta ).fill( 1.0 / muab( alpha, _beta ).states() );
192 }
193 }
194
195
196 double HAK::doGBP() {
197 if( props.verbose >= 1 )
198 cout << "Starting " << identify() << "...";
199 if( props.verbose >= 3)
200 cout << endl;
201
202 double tic = toc();
203
204 // Check whether counting numbers won't lead to problems
205 for( size_t beta = 0; beta < nrIRs(); beta++ )
206 assert( nbIR(beta).size() + IR(beta).c() != 0.0 );
207
208 // Keep old beliefs to check convergence
209 vector<Factor> old_beliefs;
210 old_beliefs.reserve( nrVars() );
211 for( size_t i = 0; i < nrVars(); i++ )
212 old_beliefs.push_back( belief( var(i) ) );
213
214 // Differences in single node beliefs
215 Diffs diffs(nrVars(), 1.0);
216
217 size_t iter = 0;
218 // do several passes over the network until maximum number of iterations has
219 // been reached or until the maximum belief difference is smaller than tolerance
220 for( iter = 0; iter < props.maxiter && diffs.maxDiff() > props.tol; iter++ ) {
221 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
222 foreach( const Neighbor &alpha, nbIR(beta) ) {
223 size_t _beta = alpha.dual;
224 muab( alpha, _beta ) = _Qa[alpha].marginal(IR(beta)).divided_by( muba(alpha,_beta) );
225 }
226
227 Factor Qb_new;
228 foreach( const Neighbor &alpha, nbIR(beta) ) {
229 size_t _beta = alpha.dual;
230 Qb_new *= muab(alpha,_beta) ^ (1 / (nbIR(beta).size() + IR(beta).c()));
231 }
232
233 Qb_new.normalize( Prob::NORMPROB );
234 if( Qb_new.hasNaNs() ) {
235 cout << "HAK::doGBP: Qb_new has NaNs!" << endl;
236 return 1.0;
237 }
238 // _Qb[beta] = Qb_new.makeZero(1e-100); // damping?
239 _Qb[beta] = Qb_new;
240
241 foreach( const Neighbor &alpha, nbIR(beta) ) {
242 size_t _beta = alpha.dual;
243
244 muba(alpha,_beta) = _Qb[beta].divided_by( muab(alpha,_beta) );
245
246 Factor Qa_new = OR(alpha);
247 foreach( const Neighbor &gamma, nbOR(alpha) )
248 Qa_new *= muba(alpha,gamma.iter);
249 Qa_new ^= (1.0 / OR(alpha).c());
250 Qa_new.normalize( Prob::NORMPROB );
251 if( Qa_new.hasNaNs() ) {
252 cout << "HAK::doGBP: Qa_new has NaNs!" << endl;
253 return 1.0;
254 }
255 // _Qa[alpha] = Qa_new.makeZero(1e-100); // damping?
256 _Qa[alpha] = Qa_new;
257 }
258 }
259
260 // Calculate new single variable beliefs and compare with old ones
261 for( size_t i = 0; i < nrVars(); i++ ) {
262 Factor new_belief = belief( var( i ) );
263 diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
264 old_beliefs[i] = new_belief;
265 }
266
267 if( props.verbose >= 3 )
268 cout << "HAK::doGBP: maxdiff " << diffs.maxDiff() << " after " << iter+1 << " passes" << endl;
269 }
270
271 if( diffs.maxDiff() > maxdiff )
272 maxdiff = diffs.maxDiff();
273
274 if( props.verbose >= 1 ) {
275 if( diffs.maxDiff() > props.tol ) {
276 if( props.verbose == 1 )
277 cout << endl;
278 cout << "HAK::doGBP: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
279 } else {
280 if( props.verbose >= 2 )
281 cout << "HAK::doGBP: ";
282 cout << "converged in " << iter << " passes (" << toc() - tic << " clocks)." << endl;
283 }
284 }
285
286 return diffs.maxDiff();
287 }
288
289
290 double HAK::doDoubleLoop() {
291 if( props.verbose >= 1 )
292 cout << "Starting " << identify() << "...";
293 if( props.verbose >= 3)
294 cout << endl;
295
296 double tic = toc();
297
298 // Save original outer regions
299 vector<FRegion> org_ORs = ORs;
300
301 // Save original inner counting numbers and set negative counting numbers to zero
302 vector<double> org_IR_cs( nrIRs(), 0.0 );
303 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
304 org_IR_cs[beta] = IR(beta).c();
305 if( IR(beta).c() < 0.0 )
306 IR(beta).c() = 0.0;
307 }
308
309 // Keep old beliefs to check convergence
310 vector<Factor> old_beliefs;
311 old_beliefs.reserve( nrVars() );
312 for( size_t i = 0; i < nrVars(); i++ )
313 old_beliefs.push_back( belief( var(i) ) );
314
315 // Differences in single node beliefs
316 Diffs diffs(nrVars(), 1.0);
317
318 size_t outer_maxiter = props.maxiter;
319 double outer_tol = props.tol;
320 size_t outer_verbose = props.verbose;
321 double org_maxdiff = maxdiff;
322
323 // Set parameters for inner loop
324 props.maxiter = 5;
325 props.verbose = outer_verbose ? outer_verbose - 1 : 0;
326
327 size_t outer_iter = 0;
328 for( outer_iter = 0; outer_iter < outer_maxiter && diffs.maxDiff() > outer_tol; outer_iter++ ) {
329 // Calculate new outer regions
330 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
331 OR(alpha) = org_ORs[alpha];
332 foreach( const Neighbor &beta, nbOR(alpha) )
333 OR(alpha) *= _Qb[beta] ^ ((IR(beta).c() - org_IR_cs[beta]) / nbIR(beta).size());
334 }
335
336 // Inner loop
337 if( isnan( doGBP() ) )
338 return 1.0;
339
340 // Calculate new single variable beliefs and compare with old ones
341 for( size_t i = 0; i < nrVars(); ++i ) {
342 Factor new_belief = belief( var( i ) );
343 diffs.push( dist( new_belief, old_beliefs[i], Prob::DISTLINF ) );
344 old_beliefs[i] = new_belief;
345 }
346
347 if( props.verbose >= 3 )
348 cout << "HAK::doDoubleLoop: maxdiff " << diffs.maxDiff() << " after " << outer_iter+1 << " passes" << endl;
349 }
350
351 // restore _maxiter, _verbose and _maxdiff
352 props.maxiter = outer_maxiter;
353 props.verbose = outer_verbose;
354 maxdiff = org_maxdiff;
355
356 if( diffs.maxDiff() > maxdiff )
357 maxdiff = diffs.maxDiff();
358
359 // Restore original outer regions
360 ORs = org_ORs;
361
362 // Restore original inner counting numbers
363 for( size_t beta = 0; beta < nrIRs(); ++beta )
364 IR(beta).c() = org_IR_cs[beta];
365
366 if( props.verbose >= 1 ) {
367 if( diffs.maxDiff() > props.tol ) {
368 if( props.verbose == 1 )
369 cout << endl;
370 cout << "HAK::doDoubleLoop: WARNING: not converged within " << outer_maxiter << " passes (" << toc() - tic << " clocks)...final maxdiff:" << diffs.maxDiff() << endl;
371 } else {
372 if( props.verbose >= 3 )
373 cout << "HAK::doDoubleLoop: ";
374 cout << "converged in " << outer_iter << " passes (" << toc() - tic << " clocks)." << endl;
375 }
376 }
377
378 return diffs.maxDiff();
379 }
380
381
382 double HAK::run() {
383 if( props.doubleloop )
384 return doDoubleLoop();
385 else
386 return doGBP();
387 }
388
389
390 Factor HAK::belief( const VarSet &ns ) const {
391 vector<Factor>::const_iterator beta;
392 for( beta = _Qb.begin(); beta != _Qb.end(); beta++ )
393 if( beta->vars() >> ns )
394 break;
395 if( beta != _Qb.end() )
396 return( beta->marginal(ns) );
397 else {
398 vector<Factor>::const_iterator alpha;
399 for( alpha = _Qa.begin(); alpha != _Qa.end(); alpha++ )
400 if( alpha->vars() >> ns )
401 break;
402 assert( alpha != _Qa.end() );
403 return( alpha->marginal(ns) );
404 }
405 }
406
407
408 Factor HAK::belief( const Var &n ) const {
409 return belief( (VarSet)n );
410 }
411
412
413 vector<Factor> HAK::beliefs() const {
414 vector<Factor> result;
415 for( size_t beta = 0; beta < nrIRs(); beta++ )
416 result.push_back( Qb(beta) );
417 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
418 result.push_back( Qa(alpha) );
419 return result;
420 }
421
422
423 Real HAK::logZ() const {
424 Real sum = 0.0;
425 for( size_t beta = 0; beta < nrIRs(); beta++ )
426 sum += IR(beta).c() * Qb(beta).entropy();
427 for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
428 sum += OR(alpha).c() * Qa(alpha).entropy();
429 sum += (OR(alpha).log0() * Qa(alpha)).totalSum();
430 }
431 return sum;
432 }
433
434
435 } // end of namespace dai