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