Multiple changes: changes in build system, one workaround and one bug fix
[libdai.git] / src / mr.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <dai/dai_config.h>
10 #ifdef DAI_WITH_MR
11
12
13 #include <cstdio>
14 #include <ctime>
15 #include <cmath>
16 #include <cstdlib>
17 #include <dai/mr.h>
18 #include <dai/bp.h>
19 #include <dai/jtree.h>
20 #include <dai/util.h>
21 #include <dai/bbp.h>
22
23
24 namespace dai {
25
26
27 using namespace std;
28
29
30 void MR::setProperties( const PropertySet &opts ) {
31 DAI_ASSERT( opts.hasKey("tol") );
32 DAI_ASSERT( opts.hasKey("updates") );
33 DAI_ASSERT( opts.hasKey("inits") );
34
35 props.tol = opts.getStringAs<Real>("tol");
36 props.updates = opts.getStringAs<Properties::UpdateType>("updates");
37 props.inits = opts.getStringAs<Properties::InitType>("inits");
38 if( opts.hasKey("verbose") )
39 props.verbose = opts.getStringAs<size_t>("verbose");
40 else
41 props.verbose = 0;
42 }
43
44
45 PropertySet MR::getProperties() const {
46 PropertySet opts;
47 opts.set( "tol", props.tol );
48 opts.set( "verbose", props.verbose );
49 opts.set( "updates", props.updates );
50 opts.set( "inits", props.inits );
51 return opts;
52 }
53
54
55 string MR::printProperties() const {
56 stringstream s( stringstream::out );
57 s << "[";
58 s << "tol=" << props.tol << ",";
59 s << "verbose=" << props.verbose << ",";
60 s << "updates=" << props.updates << ",";
61 s << "inits=" << props.inits << "]";
62 return s.str();
63 }
64
65
66 Real MR::T(size_t i, sub_nb A) {
67 sub_nb _nbi_min_A(G.nb(i).size());
68 _nbi_min_A.set();
69 _nbi_min_A &= ~A;
70
71 Real res = theta[i];
72 for( size_t _j = 0; _j < _nbi_min_A.size(); _j++ )
73 if( _nbi_min_A.test(_j) )
74 res += atanh(tJ[i][_j] * M[i][_j]);
75 return tanh(res);
76 }
77
78
79 Real MR::T(size_t i, size_t _j) {
80 sub_nb j(G.nb(i).size());
81 j.set(_j);
82 return T(i,j);
83 }
84
85
86 Real MR::Omega(size_t i, size_t _j, size_t _l) {
87 sub_nb jl(G.nb(i).size());
88 jl.set(_j);
89 jl.set(_l);
90 Real Tijl = T(i,jl);
91 return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
92 }
93
94
95 Real MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
96 sub_nb jll(G.nb(i).size());
97 jll.set(_j);
98 Real Tij = T(i,jll);
99 jll.set(_l1);
100 jll.set(_l2);
101 Real Tijll = T(i,jll);
102
103 return (Tijll - Tij) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Tijll + tJ[i][_l2] * M[i][_l2] * Tijll);
104 }
105
106
107 Real MR::Gamma(size_t i, size_t _l1, size_t _l2) {
108 sub_nb ll(G.nb(i).size());
109 Real Ti = T(i,ll);
110 ll.set(_l1);
111 ll.set(_l2);
112 Real Till = T(i,ll);
113
114 return (Till - Ti) / (1.0 + tJ[i][_l1] * tJ[i][_l2] * M[i][_l1] * M[i][_l2] + tJ[i][_l1] * M[i][_l1] * Till + tJ[i][_l2] * M[i][_l2] * Till);
115 }
116
117
118 Real MR::_tJ(size_t i, sub_nb A) {
119 sub_nb::size_type _j = A.find_first();
120 if( _j == sub_nb::npos )
121 return 1.0;
122 else
123 return tJ[i][_j] * _tJ(i, A.reset(_j));
124 }
125
126
127 Real MR::appM(size_t i, sub_nb A) {
128 sub_nb::size_type _j = A.find_first();
129 if( _j == sub_nb::npos )
130 return 1.0;
131 else {
132 sub_nb A_j(A); A_j.reset(_j);
133
134 Real result = M[i][_j] * appM(i, A_j);
135 for( size_t _k = 0; _k < A_j.size(); _k++ )
136 if( A_j.test(_k) ) {
137 sub_nb A_jk(A_j); A_jk.reset(_k);
138 result += cors[i][_j][_k] * appM(i,A_jk);
139 }
140
141 return result;
142 }
143 }
144
145
146 void MR::sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd) {
147 *sum_even = 0.0;
148 *sum_odd = 0.0;
149
150 sub_nb B(A.size());
151 do {
152 if( B.count() % 2 )
153 *sum_odd += _tJ(j,B) * appM(j,B);
154 else
155 *sum_even += _tJ(j,B) * appM(j,B);
156
157 // calc next subset B
158 size_t bit = 0;
159 for( ; bit < A.size(); bit++ )
160 if( A.test(bit) ) {
161 if( B.test(bit) )
162 B.reset(bit);
163 else {
164 B.set(bit);
165 break;
166 }
167 }
168 } while (!B.none());
169 }
170
171
172 void MR::propagateCavityFields() {
173 Real sum_even, sum_odd;
174 Real maxdev;
175 size_t maxruns = 1000;
176
177 for( size_t i = 0; i < G.nrNodes(); i++ )
178 bforeach( const Neighbor &j, G.nb(i) )
179 M[i][j.iter] = 0.1;
180
181 size_t run=0;
182 do {
183 maxdev=0.0;
184 run++;
185 for( size_t i = 0; i < G.nrNodes(); i++ ) {
186 bforeach( const Neighbor &j, G.nb(i) ) {
187 size_t _j = j.iter;
188 size_t _i = G.findNb(j,i);
189 DAI_ASSERT( G.nb(j,_i) == i );
190
191 Real newM = 0.0;
192 if( props.updates == Properties::UpdateType::FULL ) {
193 // find indices in nb(j) that do not correspond with i
194 sub_nb _nbj_min_i(G.nb(j).size());
195 _nbj_min_i.set();
196 _nbj_min_i.reset(_i);
197
198 // find indices in nb(i) that do not correspond with j
199 sub_nb _nbi_min_j(G.nb(i).size());
200 _nbi_min_j.set();
201 _nbi_min_j.reset(_j);
202
203 sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
204 newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
205
206 sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
207 Real denom = sum_even + tanh(theta[i]) * sum_odd;
208 Real numer = 0.0;
209 for(size_t _k=0; _k < G.nb(i).size(); _k++) if(_k != _j) {
210 sub_nb _nbi_min_jk(_nbi_min_j);
211 _nbi_min_jk.reset(_k);
212 sum_subs(i, _nbi_min_jk, &sum_even, &sum_odd);
213 numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
214 }
215 newM -= numer / denom;
216 } else if( props.updates == Properties::UpdateType::LINEAR ) {
217 newM = T(j,_i);
218 for(size_t _l=0; _l<G.nb(i).size(); _l++) if( _l != _j )
219 newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
220 for(size_t _l1=0; _l1<G.nb(j).size(); _l1++) if( _l1 != _i )
221 for( size_t _l2=_l1+1; _l2<G.nb(j).size(); _l2++) if( _l2 != _i)
222 newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
223 }
224
225 Real dev = newM - M[i][_j];
226 // dev *= 0.02;
227 if( abs(dev) >= maxdev )
228 maxdev = abs(dev);
229
230 newM = M[i][_j] + dev;
231 if( abs(newM) > 1.0 )
232 newM = (newM > 0.0) ? 1.0 : -1.0;
233 M[i][_j] = newM;
234 }
235 }
236 } while((maxdev>props.tol)&&(run<maxruns));
237
238 _iters = run;
239 if( maxdev > _maxdiff )
240 _maxdiff = maxdev;
241
242 if(run==maxruns){
243 if( props.verbose >= 1 )
244 cerr << "MR::propagateCavityFields: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
245 }
246 }
247
248
249 void MR::calcMagnetizations() {
250 for( size_t i = 0; i < G.nrNodes(); i++ ) {
251 if( props.updates == Properties::UpdateType::FULL ) {
252 // find indices in nb(i)
253 sub_nb _nbi( G.nb(i).size() );
254 _nbi.set();
255
256 // calc numerator1 and denominator1
257 Real sum_even, sum_odd;
258 sum_subs(i, _nbi, &sum_even, &sum_odd);
259
260 Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
261
262 } else if( props.updates == Properties::UpdateType::LINEAR ) {
263 sub_nb empty( G.nb(i).size() );
264 Mag[i] = T(i,empty);
265
266 for( size_t _l1 = 0; _l1 < G.nb(i).size(); _l1++ )
267 for( size_t _l2 = _l1 + 1; _l2 < G.nb(i).size(); _l2++ )
268 Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
269 }
270 if( abs( Mag[i] ) > 1.0 )
271 Mag[i] = (Mag[i] > 0.0) ? 1.0 : -1.0;
272 }
273 }
274
275
276 Real MR::calcCavityCorrelations() {
277 Real md = 0.0;
278 for( size_t i = 0; i < nrVars(); i++ ) {
279 vector<Factor> pairq;
280 if( props.inits == Properties::InitType::EXACT ) {
281 JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
282 jtcav.makeCavity( i );
283 pairq = calcPairBeliefs( jtcav, delta(i), false, true );
284 } else if( props.inits == Properties::InitType::CLAMPING ) {
285 BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
286 bpcav.makeCavity( i );
287
288 pairq = calcPairBeliefs( bpcav, delta(i), false, true );
289 md = std::max( md, bpcav.maxDiff() );
290 } else if( props.inits == Properties::InitType::RESPPROP ) {
291 BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
292 bpcav.makeCavity( i );
293 bpcav.makeCavity( i );
294 bpcav.init();
295 bpcav.run();
296
297 BBP bbp( &bpcav, PropertySet()("verbose",(size_t)0)("tol",(Real)1.0e-9)("maxiter",(size_t)10000)("damping",(Real)0.0)("updates",string("SEQ_MAX")) );
298 bforeach( const Neighbor &j, G.nb(i) ) {
299 // Create weights for magnetization of some spin
300 Prob p( 2, 0.0 );
301 p.set( 0, -1.0 );
302 p.set( 1, 1.0 );
303
304 // BBP cost function would be the magnetization of spin j
305 vector<Prob> b1_adj;
306 b1_adj.reserve( nrVars() );
307 for( size_t l = 0; l < nrVars(); l++ )
308 if( l == j )
309 b1_adj.push_back( p );
310 else
311 b1_adj.push_back( Prob( 2, 0.0 ) );
312 bbp.init_V( b1_adj );
313
314 // run BBP to estimate adjoints
315 bbp.run();
316
317 bforeach( const Neighbor &k, G.nb(i) ) {
318 if( k != j )
319 cors[i][j.iter][k.iter] = (bbp.adj_psi_V(k)[1] - bbp.adj_psi_V(k)[0]);
320 else
321 cors[i][j.iter][k.iter] = 0.0;
322 }
323 }
324 }
325
326 if( props.inits != Properties::InitType::RESPPROP ) {
327 for( size_t jk = 0; jk < pairq.size(); jk++ ) {
328 VarSet::const_iterator kit = pairq[jk].vars().begin();
329 size_t j = findVar( *(kit) );
330 size_t k = findVar( *(++kit) );
331 pairq[jk].normalize();
332 Real cor = (pairq[jk][3] - pairq[jk][2] - pairq[jk][1] + pairq[jk][0]) - (pairq[jk][3] + pairq[jk][2] - pairq[jk][1] - pairq[jk][0]) * (pairq[jk][3] - pairq[jk][2] + pairq[jk][1] - pairq[jk][0]);
333
334 size_t _j = G.findNb(i,j);
335 size_t _k = G.findNb(i,k);
336 cors[i][_j][_k] = cor;
337 cors[i][_k][_j] = cor;
338 }
339 }
340
341 }
342 return md;
343 }
344
345
346 Real MR::run() {
347 if( supported ) {
348 if( props.verbose >= 1 )
349 cerr << "Starting " << identify() << "...";
350
351 double tic = toc();
352
353 // approximate correlations of cavity spins
354 Real md = calcCavityCorrelations();
355 if( md > _maxdiff )
356 _maxdiff = md;
357
358 // solve messages
359 propagateCavityFields();
360
361 // calculate magnetizations
362 calcMagnetizations();
363
364 if( props.verbose >= 1 )
365 cerr << name() << " needed " << toc() - tic << " seconds." << endl;
366
367 return _maxdiff;
368 } else
369 return 1.0;
370 }
371
372
373 Factor MR::beliefV( size_t i ) const {
374 if( supported ) {
375 Real x[2];
376 x[0] = 0.5 - Mag[i] / 2.0;
377 x[1] = 0.5 + Mag[i] / 2.0;
378
379 return Factor( var(i), x );
380 } else
381 return Factor();
382 }
383
384
385 Factor MR::belief (const VarSet &ns) const {
386 if( ns.size() == 0 )
387 return Factor();
388 else if( ns.size() == 1 )
389 return beliefV( findVar( *(ns.begin()) ) );
390 else {
391 DAI_THROW(BELIEF_NOT_AVAILABLE);
392 return Factor();
393 }
394 }
395
396
397 vector<Factor> MR::beliefs() const {
398 vector<Factor> result;
399 for( size_t i = 0; i < nrVars(); i++ )
400 result.push_back( beliefV( i ) );
401 return result;
402 }
403
404
405 MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
406 setProperties( opts );
407
408 size_t N = fg.nrVars();
409
410 // check whether all vars in fg are binary
411 for( size_t i = 0; i < N; i++ )
412 if( (fg.var(i).states() > 2) ) {
413 supported = false;
414 break;
415 }
416 if( !supported )
417 DAI_THROWE(NOT_IMPLEMENTED,"MR only supports binary variables");
418
419 // check whether all interactions are pairwise or single
420 // and construct Markov graph
421 G = GraphAL(N);
422 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
423 const Factor &psi = fg.factor(I);
424 if( psi.vars().size() > 2 ) {
425 supported = false;
426 break;
427 } else if( psi.vars().size() == 2 ) {
428 VarSet::const_iterator jit = psi.vars().begin();
429 size_t i = fg.findVar( *(jit) );
430 size_t j = fg.findVar( *(++jit) );
431 G.addEdge( i, j, false );
432 }
433 }
434 if( !supported )
435 DAI_THROWE(NOT_IMPLEMENTED,"MR does not support higher order interactions (only single and pairwise are supported)");
436
437 // construct theta
438 theta.clear();
439 theta.resize( N, 0.0 );
440
441 // construct tJ
442 tJ.resize( N );
443 for( size_t i = 0; i < N; i++ )
444 tJ[i].resize( G.nb(i).size(), 0.0 );
445
446 // initialize theta and tJ
447 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
448 const Factor &psi = fg.factor(I);
449 if( psi.vars().size() == 1 ) {
450 size_t i = fg.findVar( *(psi.vars().begin()) );
451 theta[i] += 0.5 * log(psi[1] / psi[0]);
452 } else if( psi.vars().size() == 2 ) {
453 VarSet::const_iterator jit = psi.vars().begin();
454 size_t i = fg.findVar( *(jit) );
455 size_t j = fg.findVar( *(++jit) );
456
457 Real w_ij = 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
458 tJ[i][G.findNb(i,j)] += w_ij;
459 tJ[j][G.findNb(j,i)] += w_ij;
460
461 theta[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
462 theta[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
463 }
464 }
465 for( size_t i = 0; i < N; i++ )
466 bforeach( const Neighbor &j, G.nb(i) )
467 tJ[i][j.iter] = tanh( tJ[i][j.iter] );
468
469 // construct M
470 M.resize( N );
471 for( size_t i = 0; i < N; i++ )
472 M[i].resize( G.nb(i).size() );
473
474 // construct cors
475 cors.resize( N );
476 for( size_t i = 0; i < N; i++ )
477 cors[i].resize( G.nb(i).size() );
478 for( size_t i = 0; i < N; i++ )
479 for( size_t _j = 0; _j < cors[i].size(); _j++ )
480 cors[i][_j].resize( G.nb(i).size() );
481
482 // construct Mag
483 Mag.resize( N );
484 }
485
486
487 } // end of namespace dai
488
489
490 #endif