c23a6a801dfc9bf63d7b5e41037f565861c31531
1 /* This file is part of libDAI - http://www.libdai.org/
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.
7 * Copyright (C) 2007 Bastian Wemmenhove
8 * Copyright (C) 2007-2010 Joris Mooij [joris dot mooij at libdai dot org]
9 * Copyright (C) 2007 Radboud University Nijmegen, The Netherlands
19 #include <dai/jtree.h>
29 const char *MR::Name
= "MR";
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") );
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");
45 PropertySet
MR::getProperties() const {
47 opts
.Set( "tol", props
.tol
);
48 opts
.Set( "verbose", props
.verbose
);
49 opts
.Set( "updates", props
.updates
);
50 opts
.Set( "inits", props
.inits
);
55 string
MR::printProperties() const {
56 stringstream
s( stringstream::out
);
58 s
<< "tol=" << props
.tol
<< ",";
59 s
<< "verbose=" << props
.verbose
<< ",";
60 s
<< "updates=" << props
.updates
<< ",";
61 s
<< "inits=" << props
.inits
<< "]";
66 void MR::init(size_t Nin
, Real
*_w
, Real
*_th
) {
79 if( _w
[i
*N
+j
] != 0.0 ) {
81 tJ
[i
][con
[i
]] = tanh(_w
[i
*N
+j
]);
92 Real
MR::init_cor_resp() {
93 size_t j
,k
,l
, runx
,i2
;
94 Real variab1
, variab2
;
107 vector
<vector
<Real
> > tJ_org
;
108 vector
<vector
<size_t> > nb_org
;
109 vector
<size_t> con_org
;
110 vector
<Real
> theta_org
;
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);
120 // save original tJ, nb
127 for(cavity
=0; cavity
<N
; cavity
++){ // for each spin to be removed
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
];
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];
151 for(size_t i
=0; i
<con
[cavity
]; i
++)
152 con
[nb
[cavity
][i
]]--;
156 // Do everything starting from the new graph********
161 for(size_t i
=0; i
<kmax
*N
; i
++)
162 xfield
[i
] = 3.0*(2*rnd_uniform()-1.);
164 for(i2
=0; i2
<concav
; i2
++){ // Subsequently apply a field to each cavity spin ****************
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.;
171 do { // From here start the response and belief propagation
176 for(size_t i
=0; i
<con
[k
]; i
++)
177 thbJsite
[i
] = tJ
[k
][i
];
178 for(l
=0; l
<con
[k
]; l
++){
181 if(k
==s2
) rinter
+= 1.;
182 for(j
=0; j
<con
[k
]; j
++)
184 variab2
= tanh(xfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]);
185 variab1
= thbJsite
[j
]*variab2
;
186 xinter
*= (1.+variab1
)/(1.-variab1
);
188 rinter
+= thbJsite
[j
]*rfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
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
]);
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
]);
205 } while((md
> props
.tol
)&&(runx
<runs
)); // Precision condition reached -> BP and RP finished
207 if( props
.verbose
>= 2 )
208 cerr
<< "init_cor_resp: Convergence not reached (md=" << md
<< ")..." << endl
;
212 // compute the observables (i.e. magnetizations and responses)******
214 for(size_t i
=0; i
<concav
; i
++){
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
);
224 xinter
= tanh(0.5*log(xinter
)+theta
[nb
[cavity
][i
]]);
225 res
[i
] = rinter
*(1-xinter
*xinter
);
228 // *******************
230 for(size_t i
=0; i
<concav
; i
++)
231 if(nb
[cavity
][i
]!=s2
)
233 cors
[cavity
][i2
][i
] = res
[i
];
235 cors
[cavity
][i2
][i
] = 0;
236 } // close for i2 = 0...concav
239 // restore nb, tJ, con
249 Real
MR::T(size_t i
, sub_nb A
) {
250 sub_nb
_nbi_min_A(con
[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
]);
262 Real
MR::T(size_t i
, size_t _j
) {
269 Real
MR::Omega(size_t i
, size_t _j
, size_t _l
) {
274 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
278 Real
MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
284 Real Tijll
= T(i
,jll
);
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
);
290 Real
MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
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
);
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
)
306 return tJ
[i
][_j
] * _tJ(i
, A
.reset(_j
));
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
)
315 sub_nb
A_j(A
); A_j
.reset(_j
);
317 Real result
= M
[i
][_j
] * appM(i
, A_j
);
318 for( size_t _k
= 0; _k
< A_j
.size(); _k
++ )
320 sub_nb
A_jk(A_j
); A_jk
.reset(_k
);
321 result
+= cors
[i
][_j
][_k
] * appM(i
,A_jk
);
329 void MR::sum_subs(size_t j
, sub_nb A
, Real
*sum_even
, Real
*sum_odd
) {
336 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
338 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
340 // calc next subset B
342 for( ; bit
< A
.size(); bit
++ )
355 void MR::solvemcav() {
356 Real sum_even
, sum_odd
;
358 size_t maxruns
= 1000;
361 for(size_t i
=0; i
<N
; i
++)
362 for(size_t _j
=0; _j
<con
[i
]; _j
++)
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
);
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
]);
380 _nbj_min_i
.reset(kindex
[i
][_j
]);
382 // find indices in nb[i] that do not correspond with j
383 sub_nb
_nbi_min_j(con
[i
]);
385 _nbi_min_j
.reset(_j
);
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
);
390 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
391 Real denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
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
);
399 newM
-= numer
/ denom
;
400 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
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
];
409 Real dev
= newM
- M
[i
][_j
];
411 if( fabs(dev
) >= maxdev
)
414 newM
= M
[i
][_j
] + dev
;
415 if( fabs(newM
) > 1.0 )
420 } while((maxdev
>props
.tol
)&&(run
<maxruns
));
423 if( maxdev
> _maxdiff
)
427 if( props
.verbose
>= 1 )
428 cerr
<< "solve_mcav: Convergence not reached (maxdev=" << maxdev
<< ")..." << endl
;
434 for(size_t i
=0; i
<N
; i
++) {
435 if( props
.updates
== Properties::UpdateType::FULL
) {
436 // find indices in nb[i]
440 // calc numerator1 and denominator1
441 Real sum_even
, sum_odd
;
442 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
444 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
446 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
447 sub_nb
empty(con
[i
]);
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
];
455 Mag
[i
] = sign(Mag
[i
]);
460 Real
MR::init_cor() {
462 for( size_t i
= 0; i
< nrVars(); i
++ ) {
463 vector
<Factor
> pairq
;
464 if( props
.inits
== Properties::InitType::CLAMPING
) {
465 BP
bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real
)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
466 bpcav
.makeCavity( i
);
467 pairq
= calcPairBeliefs( bpcav
, delta(i
), false, true );
468 md
= std::max( md
, bpcav
.maxDiff() );
469 } else if( props
.inits
== Properties::InitType::EXACT
) {
470 JTree
jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
471 jtcav
.makeCavity( i
);
472 pairq
= calcPairBeliefs( jtcav
, delta(i
), false, true );
474 for( size_t jk
= 0; jk
< pairq
.size(); jk
++ ) {
475 VarSet::const_iterator kit
= pairq
[jk
].vars().begin();
476 size_t j
= findVar( *(kit
) );
477 size_t k
= findVar( *(++kit
) );
478 pairq
[jk
].normalize();
479 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]);
480 for( size_t _j
= 0; _j
< con
[i
]; _j
++ ) if( nb
[i
][_j
] == j
)
481 for( size_t _k
= 0; _k
< con
[i
]; _k
++ ) if( nb
[i
][_k
] == k
) {
482 cors
[i
][_j
][_k
] = cor
;
483 cors
[i
][_k
][_j
] = cor
;
491 string
MR::identify() const {
492 return string(Name
) + printProperties();
498 if( props
.verbose
>= 1 )
499 cerr
<< "Starting " << identify() << "...";
504 for(size_t i
=0; i
<N
; i
++)
508 for(size_t i
=0; i
<N
; i
++)
509 cors
[i
].resize(kmax
);
510 for(size_t i
=0; i
<N
; i
++)
511 for(size_t j
=0; j
<kmax
; j
++)
512 cors
[i
][j
].resize(kmax
);
515 for(size_t i
=0; i
<N
; i
++)
516 kindex
[i
].resize(kmax
);
518 if( props
.inits
== Properties::InitType::RESPPROP
) {
519 Real md
= init_cor_resp();
522 } else if( props
.inits
== Properties::InitType::EXACT
) {
523 Real md
= init_cor();
527 else if( props
.inits
== Properties::InitType::CLAMPING
) {
528 Real md
= init_cor();
538 if( props
.verbose
>= 1 )
539 cerr
<< Name
<< " needed " << toc() - tic
<< " seconds." << endl
;
547 void MR::makekindex() {
548 for(size_t i
=0; i
<N
; i
++)
549 for(size_t j
=0; j
<con
[i
]; j
++) {
550 size_t ij
= nb
[i
][j
]; // ij is the j'th neighbour of spin i
552 while( nb
[ij
][k
] != i
)
554 kindex
[i
][j
] = k
; // the j'th neighbour of spin i has spin i as its k'th neighbour
559 Factor
MR::beliefV( size_t i
) const {
562 x
[0] = 0.5 - Mag
[i
] / 2.0;
563 x
[1] = 0.5 + Mag
[i
] / 2.0;
565 return Factor( var(i
), x
);
571 vector
<Factor
> MR::beliefs() const {
572 vector
<Factor
> result
;
573 for( size_t i
= 0; i
< nrVars(); i
++ )
574 result
.push_back( belief( var(i
) ) );
580 MR::MR( const FactorGraph
&fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), supported(true), _maxdiff(0.0), _iters(0) {
581 setProperties( opts
);
583 // check whether all vars in fg are binary
584 // check whether connectivity is <= kmax
585 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
586 if( (fg
.var(i
).states() > 2) || (fg
.delta(i
).size() > kmax
) ) {
592 DAI_THROWE(NOT_IMPLEMENTED
,"MR only supports binary variables with low connectivity");
594 // check whether all interactions are pairwise or single
595 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
596 if( fg
.factor(I
).vars().size() > 2 ) {
602 DAI_THROWE(NOT_IMPLEMENTED
,"MR does not support higher order interactions (only single and pairwise are supported)");
605 size_t Nin
= fg
.nrVars();
607 Real
*w
= new Real
[Nin
*Nin
];
608 Real
*th
= new Real
[Nin
];
610 for( size_t i
= 0; i
< Nin
; i
++ ) {
612 for( size_t j
= 0; j
< Nin
; j
++ )
616 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
617 const Factor
&psi
= fg
.factor(I
);
618 if( psi
.vars().size() == 1 ) {
619 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
620 th
[i
] += 0.5 * log(psi
[1] / psi
[0]);
621 } else if( psi
.vars().size() == 2 ) {
622 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
623 VarSet::const_iterator jit
= psi
.vars().begin();
624 size_t j
= fg
.findVar( *(++jit
) );
626 w
[i
*Nin
+j
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
627 w
[j
*Nin
+i
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
629 th
[i
] += 0.25 * log(psi
[3] / psi
[2] * psi
[1] / psi
[0]);
630 th
[j
] += 0.25 * log(psi
[3] / psi
[1] * psi
[2] / psi
[0]);
641 } // end of namespace dai