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