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 bool MR::checkProperties() {
43 if( !HasProperty("updates") )
45 if( !HasProperty("inits") )
47 if( !HasProperty("verbose") )
49 if( !HasProperty("tol") )
52 ConvertPropertyTo
<UpdateType
>("updates");
53 ConvertPropertyTo
<InitType
>("inits");
54 ConvertPropertyTo
<size_t>("verbose");
55 ConvertPropertyTo
<double>("tol");
61 // init N, con, nb, tJ, theta
62 void MR::init(size_t Nin
, double *_w
, double *_th
) {
75 if( _w
[i
*N
+j
] != 0.0 ) {
77 tJ
[i
][con
[i
]] = tanh(_w
[i
*N
+j
]);
89 double MR::init_cor_resp() {
90 size_t j
,k
,l
, runx
,i2
;
91 double variab1
, variab2
;
93 double thbJsite
[kmax
];
104 vector
<vector
<double> > tJ_org
;
105 vector
<vector
<size_t> > nb_org
;
106 vector
<size_t> con_org
;
107 vector
<double> theta_org
;
109 vector
<double> xfield(N
*kmax
,0.0);
110 vector
<double> rfield(N
*kmax
,0.0);
111 vector
<double> Hfield(N
,0.0);
112 vector
<double> devs(N
*kmax
,0.0);
113 vector
<double> devs2(N
*kmax
,0.0);
114 vector
<double> dev(N
,0.0);
115 vector
<double> avmag(N
,0.0);
117 // save original tJ, nb
124 for(cavity
=0; cavity
<N
; cavity
++){ // for each spin to be removed
131 // Adapt the graph variables nb[], tJ[] and con[]
132 for(size_t i
=0; i
<con
[cavity
]; i
++) {
133 size_t ij
= nb
[cavity
][i
];
137 if(nb
[ij
][j
]==cavity
){
138 while(j
<(con
[ij
]-1)){
139 nb
[ij
][j
]=nb
[ij
][j
+1];
140 tJ
[ij
][j
] = tJ
[ij
][j
+1];
148 for(size_t i
=0; i
<con
[cavity
]; i
++)
149 con
[nb
[cavity
][i
]]--;
153 // Do everything starting from the new graph********
158 for(size_t i
=0; i
<kmax
*N
; i
++)
159 xfield
[i
] = 3.0*(2*rnd_uniform()-1.);
161 for(i2
=0; i2
<concav
; i2
++){ // Subsequently apply a field to each cavity spin ****************
163 s2
= nb
[cavity
][i2
]; // identify the index of the cavity spin
164 for(size_t i
=0; i
<con
[s2
]; i
++)
165 rfield
[kmax
*s2
+i
] = 1.;
168 do { // From here start the response and belief propagation
173 for(size_t i
=0; i
<con
[k
]; i
++)
174 thbJsite
[i
] = tJ
[k
][i
];
175 for(l
=0; l
<con
[k
]; l
++){
178 if(k
==s2
) rinter
+= 1.;
179 for(j
=0; j
<con
[k
]; j
++)
181 variab2
= tanh(xfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]);
182 variab1
= thbJsite
[j
]*variab2
;
183 xinter
*= (1.+variab1
)/(1.-variab1
);
185 rinter
+= thbJsite
[j
]*rfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
188 variab1
= 0.5*log(xinter
);
189 xinter
= variab1
+ theta
[k
];
190 devs
[kmax
*k
+l
] = xinter
-xfield
[kmax
*k
+l
];
191 xfield
[kmax
*k
+l
] = xfield
[kmax
*k
+l
]+devs
[kmax
*k
+l
]*eps
;
192 if( fabs(devs
[kmax
*k
+l
]) > md
)
193 md
= fabs(devs
[kmax
*k
+l
]);
195 devs2
[kmax
*k
+l
] = rinter
-rfield
[kmax
*k
+l
];
196 rfield
[kmax
*k
+l
] = rfield
[kmax
*k
+l
]+devs2
[kmax
*k
+l
]*eps
;
197 if( fabs(devs2
[kmax
*k
+l
]) > md
)
198 md
= fabs(devs2
[kmax
*k
+l
]);
202 } while((md
> Tol())&&(runx
<runs
)); // Precision condition reached -> BP and RP finished
205 cout
<< "init_cor_resp: Convergence not reached (md=" << md
<< ")..." << endl
;
209 // compute the observables (i.e. magnetizations and responses)******
211 for(size_t i
=0; i
<concav
; i
++){
215 for(j
=0; j
<con
[nb
[cavity
][i
]]; j
++){
216 variab2
= tanh(xfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]);
217 variab1
= tJ
[nb
[cavity
][i
]][j
]*variab2
;
218 rinter
+= tJ
[nb
[cavity
][i
]][j
]*rfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
219 xinter
*= (1.+variab1
)/(1.-variab1
);
221 xinter
= tanh(0.5*log(xinter
)+theta
[nb
[cavity
][i
]]);
222 res
[i
] = rinter
*(1-xinter
*xinter
);
225 // *******************
227 for(size_t i
=0; i
<concav
; i
++)
228 if(nb
[cavity
][i
]!=s2
)
230 cors
[cavity
][i2
][i
] = res
[i
];
232 cors
[cavity
][i2
][i
] = 0;
233 } // close for i2 = 0...concav
236 // restore nb, tJ, con
246 double MR::T(size_t i
, sub_nb A
) {
247 // i is a variable index
248 // A is a subset of nb[i]
250 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
252 sub_nb
_nbi_min_A(con
[i
]);
253 for( size_t __j
= 0; __j
< A
.size(); __j
++ )
254 _nbi_min_A
-= A
[__j
];
256 double res
= theta
[i
];
257 for( size_t __j
= 0; __j
< _nbi_min_A
.size(); __j
++ ) {
258 size_t _j
= _nbi_min_A
[__j
];
259 res
+= atanh(tJ
[i
][_j
] * M
[i
][_j
]);
265 double MR::T(size_t i
, size_t _j
) {
273 double MR::Omega(size_t i
, size_t _j
, size_t _l
) {
278 double Tijl
= T(i
,jl
);
279 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
283 double MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
287 double Tij
= T(i
,jll
);
290 double Tijll
= T(i
,jll
);
292 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
);
296 double MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
302 double Till
= T(i
,ll
);
304 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
);
308 double MR::_tJ(size_t i
, sub_nb A
) {
309 // i is a variable index
310 // A is a subset of nb[i]
312 // calculate the product of all tJ[i][_j] for _j in A
314 size_t Asize
= A
.size();
317 // case 1: return tJ[i][A[0]];
318 // case 2: return tJ[i][A[0]] * tJ[i][A[1]];
319 // case 3: return tJ[i][A[0]] * tJ[i][A[1]] * tJ[i][A[2]];
321 size_t __j
= Asize
- 1;
324 return tJ
[i
][_j
] * _tJ(i
, A_j
);
330 double MR::appM(size_t i
, sub_nb A
) {
331 // i is a variable index
332 // A is a subset of nb[i]
334 // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
335 // defined as the sum over all partitions of A into subsets of cardinality two at most of the
336 // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
337 // entries of the partitions
339 size_t Asize
= A
.size();
342 // case 1: return M[i][A[0]];
343 // case 2: return M[i][A[0]] * M[i][A[1]] + cors[i][A[0]][A[1]];
344 // case 3: return M[i][A[0]] * M[i][A[1]] * M[i][A[2]] + M[i][A[0]] * cors[i][A[1]][A[2]] + M[i][A[1]] * cors[i][A[0]][A[2]] + M[i][A[2]] * cors[i][A[0]][A[1]];
346 size_t __j
= Asize
- 1;
350 double result
= M
[i
][_j
] * appM(i
, A_j
);
351 for( size_t __k
= 0; __k
< __j
; __k
++ ) {
353 result
+= cors
[i
][_j
][_k
] * appM(i
,A_j
- _k
);
361 void MR::sum_subs(size_t j
, sub_nb A
, double *sum_even
, double *sum_odd
) {
362 // j is a variable index
363 // A is a subset of nb[j]
365 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
372 do { // for all subsets of A
374 // construct subset B of A corresponding to S
377 size_t Ssize
= S
.size();
378 for( size_t bit
= 0; bit
< Ssize
; bit
++ )
382 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
384 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
387 } while( !S
.empty() );
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( Updates() == 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
-= kindex
[i
][_j
];
417 // find indices in nb[i] that do not correspond with j
418 sub_nb
_nbi_min_j(con
[i
]);
421 sum_subs(j
, _nbj_min_i
, &sum_even
, &sum_odd
);
422 newM
= (tanh(theta
[j
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[j
]) * sum_odd
);
424 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
425 double denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
427 for(size_t _k
=0; _k
<con
[i
]; _k
++) if(_k
!= _j
) {
428 sum_subs(i
, _nbi_min_j
- _k
, &sum_even
, &sum_odd
);
429 numer
+= tJ
[i
][_k
] * cors
[i
][_j
][_k
] * (tanh(theta
[i
]) * sum_even
+ sum_odd
);
431 newM
-= numer
/ denom
;
432 } else if( Updates() == UpdateType::LINEAR
) {
434 for(size_t _l
=0; _l
<con
[i
]; _l
++) if( _l
!= _j
)
435 newM
-= Omega(i
,_j
,_l
) * tJ
[i
][_l
] * cors
[i
][_j
][_l
];
436 for(size_t _l1
=0; _l1
<con
[j
]; _l1
++) if( _l1
!= _i
)
437 for( size_t _l2
=_l1
+1; _l2
<con
[j
]; _l2
++) if( _l2
!= _i
)
438 newM
+= Gamma(j
,_i
,_l1
,_l2
) * tJ
[j
][_l1
] * tJ
[j
][_l2
] * cors
[j
][_l1
][_l2
];
441 double dev
= newM
- M
[i
][_j
];
443 if( fabs(dev
) >= maxdev
)
446 newM
= M
[i
][_j
] + dev
;
447 if( fabs(newM
) > 1.0 )
452 } while((maxdev
>Tol())&&(run
<maxruns
));
454 updateMaxDiff( maxdev
);
458 cout
<< "solve_mcav: Convergence not reached (maxdev=" << maxdev
<< ")..." << endl
;
464 for(size_t i
=0; i
<N
; i
++) {
465 if( Updates() == UpdateType::FULL
) {
466 // find indices in nb[i]
469 // calc numerator1 and denominator1
470 double sum_even
, sum_odd
;
471 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
473 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
475 } else if( Updates() == UpdateType::LINEAR
) {
476 sub_nb
empty(con
[i
]);
480 for(size_t _l1
=0; _l1
<con
[i
]; _l1
++)
481 for( size_t _l2
=_l1
+1; _l2
<con
[i
]; _l2
++)
482 Mag
[i
] += Gamma(i
,_l1
,_l2
) * tJ
[i
][_l1
] * tJ
[i
][_l2
] * cors
[i
][_l1
][_l2
];
485 Mag
[i
] = sign(Mag
[i
]);
490 void MR::init_cor() {
491 for( size_t i
= 0; i
< nrVars(); i
++ ) {
492 vector
<Factor
> pairq
;
493 if( Inits() == InitType::CLAMPING
) {
494 BP
bpcav(*this, Properties()("updates",string("SEQMAX"))("tol", string("1e-9"))("maxiter", string("1000UL"))("verbose", string("0UL"))("logdomain", string("0")));
495 bpcav
.makeCavity( i
);
496 pairq
= calcPairBeliefs( bpcav
, delta(i
), false );
497 } else if( Inits() == InitType::EXACT
) {
498 JTree
jtcav(*this, Properties()("updates",string("HUGIN"))("verbose", string("0UL")) );
499 jtcav
.makeCavity( i
);
500 pairq
= calcPairBeliefs( jtcav
, delta(i
), false );
502 for( size_t jk
= 0; jk
< pairq
.size(); jk
++ ) {
503 VarSet::const_iterator kit
= pairq
[jk
].vars().begin();
504 size_t j
= findVar( *(kit
) );
505 size_t k
= findVar( *(++kit
) );
506 pairq
[jk
].normalize(Prob::NORMPROB
);
507 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]);
508 for( size_t _j
= 0; _j
< con
[i
]; _j
++ ) if( nb
[i
][_j
] == j
)
509 for( size_t _k
= 0; _k
< con
[i
]; _k
++ ) if( nb
[i
][_k
] == k
) {
510 cors
[i
][_j
][_k
] = cor
;
511 cors
[i
][_k
][_j
] = cor
;
518 string
MR::identify() const {
519 stringstream
result (stringstream::out
);
520 result
<< Name
<< GetProperties();
528 cout
<< "Starting " << identify() << "...";
531 // Diffs diffs(nrVars(), 1.0);
534 for(size_t i
=0; i
<N
; i
++)
538 for(size_t i
=0; i
<N
; i
++)
539 cors
[i
].resize(kmax
);
540 for(size_t i
=0; i
<N
; i
++)
541 for(size_t j
=0; j
<kmax
; j
++)
542 cors
[i
][j
].resize(kmax
);
545 for(size_t i
=0; i
<N
; i
++)
546 kindex
[i
].resize(kmax
);
548 if( Inits() == InitType::RESPPROP
)
549 updateMaxDiff( init_cor_resp() );
550 else if( Inits() == InitType::EXACT
)
551 init_cor(); // FIXME no MaxDiff() calculation
552 else if( Inits() == InitType::CLAMPING
)
553 init_cor(); // FIXME no MaxDiff() calculation
561 cout
<< "MR needed " << toc() - tic
<< " clocks." << endl
;
569 void MR::makekindex() {
570 for(size_t i
=0; i
<N
; i
++)
571 for(size_t j
=0; j
<con
[i
]; j
++) {
572 size_t ij
= nb
[i
][j
]; // ij is the j'th neighbour of spin i
574 while( nb
[ij
][k
] != i
)
576 kindex
[i
][j
] = k
; // the j'th neighbour of spin i has spin i as its k'th neighbour
581 Factor
MR::belief( const Var
&n
) const {
583 size_t i
= findVar( n
);
586 x
[0] = 0.5 - Mag
[i
] / 2.0;
587 x
[1] = 0.5 + Mag
[i
] / 2.0;
589 return Factor( n
, x
);
595 vector
<Factor
> MR::beliefs() const {
596 vector
<Factor
> result
;
597 for( size_t i
= 0; i
< nrVars(); i
++ )
598 result
.push_back( belief( var(i
) ) );
604 MR::MR( const FactorGraph
&fg
, const Properties
&opts
) : DAIAlgFG(fg
, opts
), supported(true) {
605 // check whether all vars in fg are binary
606 // check whether connectivity is <= kmax
607 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
608 if( (fg
.var(i
).states() > 2) || (fg
.delta(i
).size() > kmax
) ) {
616 // check whether all interactions are pairwise or single
617 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
618 if( fg
.factor(I
).vars().size() > 2 ) {
627 size_t Nin
= fg
.nrVars();
629 double *w
= new double[Nin
*Nin
];
630 double *th
= new double[Nin
];
632 for( size_t i
= 0; i
< Nin
; i
++ ) {
634 for( size_t j
= 0; j
< Nin
; j
++ )
638 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
639 const Factor
&psi
= fg
.factor(I
);
640 if( psi
.vars().size() == 1 ) {
641 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
642 th
[i
] += 0.5 * log(psi
[1] / psi
[0]);
643 } else if( psi
.vars().size() == 2 ) {
644 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
645 VarSet::const_iterator jit
= psi
.vars().begin();
646 size_t j
= fg
.findVar( *(++jit
) );
648 w
[i
*Nin
+j
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
649 w
[j
*Nin
+i
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
651 th
[i
] += 0.25 * log(psi
[3] / psi
[2] * psi
[1] / psi
[0]);
652 th
[j
] += 0.25 * log(psi
[3] / psi
[1] * psi
[2] / psi
[0]);
663 } // end of namespace dai