Added source code for one of the winning solvers of the UAI 2010 Approximate Inferenc...
[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/util.h>
17 #include <dai/index.h>
18 #include <dai/jtree.h>
19
20
21 using namespace std;
22 using namespace dai;
23
24
25 // Type for storing a joint state of all variables
26 typedef std::vector<size_t> stateVec;
27
28
29 struct PRbest {
30 Real value;
31 Real maxdiff;
32 bool ready;
33 PRbest() : value(0.0), maxdiff(INFINITY), ready(false) {}
34 };
35
36 struct MARbest {
37 vector<Factor> beliefs;
38 Real maxdiff;
39 bool ready;
40 MARbest() : beliefs(), maxdiff(INFINITY), ready(false) {}
41 };
42
43 struct MPEbest {
44 Real value;
45 stateVec state;
46 bool ready;
47 MPEbest() : value(-INFINITY), state(), ready(false) {}
48 };
49
50
51 /// Reads "evidence" (a mapping from observed variable labels to the observed values) from a UAI evidence file
52 vector<map<size_t, size_t> > ReadUAIEvidenceFile( char* filename, size_t verbose ) {
53 vector<map<size_t, size_t> > evid;
54
55 // open file
56 ifstream is;
57 is.open( filename );
58 if( is.is_open() ) {
59 // read number of evidence cases
60 size_t nr_evid;
61 is >> nr_evid;
62 if( is.fail() )
63 DAI_THROWE(INVALID_EVIDENCE_FILE,"Cannot read number of evidence cases");
64 if( verbose >= 2 )
65 cout << "Reading " << nr_evid << " evidence cases..." << endl;
66
67 evid.resize( nr_evid );
68 for( size_t i = 0; i < nr_evid; i++ ) {
69 // read number of variables
70 size_t nr_obs;
71 is >> nr_obs;
72 if( is.fail() )
73 DAI_THROWE(INVALID_EVIDENCE_FILE,"Evidence case " + toString(i) + ": Cannot read number of observations");
74 if( verbose >= 2 )
75 cout << "Evidence case " << i << ": reading " << nr_obs << " observations..." << endl;
76
77 // for each observation, read the variable label and the observed value
78 for( size_t j = 0; j < nr_obs; j++ ) {
79 size_t label, val;
80 is >> label;
81 if( is.fail() )
82 DAI_THROWE(INVALID_EVIDENCE_FILE,"Evidence case " + toString(i) + ": Cannot read label for " + toString(j) + "'th observed variable");
83 is >> val;
84 if( is.fail() )
85 DAI_THROWE(INVALID_EVIDENCE_FILE,"Evidence case " + toString(i) + ": Cannot read value of " + toString(j) + "'th observed variable");
86 if( verbose >= 3 )
87 cout << " variable: " << label << ", value: " << val << endl;
88 evid[i][label] = val;
89 }
90 }
91
92 // close file
93 is.close();
94 } else
95 DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
96
97 if( evid.size() == 0 )
98 evid.resize( 1 );
99
100 return evid;
101 }
102
103
104 /// Reads factor graph (as a pair of a variable vector and factor vector) from a UAI factor graph file
105 void ReadUAIFGFile( const char *filename, size_t verbose, vector<Var>& vars, vector<Factor>& factors, vector<Permute>& permutations ) {
106 vars.clear();
107 factors.clear();
108 permutations.clear();
109
110 // open file
111 ifstream is;
112 is.open( filename );
113 if( is.is_open() ) {
114 size_t nrFacs, nrVars;
115 string line;
116
117 // read header line
118 getline(is,line);
119 if( is.fail() || (line != "BAYES" && line != "MARKOV" && line != "BAYES\r" && line != "MARKOV\r") )
120 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"UAI factor graph file should start with \"BAYES\" or \"MARKOV\"");
121 if( verbose >= 2 )
122 cout << "Reading " << line << " network..." << endl;
123
124 // read number of variables
125 is >> nrVars;
126 if( is.fail() )
127 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of variables");
128 if( verbose >= 2 )
129 cout << "Reading " << nrVars << " variables..." << endl;
130
131 // for each variable, read its number of states
132 vars.reserve( nrVars );
133 for( size_t i = 0; i < nrVars; i++ ) {
134 size_t dim;
135 is >> dim;
136 if( is.fail() )
137 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of states of " + toString(i) + "'th variable");
138 vars.push_back( Var( i, dim ) );
139 }
140
141 // read number of factors
142 is >> nrFacs;
143 if( is.fail() )
144 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
145 if( verbose >= 2 )
146 cout << "Reading " << nrFacs << " factors..." << endl;
147
148 // for each factor, read the variables on which it depends
149 vector<vector<Var> > factorVars;
150 factors.reserve( nrFacs );
151 factorVars.reserve( nrFacs );
152 for( size_t I = 0; I < nrFacs; I++ ) {
153 if( verbose >= 3 )
154 cout << "Reading factor " << I << "..." << endl;
155
156 // read number of variables for factor I
157 size_t I_nrVars;
158 is >> I_nrVars;
159 if( is.fail() )
160 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of variables for " + toString(I) + "'th factor");
161 if( verbose >= 3 )
162 cout << " which depends on " << I_nrVars << " variables" << endl;
163
164 // read the variable labels
165 vector<long> I_labels;
166 vector<size_t> I_dims;
167 I_labels.reserve( I_nrVars );
168 I_dims.reserve( I_nrVars );
169 factorVars[I].reserve( I_nrVars );
170 for( size_t _i = 0; _i < I_nrVars; _i++ ) {
171 long label;
172 is >> label;
173 if( is.fail() )
174 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read variable labels for " + toString(I) + "'th factor");
175 I_labels.push_back( label );
176 I_dims.push_back( vars[label].states() );
177 factorVars[I].push_back( vars[label] );
178 }
179 if( verbose >= 3 )
180 cout << " labels: " << I_labels << ", dimensions " << I_dims << endl;
181
182 // add the factor and the labels
183 factors.push_back( Factor( VarSet( factorVars[I].begin(), factorVars[I].end(), factorVars[I].size() ), (Real)0 ) );
184 }
185
186 // for each factor, read its values
187 permutations.reserve( nrFacs );
188 for( size_t I = 0; I < nrFacs; I++ ) {
189 if( verbose >= 3 )
190 cout << "Reading factor " << I << "..." << endl;
191
192 // calculate permutation object, reversing the indexing in factorVars[I] first
193 Permute permindex( factorVars[I], true );
194 permutations.push_back( permindex );
195
196 // read factor values
197 size_t nrNonZeros;
198 is >> nrNonZeros;
199 if( is.fail() )
200 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of nonzero factor values for " + toString(I) + "'th factor");
201 if( verbose >= 3 )
202 cout << " number of nonzero values: " << nrNonZeros << endl;
203 DAI_ASSERT( nrNonZeros == factors[I].nrStates() );
204 for( size_t li = 0; li < nrNonZeros; li++ ) {
205 Real val;
206 is >> val;
207 if( is.fail() )
208 DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read factor values of " + toString(I) + "'th factor");
209 // assign value after calculating its linear index corresponding to the permutation
210 if( verbose >= 4 )
211 cout << " " << li << "'th value " << val << " corresponds with index " << permindex.convertLinearIndex(li) << endl;
212 factors[I].set( permindex.convertLinearIndex( li ), val );
213 }
214 }
215 if( verbose >= 3 )
216 cout << "variables:" << vars << endl;
217 if( verbose >= 3 )
218 cout << "factors:" << factors << endl;
219
220 // close file
221 is.close();
222 } else
223 DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
224 }
225
226
227 int main( int argc, char *argv[] ) {
228 if ( argc != 5 ) {
229 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
230 cout << "It is one of the winning solvers that participated in the" << endl;
231 cout << "UAI 2010 Approximate Inference Challenge" << endl;
232 cout << "(see http://www.cs.huji.ac.il/project/UAI10/)" << endl << endl;
233 cout << "Usage: " << argv[0] << " <filename.uai> <filename.uai.evid> <seed> <task>" << endl << endl;
234 return 1;
235 } else {
236 double starttic = toc();
237
238 size_t verbose = 1;
239 size_t ia_verbose = 0;
240 bool surgery = true;
241 if( verbose )
242 cout << "Solver: " << argv[0] << endl;
243
244 // set random seed
245 size_t seed = fromString<size_t>( argv[3] );
246 rnd_seed( seed );
247 if( verbose )
248 cout << "Seed: " << seed << endl;
249
250 // check whether the task is valid
251 string task( argv[4] );
252 if( task != string("PR") && task != string("MPE") && task != string("MAR") )
253 DAI_THROWE(RUNTIME_ERROR,"Unknown task");
254 if( verbose )
255 cout << "Task: " << task << endl;
256
257 // output other command line options
258 if( verbose ) {
259 cout << "Factorgraph filename: " << argv[1] << endl;
260 cout << "Evidence filename: " << argv[2] << endl;
261 }
262
263 // get time and memory limits
264 char *buf = getenv( "UAI_TIME" );
265 double UAI_time = INFINITY;
266 if( buf != NULL )
267 UAI_time = fromString<double>( buf );
268 buf = getenv( "UAI_MEMORY" );
269 size_t UAI_memory = 0;
270 if( buf != NULL ) {
271 UAI_memory = fromString<double>( buf ) * 1024 * 1024 * 1024;
272 }
273 if( verbose ) {
274 cout << "Time limit: " << UAI_time << endl;
275 cout << "Memory limit: " << UAI_memory << endl;
276 }
277
278 // build output file name
279 vector<string> pathComponents;
280 tokenizeString( string(argv[1]), pathComponents, "/" );
281 string outfile = pathComponents.back() + "." + task;
282 if( verbose )
283 cout << "Output filename: " << outfile << endl;
284
285 // open output stream
286 ofstream os;
287 os.open( outfile.c_str() );
288 if( !os.is_open() )
289 DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + outfile);
290 if( verbose )
291 cout << "Opened output stream" << endl;
292
293 // read factor graph
294 vector<Var> vars;
295 vector<Factor> facs0;
296 vector<Permute> permutations;
297 if( verbose )
298 cout << "Reading factor graph..." << endl;
299 ReadUAIFGFile( argv[1], verbose, vars, facs0, permutations );
300 if( verbose )
301 cout << "Successfully read factor graph" << endl;
302
303 // check if it could be a grid
304 bool couldBeGrid = true;
305 FactorGraph fg0( facs0.begin(), facs0.end(), vars.begin(), vars.end(), facs0.size(), vars.size() );
306 for( size_t i = 0; i < fg0.nrVars(); i++ )
307 if( fg0.delta(i).size() > 4 ) {
308 couldBeGrid = false;
309 break;
310 }
311 if( couldBeGrid )
312 for( size_t I = 0; I < fg0.nrFactors(); I++ )
313 if( fg0.factor(I).vars().size() > 2 ) {
314 couldBeGrid = false;
315 break;
316 }
317 if( verbose ) {
318 if( couldBeGrid )
319 cout << "This could be a grid!" << endl;
320 else
321 cout << "This cannot be a grid!" << endl;
322 }
323
324 // read evidence
325 if( verbose )
326 cout << "Reading evidence..." << endl;
327 vector<map<size_t,size_t> > evid = ReadUAIEvidenceFile( argv[2], verbose );
328 if( verbose )
329 cout << "Successfully read " << evid.size() << " evidence cases" << endl;
330
331 // write output header
332 if( verbose )
333 cout << " Writing header to file..." << endl;
334 os << task << endl;
335
336 // construct clamped factor graphs
337 if( verbose )
338 cout << "Constructing clamped factor graphs..." << endl;
339 vector<FactorGraph> fgs;
340 fgs.reserve( evid.size() );
341 for( size_t ev = 0; ev < evid.size(); ev++ ) {
342 if( verbose )
343 cout << " Evidence case " << ev << "..." << endl;
344 // copy vector of factors
345 vector<Factor> facs( facs0 );
346
347 // change factor graph to reflect observed evidence
348 if( verbose )
349 cout << " Applying evidence..." << endl;
350 if( surgery ) {
351 // replace factors with clamped variables with slices
352 for( size_t I = 0; I < facs.size(); I++ ) {
353 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ ) {
354 if( facs[I].vars() >> vars[e->first] ) {
355 if( verbose >= 2 )
356 cout << " Clamping " << e->first << " to value " << e->second << " in factor " << I << " = " << facs[I].vars() << endl;
357 facs[I] = facs[I].slice( vars[e->first], e->second );
358 if( verbose >= 2 )
359 cout << " ...remaining vars: " << facs[I].vars() << endl;
360 }
361 }
362 }
363 // remove empty factors
364 Real logZcorr = 0.0;
365 for( vector<Factor>::iterator I = facs.begin(); I != facs.end(); )
366 if( I->vars().size() == 0 ) {
367 logZcorr += std::log( (Real)(*I)[0] );
368 I = facs.erase( I );
369 } else
370 I++;
371 // multiply with logZcorr constant
372 if( facs.size() == 0 )
373 facs.push_back( Factor( VarSet(), std::exp(logZcorr) ) );
374 else
375 facs.front() *= std::exp(logZcorr);
376 }
377 // add delta factors corresponding to observed variable values
378 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ )
379 facs.push_back( createFactorDelta( vars[e->first], e->second ) );
380
381 // construct clamped factor graph
382 if( verbose )
383 cout << " Constructing factor graph..." << endl;
384 fgs.push_back( FactorGraph( facs.begin(), facs.end(), vars.begin(), vars.end(), facs.size(), vars.size() ) );
385 }
386
387 // variables for storing best results so far
388 vector<PRbest> bestPR( evid.size() );
389 vector<MARbest> bestMAR( evid.size() );
390 vector<MPEbest> bestMPE( evid.size() );
391 for( size_t ev = 0; ev < evid.size(); ev++ )
392 bestMPE[ev].state = stateVec( fgs[ev].nrVars(), 0 );
393 vector<size_t> ev2go;
394 ev2go.reserve( evid.size() );
395 for( size_t ev = 0; ev < evid.size(); ev++ )
396 ev2go.push_back( ev );
397
398 // solve inference problems
399 if( verbose )
400 cout << "Solving inference problems..." << endl;
401 bool first = true;
402 size_t nrsolvers = 3;
403 vector<size_t> nrsubsolvers( nrsolvers );
404 nrsubsolvers[0] = 2;
405 nrsubsolvers[1] = 1;
406 nrsubsolvers[2] = 1;
407 double MPEdamping = 0.49;
408 // for each (sub)solver
409 for( size_t solver = 0; solver < nrsolvers; solver++ ) {
410 if( verbose )
411 cout << " Solver " << solver << endl;
412
413 // for each evidence case
414 size_t subsolver = 0;
415 for( long _ev = 0; _ev < (long)ev2go.size(); ) {
416 bool improved = false;
417 size_t ev = ev2go[_ev];
418 if( verbose )
419 cout << " Evidence case " << ev << ", subsolver = " << subsolver << "..." << endl;
420
421 // construct inference algorithm on clamped factor graph
422 if( verbose )
423 cout << " Constructing inference algorithm..." << endl;
424 InfAlg *ia;
425 double tic = toc();
426 bool failed = false;
427 try {
428 // construct
429 if( solver == 0 ) { // the quick one
430 double remtime = (UAI_time - (toc() - starttic)) * 0.9;
431 if( remtime < 1.0 )
432 remtime = 1.0 ;
433 double maxtime = remtime / (ev2go.size() - _ev);
434 if( verbose ) {
435 cout << " Past time: " << (toc() - starttic) << endl;
436 cout << " Remaining time:" << remtime << endl;
437 cout << " Allotted time: " << maxtime << endl;
438 }
439 string maxtimestr;
440 if( maxtime != INFINITY )
441 maxtimestr = ",maxtime=" + toString(maxtime);
442 // quick and dirty...
443 if( task == "MPE" )
444 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] );
445 else {
446 if( couldBeGrid )
447 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9,maxiter=10000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
448 else
449 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] );
450 }
451 } else if( solver == 1 ) { // the exact one
452 string maxmemstr;
453 if( UAI_memory != 0 )
454 maxmemstr = ",maxmem=" + toString(UAI_memory);
455 if( task == "MPE" )
456 ia = newInfAlgFromString( "JTREE[inference=MAXPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
457 else
458 ia = newInfAlgFromString( "JTREE[inference=SUMPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
459 } else if( solver == 2 ) { // the decent one
460 double remtime = (UAI_time - (toc() - starttic));
461 if( remtime < 1.0 )
462 remtime = 1.0;
463 double maxtime = 0.95 * remtime / (ev2go.size() - _ev);
464 if( verbose ) {
465 cout << " Past time: " << (toc() - starttic) << endl;
466 cout << " Remaining time:" << remtime << endl;
467 cout << " Allotted time: " << maxtime << endl;
468 }
469 if( task == "MPE" )
470 maxtime /= fgs[ev].nrVars();
471 string maxtimestr;
472 if( maxtime != INFINITY )
473 maxtimestr = ",maxtime=" + toString(maxtime);
474 if( task == "MPE" )
475 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] );
476 else {
477 if( couldBeGrid )
478 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
479 else {
480 if( task == "PR" )
481 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=MIN,init=UNIFORM,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
482 else
483 ia = newInfAlgFromString( "GIBBS[maxiter=1000000000,burnin=1000,restart=10000000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
484 }
485 }
486 }
487
488 // initialize
489 if( verbose )
490 cout << " Initializing inference algorithm..." << endl;
491 ia->init();
492 // run
493 if( verbose )
494 cout << " Running inference algorithm..." << endl;
495 ia->run();
496 if( verbose )
497 cout << " Inference algorithm finished..." << endl;
498 } catch( Exception &e ) {
499 failed = true;
500 if( verbose ) {
501 cout << " Inference algorithm failed...!" << endl;
502 cout << " Exception: " << e.what() << endl;
503 }
504 }
505
506 if( verbose )
507 cout << " Used time: " << toc() - tic << endl;
508 if( !failed && verbose ) {
509 try {
510 cout << " Number of iterations: " << ia->Iterations() << endl;
511 } catch( Exception &e ) {
512 cout << " Number of iterations: N/A" << endl;
513 }
514 try {
515 cout << " Final maxdiff: " << ia->maxDiff() << endl;
516 } catch( Exception &e ) {
517 cout << " Final maxdiff: N/A" << endl;
518 }
519 }
520
521 // update results for inference task
522 if( !failed ) {
523 if( task == "PR" ) {
524 PRbest cur;
525
526 // calculate PR value
527 cur.value = ia->logZ() / dai::log((Real)10.0);
528
529 // get maxdiff
530 try {
531 cur.maxdiff = ia->maxDiff();
532 } catch( Exception &e ) {
533 cur.maxdiff = 1e-9;
534 }
535
536 // only update if this run has converged
537 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestPR[ev].maxdiff)) && !dai::isnan(cur.value) ) {
538 // if this was exact inference, we are ready
539 if( solver == 1 ) {
540 ev2go.erase( ev2go.begin() + _ev );
541 _ev--;
542 cur.ready = true;
543 }
544
545 if( verbose )
546 cout << " Replacing best PR value so far (" << bestPR[ev].value << ") with new value " << cur.value << endl;
547 bestPR[ev] = cur;
548 improved = true;
549 } else {
550 if( verbose )
551 cout << " Discarding PR value " << cur.value << endl;
552 }
553 } else if( task == "MAR" ) {
554 MARbest cur;
555
556 // get variable beliefs
557 bool hasnans = false;
558 cur.beliefs.reserve( fgs[ev].nrVars() );
559 for( size_t i = 0; i < fgs[ev].nrVars(); i++ ) {
560 cur.beliefs.push_back( ia->beliefV(i) );
561 if( cur.beliefs.back().hasNaNs() )
562 hasnans = true;
563 }
564
565 // get maxdiff
566 try {
567 cur.maxdiff = ia->maxDiff();
568 } catch( Exception &e ) {
569 cur.maxdiff = 1e-9;
570 }
571
572 // only update if this run has converged
573 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestMAR[ev].maxdiff)) && !hasnans ) {
574 // if this was exact inference, we are ready
575 if( solver == 1 ) {
576 ev2go.erase( ev2go.begin() + _ev );
577 _ev--;
578 cur.ready = true;
579 }
580
581 if( verbose )
582 cout << " Replacing best beliefs so far with new beliefs" << endl;
583 bestMAR[ev] = cur;
584 improved = true;
585 } else {
586 if( verbose )
587 cout << " Discarding beliefs" << endl;
588 }
589 } else if( task == "MPE" ) {
590 MPEbest cur;
591
592 // calculate MPE state
593 cur.state = ia->findMaximum();
594
595 // calculate MPE value
596 cur.value = fgs[ev].logScore( cur.state );
597
598 // update best MPE state and value
599 if( cur.value > bestMPE[ev].value && !dai::isnan(cur.value) ) {
600 // if this was exact inference, we are ready
601 if( solver == 1 ) {
602 ev2go.erase( ev2go.begin() + _ev );
603 _ev--;
604 cur.ready = true;
605 }
606
607 if( verbose )
608 cout << " Replacing best MPE value so far (" << bestMPE[ev].value << ") with new value " << cur.value << endl;
609 bestMPE[ev] = cur;
610 improved = true;
611 } else {
612 if( verbose )
613 cout << " New MPE value " << cur.value << " not better than best one so far " << bestMPE[ev].value << endl;
614 }
615 }
616 }
617
618 // remove inference algorithm
619 if( verbose )
620 cout << " Cleaning up..." << endl;
621 if( !failed )
622 delete ia;
623
624 // write current best output to stream
625 if( improved ) {
626 if( verbose )
627 cout << " Writing output..." << endl;
628 if( first )
629 first = false;
630 else
631 os << "-BEGIN-" << endl;
632 os << evid.size() << endl;
633 for( size_t ev = 0; ev < evid.size(); ev++ ) {
634 if( task == "PR" ) {
635 // output probability of evidence
636 os << bestPR[ev].value << endl;
637 } else if( task == "MAR" ) {
638 // output variable marginals
639 os << bestMAR[ev].beliefs.size() << " ";
640 for( size_t i = 0; i < bestMAR[ev].beliefs.size(); i++ ) {
641 os << bestMAR[ev].beliefs[i].nrStates() << " ";
642 for( size_t s = 0; s < bestMAR[ev].beliefs[i].nrStates(); s++ )
643 os << bestMAR[ev].beliefs[i][s] << " ";
644 }
645 os << endl;
646 } else if( task == "MPE" ) {
647 // output MPE state
648 os << fgs[ev].nrVars() << " ";
649 for( size_t i = 0; i < fgs[ev].nrVars(); i++ )
650 os << bestMPE[ev].state[i] << " ";
651 os << endl;
652 }
653 }
654 os.flush();
655 }
656
657 if( verbose )
658 cout << " Done..." << endl;
659
660 if( !improved )
661 subsolver++;
662 if( improved || subsolver >= nrsubsolvers[solver] || couldBeGrid ) {
663 subsolver = 0;
664 _ev++;
665 }
666 }
667
668 if( task == "MPE" && solver == 2 && (toc() - starttic) < UAI_time && MPEdamping != 0.0 ) {
669 MPEdamping /= 2.0;
670 solver--; // repeat this one
671 }
672 if( ev2go.size() == 0 )
673 break;
674 }
675
676 // close output file
677 if( verbose )
678 cout << "Closing output file..." << endl;
679 os.close();
680
681 if( verbose )
682 cout << "Done!" << endl;
683 }
684
685 return 0;
686 }