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