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