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