Updated copyrights
[libdai.git] / src / mr.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
2 Radboud University Nijmegen, The Netherlands /
3 Max Planck Institute for Biological Cybernetics, Germany
4
5 This file is part of libDAI.
6
7 libDAI is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 libDAI is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with libDAI; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22
23 #include <cstdio>
24 #include <ctime>
25 #include <cmath>
26 #include <cstdlib>
27 #include <dai/mr.h>
28 #include <dai/bp.h>
29 #include <dai/jtree.h>
30 #include <dai/util.h>
31 #include <dai/diffs.h>
32
33
34 namespace dai {
35
36
37 using namespace std;
38
39
40 const char *MR::Name = "MR";
41
42
43 void MR::setProperties( const PropertySet &opts ) {
44 assert( opts.hasKey("tol") );
45 assert( opts.hasKey("verbose") );
46 assert( opts.hasKey("updates") );
47 assert( opts.hasKey("inits") );
48
49 props.tol = opts.getStringAs<double>("tol");
50 props.verbose = opts.getStringAs<size_t>("verbose");
51 props.updates = opts.getStringAs<Properties::UpdateType>("updates");
52 props.inits = opts.getStringAs<Properties::InitType>("inits");
53 }
54
55
56 PropertySet MR::getProperties() const {
57 PropertySet opts;
58 opts.Set( "tol", props.tol );
59 opts.Set( "verbose", props.verbose );
60 opts.Set( "updates", props.updates );
61 opts.Set( "inits", props.inits );
62 return opts;
63 }
64
65
66 string MR::printProperties() const {
67 stringstream s( stringstream::out );
68 s << "[";
69 s << "tol=" << props.tol << ",";
70 s << "verbose=" << props.verbose << ",";
71 s << "updates=" << props.updates << ",";
72 s << "inits=" << props.inits << "]";
73 return s.str();
74 }
75
76
77 // init N, con, nb, tJ, theta
78 void MR::init(size_t Nin, double *_w, double *_th) {
79 size_t i,j;
80
81 N = Nin;
82
83 con.resize(N);
84 nb.resize(N);
85 tJ.resize(N);
86 for(i=0; i<N; i++ ) {
87 nb[i].resize(kmax);
88 tJ[i].resize(kmax);
89 con[i]=0;
90 for(j=0; j<N; j++ )
91 if( _w[i*N+j] != 0.0 ) {
92 nb[i][con[i]] = j;
93 tJ[i][con[i]] = tanh(_w[i*N+j]);
94 con[i]++;
95 }
96 }
97
98 theta.resize(N);
99 for(i=0; i<N; i++)
100 theta[i] = _th[i];
101 }
102
103
104 // calculate cors
105 double MR::init_cor_resp() {
106 size_t j,k,l, runx,i2;
107 double variab1, variab2;
108 double md, maxdev;
109 double thbJsite[kmax];
110 double xinter;
111 double rinter;
112 double res[kmax];
113 size_t s2;
114 size_t flag;
115 size_t concav;
116 size_t runs = 3000;
117 double eps = 0.2;
118 size_t cavity;
119
120 vector<vector<double> > tJ_org;
121 vector<vector<size_t> > nb_org;
122 vector<size_t> con_org;
123 vector<double> theta_org;
124
125 vector<double> xfield(N*kmax,0.0);
126 vector<double> rfield(N*kmax,0.0);
127 vector<double> Hfield(N,0.0);
128 vector<double> devs(N*kmax,0.0);
129 vector<double> devs2(N*kmax,0.0);
130 vector<double> dev(N,0.0);
131 vector<double> avmag(N,0.0);
132
133 // save original tJ, nb
134 nb_org = nb;
135 tJ_org = tJ;
136 con_org = con;
137 theta_org = theta;
138
139 maxdev = 0.0;
140 for(cavity=0; cavity<N; cavity++){ // for each spin to be removed
141 con = con_org;
142 concav=con[cavity];
143
144 nb = nb_org;
145 tJ = tJ_org;
146
147 // Adapt the graph variables nb[], tJ[] and con[]
148 for(size_t i=0; i<con[cavity]; i++) {
149 size_t ij = nb[cavity][i];
150 flag=0;
151 j=0;
152 do{
153 if(nb[ij][j]==cavity){
154 while(j<(con[ij]-1)){
155 nb[ij][j]=nb[ij][j+1];
156 tJ[ij][j] = tJ[ij][j+1];
157 j++;
158 }
159 flag=1;
160 }
161 j++;
162 } while(flag==0);
163 }
164 for(size_t i=0; i<con[cavity]; i++)
165 con[nb[cavity][i]]--;
166 con[cavity] = 0;
167
168
169 // Do everything starting from the new graph********
170
171 makekindex();
172 theta = theta_org;
173
174 for(size_t i=0; i<kmax*N; i++)
175 xfield[i] = 3.0*(2*rnd_uniform()-1.);
176
177 for(i2=0; i2<concav; i2++){ // Subsequently apply a field to each cavity spin ****************
178
179 s2 = nb[cavity][i2]; // identify the index of the cavity spin
180 for(size_t i=0; i<con[s2]; i++)
181 rfield[kmax*s2+i] = 1.;
182
183 runx=0;
184 do { // From here start the response and belief propagation
185 runx++;
186 md=0.0;
187 for(k=0; k<N; k++){
188 if(k!=cavity) {
189 for(size_t i=0; i<con[k]; i++)
190 thbJsite[i] = tJ[k][i];
191 for(l=0; l<con[k]; l++){
192 xinter = 1.;
193 rinter = 0.;
194 if(k==s2) rinter += 1.;
195 for(j=0; j<con[k]; j++)
196 if(j!=l){
197 variab2 = tanh(xfield[kmax*nb[k][j]+kindex[k][j]]);
198 variab1 = thbJsite[j]*variab2;
199 xinter *= (1.+variab1)/(1.-variab1);
200
201 rinter += thbJsite[j]*rfield[kmax*nb[k][j]+kindex[k][j]]*(1-variab2*variab2)/(1-variab1*variab1);
202 }
203
204 variab1 = 0.5*log(xinter);
205 xinter = variab1 + theta[k];
206 devs[kmax*k+l] = xinter-xfield[kmax*k+l];
207 xfield[kmax*k+l] = xfield[kmax*k+l]+devs[kmax*k+l]*eps;
208 if( fabs(devs[kmax*k+l]) > md )
209 md = fabs(devs[kmax*k+l]);
210
211 devs2[kmax*k+l] = rinter-rfield[kmax*k+l];
212 rfield[kmax*k+l] = rfield[kmax*k+l]+devs2[kmax*k+l]*eps;
213 if( fabs(devs2[kmax*k+l]) > md )
214 md = fabs(devs2[kmax*k+l]);
215 }
216 }
217 }
218 } while((md > props.tol)&&(runx<runs)); // Precision condition reached -> BP and RP finished
219 if(runx==runs)
220 if( props.verbose >= 2 )
221 cout << "init_cor_resp: Convergence not reached (md=" << md << ")..." << endl;
222 if(md > maxdev)
223 maxdev = md;
224
225 // compute the observables (i.e. magnetizations and responses)******
226
227 for(size_t i=0; i<concav; i++){
228 rinter = 0.;
229 xinter = 1.;
230 if(i!=i2)
231 for(j=0; j<con[nb[cavity][i]]; j++){
232 variab2 = tanh(xfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]);
233 variab1 = tJ[nb[cavity][i]][j]*variab2;
234 rinter += tJ[nb[cavity][i]][j]*rfield[kmax*nb[nb[cavity][i]][j]+kindex[nb[cavity][i]][j]]*(1-variab2*variab2)/(1-variab1*variab1);
235 xinter *= (1.+variab1)/(1.-variab1);
236 }
237 xinter = tanh(0.5*log(xinter)+theta[nb[cavity][i]]);
238 res[i] = rinter*(1-xinter*xinter);
239 }
240
241 // *******************
242
243 for(size_t i=0; i<concav; i++)
244 if(nb[cavity][i]!=s2)
245 // if(i!=i2)
246 cors[cavity][i2][i] = res[i];
247 else
248 cors[cavity][i2][i] = 0;
249 } // close for i2 = 0...concav
250 }
251
252 // restore nb, tJ, con
253 tJ = tJ_org;
254 nb = nb_org;
255 con = con_org;
256 theta = theta_org;
257
258 return maxdev;
259 }
260
261
262 double MR::T(size_t i, sub_nb A) {
263 // i is a variable index
264 // A is a subset of nb[i]
265 //
266 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
267
268 sub_nb _nbi_min_A(con[i]);
269 _nbi_min_A.set();
270 _nbi_min_A &= ~A;
271
272 double res = theta[i];
273 for( size_t _j = 0; _j < _nbi_min_A.size(); _j++ )
274 if( _nbi_min_A.test(_j) )
275 res += atanh(tJ[i][_j] * M[i][_j]);
276 return tanh(res);
277 }
278
279
280 double MR::T(size_t i, size_t _j) {
281 sub_nb j(con[i]);
282 j.set(_j);
283 return T(i,j);
284 }
285
286
287 double MR::Omega(size_t i, size_t _j, size_t _l) {
288 sub_nb jl(con[i]);
289 jl.set(_j);
290 jl.set(_l);
291 double Tijl = T(i,jl);
292 return Tijl / (1.0 + tJ[i][_l] * M[i][_l] * Tijl);
293 }
294
295
296 double MR::Gamma(size_t i, size_t _j, size_t _l1, size_t _l2) {
297 sub_nb jll(con[i]);
298 jll.set(_j);
299 double Tij = T(i,jll);
300 jll.set(_l1);
301 jll.set(_l2);
302 double Tijll = T(i,jll);
303
304 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);
305 }
306
307
308 double MR::Gamma(size_t i, size_t _l1, size_t _l2) {
309 sub_nb ll(con[i]);
310 double Ti = T(i,ll);
311 ll.set(_l1);
312 ll.set(_l2);
313 double Till = T(i,ll);
314
315 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);
316 }
317
318
319 double MR::_tJ(size_t i, sub_nb A) {
320 // i is a variable index
321 // A is a subset of nb[i]
322 //
323 // calculate the product of all tJ[i][_j] for _j in A
324
325 sub_nb::size_type _j = A.find_first();
326 if( _j == sub_nb::npos )
327 return 1.0;
328 else
329 return tJ[i][_j] * _tJ(i, A.reset(_j));
330 }
331
332
333 double MR::appM(size_t i, sub_nb A) {
334 // i is a variable index
335 // A is a subset of nb[i]
336 //
337 // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
338 // defined as the sum over all partitions of A into subsets of cardinality two at most of the
339 // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
340 // entries of the partitions
341
342 sub_nb::size_type _j = A.find_first();
343 if( _j == sub_nb::npos )
344 return 1.0;
345 else {
346 sub_nb A_j(A); A_j.reset(_j);
347
348 double result = M[i][_j] * appM(i, A_j);
349 for( size_t _k = 0; _k < A_j.size(); _k++ )
350 if( A_j.test(_k) ) {
351 sub_nb A_jk(A_j); A_jk.reset(_k);
352 result += cors[i][_j][_k] * appM(i,A_jk);
353 }
354
355 return result;
356 }
357 }
358
359
360 void MR::sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd) {
361 // j is a variable index
362 // A is a subset of nb[j]
363
364 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
365
366 *sum_even = 0.0;
367 *sum_odd = 0.0;
368
369 sub_nb B(A.size());
370 do {
371 if( B.count() % 2 )
372 *sum_odd += _tJ(j,B) * appM(j,B);
373 else
374 *sum_even += _tJ(j,B) * appM(j,B);
375
376 // calc next subset B
377 size_t bit = 0;
378 for( ; bit < A.size(); bit++ )
379 if( A.test(bit) ) {
380 if( B.test(bit) )
381 B.reset(bit);
382 else {
383 B.set(bit);
384 break;
385 }
386 }
387 } while (!B.none());
388 }
389
390
391 void MR::solvemcav() {
392 double sum_even, sum_odd;
393 double maxdev;
394 size_t maxruns = 1000;
395
396 makekindex();
397 for(size_t i=0; i<N; i++)
398 for(size_t _j=0; _j<con[i]; _j++)
399 M[i][_j]=0.1;
400
401 size_t run=0;
402 do {
403 maxdev=0.0;
404 run++;
405 for(size_t i=0; i<N; i++){ // for all i
406 for(size_t _j=0; _j<con[i]; _j++){ // for all j in N_i
407 size_t _i = kindex[i][_j];
408 size_t j = nb[i][_j];
409 assert( nb[j][_i] == i );
410
411 double newM = 0.0;
412 if( props.updates == Properties::UpdateType::FULL ) {
413 // find indices in nb[j] that do not correspond with i
414 sub_nb _nbj_min_i(con[j]);
415 _nbj_min_i.set();
416 _nbj_min_i.reset(kindex[i][_j]);
417
418 // find indices in nb[i] that do not correspond with j
419 sub_nb _nbi_min_j(con[i]);
420 _nbi_min_j.set();
421 _nbi_min_j.reset(_j);
422
423 sum_subs(j, _nbj_min_i, &sum_even, &sum_odd);
424 newM = (tanh(theta[j]) * sum_even + sum_odd) / (sum_even + tanh(theta[j]) * sum_odd);
425
426 sum_subs(i, _nbi_min_j, &sum_even, &sum_odd);
427 double denom = sum_even + tanh(theta[i]) * sum_odd;
428 double numer = 0.0;
429 for(size_t _k=0; _k<con[i]; _k++) if(_k != _j) {
430 sub_nb _nbi_min_jk(_nbi_min_j);
431 _nbi_min_jk.reset(_k);
432 sum_subs(i, _nbi_min_jk, &sum_even, &sum_odd);
433 numer += tJ[i][_k] * cors[i][_j][_k] * (tanh(theta[i]) * sum_even + sum_odd);
434 }
435 newM -= numer / denom;
436 } else if( props.updates == Properties::UpdateType::LINEAR ) {
437 newM = T(j,_i);
438 for(size_t _l=0; _l<con[i]; _l++) if( _l != _j )
439 newM -= Omega(i,_j,_l) * tJ[i][_l] * cors[i][_j][_l];
440 for(size_t _l1=0; _l1<con[j]; _l1++) if( _l1 != _i )
441 for( size_t _l2=_l1+1; _l2<con[j]; _l2++) if( _l2 != _i)
442 newM += Gamma(j,_i,_l1,_l2) * tJ[j][_l1] * tJ[j][_l2] * cors[j][_l1][_l2];
443 }
444
445 double dev = newM - M[i][_j];
446 // dev *= 0.02;
447 if( fabs(dev) >= maxdev )
448 maxdev = fabs(dev);
449
450 newM = M[i][_j] + dev;
451 if( fabs(newM) > 1.0 )
452 newM = sign(newM);
453 M[i][_j] = newM;
454 }
455 }
456 } while((maxdev>props.tol)&&(run<maxruns));
457
458 _iters = run;
459 if( maxdev > _maxdiff )
460 _maxdiff = maxdev;
461
462 if(run==maxruns){
463 if( props.verbose >= 1 )
464 cout << "solve_mcav: Convergence not reached (maxdev=" << maxdev << ")..." << endl;
465 }
466 }
467
468
469 void MR::solveM() {
470 for(size_t i=0; i<N; i++) {
471 if( props.updates == Properties::UpdateType::FULL ) {
472 // find indices in nb[i]
473 sub_nb _nbi(con[i]);
474 _nbi.set();
475
476 // calc numerator1 and denominator1
477 double sum_even, sum_odd;
478 sum_subs(i, _nbi, &sum_even, &sum_odd);
479
480 Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);
481
482 } else if( props.updates == Properties::UpdateType::LINEAR ) {
483 sub_nb empty(con[i]);
484 Mag[i] = T(i,empty);
485
486 for(size_t _l1=0; _l1<con[i]; _l1++)
487 for( size_t _l2=_l1+1; _l2<con[i]; _l2++)
488 Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];
489 }
490 if(fabs(Mag[i])>1.)
491 Mag[i] = sign(Mag[i]);
492 }
493 }
494
495
496 void MR::init_cor() {
497 for( size_t i = 0; i < nrVars(); i++ ) {
498 vector<Factor> pairq;
499 if( props.inits == Properties::InitType::CLAMPING ) {
500 BP bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("10000UL"))("verbose", string("0UL"))("logdomain", string("0")));
501 bpcav.makeCavity( i );
502 pairq = calcPairBeliefs( bpcav, delta(i), false );
503 } else if( props.inits == Properties::InitType::EXACT ) {
504 JTree jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", string("0UL")) );
505 jtcav.makeCavity( i );
506 pairq = calcPairBeliefs( jtcav, delta(i), false );
507 }
508 for( size_t jk = 0; jk < pairq.size(); jk++ ) {
509 VarSet::const_iterator kit = pairq[jk].vars().begin();
510 size_t j = findVar( *(kit) );
511 size_t k = findVar( *(++kit) );
512 pairq[jk].normalize();
513 double 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]);
514 for( size_t _j = 0; _j < con[i]; _j++ ) if( nb[i][_j] == j )
515 for( size_t _k = 0; _k < con[i]; _k++ ) if( nb[i][_k] == k ) {
516 cors[i][_j][_k] = cor;
517 cors[i][_k][_j] = cor;
518 }
519 }
520 }
521 }
522
523
524 string MR::identify() const {
525 return string(Name) + printProperties();
526 }
527
528
529 double MR::run() {
530 if( supported ) {
531 if( props.verbose >= 1 )
532 cout << "Starting " << identify() << "...";
533
534 double tic = toc();
535 // Diffs diffs(nrVars(), 1.0);
536
537 M.resize(N);
538 for(size_t i=0; i<N; i++)
539 M[i].resize(kmax);
540
541 cors.resize(N);
542 for(size_t i=0; i<N; i++)
543 cors[i].resize(kmax);
544 for(size_t i=0; i<N; i++)
545 for(size_t j=0; j<kmax; j++)
546 cors[i][j].resize(kmax);
547
548 kindex.resize(N);
549 for(size_t i=0; i<N; i++)
550 kindex[i].resize(kmax);
551
552 if( props.inits == Properties::InitType::RESPPROP ) {
553 double md = init_cor_resp();
554 if( md > _maxdiff )
555 _maxdiff = md;
556 } else if( props.inits == Properties::InitType::EXACT )
557 init_cor(); // FIXME no MaxDiff() calculation
558 else if( props.inits == Properties::InitType::CLAMPING )
559 init_cor(); // FIXME no MaxDiff() calculation
560
561 solvemcav();
562
563 Mag.resize(N);
564 solveM();
565
566 if( props.verbose >= 1 )
567 cout << Name << " needed " << toc() - tic << " seconds." << endl;
568
569 return 0.0;
570 } else
571 return 1.0;
572 }
573
574
575 void MR::makekindex() {
576 for(size_t i=0; i<N; i++)
577 for(size_t j=0; j<con[i]; j++) {
578 size_t ij = nb[i][j]; // ij is the j'th neighbour of spin i
579 size_t k=0;
580 while( nb[ij][k] != i )
581 k++;
582 kindex[i][j] = k; // the j'th neighbour of spin i has spin i as its k'th neighbour
583 }
584 }
585
586
587 Factor MR::belief( const Var &n ) const {
588 if( supported ) {
589 size_t i = findVar( n );
590
591 Real x[2];
592 x[0] = 0.5 - Mag[i] / 2.0;
593 x[1] = 0.5 + Mag[i] / 2.0;
594
595 return Factor( n, x );
596 } else
597 return Factor();
598 }
599
600
601 vector<Factor> MR::beliefs() const {
602 vector<Factor> result;
603 for( size_t i = 0; i < nrVars(); i++ )
604 result.push_back( belief( var(i) ) );
605 return result;
606 }
607
608
609
610 MR::MR( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), supported(true), _maxdiff(0.0), _iters(0) {
611 setProperties( opts );
612
613 // check whether all vars in fg are binary
614 // check whether connectivity is <= kmax
615 for( size_t i = 0; i < fg.nrVars(); i++ )
616 if( (fg.var(i).states() > 2) || (fg.delta(i).size() > kmax) ) {
617 supported = false;
618 break;
619 }
620
621 if( !supported )
622 return;
623
624 // check whether all interactions are pairwise or single
625 for( size_t I = 0; I < fg.nrFactors(); I++ )
626 if( fg.factor(I).vars().size() > 2 ) {
627 supported = false;
628 break;
629 }
630
631 if( !supported )
632 return;
633
634 // create w and th
635 size_t Nin = fg.nrVars();
636
637 double *w = new double[Nin*Nin];
638 double *th = new double[Nin];
639
640 for( size_t i = 0; i < Nin; i++ ) {
641 th[i] = 0.0;
642 for( size_t j = 0; j < Nin; j++ )
643 w[i*Nin+j] = 0.0;
644 }
645
646 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
647 const Factor &psi = fg.factor(I);
648 if( psi.vars().size() == 1 ) {
649 size_t i = fg.findVar( *(psi.vars().begin()) );
650 th[i] += 0.5 * log(psi[1] / psi[0]);
651 } else if( psi.vars().size() == 2 ) {
652 size_t i = fg.findVar( *(psi.vars().begin()) );
653 VarSet::const_iterator jit = psi.vars().begin();
654 size_t j = fg.findVar( *(++jit) );
655
656 w[i*Nin+j] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
657 w[j*Nin+i] += 0.25 * log(psi[3] * psi[0] / (psi[2] * psi[1]));
658
659 th[i] += 0.25 * log(psi[3] / psi[2] * psi[1] / psi[0]);
660 th[j] += 0.25 * log(psi[3] / psi[1] * psi[2] / psi[0]);
661 }
662 }
663
664 init(Nin, w, th);
665
666 delete th;
667 delete w;
668 }
669
670
671 } // end of namespace dai