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 // init N, con, nb, tJ, theta
66 void MR::init(size_t Nin
, double *_w
, double *_th
) {
79 if( _w
[i
*N
+j
] != 0.0 ) {
81 tJ
[i
][con
[i
]] = tanh(_w
[i
*N
+j
]);
93 double MR::init_cor_resp() {
94 size_t j
,k
,l
, runx
,i2
;
95 double variab1
, variab2
;
97 double thbJsite
[kmax
];
108 vector
<vector
<double> > tJ_org
;
109 vector
<vector
<size_t> > nb_org
;
110 vector
<size_t> con_org
;
111 vector
<double> theta_org
;
113 vector
<double> xfield(N
*kmax
,0.0);
114 vector
<double> rfield(N
*kmax
,0.0);
115 vector
<double> Hfield(N
,0.0);
116 vector
<double> devs(N
*kmax
,0.0);
117 vector
<double> devs2(N
*kmax
,0.0);
118 vector
<double> dev(N
,0.0);
119 vector
<double> avmag(N
,0.0);
121 // save original tJ, nb
128 for(cavity
=0; cavity
<N
; cavity
++){ // for each spin to be removed
135 // Adapt the graph variables nb[], tJ[] and con[]
136 for(size_t i
=0; i
<con
[cavity
]; i
++) {
137 size_t ij
= nb
[cavity
][i
];
141 if(nb
[ij
][j
]==cavity
){
142 while(j
<(con
[ij
]-1)){
143 nb
[ij
][j
]=nb
[ij
][j
+1];
144 tJ
[ij
][j
] = tJ
[ij
][j
+1];
152 for(size_t i
=0; i
<con
[cavity
]; i
++)
153 con
[nb
[cavity
][i
]]--;
157 // Do everything starting from the new graph********
162 for(size_t i
=0; i
<kmax
*N
; i
++)
163 xfield
[i
] = 3.0*(2*rnd_uniform()-1.);
165 for(i2
=0; i2
<concav
; i2
++){ // Subsequently apply a field to each cavity spin ****************
167 s2
= nb
[cavity
][i2
]; // identify the index of the cavity spin
168 for(size_t i
=0; i
<con
[s2
]; i
++)
169 rfield
[kmax
*s2
+i
] = 1.;
172 do { // From here start the response and belief propagation
177 for(size_t i
=0; i
<con
[k
]; i
++)
178 thbJsite
[i
] = tJ
[k
][i
];
179 for(l
=0; l
<con
[k
]; l
++){
182 if(k
==s2
) rinter
+= 1.;
183 for(j
=0; j
<con
[k
]; j
++)
185 variab2
= tanh(xfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]);
186 variab1
= thbJsite
[j
]*variab2
;
187 xinter
*= (1.+variab1
)/(1.-variab1
);
189 rinter
+= thbJsite
[j
]*rfield
[kmax
*nb
[k
][j
]+kindex
[k
][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
192 variab1
= 0.5*log(xinter
);
193 xinter
= variab1
+ theta
[k
];
194 devs
[kmax
*k
+l
] = xinter
-xfield
[kmax
*k
+l
];
195 xfield
[kmax
*k
+l
] = xfield
[kmax
*k
+l
]+devs
[kmax
*k
+l
]*eps
;
196 if( fabs(devs
[kmax
*k
+l
]) > md
)
197 md
= fabs(devs
[kmax
*k
+l
]);
199 devs2
[kmax
*k
+l
] = rinter
-rfield
[kmax
*k
+l
];
200 rfield
[kmax
*k
+l
] = rfield
[kmax
*k
+l
]+devs2
[kmax
*k
+l
]*eps
;
201 if( fabs(devs2
[kmax
*k
+l
]) > md
)
202 md
= fabs(devs2
[kmax
*k
+l
]);
206 } while((md
> props
.tol
)&&(runx
<runs
)); // Precision condition reached -> BP and RP finished
208 if( props
.verbose
>= 2 )
209 cout
<< "init_cor_resp: Convergence not reached (md=" << md
<< ")..." << endl
;
213 // compute the observables (i.e. magnetizations and responses)******
215 for(size_t i
=0; i
<concav
; i
++){
219 for(j
=0; j
<con
[nb
[cavity
][i
]]; j
++){
220 variab2
= tanh(xfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]);
221 variab1
= tJ
[nb
[cavity
][i
]][j
]*variab2
;
222 rinter
+= tJ
[nb
[cavity
][i
]][j
]*rfield
[kmax
*nb
[nb
[cavity
][i
]][j
]+kindex
[nb
[cavity
][i
]][j
]]*(1-variab2
*variab2
)/(1-variab1
*variab1
);
223 xinter
*= (1.+variab1
)/(1.-variab1
);
225 xinter
= tanh(0.5*log(xinter
)+theta
[nb
[cavity
][i
]]);
226 res
[i
] = rinter
*(1-xinter
*xinter
);
229 // *******************
231 for(size_t i
=0; i
<concav
; i
++)
232 if(nb
[cavity
][i
]!=s2
)
234 cors
[cavity
][i2
][i
] = res
[i
];
236 cors
[cavity
][i2
][i
] = 0;
237 } // close for i2 = 0...concav
240 // restore nb, tJ, con
250 double MR::T(size_t i
, sub_nb A
) {
251 // i is a variable index
252 // A is a subset of nb[i]
254 // calculate T{(i)}_A as defined in Rizzo&Montanari e-print (2.17)
256 sub_nb
_nbi_min_A(con
[i
]);
257 for( size_t __j
= 0; __j
< A
.size(); __j
++ )
258 _nbi_min_A
-= A
[__j
];
260 double res
= theta
[i
];
261 for( size_t __j
= 0; __j
< _nbi_min_A
.size(); __j
++ ) {
262 size_t _j
= _nbi_min_A
[__j
];
263 res
+= atanh(tJ
[i
][_j
] * M
[i
][_j
]);
269 double MR::T(size_t i
, size_t _j
) {
277 double MR::Omega(size_t i
, size_t _j
, size_t _l
) {
282 double Tijl
= T(i
,jl
);
283 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
287 double MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
291 double Tij
= T(i
,jll
);
294 double Tijll
= T(i
,jll
);
296 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
);
300 double MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
306 double Till
= T(i
,ll
);
308 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
);
312 double MR::_tJ(size_t i
, sub_nb A
) {
313 // i is a variable index
314 // A is a subset of nb[i]
316 // calculate the product of all tJ[i][_j] for _j in A
318 size_t Asize
= A
.size();
321 // case 1: return tJ[i][A[0]];
322 // case 2: return tJ[i][A[0]] * tJ[i][A[1]];
323 // case 3: return tJ[i][A[0]] * tJ[i][A[1]] * tJ[i][A[2]];
325 size_t __j
= Asize
- 1;
328 return tJ
[i
][_j
] * _tJ(i
, A_j
);
334 double MR::appM(size_t i
, sub_nb A
) {
335 // i is a variable index
336 // A is a subset of nb[i]
338 // calculate the moment of variables in A from M and cors, neglecting higher order cumulants,
339 // defined as the sum over all partitions of A into subsets of cardinality two at most of the
340 // product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
341 // entries of the partitions
343 size_t Asize
= A
.size();
346 // case 1: return M[i][A[0]];
347 // case 2: return M[i][A[0]] * M[i][A[1]] + cors[i][A[0]][A[1]];
348 // 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]];
350 size_t __j
= Asize
- 1;
354 double result
= M
[i
][_j
] * appM(i
, A_j
);
355 for( size_t __k
= 0; __k
< __j
; __k
++ ) {
357 result
+= cors
[i
][_j
][_k
] * appM(i
,A_j
- _k
);
365 void MR::sum_subs(size_t j
, sub_nb A
, double *sum_even
, double *sum_odd
) {
366 // j is a variable index
367 // A is a subset of nb[j]
369 // calculate sum over all even/odd subsets B of A of _tJ(j,B) appM(j,B)
376 do { // for all subsets of A
378 // construct subset B of A corresponding to S
381 size_t Ssize
= S
.size();
382 for( size_t bit
= 0; bit
< Ssize
; bit
++ )
386 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
388 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
391 } while( !S
.empty() );
395 void MR::solvemcav() {
396 double sum_even
, sum_odd
;
398 size_t maxruns
= 1000;
401 for(size_t i
=0; i
<N
; i
++)
402 for(size_t _j
=0; _j
<con
[i
]; _j
++)
409 for(size_t i
=0; i
<N
; i
++){ // for all i
410 for(size_t _j
=0; _j
<con
[i
]; _j
++){ // for all j in N_i
411 size_t _i
= kindex
[i
][_j
];
412 size_t j
= nb
[i
][_j
];
413 assert( nb
[j
][_i
] == i
);
416 if( props
.updates
== Properties::UpdateType::FULL
) {
417 // find indices in nb[j] that do not correspond with i
418 sub_nb
_nbj_min_i(con
[j
]);
419 _nbj_min_i
-= kindex
[i
][_j
];
421 // find indices in nb[i] that do not correspond with j
422 sub_nb
_nbi_min_j(con
[i
]);
425 sum_subs(j
, _nbj_min_i
, &sum_even
, &sum_odd
);
426 newM
= (tanh(theta
[j
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[j
]) * sum_odd
);
428 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
429 double denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
431 for(size_t _k
=0; _k
<con
[i
]; _k
++) if(_k
!= _j
) {
432 sum_subs(i
, _nbi_min_j
- _k
, &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
));
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]
474 // calc numerator1 and denominator1
475 double sum_even
, sum_odd
;
476 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
478 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
480 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
481 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("1000UL"))("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(Prob::NORMPROB
);
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 stringstream
result (stringstream::out
);
525 result
<< Name
<< getProperties();
532 if( props
.verbose
>= 1 )
533 cout
<< "Starting " << identify() << "...";
536 // Diffs diffs(nrVars(), 1.0);
539 for(size_t i
=0; i
<N
; i
++)
543 for(size_t i
=0; i
<N
; i
++)
544 cors
[i
].resize(kmax
);
545 for(size_t i
=0; i
<N
; i
++)
546 for(size_t j
=0; j
<kmax
; j
++)
547 cors
[i
][j
].resize(kmax
);
550 for(size_t i
=0; i
<N
; i
++)
551 kindex
[i
].resize(kmax
);
553 if( props
.inits
== Properties::InitType::RESPPROP
) {
554 double md
= init_cor_resp();
557 } else if( props
.inits
== Properties::InitType::EXACT
)
558 init_cor(); // FIXME no MaxDiff() calculation
559 else if( props
.inits
== Properties::InitType::CLAMPING
)
560 init_cor(); // FIXME no MaxDiff() calculation
567 if( props
.verbose
>= 1 )
568 cout
<< "MR needed " << toc() - tic
<< " clocks." << endl
;
576 void MR::makekindex() {
577 for(size_t i
=0; i
<N
; i
++)
578 for(size_t j
=0; j
<con
[i
]; j
++) {
579 size_t ij
= nb
[i
][j
]; // ij is the j'th neighbour of spin i
581 while( nb
[ij
][k
] != i
)
583 kindex
[i
][j
] = k
; // the j'th neighbour of spin i has spin i as its k'th neighbour
588 Factor
MR::belief( const Var
&n
) const {
590 size_t i
= findVar( n
);
593 x
[0] = 0.5 - Mag
[i
] / 2.0;
594 x
[1] = 0.5 + Mag
[i
] / 2.0;
596 return Factor( n
, x
);
602 vector
<Factor
> MR::beliefs() const {
603 vector
<Factor
> result
;
604 for( size_t i
= 0; i
< nrVars(); i
++ )
605 result
.push_back( belief( var(i
) ) );
611 MR::MR( const FactorGraph
&fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), supported(true), maxdiff(0.0) {
612 setProperties( opts
);
614 // check whether all vars in fg are binary
615 // check whether connectivity is <= kmax
616 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
617 if( (fg
.var(i
).states() > 2) || (fg
.delta(i
).size() > kmax
) ) {
625 // check whether all interactions are pairwise or single
626 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
627 if( fg
.factor(I
).vars().size() > 2 ) {
636 size_t Nin
= fg
.nrVars();
638 double *w
= new double[Nin
*Nin
];
639 double *th
= new double[Nin
];
641 for( size_t i
= 0; i
< Nin
; i
++ ) {
643 for( size_t j
= 0; j
< Nin
; j
++ )
647 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
648 const Factor
&psi
= fg
.factor(I
);
649 if( psi
.vars().size() == 1 ) {
650 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
651 th
[i
] += 0.5 * log(psi
[1] / psi
[0]);
652 } else if( psi
.vars().size() == 2 ) {
653 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
654 VarSet::const_iterator jit
= psi
.vars().begin();
655 size_t j
= fg
.findVar( *(++jit
) );
657 w
[i
*Nin
+j
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
658 w
[j
*Nin
+i
] += 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
660 th
[i
] += 0.25 * log(psi
[3] / psi
[2] * psi
[1] / psi
[0]);
661 th
[j
] += 0.25 * log(psi
[3] / psi
[1] * psi
[2] / psi
[0]);
672 } // end of namespace dai