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