Changed testregression such that it doesn't fail on platform-dependent numerical...
[libdai.git] / examples / uai2010-aie-solver.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <iostream>
10 #include <fstream>
11 #include <ostream>
12 #include <cstdlib>
13 #include <dai/alldai.h>
14 #include <dai/io.h>
15
16
17 using namespace std;
18 using namespace dai;
19
20
21 // Type for storing a joint state of all variables
22 typedef std::vector<size_t> stateVec;
23
24
25 struct PRbest {
26 Real value;
27 Real maxdiff;
28 bool ready;
29 PRbest() : value(0.0), maxdiff(INFINITY), ready(false) {}
30 };
31
32 struct MARbest {
33 vector<Factor> beliefs;
34 Real maxdiff;
35 bool ready;
36 MARbest() : beliefs(), maxdiff(INFINITY), ready(false) {}
37 };
38
39 struct MPEbest {
40 Real value;
41 stateVec state;
42 bool ready;
43 MPEbest() : value(-INFINITY), state(), ready(false) {}
44 };
45
46
47 int main( int argc, char *argv[] ) {
48 if ( argc != 5 ) {
49 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
50 cout << "It is one of the winning solvers that participated in the" << endl;
51 cout << "UAI 2010 Approximate Inference Challenge" << endl;
52 cout << "(see http://www.cs.huji.ac.il/project/UAI10/)" << endl << endl;
53 cout << "Usage: " << argv[0] << " <filename.uai> <filename.uai.evid> <seed> <task>" << endl << endl;
54 return 1;
55 } else {
56 double starttic = toc();
57
58 size_t verbose = 1; // verbosity
59 size_t ia_verbose = 0; // verbosity of inference algorithms
60 bool surgery = true; // change factor graph structure based on evidence
61 if( verbose )
62 cout << "Solver: " << argv[0] << endl;
63
64 // set random seed
65 size_t seed = fromString<size_t>( argv[3] );
66 rnd_seed( seed );
67 if( verbose )
68 cout << "Seed: " << seed << endl;
69
70 // check whether the task is valid
71 string task( argv[4] );
72 if( task != string("PR") && task != string("MPE") && task != string("MAR") )
73 DAI_THROWE(RUNTIME_ERROR,"Unknown task");
74 if( verbose )
75 cout << "Task: " << task << endl;
76
77 // output other command line options
78 if( verbose ) {
79 cout << "Factorgraph filename: " << argv[1] << endl;
80 cout << "Evidence filename: " << argv[2] << endl;
81 }
82
83 // get time and memory limits
84 char *buf = getenv( "UAI_TIME" );
85 double UAI_time = INFINITY;
86 if( buf != NULL )
87 UAI_time = fromString<double>( buf );
88 buf = getenv( "UAI_MEMORY" );
89 size_t UAI_memory = 0;
90 if( buf != NULL ) {
91 UAI_memory = fromString<double>( buf ) * 1024 * 1024 * 1024;
92 }
93 if( verbose ) {
94 cout << "Time limit: " << UAI_time << endl;
95 cout << "Memory limit: " << UAI_memory << endl;
96 }
97
98 // build output file name
99 vector<string> pathComponents = tokenizeString( string(argv[1]), true, "/" );
100 string outfile = pathComponents.back() + "." + task;
101 if( verbose )
102 cout << "Output filename: " << outfile << endl;
103
104 // open output stream
105 ofstream os;
106 os.open( outfile.c_str() );
107 if( !os.is_open() )
108 DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + outfile);
109 if( verbose )
110 cout << "Opened output stream" << endl;
111
112 // read factor graph
113 vector<Var> vars;
114 vector<Factor> facs0;
115 vector<Permute> permutations;
116 if( verbose )
117 cout << "Reading factor graph..." << endl;
118 ReadUaiAieFactorGraphFile( argv[1], verbose, vars, facs0, permutations );
119 if( verbose )
120 cout << "Successfully read factor graph" << endl;
121
122 // check if it could be a grid
123 bool couldBeGrid = true;
124 FactorGraph fg0( facs0.begin(), facs0.end(), vars.begin(), vars.end(), facs0.size(), vars.size() );
125 for( size_t i = 0; i < fg0.nrVars(); i++ )
126 if( fg0.delta(i).size() > 4 ) {
127 couldBeGrid = false;
128 break;
129 }
130 if( couldBeGrid )
131 for( size_t I = 0; I < fg0.nrFactors(); I++ )
132 if( fg0.factor(I).vars().size() > 2 ) {
133 couldBeGrid = false;
134 break;
135 }
136 if( verbose ) {
137 if( couldBeGrid )
138 cout << "This could be a grid!" << endl;
139 else
140 cout << "This cannot be a grid!" << endl;
141 }
142
143 // read evidence
144 if( verbose )
145 cout << "Reading evidence..." << endl;
146 vector<map<size_t,size_t> > evid = ReadUaiAieEvidenceFile( argv[2], verbose );
147 if( verbose )
148 cout << "Successfully read " << evid.size() << " evidence cases" << endl;
149
150 // write output header
151 if( verbose )
152 cout << " Writing header to file..." << endl;
153 os << task << endl;
154
155 // construct clamped factor graphs
156 if( verbose )
157 cout << "Constructing clamped factor graphs..." << endl;
158 vector<FactorGraph> fgs;
159 fgs.reserve( evid.size() );
160 for( size_t ev = 0; ev < evid.size(); ev++ ) {
161 if( verbose )
162 cout << " Evidence case " << ev << "..." << endl;
163 // copy vector of factors
164 vector<Factor> facs( facs0 );
165
166 // change factor graph to reflect observed evidence
167 if( verbose )
168 cout << " Applying evidence..." << endl;
169 if( surgery ) {
170 // replace factors with clamped variables with slices
171 for( size_t I = 0; I < facs.size(); I++ ) {
172 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ ) {
173 if( facs[I].vars() >> vars[e->first] ) {
174 if( verbose >= 2 )
175 cout << " Clamping " << e->first << " to value " << e->second << " in factor " << I << " = " << facs[I].vars() << endl;
176 facs[I] = facs[I].slice( vars[e->first], e->second );
177 if( verbose >= 2 )
178 cout << " ...remaining vars: " << facs[I].vars() << endl;
179 }
180 }
181 }
182 // remove empty factors
183 Real logZcorr = 0.0;
184 for( vector<Factor>::iterator I = facs.begin(); I != facs.end(); )
185 if( I->vars().size() == 0 ) {
186 logZcorr += std::log( (Real)(*I)[0] );
187 I = facs.erase( I );
188 } else
189 I++;
190 // multiply with logZcorr constant
191 if( facs.size() == 0 )
192 facs.push_back( Factor( VarSet(), std::exp(logZcorr) ) );
193 else
194 facs.front() *= std::exp(logZcorr);
195 }
196 // add delta factors corresponding to observed variable values
197 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ )
198 facs.push_back( createFactorDelta( vars[e->first], e->second ) );
199
200 // construct clamped factor graph
201 if( verbose )
202 cout << " Constructing factor graph..." << endl;
203 fgs.push_back( FactorGraph( facs.begin(), facs.end(), vars.begin(), vars.end(), facs.size(), vars.size() ) );
204 }
205
206 // variables for storing best results so far
207 vector<PRbest> bestPR( evid.size() );
208 vector<MARbest> bestMAR( evid.size() );
209 vector<MPEbest> bestMPE( evid.size() );
210 for( size_t ev = 0; ev < evid.size(); ev++ )
211 bestMPE[ev].state = stateVec( fgs[ev].nrVars(), 0 );
212 vector<size_t> ev2go;
213 ev2go.reserve( evid.size() );
214 for( size_t ev = 0; ev < evid.size(); ev++ )
215 ev2go.push_back( ev );
216
217 // solve inference problems
218 if( verbose )
219 cout << "Solving inference problems..." << endl;
220 bool first = true;
221 size_t nrsolvers = 3;
222 vector<size_t> nrsubsolvers( nrsolvers );
223 nrsubsolvers[0] = 2;
224 nrsubsolvers[1] = 1;
225 nrsubsolvers[2] = 1;
226 double MPEdamping = 0.49;
227 // for each (sub)solver
228 for( size_t solver = 0; solver < nrsolvers; solver++ ) {
229 if( verbose )
230 cout << " Solver " << solver << endl;
231
232 // for each evidence case
233 size_t subsolver = 0;
234 for( long _ev = 0; _ev < (long)ev2go.size(); ) {
235 bool improved = false;
236 size_t ev = ev2go[_ev];
237 if( verbose )
238 cout << " Evidence case " << ev << ", subsolver = " << subsolver << "..." << endl;
239
240 // construct inference algorithm on clamped factor graph
241 if( verbose )
242 cout << " Constructing inference algorithm..." << endl;
243 InfAlg *ia = NULL;
244 double tic = toc();
245 bool failed = false;
246 try {
247 // construct
248 if( solver == 0 ) { // the quick one
249 double remtime = (UAI_time - (toc() - starttic)) * 0.9;
250 if( remtime < 1.0 )
251 remtime = 1.0 ;
252 double maxtime = remtime / (ev2go.size() - _ev);
253 if( verbose ) {
254 cout << " Past time: " << (toc() - starttic) << endl;
255 cout << " Remaining time:" << remtime << endl;
256 cout << " Allotted time: " << maxtime << endl;
257 }
258 string maxtimestr;
259 if( maxtime != INFINITY )
260 maxtimestr = ",maxtime=" + toString(maxtime);
261 // quick and dirty...
262 if( task == "MPE" )
263 ia = newInfAlgFromString( "BP[inference=MAXPROD,updates=SEQRND,logdomain=" + toString(subsolver) + ",tol=1e-9,maxiter=10000" + maxtimestr + ",damping=0.1,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
264 else {
265 if( couldBeGrid )
266 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9,maxiter=10000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
267 else
268 ia = newInfAlgFromString( "BP[inference=SUMPROD,updates=SEQRND,logdomain=" + toString(subsolver) + ",tol=1e-9,maxiter=10000" + maxtimestr + ",damping=0.0,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
269 }
270 } else if( solver == 1 ) { // the exact one
271 string maxmemstr;
272 if( UAI_memory != 0 )
273 maxmemstr = ",maxmem=" + toString(UAI_memory);
274 if( task == "MPE" )
275 ia = newInfAlgFromString( "JTREE[inference=MAXPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
276 else
277 ia = newInfAlgFromString( "JTREE[inference=SUMPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
278 } else if( solver == 2 ) { // the decent one
279 double remtime = (UAI_time - (toc() - starttic));
280 if( remtime < 1.0 )
281 remtime = 1.0;
282 double maxtime = 0.95 * remtime / (ev2go.size() - _ev);
283 if( verbose ) {
284 cout << " Past time: " << (toc() - starttic) << endl;
285 cout << " Remaining time:" << remtime << endl;
286 cout << " Allotted time: " << maxtime << endl;
287 }
288 if( task == "MPE" )
289 maxtime /= fgs[ev].nrVars();
290 string maxtimestr;
291 if( maxtime != INFINITY )
292 maxtimestr = ",maxtime=" + toString(maxtime);
293 if( task == "MPE" )
294 ia = newInfAlgFromString( "DECMAP[ianame=BP,iaopts=[inference=MAXPROD,updates=SEQRND,logdomain=1,tol=1e-9,maxiter=10000" + maxtimestr + ",damping=" + toString(MPEdamping) + ",verbose=0],reinit=1,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
295 else {
296 if( couldBeGrid )
297 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
298 else {
299 if( task == "PR" )
300 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=MIN,init=UNIFORM,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
301 else
302 ia = newInfAlgFromString( "GIBBS[maxiter=1000000000,burnin=1000,restart=10000000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
303 }
304 }
305 }
306
307 // initialize
308 if( verbose )
309 cout << " Initializing inference algorithm..." << endl;
310 ia->init();
311 // run
312 if( verbose )
313 cout << " Running inference algorithm..." << endl;
314 ia->run();
315 if( verbose )
316 cout << " Inference algorithm finished..." << endl;
317 } catch( Exception &e ) {
318 failed = true;
319 if( verbose ) {
320 cout << " Inference algorithm failed...!" << endl;
321 cout << " Exception: " << e.what() << endl;
322 }
323 }
324
325 if( verbose )
326 cout << " Used time: " << toc() - tic << endl;
327 if( !failed && verbose ) {
328 try {
329 cout << " Number of iterations: " << ia->Iterations() << endl;
330 } catch( Exception &e ) {
331 cout << " Number of iterations: N/A" << endl;
332 }
333 try {
334 cout << " Final maxdiff: " << ia->maxDiff() << endl;
335 } catch( Exception &e ) {
336 cout << " Final maxdiff: N/A" << endl;
337 }
338 }
339
340 // update results for inference task
341 if( !failed ) {
342 if( task == "PR" ) {
343 PRbest cur;
344
345 // calculate PR value
346 cur.value = ia->logZ() / dai::log((Real)10.0);
347
348 // get maxdiff
349 try {
350 cur.maxdiff = ia->maxDiff();
351 } catch( Exception &e ) {
352 cur.maxdiff = 1e-9;
353 }
354
355 // only update if this run has converged
356 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestPR[ev].maxdiff)) && !dai::isnan(cur.value) ) {
357 // if this was exact inference, we are ready
358 if( solver == 1 ) {
359 ev2go.erase( ev2go.begin() + _ev );
360 _ev--;
361 cur.ready = true;
362 }
363
364 if( verbose )
365 cout << " Replacing best PR value so far (" << bestPR[ev].value << ") with new value " << cur.value << endl;
366 bestPR[ev] = cur;
367 improved = true;
368 } else {
369 if( verbose )
370 cout << " Discarding PR value " << cur.value << endl;
371 }
372 } else if( task == "MAR" ) {
373 MARbest cur;
374
375 // get variable beliefs
376 bool hasnans = false;
377 cur.beliefs.reserve( fgs[ev].nrVars() );
378 for( size_t i = 0; i < fgs[ev].nrVars(); i++ ) {
379 cur.beliefs.push_back( ia->beliefV(i) );
380 if( cur.beliefs.back().hasNaNs() )
381 hasnans = true;
382 }
383
384 // get maxdiff
385 try {
386 cur.maxdiff = ia->maxDiff();
387 } catch( Exception &e ) {
388 cur.maxdiff = 1e-9;
389 }
390
391 // only update if this run has converged
392 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestMAR[ev].maxdiff)) && !hasnans ) {
393 // if this was exact inference, we are ready
394 if( solver == 1 ) {
395 ev2go.erase( ev2go.begin() + _ev );
396 _ev--;
397 cur.ready = true;
398 }
399
400 if( verbose )
401 cout << " Replacing best beliefs so far with new beliefs" << endl;
402 bestMAR[ev] = cur;
403 improved = true;
404 } else {
405 if( verbose )
406 cout << " Discarding beliefs" << endl;
407 }
408 } else if( task == "MPE" ) {
409 MPEbest cur;
410
411 // calculate MPE state
412 cur.state = ia->findMaximum();
413
414 // calculate MPE value
415 cur.value = fgs[ev].logScore( cur.state );
416
417 // update best MPE state and value
418 if( cur.value > bestMPE[ev].value && !dai::isnan(cur.value) ) {
419 // if this was exact inference, we are ready
420 if( solver == 1 ) {
421 ev2go.erase( ev2go.begin() + _ev );
422 _ev--;
423 cur.ready = true;
424 }
425
426 if( verbose )
427 cout << " Replacing best MPE value so far (" << bestMPE[ev].value << ") with new value " << cur.value << endl;
428 bestMPE[ev] = cur;
429 improved = true;
430 } else {
431 if( verbose )
432 cout << " New MPE value " << cur.value << " not better than best one so far " << bestMPE[ev].value << endl;
433 }
434 }
435 }
436
437 // remove inference algorithm
438 if( verbose )
439 cout << " Cleaning up..." << endl;
440 if( !failed )
441 delete ia;
442
443 // write current best output to stream
444 if( improved ) {
445 if( verbose )
446 cout << " Writing output..." << endl;
447 if( first )
448 first = false;
449 else
450 os << "-BEGIN-" << endl;
451 os << evid.size() << endl;
452 for( size_t ev = 0; ev < evid.size(); ev++ ) {
453 if( task == "PR" ) {
454 // output probability of evidence
455 os << bestPR[ev].value << endl;
456 } else if( task == "MAR" ) {
457 // output variable marginals
458 os << bestMAR[ev].beliefs.size() << " ";
459 for( size_t i = 0; i < bestMAR[ev].beliefs.size(); i++ ) {
460 os << bestMAR[ev].beliefs[i].nrStates() << " ";
461 for( size_t s = 0; s < bestMAR[ev].beliefs[i].nrStates(); s++ )
462 os << bestMAR[ev].beliefs[i][s] << " ";
463 }
464 os << endl;
465 } else if( task == "MPE" ) {
466 // output MPE state
467 os << fgs[ev].nrVars() << " ";
468 for( size_t i = 0; i < fgs[ev].nrVars(); i++ )
469 os << bestMPE[ev].state[i] << " ";
470 os << endl;
471 }
472 }
473 os.flush();
474 }
475
476 if( verbose )
477 cout << " Done..." << endl;
478
479 if( !improved )
480 subsolver++;
481 if( improved || subsolver >= nrsubsolvers[solver] || couldBeGrid ) {
482 subsolver = 0;
483 _ev++;
484 }
485 }
486
487 if( task == "MPE" && solver == 2 && (toc() - starttic) < UAI_time && MPEdamping != 0.0 ) {
488 MPEdamping /= 2.0;
489 solver--; // repeat this one
490 }
491 if( ev2go.size() == 0 )
492 break;
493 }
494
495 // close output file
496 if( verbose )
497 cout << "Closing output file..." << endl;
498 os.close();
499
500 if( verbose )
501 cout << "Done!" << endl;
502 }
503
504 return 0;
505 }