1 /* This file is part of libDAI - http://www.libdai.org/
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
9 #include <dai/dai_config.h>
19 #include <dai/jtree.h>
30 void MR::setProperties( const PropertySet
&opts
) {
31 DAI_ASSERT( opts
.hasKey("tol") );
32 DAI_ASSERT( opts
.hasKey("updates") );
33 DAI_ASSERT( opts
.hasKey("inits") );
35 props
.tol
= opts
.getStringAs
<Real
>("tol");
36 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
37 props
.inits
= opts
.getStringAs
<Properties::InitType
>("inits");
38 if( opts
.hasKey("verbose") )
39 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
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 Real
MR::T(size_t i
, sub_nb A
) {
67 sub_nb
_nbi_min_A(G
.nb(i
).size());
72 for( size_t _j
= 0; _j
< _nbi_min_A
.size(); _j
++ )
73 if( _nbi_min_A
.test(_j
) )
74 res
+= atanh(tJ
[i
][_j
] * M
[i
][_j
]);
79 Real
MR::T(size_t i
, size_t _j
) {
80 sub_nb
j(G
.nb(i
).size());
86 Real
MR::Omega(size_t i
, size_t _j
, size_t _l
) {
87 sub_nb
jl(G
.nb(i
).size());
91 return Tijl
/ (1.0 + tJ
[i
][_l
] * M
[i
][_l
] * Tijl
);
95 Real
MR::Gamma(size_t i
, size_t _j
, size_t _l1
, size_t _l2
) {
96 sub_nb
jll(G
.nb(i
).size());
101 Real Tijll
= T(i
,jll
);
103 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
);
107 Real
MR::Gamma(size_t i
, size_t _l1
, size_t _l2
) {
108 sub_nb
ll(G
.nb(i
).size());
114 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
);
118 Real
MR::_tJ(size_t i
, sub_nb A
) {
119 sub_nb::size_type _j
= A
.find_first();
120 if( _j
== sub_nb::npos
)
123 return tJ
[i
][_j
] * _tJ(i
, A
.reset(_j
));
127 Real
MR::appM(size_t i
, sub_nb A
) {
128 sub_nb::size_type _j
= A
.find_first();
129 if( _j
== sub_nb::npos
)
132 sub_nb
A_j(A
); A_j
.reset(_j
);
134 Real result
= M
[i
][_j
] * appM(i
, A_j
);
135 for( size_t _k
= 0; _k
< A_j
.size(); _k
++ )
137 sub_nb
A_jk(A_j
); A_jk
.reset(_k
);
138 result
+= cors
[i
][_j
][_k
] * appM(i
,A_jk
);
146 void MR::sum_subs(size_t j
, sub_nb A
, Real
*sum_even
, Real
*sum_odd
) {
153 *sum_odd
+= _tJ(j
,B
) * appM(j
,B
);
155 *sum_even
+= _tJ(j
,B
) * appM(j
,B
);
157 // calc next subset B
159 for( ; bit
< A
.size(); bit
++ )
172 void MR::propagateCavityFields() {
173 Real sum_even
, sum_odd
;
175 size_t maxruns
= 1000;
177 for( size_t i
= 0; i
< G
.nrNodes(); i
++ )
178 bforeach( const Neighbor
&j
, G
.nb(i
) )
185 for( size_t i
= 0; i
< G
.nrNodes(); i
++ ) {
186 bforeach( const Neighbor
&j
, G
.nb(i
) ) {
188 size_t _i
= G
.findNb(j
,i
);
189 DAI_ASSERT( G
.nb(j
,_i
) == i
);
192 if( props
.updates
== Properties::UpdateType::FULL
) {
193 // find indices in nb(j) that do not correspond with i
194 sub_nb
_nbj_min_i(G
.nb(j
).size());
196 _nbj_min_i
.reset(_i
);
198 // find indices in nb(i) that do not correspond with j
199 sub_nb
_nbi_min_j(G
.nb(i
).size());
201 _nbi_min_j
.reset(_j
);
203 sum_subs(j
, _nbj_min_i
, &sum_even
, &sum_odd
);
204 newM
= (tanh(theta
[j
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[j
]) * sum_odd
);
206 sum_subs(i
, _nbi_min_j
, &sum_even
, &sum_odd
);
207 Real denom
= sum_even
+ tanh(theta
[i
]) * sum_odd
;
209 for(size_t _k
=0; _k
< G
.nb(i
).size(); _k
++) if(_k
!= _j
) {
210 sub_nb
_nbi_min_jk(_nbi_min_j
);
211 _nbi_min_jk
.reset(_k
);
212 sum_subs(i
, _nbi_min_jk
, &sum_even
, &sum_odd
);
213 numer
+= tJ
[i
][_k
] * cors
[i
][_j
][_k
] * (tanh(theta
[i
]) * sum_even
+ sum_odd
);
215 newM
-= numer
/ denom
;
216 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
218 for(size_t _l
=0; _l
<G
.nb(i
).size(); _l
++) if( _l
!= _j
)
219 newM
-= Omega(i
,_j
,_l
) * tJ
[i
][_l
] * cors
[i
][_j
][_l
];
220 for(size_t _l1
=0; _l1
<G
.nb(j
).size(); _l1
++) if( _l1
!= _i
)
221 for( size_t _l2
=_l1
+1; _l2
<G
.nb(j
).size(); _l2
++) if( _l2
!= _i
)
222 newM
+= Gamma(j
,_i
,_l1
,_l2
) * tJ
[j
][_l1
] * tJ
[j
][_l2
] * cors
[j
][_l1
][_l2
];
225 Real dev
= newM
- M
[i
][_j
];
227 if( abs(dev
) >= maxdev
)
230 newM
= M
[i
][_j
] + dev
;
231 if( abs(newM
) > 1.0 )
232 newM
= (newM
> 0.0) ? 1.0 : -1.0;
236 } while((maxdev
>props
.tol
)&&(run
<maxruns
));
239 if( maxdev
> _maxdiff
)
243 if( props
.verbose
>= 1 )
244 cerr
<< "MR::propagateCavityFields: Convergence not reached (maxdev=" << maxdev
<< ")..." << endl
;
249 void MR::calcMagnetizations() {
250 for( size_t i
= 0; i
< G
.nrNodes(); i
++ ) {
251 if( props
.updates
== Properties::UpdateType::FULL
) {
252 // find indices in nb(i)
253 sub_nb
_nbi( G
.nb(i
).size() );
256 // calc numerator1 and denominator1
257 Real sum_even
, sum_odd
;
258 sum_subs(i
, _nbi
, &sum_even
, &sum_odd
);
260 Mag
[i
] = (tanh(theta
[i
]) * sum_even
+ sum_odd
) / (sum_even
+ tanh(theta
[i
]) * sum_odd
);
262 } else if( props
.updates
== Properties::UpdateType::LINEAR
) {
263 sub_nb
empty( G
.nb(i
).size() );
266 for( size_t _l1
= 0; _l1
< G
.nb(i
).size(); _l1
++ )
267 for( size_t _l2
= _l1
+ 1; _l2
< G
.nb(i
).size(); _l2
++ )
268 Mag
[i
] += Gamma(i
,_l1
,_l2
) * tJ
[i
][_l1
] * tJ
[i
][_l2
] * cors
[i
][_l1
][_l2
];
270 if( abs( Mag
[i
] ) > 1.0 )
271 Mag
[i
] = (Mag
[i
] > 0.0) ? 1.0 : -1.0;
276 Real
MR::calcCavityCorrelations() {
278 for( size_t i
= 0; i
< nrVars(); i
++ ) {
279 vector
<Factor
> pairq
;
280 if( props
.inits
== Properties::InitType::EXACT
) {
281 JTree
jtcav(*this, PropertySet()("updates", string("HUGIN"))("verbose", (size_t)0) );
282 jtcav
.makeCavity( i
);
283 pairq
= calcPairBeliefs( jtcav
, delta(i
), false, true );
284 } else if( props
.inits
== Properties::InitType::CLAMPING
) {
285 BP
bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real
)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
286 bpcav
.makeCavity( i
);
288 pairq
= calcPairBeliefs( bpcav
, delta(i
), false, true );
289 md
= std::max( md
, bpcav
.maxDiff() );
290 } else if( props
.inits
== Properties::InitType::RESPPROP
) {
291 BP
bpcav(*this, PropertySet()("updates", string("SEQMAX"))("tol", (Real
)1.0e-9)("maxiter", (size_t)10000)("verbose", (size_t)0)("logdomain", false));
292 bpcav
.makeCavity( i
);
293 bpcav
.makeCavity( i
);
297 BBP
bbp( &bpcav
, PropertySet()("verbose",(size_t)0)("tol",(Real
)1.0e-9)("maxiter",(size_t)10000)("damping",(Real
)0.0)("updates",string("SEQ_MAX")) );
298 bforeach( const Neighbor
&j
, G
.nb(i
) ) {
299 // Create weights for magnetization of some spin
304 // BBP cost function would be the magnetization of spin j
306 b1_adj
.reserve( nrVars() );
307 for( size_t l
= 0; l
< nrVars(); l
++ )
309 b1_adj
.push_back( p
);
311 b1_adj
.push_back( Prob( 2, 0.0 ) );
312 bbp
.init_V( b1_adj
);
314 // run BBP to estimate adjoints
317 bforeach( const Neighbor
&k
, G
.nb(i
) ) {
319 cors
[i
][j
.iter
][k
.iter
] = (bbp
.adj_psi_V(k
)[1] - bbp
.adj_psi_V(k
)[0]);
321 cors
[i
][j
.iter
][k
.iter
] = 0.0;
326 if( props
.inits
!= Properties::InitType::RESPPROP
) {
327 for( size_t jk
= 0; jk
< pairq
.size(); jk
++ ) {
328 VarSet::const_iterator kit
= pairq
[jk
].vars().begin();
329 size_t j
= findVar( *(kit
) );
330 size_t k
= findVar( *(++kit
) );
331 pairq
[jk
].normalize();
332 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]);
334 size_t _j
= G
.findNb(i
,j
);
335 size_t _k
= G
.findNb(i
,k
);
336 cors
[i
][_j
][_k
] = cor
;
337 cors
[i
][_k
][_j
] = cor
;
348 if( props
.verbose
>= 1 )
349 cerr
<< "Starting " << identify() << "...";
353 // approximate correlations of cavity spins
354 Real md
= calcCavityCorrelations();
359 propagateCavityFields();
361 // calculate magnetizations
362 calcMagnetizations();
364 if( props
.verbose
>= 1 )
365 cerr
<< name() << " needed " << toc() - tic
<< " seconds." << endl
;
373 Factor
MR::beliefV( size_t i
) const {
376 x
[0] = 0.5 - Mag
[i
] / 2.0;
377 x
[1] = 0.5 + Mag
[i
] / 2.0;
379 return Factor( var(i
), x
);
385 Factor
MR::belief (const VarSet
&ns
) const {
388 else if( ns
.size() == 1 )
389 return beliefV( findVar( *(ns
.begin()) ) );
391 DAI_THROW(BELIEF_NOT_AVAILABLE
);
397 vector
<Factor
> MR::beliefs() const {
398 vector
<Factor
> result
;
399 for( size_t i
= 0; i
< nrVars(); i
++ )
400 result
.push_back( beliefV( i
) );
405 MR::MR( const FactorGraph
&fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), supported(true), _maxdiff(0.0), _iters(0) {
406 setProperties( opts
);
408 size_t N
= fg
.nrVars();
410 // check whether all vars in fg are binary
411 for( size_t i
= 0; i
< N
; i
++ )
412 if( (fg
.var(i
).states() > 2) ) {
417 DAI_THROWE(NOT_IMPLEMENTED
,"MR only supports binary variables");
419 // check whether all interactions are pairwise or single
420 // and construct Markov graph
422 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
423 const Factor
&psi
= fg
.factor(I
);
424 if( psi
.vars().size() > 2 ) {
427 } else if( psi
.vars().size() == 2 ) {
428 VarSet::const_iterator jit
= psi
.vars().begin();
429 size_t i
= fg
.findVar( *(jit
) );
430 size_t j
= fg
.findVar( *(++jit
) );
431 G
.addEdge( i
, j
, false );
435 DAI_THROWE(NOT_IMPLEMENTED
,"MR does not support higher order interactions (only single and pairwise are supported)");
439 theta
.resize( N
, 0.0 );
443 for( size_t i
= 0; i
< N
; i
++ )
444 tJ
[i
].resize( G
.nb(i
).size(), 0.0 );
446 // initialize theta and tJ
447 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
448 const Factor
&psi
= fg
.factor(I
);
449 if( psi
.vars().size() == 1 ) {
450 size_t i
= fg
.findVar( *(psi
.vars().begin()) );
451 theta
[i
] += 0.5 * log(psi
[1] / psi
[0]);
452 } else if( psi
.vars().size() == 2 ) {
453 VarSet::const_iterator jit
= psi
.vars().begin();
454 size_t i
= fg
.findVar( *(jit
) );
455 size_t j
= fg
.findVar( *(++jit
) );
457 Real w_ij
= 0.25 * log(psi
[3] * psi
[0] / (psi
[2] * psi
[1]));
458 tJ
[i
][G
.findNb(i
,j
)] += w_ij
;
459 tJ
[j
][G
.findNb(j
,i
)] += w_ij
;
461 theta
[i
] += 0.25 * log(psi
[3] / psi
[2] * psi
[1] / psi
[0]);
462 theta
[j
] += 0.25 * log(psi
[3] / psi
[1] * psi
[2] / psi
[0]);
465 for( size_t i
= 0; i
< N
; i
++ )
466 bforeach( const Neighbor
&j
, G
.nb(i
) )
467 tJ
[i
][j
.iter
] = tanh( tJ
[i
][j
.iter
] );
471 for( size_t i
= 0; i
< N
; i
++ )
472 M
[i
].resize( G
.nb(i
).size() );
476 for( size_t i
= 0; i
< N
; i
++ )
477 cors
[i
].resize( G
.nb(i
).size() );
478 for( size_t i
= 0; i
< N
; i
++ )
479 for( size_t _j
= 0; _j
< cors
[i
].size(); _j
++ )
480 cors
[i
][_j
].resize( G
.nb(i
).size() );
487 } // end of namespace dai