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