1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
4 This file is part of libDAI.
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
28 #include <dai/jtree.h>
30 #include <dai/diffs.h>
39 const char *MR::Name
= "MR";
42 void MR::setProperties( const PropertySet
&opts
) {
43 assert( opts
.hasKey("tol") );
44 assert( opts
.hasKey("verbose") );
45 assert( opts
.hasKey("updates") );
46 assert( opts
.hasKey("inits") );
48 props
.tol
= opts
.getStringAs
<double>("tol");
49 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
50 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
51 props
.inits
= opts
.getStringAs
<Properties::InitType
>("inits");
55 PropertySet
MR::getProperties() const {
57 opts
.Set( "tol", props
.tol
);
58 opts
.Set( "verbose", props
.verbose
);
59 opts
.Set( "updates", props
.updates
);
60 opts
.Set( "inits", props
.inits
);
65 string
MR::printProperties() const {
66 stringstream
s( stringstream::out
);
68 s
<< "tol=" << props
.tol
<< ",";
69 s
<< "verbose=" << props
.verbose
<< ",";
70 s
<< "updates=" << props
.updates
<< ",";
71 s
<< "inits=" << props
.inits
<< "]";
76 // init N, con, nb, tJ, theta
77 void MR::init(size_t Nin
, double *_w
, double *_th
) {
90 if( _w
[i
*N
+j
] != 0.0 ) {
92 tJ
[i
][con
[i
]] = tanh(_w
[i
*N
+j
]);
104 double MR::init_cor_resp() {
105 size_t j
,k
,l
, runx
,i2
;
106 double variab1
, variab2
;
108 double thbJsite
[kmax
];
119 vector
<vector
<double> > tJ_org
;
120 vector
<vector
<size_t> > nb_org
;
121 vector
<size_t> con_org
;
122 vector
<double> theta_org
;
124 vector
<double> xfield(N
*kmax
,0.0);
125 vector
<double> rfield(N
*kmax
,0.0);
126 vector
<double> Hfield(N
,0.0);
127 vector
<double> devs(N
*kmax
,0.0);
128 vector
<double> devs2(N
*kmax
,0.0);
129 vector
<double> dev(N
,0.0);
130 vector
<double> avmag(N
,0.0);
132 // save original tJ, nb
139 for(cavity
=0; cavity
<N
; cavity
++){ // for each spin to be removed
146 // Adapt the graph variables nb[], tJ[] and con[]
147 for(size_t i
=0; i
<con
[cavity
]; i
++) {
148 size_t ij
= nb
[cavity
][i
];
152 if(nb
[ij
][j
]==cavity
){
153 while(j
<(con
[ij
]-1)){
154 nb
[ij
][j
]=nb
[ij
][j
+1];
155 tJ
[ij
][j
] = tJ
[ij
][j
+1];
163 for(size_t i
=0; i
<con
[cavity
]; i
++)
164 con
[nb
[cavity
][i
]]--;
168 // Do everything starting from the new graph********
173 for(size_t i
=0; i
<kmax
*N
; i
++)
174 xfield
[i
] = 3.0*(2*rnd_uniform()-1.);
176 for(i2
=0; i2
<concav
; i2
++){ // Subsequently apply a field to each cavity spin ****************
178 s2
= nb
[cavity
][i2
]; // identify the index of the cavity spin
179 for(size_t i
=0; i
<con
[s2
]; i
++)
180 rfield
[kmax
*s2
+i
] = 1.;
183 do { // From here start the response and belief propagation
188 for(size_t i
=0; i
<con
[k
]; i
++)
189 thbJsite
[i
] = tJ
[k
][i
];
190 for(l
=0; l
<con
[k
]; l
++){
193 if(k
==s2
) rinter
+= 1.;
194 for(j
=0; j
<con
[k
]; j
++)
196 variab2
= tanh(xfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]);
197 variab1
= thbJsite
[j
]*variab2
;
198 xinter
*= (1.+variab1
)/(1.-variab1
);
200 rinter
+= thbJsite
[j
]*rfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
203 variab1
= 0.5*log(xinter
);
204 xinter
= variab1
+ theta
[k
];
205 devs
[kmax
*k
+l
] = xinter
-xfield
[kmax
*k
+l
];
206 xfield
[kmax
*k
+l
] = xfield
[kmax
*k
+l
]+devs
[kmax
*k
+l
]*eps
;
207 if( fabs(devs
[kmax
*k
+l
]) > md
)
208 md
= fabs(devs
[kmax
*k
+l
]);
210 devs2
[kmax
*k
+l
] = rinter
-rfield
[kmax
*k
+l
];
211 rfield
[kmax
*k
+l
] = rfield
[kmax
*k
+l
]+devs2
[kmax
*k
+l
]*eps
;
212 if( fabs(devs2
[kmax
*k
+l
]) > md
)
213 md
= fabs(devs2
[kmax
*k
+l
]);
217 } while((md
> props
.tol
)&&(runx
<runs
)); // Precision condition reached -> BP and RP finished
219 if( props
.verbose
>= 2 )
220 cout
<< "init_cor_resp: Convergence not reached (md=" << md
<< ")..." << endl
;
224 // compute the observables (i.e. magnetizations and responses)******
226 for(size_t i
=0; i
<concav
; i
++){
230 for(j
=0; j
<con
[nb
[cavity
][i
]]; j
++){
231 variab2
= tanh(xfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]);
232 variab1
= tJ
[nb
[cavity
][i
]][j
]*variab2
;
233 rinter
+= tJ
[nb
[cavity
][i
]][j
]*rfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
234 xinter
*= (1.+variab1
)/(1.-variab1
);
236 xinter
= tanh(0.5*log(xinter
)+theta
[nb
[cavity
][i
]]);
237 res
[i
] = rinter
*(1-xinter
*xinter
);
240 // *******************
242 for(size_t i
=0; i
<concav
; i
++)
243 if(nb
[cavity
][i
]!=s2
)
245 cors
[cavity
][i2
][i
] = res
[i
];
247 cors
[cavity
][i2
][i
] = 0;
248 } // close for i2 = 0...concav
251 // restore nb, tJ, con
261 double MR::T(size_t i
, sub_nb A
) {
262 // i is a variable index
263 // A is a subset of nb[i]
265 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
267 sub_nb
_nbi_min_A(con
[i
]);
271 double res
= theta
[i
];
272 for( size_t _j
= 0; _j
< _nbi_min_A
.size(); _j
++ )
273 if( _nbi_min_A
.test(_j
) )
274 res
+= atanh(tJ
[i
][_j
] * M
[i
][_j
]);
279 double MR::T(size_t i
, size_t _j
) {
286 double MR::Omega(size_t i
, size_t _j
, size_t _l
) {
290 double Tijl
= T(i
,jl
);
291 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
295 double MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
298 double Tij
= T(i
,jll
);
301 double Tijll
= T(i
,jll
);
303 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
);
307 double MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
312 double Till
= T(i
,ll
);
314 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
);
318 double MR::_tJ(size_t i
, sub_nb A
) {
319 // i is a variable index
320 // A is a subset of nb[i]
322 // calculate the product of all tJ[i][_j] for _j in A
324 sub_nb::size_type _j
= A
.find_first();
325 if( _j
== sub_nb::npos
)
328 return tJ
[i
][_j
] * _tJ(i
, A
.reset(_j
));
332 double MR::appM(size_t i
, sub_nb A
) {
333 // i is a variable index
334 // A is a subset of nb[i]
336 // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
337 // defined as the sum over all partitions of A into subsets of cardinality two at most of the
338 // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
339 // entries of the partitions
341 sub_nb::size_type _j
= A
.find_first();
342 if( _j
== sub_nb::npos
)
345 sub_nb
A_j(A
); A_j
.reset(_j
);
347 double result
= M
[i
][_j
] * appM(i
, A_j
);
348 for( size_t _k
= 0; _k
< A_j
.size(); _k
++ )
350 sub_nb
A_jk(A_j
); A_jk
.reset(_k
);
351 result
+= cors
[i
][_j
][_k
] * appM(i
,A_jk
);
359 void MR::sum_subs(size_t j
, sub_nb A
, double *sum_even
, double *sum_odd
) {
360 // j is a variable index
361 // A is a subset of nb[j]
363 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
371 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
373 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
375 // calc next subset B
377 for( ; bit
< A
.size(); bit
++ )
390 void MR::solvemcav() {
391 double sum_even
, sum_odd
;
393 size_t maxruns
= 1000;
396 for(size_t i
=0; i
<N
; i
++)
397 for(size_t _j
=0; _j
<con
[i
]; _j
++)
404 for(size_t i
=0; i
<N
; i
++){ // for all i
405 for(size_t _j
=0; _j
<con
[i
]; _j
++){ // for all j in N_i
406 size_t _i
= kindex
[i
][_j
];
407 size_t j
= nb
[i
][_j
];
408 assert( nb
[j
][_i
] == i
);
411 if( props
.updates
== Properties::UpdateType::FULL
) {
412 // find indices in nb[j] that do not correspond with i
413 sub_nb
_nbj_min_i(con
[j
]);
415 _nbj_min_i
.reset(kindex
[i
][_j
]);
417 // find indices in nb[i] that do not correspond with j
418 sub_nb
_nbi_min_j(con
[i
]);
420 _nbi_min_j
.reset(_j
);
422 sum_subs(j
, _nbj_min_i
, &sum_even
, &sum_odd
);
423 newM
= (tanh(theta
[j
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[j
]) * sum_odd
);
425 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
426 double denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
428 for(size_t _k
=0; _k
<con
[i
]; _k
++) if(_k
!= _j
) {
429 sub_nb
_nbi_min_jk(_nbi_min_j
);
430 _nbi_min_jk
.reset(_k
);
431 sum_subs(i
, _nbi_min_jk
, &sum_even
, &sum_odd
);
432 numer
+= tJ
[i
][_k
] * cors
[i
][_j
][_k
] * (tanh(theta
[i
]) * sum_even
+ sum_odd
);
434 newM
-= numer
/ denom
;
435 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
437 for(size_t _l
=0; _l
<con
[i
]; _l
++) if( _l
!= _j
)
438 newM
-= Omega(i
,_j
,_l
) * tJ
[i
][_l
] * cors
[i
][_j
][_l
];
439 for(size_t _l1
=0; _l1
<con
[j
]; _l1
++) if( _l1
!= _i
)
440 for( size_t _l2
=_l1
+1; _l2
<con
[j
]; _l2
++) if( _l2
!= _i
)
441 newM
+= Gamma(j
,_i
,_l1
,_l2
) * tJ
[j
][_l1
] * tJ
[j
][_l2
] * cors
[j
][_l1
][_l2
];
444 double dev
= newM
- M
[i
][_j
];
446 if( fabs(dev
) >= maxdev
)
449 newM
= M
[i
][_j
] + dev
;
450 if( fabs(newM
) > 1.0 )
455 } while((maxdev
>props
.tol
)&&(run
<maxruns
));
458 if( maxdev
> _maxdiff
)
462 if( props
.verbose
>= 1 )
463 cout
<< "solve_mcav: Convergence not reached (maxdev=" << maxdev
<< ")..." << endl
;
469 for(size_t i
=0; i
<N
; i
++) {
470 if( props
.updates
== Properties::UpdateType::FULL
) {
471 // find indices in nb[i]
475 // calc numerator1 and denominator1
476 double sum_even
, sum_odd
;
477 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
479 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
481 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
482 sub_nb
empty(con
[i
]);
485 for(size_t _l1
=0; _l1
<con
[i
]; _l1
++)
486 for( size_t _l2
=_l1
+1; _l2
<con
[i
]; _l2
++)
487 Mag
[i
] += Gamma(i
,_l1
,_l2
) * tJ
[i
][_l1
] * tJ
[i
][_l2
] * cors
[i
][_l1
][_l2
];
490 Mag
[i
] = sign(Mag
[i
]);
495 void MR::init_cor() {
496 for( size_t i
= 0; i
< nrVars(); i
++ ) {
497 vector
<Factor
> pairq
;
498 if( props
.inits
== Properties::InitType::CLAMPING
) {
499 BP
bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("10000UL"))("verbose", string("0UL"))("logdomain", string("0")));
500 bpcav
.makeCavity( i
);
501 pairq
= calcPairBeliefs( bpcav
, delta(i
), false );
502 } else if( props
.inits
== Properties::InitType::EXACT
) {
503 JTree
jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", string("0UL")) );
504 jtcav
.makeCavity( i
);
505 pairq
= calcPairBeliefs( jtcav
, delta(i
), false );
507 for( size_t jk
= 0; jk
< pairq
.size(); jk
++ ) {
508 VarSet::const_iterator kit
= pairq
[jk
].vars().begin();
509 size_t j
= findVar( *(kit
) );
510 size_t k
= findVar( *(++kit
) );
511 pairq
[jk
].normalize();
512 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]);
513 for( size_t _j
= 0; _j
< con
[i
]; _j
++ ) if( nb
[i
][_j
] == j
)
514 for( size_t _k
= 0; _k
< con
[i
]; _k
++ ) if( nb
[i
][_k
] == k
) {
515 cors
[i
][_j
][_k
] = cor
;
516 cors
[i
][_k
][_j
] = cor
;
523 string
MR::identify() const {
524 return string(Name
) + printProperties();
530 if( props
.verbose
>= 1 )
531 cout
<< "Starting " << identify() << "...";
534 // Diffs diffs(nrVars(), 1.0);
537 for(size_t i
=0; i
<N
; i
++)
541 for(size_t i
=0; i
<N
; i
++)
542 cors
[i
].resize(kmax
);
543 for(size_t i
=0; i
<N
; i
++)
544 for(size_t j
=0; j
<kmax
; j
++)
545 cors
[i
][j
].resize(kmax
);
548 for(size_t i
=0; i
<N
; i
++)
549 kindex
[i
].resize(kmax
);
551 if( props
.inits
== Properties::InitType::RESPPROP
) {
552 double md
= init_cor_resp();
555 } else if( props
.inits
== Properties::InitType::EXACT
)
556 init_cor(); // FIXME no MaxDiff() calculation
557 else if( props
.inits
== Properties::InitType::CLAMPING
)
558 init_cor(); // FIXME no MaxDiff() calculation
565 if( props
.verbose
>= 1 )
566 cout
<< Name
<< " needed " << toc() - tic
<< " seconds." << endl
;
574 void MR::makekindex() {
575 for(size_t i
=0; i
<N
; i
++)
576 for(size_t j
=0; j
<con
[i
]; j
++) {
577 size_t ij
= nb
[i
][j
]; // ij is the j'th neighbour of spin i
579 while( nb
[ij
][k
] != i
)
581 kindex
[i
][j
] = k
; // the j'th neighbour of spin i has spin i as its k'th neighbour
586 Factor
MR::belief( const Var
&n
) const {
588 size_t i
= findVar( n
);
591 x
[0] = 0.5 - Mag
[i
] / 2.0;
592 x
[1] = 0.5 + Mag
[i
] / 2.0;
594 return Factor( n
, x
);
600 vector
<Factor
> MR::beliefs() const {
601 vector
<Factor
> result
;
602 for( size_t i
= 0; i
< nrVars(); i
++ )
603 result
.push_back( belief( var(i
) ) );
609 MR::MR( const FactorGraph
&fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), supported(true), _maxdiff(0.0), _iters(0) {
610 setProperties( opts
);
612 // check whether all vars in fg are binary
613 // check whether connectivity is <= kmax
614 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
615 if( (fg
.var(i
).states() > 2) || (fg
.delta(i
).size() > kmax
) ) {
623 // check whether all interactions are pairwise or single
624 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
625 if( fg
.factor(I
).vars().size() > 2 ) {
634 size_t Nin
= fg
.nrVars();
636 double *w
= new double[Nin
*Nin
];
637 double *th
= new double[Nin
];
639 for( size_t i
= 0; i
< Nin
; i
++ ) {
641 for( size_t j
= 0; j
< Nin
; j
++ )
645 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
646 const Factor
&psi
= fg
.factor(I
);
647 if( psi
.vars().size() == 1 ) {
648 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
649 th
[i
] += 0.5 * log(psi
[1] / psi
[0]);
650 } else if( psi
.vars().size() == 2 ) {
651 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
652 VarSet::const_iterator jit
= psi
.vars().begin();
653 size_t j
= fg
.findVar( *(++jit
) );
655 w
[i
*Nin
+j
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
656 w
[j
*Nin
+i
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
658 th
[i
] += 0.25 * log(psi
[3] / psi
[2] * psi
[1] / psi
[0]);
659 th
[j
] += 0.25 * log(psi
[3] / psi
[1] * psi
[2] / psi
[0]);
670 } // end of namespace dai