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