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