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