1 /* This file is part of libDAI - http://www.libdai.org/
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
7 * Copyright (C) 2006-2009 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
18 #include <dai/alldai.h>
27 const char *LC::Name
= "LC";
30 void LC::setProperties( const PropertySet
&opts
) {
31 DAI_ASSERT( opts
.hasKey("tol") );
32 DAI_ASSERT( opts
.hasKey("maxiter") );
33 DAI_ASSERT( opts
.hasKey("verbose") );
34 DAI_ASSERT( opts
.hasKey("cavity") );
35 DAI_ASSERT( opts
.hasKey("updates") );
37 props
.tol
= opts
.getStringAs
<Real
>("tol");
38 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
39 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
40 props
.cavity
= opts
.getStringAs
<Properties::CavityType
>("cavity");
41 props
.updates
= opts
.getStringAs
<Properties::UpdateType
>("updates");
42 if( opts
.hasKey("cavainame") )
43 props
.cavainame
= opts
.getStringAs
<string
>("cavainame");
44 if( opts
.hasKey("cavaiopts") )
45 props
.cavaiopts
= opts
.getStringAs
<PropertySet
>("cavaiopts");
46 if( opts
.hasKey("reinit") )
47 props
.reinit
= opts
.getStringAs
<bool>("reinit");
48 if( opts
.hasKey("damping") )
49 props
.damping
= opts
.getStringAs
<Real
>("damping");
55 PropertySet
LC::getProperties() const {
57 opts
.Set( "tol", props
.tol
);
58 opts
.Set( "maxiter", props
.maxiter
);
59 opts
.Set( "verbose", props
.verbose
);
60 opts
.Set( "cavity", props
.cavity
);
61 opts
.Set( "updates", props
.updates
);
62 opts
.Set( "cavainame", props
.cavainame
);
63 opts
.Set( "cavaiopts", props
.cavaiopts
);
64 opts
.Set( "reinit", props
.reinit
);
65 opts
.Set( "damping", props
.damping
);
70 string
LC::printProperties() const {
71 stringstream
s( stringstream::out
);
73 s
<< "tol=" << props
.tol
<< ",";
74 s
<< "maxiter=" << props
.maxiter
<< ",";
75 s
<< "verbose=" << props
.verbose
<< ",";
76 s
<< "cavity=" << props
.cavity
<< ",";
77 s
<< "updates=" << props
.updates
<< ",";
78 s
<< "cavainame=" << props
.cavainame
<< ",";
79 s
<< "cavaiopts=" << props
.cavaiopts
<< ",";
80 s
<< "reinit=" << props
.reinit
<< ",";
81 s
<< "damping=" << props
.damping
<< "]";
86 LC::LC( const FactorGraph
& fg
, const PropertySet
&opts
) : DAIAlgFG(fg
), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(0.0), _iters(0), props() {
87 setProperties( opts
);
90 _pancakes
.resize( nrVars() );
93 for( size_t i
=0; i
< nrVars(); i
++ )
94 _cavitydists
.push_back(Factor( delta(i
) ));
97 _phis
.reserve( nrVars() );
98 for( size_t i
= 0; i
< nrVars(); i
++ ) {
99 _phis
.push_back( vector
<Factor
>() );
100 _phis
[i
].reserve( nbV(i
).size() );
101 foreach( const Neighbor
&I
, nbV(i
) )
102 _phis
[i
].push_back( Factor( factor(I
).vars() / var(i
) ) );
106 _beliefs
.reserve( nrVars() );
107 for( size_t i
=0; i
< nrVars(); i
++ )
108 _beliefs
.push_back(Factor(var(i
)));
112 string
LC::identify() const {
113 return string(Name
) + printProperties();
117 void LC::CalcBelief (size_t i
) {
118 _beliefs
[i
] = _pancakes
[i
].marginal(var(i
));
122 Real
LC::CalcCavityDist (size_t i
, const std::string
&name
, const PropertySet
&opts
) {
126 if( props
.verbose
>= 2 )
127 cerr
<< "Initing cavity " << var(i
) << "(" << delta(i
).size() << " vars, " << delta(i
).nrStates() << " states)" << endl
;
129 if( props
.cavity
== Properties::CavityType::UNIFORM
)
130 Bi
= Factor(delta(i
));
132 InfAlg
*cav
= newInfAlg( name
, *this, opts
);
133 cav
->makeCavity( i
);
135 if( props
.cavity
== Properties::CavityType::FULL
)
136 Bi
= calcMarginal( *cav
, cav
->fg().delta(i
), props
.reinit
);
137 else if( props
.cavity
== Properties::CavityType::PAIR
) {
138 vector
<Factor
> pairbeliefs
= calcPairBeliefs( *cav
, cav
->fg().delta(i
), props
.reinit
, false );
139 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
140 Bi
*= pairbeliefs
[ij
];
141 } else if( props
.cavity
== Properties::CavityType::PAIR2
) {
142 vector
<Factor
> pairbeliefs
= calcPairBeliefs( *cav
, cav
->fg().delta(i
), props
.reinit
, true );
143 for( size_t ij
= 0; ij
< pairbeliefs
.size(); ij
++ )
144 Bi
*= pairbeliefs
[ij
];
146 maxdiff
= cav
->maxDiff();
150 _cavitydists
[i
] = Bi
;
156 Real
LC::InitCavityDists( const std::string
&name
, const PropertySet
&opts
) {
159 if( props
.verbose
>= 1 ) {
160 cerr
<< Name
<< "::InitCavityDists: ";
161 if( props
.cavity
== Properties::CavityType::UNIFORM
)
162 cerr
<< "Using uniform initial cavity distributions" << endl
;
163 else if( props
.cavity
== Properties::CavityType::FULL
)
164 cerr
<< "Using full " << name
<< opts
<< "...";
165 else if( props
.cavity
== Properties::CavityType::PAIR
)
166 cerr
<< "Using pairwise " << name
<< opts
<< "...";
167 else if( props
.cavity
== Properties::CavityType::PAIR2
)
168 cerr
<< "Using pairwise(new) " << name
<< opts
<< "...";
172 for( size_t i
= 0; i
< nrVars(); i
++ ) {
173 Real md
= CalcCavityDist(i
, name
, opts
);
178 if( props
.verbose
>= 1 ) {
179 cerr
<< Name
<< "::InitCavityDists used " << toc() - tic
<< " seconds." << endl
;
186 long LC::SetCavityDists( std::vector
<Factor
> &Q
) {
187 if( props
.verbose
>= 1 )
188 cerr
<< Name
<< "::SetCavityDists: Setting initial cavity distributions" << endl
;
189 if( Q
.size() != nrVars() )
191 for( size_t i
= 0; i
< nrVars(); i
++ ) {
192 if( _cavitydists
[i
].vars() != Q
[i
].vars() ) {
195 _cavitydists
[i
] = Q
[i
];
202 for( size_t i
= 0; i
< nrVars(); ++i
)
203 foreach( const Neighbor
&I
, nbV(i
) )
204 if( props
.updates
== Properties::UpdateType::SEQRND
)
205 _phis
[i
][I
.iter
].randomize();
207 _phis
[i
][I
.iter
].fill(1.0);
211 Factor
LC::NewPancake (size_t i
, size_t _I
, bool & hasNaNs
) {
212 size_t I
= nbV(i
)[_I
];
213 Factor piet
= _pancakes
[i
];
215 // recalculate _pancake[i]
216 VarSet Ivars
= factor(I
).vars();
218 for( VarSet::const_iterator k
= Ivars
.begin(); k
!= Ivars
.end(); k
++ )
220 A_I
*= (_pancakes
[findVar(*k
)] * factor(I
).inverse()).marginal( Ivars
/ var(i
), false );
221 if( Ivars
.size() > 1 )
222 A_I
^= (1.0 / (Ivars
.size() - 1));
223 Factor A_Ii
= (_pancakes
[i
] * factor(I
).inverse() * _phis
[i
][_I
].inverse()).marginal( Ivars
/ var(i
), false );
224 Factor quot
= A_I
/ A_Ii
;
225 if( props
.damping
!= 0.0 )
226 quot
= (quot
^(1.0 - props
.damping
)) * (_phis
[i
][_I
]^props
.damping
);
228 piet
*= quot
/ _phis
[i
][_I
].normalized();
229 _phis
[i
][_I
] = quot
.normalized();
233 if( piet
.hasNaNs() ) {
234 cerr
<< Name
<< "::NewPancake(" << i
<< ", " << _I
<< "): has NaNs!" << endl
;
243 if( props
.verbose
>= 1 )
244 cerr
<< "Starting " << identify() << "...";
245 if( props
.verbose
>= 2 )
249 vector
<Real
> diffs( nrVars(), INFINITY
);
250 Real maxDiff
= INFINITY
;
252 Real md
= InitCavityDists( props
.cavainame
, props
.cavaiopts
);
256 for( size_t i
= 0; i
< nrVars(); i
++ ) {
257 _pancakes
[i
] = _cavitydists
[i
];
259 foreach( const Neighbor
&I
, nbV(i
) ) {
260 _pancakes
[i
] *= factor(I
);
261 if( props
.updates
== Properties::UpdateType::SEQRND
)
262 _pancakes
[i
] *= _phis
[i
][I
.iter
];
265 _pancakes
[i
].normalize();
270 vector
<Factor
> old_beliefs
;
271 for(size_t i
=0; i
< nrVars(); i
++ )
272 old_beliefs
.push_back(beliefV(i
));
274 bool hasNaNs
= false;
275 for( size_t i
=0; i
< nrVars(); i
++ )
276 if( _pancakes
[i
].hasNaNs() ) {
281 cerr
<< Name
<< "::run: initial _pancakes has NaNs!" << endl
;
285 size_t nredges
= nrEdges();
286 vector
<Edge
> update_seq
;
287 update_seq
.reserve( nredges
);
288 for( size_t i
= 0; i
< nrVars(); ++i
)
289 foreach( const Neighbor
&I
, nbV(i
) )
290 update_seq
.push_back( Edge( i
, I
.iter
) );
292 // do several passes over the network until maximum number of iterations has
293 // been reached or until the maximum belief difference is smaller than tolerance
294 for( _iters
=0; _iters
< props
.maxiter
&& maxDiff
> props
.tol
; _iters
++ ) {
295 // Sequential updates
296 if( props
.updates
== Properties::UpdateType::SEQRND
)
297 random_shuffle( update_seq
.begin(), update_seq
.end() );
299 for( size_t t
=0; t
< nredges
; t
++ ) {
300 size_t i
= update_seq
[t
].first
;
301 size_t _I
= update_seq
[t
].second
;
302 _pancakes
[i
] = NewPancake( i
, _I
, hasNaNs
);
308 // compare new beliefs with old ones
309 for(size_t i
=0; i
< nrVars(); i
++ ) {
310 diffs
[i
] = dist( beliefV(i
), old_beliefs
[i
], Prob::DISTLINF
);
311 old_beliefs
[i
] = beliefV(i
);
313 maxDiff
= max( diffs
);
315 if( props
.verbose
>= 3 )
316 cerr
<< Name
<< "::run: maxdiff " << maxDiff
<< " after " << _iters
+1 << " passes" << endl
;
319 if( maxDiff
> _maxdiff
)
322 if( props
.verbose
>= 1 ) {
323 if( maxDiff
> props
.tol
) {
324 if( props
.verbose
== 1 )
326 cerr
<< Name
<< "::run: WARNING: not converged within " << props
.maxiter
<< " passes (" << toc() - tic
<< " seconds)...final maxdiff:" << maxDiff
<< endl
;
328 if( props
.verbose
>= 2 )
329 cerr
<< Name
<< "::run: ";
330 cerr
<< "converged in " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
338 } // end of namespace dai