Multiple changes: changes in build system, one workaround and one bug fix
[libdai.git] / src / lc.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <dai/dai_config.h>
10 #ifdef DAI_WITH_LC
11
12
13 #include <iostream>
14 #include <algorithm>
15 #include <map>
16 #include <set>
17 #include <dai/lc.h>
18 #include <dai/util.h>
19 #include <dai/alldai.h>
20
21
22 namespace dai {
23
24
25 using namespace std;
26
27
28 void LC::setProperties( const PropertySet &opts ) {
29 DAI_ASSERT( opts.hasKey("tol") );
30 DAI_ASSERT( opts.hasKey("maxiter") );
31 DAI_ASSERT( opts.hasKey("cavity") );
32 DAI_ASSERT( opts.hasKey("updates") );
33
34 props.tol = opts.getStringAs<Real>("tol");
35 props.maxiter = opts.getStringAs<size_t>("maxiter");
36 props.cavity = opts.getStringAs<Properties::CavityType>("cavity");
37 props.updates = opts.getStringAs<Properties::UpdateType>("updates");
38 if( opts.hasKey("verbose") )
39 props.verbose = opts.getStringAs<size_t>("verbose");
40 else
41 props.verbose = 0;
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");
50 else
51 props.damping = 0.0;
52 }
53
54
55 PropertySet LC::getProperties() const {
56 PropertySet opts;
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 );
66 return opts;
67 }
68
69
70 string LC::printProperties() const {
71 stringstream s( stringstream::out );
72 s << "[";
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 << "]";
82 return s.str();
83 }
84
85
86 LC::LC( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _pancakes(), _cavitydists(), _phis(), _beliefs(), _maxdiff(0.0), _iters(0), props() {
87 setProperties( opts );
88
89 // create pancakes
90 _pancakes.resize( nrVars() );
91
92 // create cavitydists
93 for( size_t i=0; i < nrVars(); i++ )
94 _cavitydists.push_back(Factor( delta(i) ));
95
96 // create phis
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 bforeach( const Neighbor &I, nbV(i) )
102 _phis[i].push_back( Factor( factor(I).vars() / var(i) ) );
103 }
104
105 // create beliefs
106 _beliefs.reserve( nrVars() );
107 for( size_t i=0; i < nrVars(); i++ )
108 _beliefs.push_back(Factor(var(i)));
109 }
110
111
112 void LC::CalcBelief (size_t i) {
113 _beliefs[i] = _pancakes[i].marginal(var(i));
114 }
115
116
117 Factor LC::belief (const VarSet &ns) const {
118 if( ns.size() == 0 )
119 return Factor();
120 else if( ns.size() == 1 )
121 return beliefV( findVar( *(ns.begin()) ) );
122 else {
123 DAI_THROW(BELIEF_NOT_AVAILABLE);
124 return Factor();
125 }
126 }
127
128
129 Real LC::CalcCavityDist (size_t i, const std::string &name, const PropertySet &opts) {
130 Factor Bi;
131 Real maxdiff = 0;
132
133 if( props.verbose >= 2 )
134 cerr << "Initing cavity " << var(i) << "(" << delta(i).size() << " vars, " << delta(i).nrStates() << " states)" << endl;
135
136 if( props.cavity == Properties::CavityType::UNIFORM )
137 Bi = Factor(delta(i));
138 else {
139 InfAlg *cav = newInfAlg( name, *this, opts );
140 cav->makeCavity( i );
141
142 if( props.cavity == Properties::CavityType::FULL )
143 Bi = calcMarginal( *cav, cav->fg().delta(i), props.reinit );
144 else if( props.cavity == Properties::CavityType::PAIR ) {
145 vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, false );
146 for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
147 Bi *= pairbeliefs[ij];
148 } else if( props.cavity == Properties::CavityType::PAIR2 ) {
149 vector<Factor> pairbeliefs = calcPairBeliefs( *cav, cav->fg().delta(i), props.reinit, true );
150 for( size_t ij = 0; ij < pairbeliefs.size(); ij++ )
151 Bi *= pairbeliefs[ij];
152 }
153 maxdiff = cav->maxDiff();
154 delete cav;
155 }
156 Bi.normalize();
157 _cavitydists[i] = Bi;
158
159 return maxdiff;
160 }
161
162
163 Real LC::InitCavityDists( const std::string &name, const PropertySet &opts ) {
164 double tic = toc();
165
166 if( props.verbose >= 1 ) {
167 cerr << this->name() << "::InitCavityDists: ";
168 if( props.cavity == Properties::CavityType::UNIFORM )
169 cerr << "Using uniform initial cavity distributions" << endl;
170 else if( props.cavity == Properties::CavityType::FULL )
171 cerr << "Using full " << name << opts << "...";
172 else if( props.cavity == Properties::CavityType::PAIR )
173 cerr << "Using pairwise " << name << opts << "...";
174 else if( props.cavity == Properties::CavityType::PAIR2 )
175 cerr << "Using pairwise(new) " << name << opts << "...";
176 }
177
178 Real maxdiff = 0.0;
179 for( size_t i = 0; i < nrVars(); i++ ) {
180 Real md = CalcCavityDist(i, name, opts);
181 if( md > maxdiff )
182 maxdiff = md;
183 }
184
185 if( props.verbose >= 1 ) {
186 cerr << this->name() << "::InitCavityDists used " << toc() - tic << " seconds." << endl;
187 }
188
189 return maxdiff;
190 }
191
192
193 long LC::SetCavityDists( std::vector<Factor> &Q ) {
194 if( props.verbose >= 1 )
195 cerr << name() << "::SetCavityDists: Setting initial cavity distributions" << endl;
196 if( Q.size() != nrVars() )
197 return -1;
198 for( size_t i = 0; i < nrVars(); i++ ) {
199 if( _cavitydists[i].vars() != Q[i].vars() ) {
200 return i+1;
201 } else
202 _cavitydists[i] = Q[i];
203 }
204 return 0;
205 }
206
207
208 void LC::init() {
209 for( size_t i = 0; i < nrVars(); ++i )
210 bforeach( const Neighbor &I, nbV(i) )
211 if( props.updates == Properties::UpdateType::SEQRND )
212 _phis[i][I.iter].randomize();
213 else
214 _phis[i][I.iter].fill(1.0);
215 }
216
217
218 Factor LC::NewPancake (size_t i, size_t _I, bool & hasNaNs) {
219 size_t I = nbV(i)[_I];
220 Factor piet = _pancakes[i];
221
222 // recalculate _pancake[i]
223 VarSet Ivars = factor(I).vars();
224 Factor A_I;
225 for( VarSet::const_iterator k = Ivars.begin(); k != Ivars.end(); k++ )
226 if( var(i) != *k )
227 A_I *= (_pancakes[findVar(*k)] * factor(I).inverse()).marginal( Ivars / var(i), false );
228 if( Ivars.size() > 1 )
229 A_I ^= (1.0 / (Ivars.size() - 1));
230 Factor A_Ii = (_pancakes[i] * factor(I).inverse() * _phis[i][_I].inverse()).marginal( Ivars / var(i), false );
231 Factor quot = A_I / A_Ii;
232 if( props.damping != 0.0 )
233 quot = (quot^(1.0 - props.damping)) * (_phis[i][_I]^props.damping);
234
235 piet *= quot / _phis[i][_I].normalized();
236 _phis[i][_I] = quot.normalized();
237
238 piet.normalize();
239
240 if( piet.hasNaNs() ) {
241 cerr << name() << "::NewPancake(" << i << ", " << _I << "): has NaNs!" << endl;
242 hasNaNs = true;
243 }
244
245 return piet;
246 }
247
248
249 Real LC::run() {
250 if( props.verbose >= 1 )
251 cerr << "Starting " << identify() << "...";
252 if( props.verbose >= 2 )
253 cerr << endl;
254
255 double tic = toc();
256
257 Real md = InitCavityDists( props.cavainame, props.cavaiopts );
258 if( md > _maxdiff )
259 _maxdiff = md;
260
261 for( size_t i = 0; i < nrVars(); i++ ) {
262 _pancakes[i] = _cavitydists[i];
263
264 bforeach( const Neighbor &I, nbV(i) ) {
265 _pancakes[i] *= factor(I);
266 if( props.updates == Properties::UpdateType::SEQRND )
267 _pancakes[i] *= _phis[i][I.iter];
268 }
269
270 _pancakes[i].normalize();
271
272 CalcBelief(i);
273 }
274
275 vector<Factor> oldBeliefsV;
276 for( size_t i = 0; i < nrVars(); i++ )
277 oldBeliefsV.push_back( beliefV(i) );
278
279 bool hasNaNs = false;
280 for( size_t i=0; i < nrVars(); i++ )
281 if( _pancakes[i].hasNaNs() ) {
282 hasNaNs = true;
283 break;
284 }
285 if( hasNaNs ) {
286 cerr << name() << "::run: initial _pancakes has NaNs!" << endl;
287 return 1.0;
288 }
289
290 size_t nredges = nrEdges();
291 vector<Edge> update_seq;
292 update_seq.reserve( nredges );
293 for( size_t i = 0; i < nrVars(); ++i )
294 bforeach( const Neighbor &I, nbV(i) )
295 update_seq.push_back( Edge( i, I.iter ) );
296
297 // do several passes over the network until maximum number of iterations has
298 // been reached or until the maximum belief difference is smaller than tolerance
299 Real maxDiff = INFINITY;
300 for( _iters = 0; _iters < props.maxiter && maxDiff > props.tol; _iters++ ) {
301 // Sequential updates
302 if( props.updates == Properties::UpdateType::SEQRND )
303 random_shuffle( update_seq.begin(), update_seq.end(), rnd );
304
305 for( size_t t=0; t < nredges; t++ ) {
306 size_t i = update_seq[t].first;
307 size_t _I = update_seq[t].second;
308 _pancakes[i] = NewPancake( i, _I, hasNaNs);
309 if( hasNaNs )
310 return 1.0;
311 CalcBelief( i );
312 }
313
314 // compare new beliefs with old ones
315 maxDiff = -INFINITY;
316 for( size_t i = 0; i < nrVars(); i++ ) {
317 maxDiff = std::max( maxDiff, dist( beliefV(i), oldBeliefsV[i], DISTLINF ) );
318 oldBeliefsV[i] = beliefV(i);
319 }
320
321 if( props.verbose >= 3 )
322 cerr << name() << "::run: maxdiff " << maxDiff << " after " << _iters+1 << " passes" << endl;
323 }
324
325 if( maxDiff > _maxdiff )
326 _maxdiff = maxDiff;
327
328 if( props.verbose >= 1 ) {
329 if( maxDiff > props.tol ) {
330 if( props.verbose == 1 )
331 cerr << endl;
332 cerr << name() << "::run: WARNING: not converged within " << props.maxiter << " passes (" << toc() - tic << " seconds)...final maxdiff:" << maxDiff << endl;
333 } else {
334 if( props.verbose >= 2 )
335 cerr << name() << "::run: ";
336 cerr << "converged in " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
337 }
338 }
339
340 return maxDiff;
341 }
342
343
344 } // end of namespace dai
345
346
347 #endif