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
5 This file is part of libDAI.
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.
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.
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
29 #include <dai/jtree.h>
31 #include <dai/diffs.h>
40 const char *MR::Name
= "MR";
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") );
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");
56 PropertySet
MR::getProperties() const {
58 opts
.Set( "tol", props
.tol
);
59 opts
.Set( "verbose", props
.verbose
);
60 opts
.Set( "updates", props
.updates
);
61 opts
.Set( "inits", props
.inits
);
66 string
MR::printProperties() const {
67 stringstream
s( stringstream::out
);
69 s
<< "tol=" << props
.tol
<< ",";
70 s
<< "verbose=" << props
.verbose
<< ",";
71 s
<< "updates=" << props
.updates
<< ",";
72 s
<< "inits=" << props
.inits
<< "]";
77 // init N, con, nb, tJ, theta
78 void MR::init(size_t Nin
, double *_w
, double *_th
) {
91 if( _w
[i
*N
+j
] != 0.0 ) {
93 tJ
[i
][con
[i
]] = tanh(_w
[i
*N
+j
]);
105 double MR::init_cor_resp() {
106 size_t j
,k
,l
, runx
,i2
;
107 double variab1
, variab2
;
109 double thbJsite
[kmax
];
120 vector
<vector
<double> > tJ_org
;
121 vector
<vector
<size_t> > nb_org
;
122 vector
<size_t> con_org
;
123 vector
<double> theta_org
;
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);
133 // save original tJ, nb
140 for(cavity
=0; cavity
<N
; cavity
++){ // for each spin to be removed
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
];
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];
164 for(size_t i
=0; i
<con
[cavity
]; i
++)
165 con
[nb
[cavity
][i
]]--;
169 // Do everything starting from the new graph********
174 for(size_t i
=0; i
<kmax
*N
; i
++)
175 xfield
[i
] = 3.0*(2*rnd_uniform()-1.);
177 for(i2
=0; i2
<concav
; i2
++){ // Subsequently apply a field to each cavity spin ****************
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.;
184 do { // From here start the response and belief propagation
189 for(size_t i
=0; i
<con
[k
]; i
++)
190 thbJsite
[i
] = tJ
[k
][i
];
191 for(l
=0; l
<con
[k
]; l
++){
194 if(k
==s2
) rinter
+= 1.;
195 for(j
=0; j
<con
[k
]; j
++)
197 variab2
= tanh(xfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]);
198 variab1
= thbJsite
[j
]*variab2
;
199 xinter
*= (1.+variab1
)/(1.-variab1
);
201 rinter
+= thbJsite
[j
]*rfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
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
]);
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
]);
218 } while((md
> props
.tol
)&&(runx
<runs
)); // Precision condition reached -> BP and RP finished
220 if( props
.verbose
>= 2 )
221 cout
<< "init_cor_resp: Convergence not reached (md=" << md
<< ")..." << endl
;
225 // compute the observables (i.e. magnetizations and responses)******
227 for(size_t i
=0; i
<concav
; i
++){
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
);
237 xinter
= tanh(0.5*log(xinter
)+theta
[nb
[cavity
][i
]]);
238 res
[i
] = rinter
*(1-xinter
*xinter
);
241 // *******************
243 for(size_t i
=0; i
<concav
; i
++)
244 if(nb
[cavity
][i
]!=s2
)
246 cors
[cavity
][i2
][i
] = res
[i
];
248 cors
[cavity
][i2
][i
] = 0;
249 } // close for i2 = 0...concav
252 // restore nb, tJ, con
262 double MR::T(size_t i
, sub_nb A
) {
263 // i is a variable index
264 // A is a subset of nb[i]
266 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
268 sub_nb
_nbi_min_A(con
[i
]);
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
]);
280 double MR::T(size_t i
, size_t _j
) {
287 double MR::Omega(size_t i
, size_t _j
, size_t _l
) {
291 double Tijl
= T(i
,jl
);
292 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
296 double MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
299 double Tij
= T(i
,jll
);
302 double Tijll
= T(i
,jll
);
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
);
308 double MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
313 double Till
= T(i
,ll
);
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
);
319 double MR::_tJ(size_t i
, sub_nb A
) {
320 // i is a variable index
321 // A is a subset of nb[i]
323 // calculate the product of all tJ[i][_j] for _j in A
325 sub_nb::size_type _j
= A
.find_first();
326 if( _j
== sub_nb::npos
)
329 return tJ
[i
][_j
] * _tJ(i
, A
.reset(_j
));
333 double MR::appM(size_t i
, sub_nb A
) {
334 // i is a variable index
335 // A is a subset of nb[i]
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
342 sub_nb::size_type _j
= A
.find_first();
343 if( _j
== sub_nb::npos
)
346 sub_nb
A_j(A
); A_j
.reset(_j
);
348 double result
= M
[i
][_j
] * appM(i
, A_j
);
349 for( size_t _k
= 0; _k
< A_j
.size(); _k
++ )
351 sub_nb
A_jk(A_j
); A_jk
.reset(_k
);
352 result
+= cors
[i
][_j
][_k
] * appM(i
,A_jk
);
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]
364 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
372 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
374 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
376 // calc next subset B
378 for( ; bit
< A
.size(); bit
++ )
391 void MR::solvemcav() {
392 double sum_even
, sum_odd
;
394 size_t maxruns
= 1000;
397 for(size_t i
=0; i
<N
; i
++)
398 for(size_t _j
=0; _j
<con
[i
]; _j
++)
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
);
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
]);
416 _nbj_min_i
.reset(kindex
[i
][_j
]);
418 // find indices in nb[i] that do not correspond with j
419 sub_nb
_nbi_min_j(con
[i
]);
421 _nbi_min_j
.reset(_j
);
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
);
426 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
427 double denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
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
);
435 newM
-= numer
/ denom
;
436 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
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
];
445 double dev
= newM
- M
[i
][_j
];
447 if( fabs(dev
) >= maxdev
)
450 newM
= M
[i
][_j
] + dev
;
451 if( fabs(newM
) > 1.0 )
456 } while((maxdev
>props
.tol
)&&(run
<maxruns
));
459 if( maxdev
> _maxdiff
)
463 if( props
.verbose
>= 1 )
464 cout
<< "solve_mcav: Convergence not reached (maxdev=" << maxdev
<< ")..." << endl
;
470 for(size_t i
=0; i
<N
; i
++) {
471 if( props
.updates
== Properties::UpdateType::FULL
) {
472 // find indices in nb[i]
476 // calc numerator1 and denominator1
477 double sum_even
, sum_odd
;
478 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
480 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
482 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
483 sub_nb
empty(con
[i
]);
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
];
491 Mag
[i
] = sign(Mag
[i
]);
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 );
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
;
524 string
MR::identify() const {
525 return string(Name
) + printProperties();
531 if( props
.verbose
>= 1 )
532 cout
<< "Starting " << identify() << "...";
535 // Diffs diffs(nrVars(), 1.0);
538 for(size_t i
=0; i
<N
; i
++)
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
);
549 for(size_t i
=0; i
<N
; i
++)
550 kindex
[i
].resize(kmax
);
552 if( props
.inits
== Properties::InitType::RESPPROP
) {
553 double md
= init_cor_resp();
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
566 if( props
.verbose
>= 1 )
567 cout
<< Name
<< " needed " << toc() - tic
<< " seconds." << endl
;
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
580 while( nb
[ij
][k
] != i
)
582 kindex
[i
][j
] = k
; // the j'th neighbour of spin i has spin i as its k'th neighbour
587 Factor
MR::belief( const Var
&n
) const {
589 size_t i
= findVar( n
);
592 x
[0] = 0.5 - Mag
[i
] / 2.0;
593 x
[1] = 0.5 + Mag
[i
] / 2.0;
595 return Factor( n
, x
);
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
) ) );
610 MR::MR( const FactorGraph
&fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), supported(true), _maxdiff(0.0), _iters(0) {
611 setProperties( opts
);
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
) ) {
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 ) {
635 size_t Nin
= fg
.nrVars();
637 double *w
= new double[Nin
*Nin
];
638 double *th
= new double[Nin
];
640 for( size_t i
= 0; i
< Nin
; i
++ ) {
642 for( size_t j
= 0; j
< Nin
; j
++ )
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
) );
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]));
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]);
671 } // end of namespace dai