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