Improved tokenizeString
[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 = tokenizeString( string(argv[1]), true, "/" );
280 string outfile = pathComponents.back() + "." + task;
281 if( verbose )
282 cout << "Output filename: " << outfile << endl;
283
284 // open output stream
285 ofstream os;
286 os.open( outfile.c_str() );
287 if( !os.is_open() )
288 DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + outfile);
289 if( verbose )
290 cout << "Opened output stream" << endl;
291
292 // read factor graph
293 vector<Var> vars;
294 vector<Factor> facs0;
295 vector<Permute> permutations;
296 if( verbose )
297 cout << "Reading factor graph..." << endl;
298 ReadUAIFGFile( argv[1], verbose, vars, facs0, permutations );
299 if( verbose )
300 cout << "Successfully read factor graph" << endl;
301
302 // check if it could be a grid
303 bool couldBeGrid = true;
304 FactorGraph fg0( facs0.begin(), facs0.end(), vars.begin(), vars.end(), facs0.size(), vars.size() );
305 for( size_t i = 0; i < fg0.nrVars(); i++ )
306 if( fg0.delta(i).size() > 4 ) {
307 couldBeGrid = false;
308 break;
309 }
310 if( couldBeGrid )
311 for( size_t I = 0; I < fg0.nrFactors(); I++ )
312 if( fg0.factor(I).vars().size() > 2 ) {
313 couldBeGrid = false;
314 break;
315 }
316 if( verbose ) {
317 if( couldBeGrid )
318 cout << "This could be a grid!" << endl;
319 else
320 cout << "This cannot be a grid!" << endl;
321 }
322
323 // read evidence
324 if( verbose )
325 cout << "Reading evidence..." << endl;
326 vector<map<size_t,size_t> > evid = ReadUAIEvidenceFile( argv[2], verbose );
327 if( verbose )
328 cout << "Successfully read " << evid.size() << " evidence cases" << endl;
329
330 // write output header
331 if( verbose )
332 cout << " Writing header to file..." << endl;
333 os << task << endl;
334
335 // construct clamped factor graphs
336 if( verbose )
337 cout << "Constructing clamped factor graphs..." << endl;
338 vector<FactorGraph> fgs;
339 fgs.reserve( evid.size() );
340 for( size_t ev = 0; ev < evid.size(); ev++ ) {
341 if( verbose )
342 cout << " Evidence case " << ev << "..." << endl;
343 // copy vector of factors
344 vector<Factor> facs( facs0 );
345
346 // change factor graph to reflect observed evidence
347 if( verbose )
348 cout << " Applying evidence..." << endl;
349 if( surgery ) {
350 // replace factors with clamped variables with slices
351 for( size_t I = 0; I < facs.size(); I++ ) {
352 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ ) {
353 if( facs[I].vars() >> vars[e->first] ) {
354 if( verbose >= 2 )
355 cout << " Clamping " << e->first << " to value " << e->second << " in factor " << I << " = " << facs[I].vars() << endl;
356 facs[I] = facs[I].slice( vars[e->first], e->second );
357 if( verbose >= 2 )
358 cout << " ...remaining vars: " << facs[I].vars() << endl;
359 }
360 }
361 }
362 // remove empty factors
363 Real logZcorr = 0.0;
364 for( vector<Factor>::iterator I = facs.begin(); I != facs.end(); )
365 if( I->vars().size() == 0 ) {
366 logZcorr += std::log( (Real)(*I)[0] );
367 I = facs.erase( I );
368 } else
369 I++;
370 // multiply with logZcorr constant
371 if( facs.size() == 0 )
372 facs.push_back( Factor( VarSet(), std::exp(logZcorr) ) );
373 else
374 facs.front() *= std::exp(logZcorr);
375 }
376 // add delta factors corresponding to observed variable values
377 for( map<size_t,size_t>::const_iterator e = evid[ev].begin(); e != evid[ev].end(); e++ )
378 facs.push_back( createFactorDelta( vars[e->first], e->second ) );
379
380 // construct clamped factor graph
381 if( verbose )
382 cout << " Constructing factor graph..." << endl;
383 fgs.push_back( FactorGraph( facs.begin(), facs.end(), vars.begin(), vars.end(), facs.size(), vars.size() ) );
384 }
385
386 // variables for storing best results so far
387 vector<PRbest> bestPR( evid.size() );
388 vector<MARbest> bestMAR( evid.size() );
389 vector<MPEbest> bestMPE( evid.size() );
390 for( size_t ev = 0; ev < evid.size(); ev++ )
391 bestMPE[ev].state = stateVec( fgs[ev].nrVars(), 0 );
392 vector<size_t> ev2go;
393 ev2go.reserve( evid.size() );
394 for( size_t ev = 0; ev < evid.size(); ev++ )
395 ev2go.push_back( ev );
396
397 // solve inference problems
398 if( verbose )
399 cout << "Solving inference problems..." << endl;
400 bool first = true;
401 size_t nrsolvers = 3;
402 vector<size_t> nrsubsolvers( nrsolvers );
403 nrsubsolvers[0] = 2;
404 nrsubsolvers[1] = 1;
405 nrsubsolvers[2] = 1;
406 double MPEdamping = 0.49;
407 // for each (sub)solver
408 for( size_t solver = 0; solver < nrsolvers; solver++ ) {
409 if( verbose )
410 cout << " Solver " << solver << endl;
411
412 // for each evidence case
413 size_t subsolver = 0;
414 for( long _ev = 0; _ev < (long)ev2go.size(); ) {
415 bool improved = false;
416 size_t ev = ev2go[_ev];
417 if( verbose )
418 cout << " Evidence case " << ev << ", subsolver = " << subsolver << "..." << endl;
419
420 // construct inference algorithm on clamped factor graph
421 if( verbose )
422 cout << " Constructing inference algorithm..." << endl;
423 InfAlg *ia;
424 double tic = toc();
425 bool failed = false;
426 try {
427 // construct
428 if( solver == 0 ) { // the quick one
429 double remtime = (UAI_time - (toc() - starttic)) * 0.9;
430 if( remtime < 1.0 )
431 remtime = 1.0 ;
432 double maxtime = remtime / (ev2go.size() - _ev);
433 if( verbose ) {
434 cout << " Past time: " << (toc() - starttic) << endl;
435 cout << " Remaining time:" << remtime << endl;
436 cout << " Allotted time: " << maxtime << endl;
437 }
438 string maxtimestr;
439 if( maxtime != INFINITY )
440 maxtimestr = ",maxtime=" + toString(maxtime);
441 // quick and dirty...
442 if( task == "MPE" )
443 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] );
444 else {
445 if( couldBeGrid )
446 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9,maxiter=10000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
447 else
448 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] );
449 }
450 } else if( solver == 1 ) { // the exact one
451 string maxmemstr;
452 if( UAI_memory != 0 )
453 maxmemstr = ",maxmem=" + toString(UAI_memory);
454 if( task == "MPE" )
455 ia = newInfAlgFromString( "JTREE[inference=MAXPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
456 else
457 ia = newInfAlgFromString( "JTREE[inference=SUMPROD,updates=HUGIN" + maxmemstr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
458 } else if( solver == 2 ) { // the decent one
459 double remtime = (UAI_time - (toc() - starttic));
460 if( remtime < 1.0 )
461 remtime = 1.0;
462 double maxtime = 0.95 * remtime / (ev2go.size() - _ev);
463 if( verbose ) {
464 cout << " Past time: " << (toc() - starttic) << endl;
465 cout << " Remaining time:" << remtime << endl;
466 cout << " Allotted time: " << maxtime << endl;
467 }
468 if( task == "MPE" )
469 maxtime /= fgs[ev].nrVars();
470 string maxtimestr;
471 if( maxtime != INFINITY )
472 maxtimestr = ",maxtime=" + toString(maxtime);
473 if( task == "MPE" )
474 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] );
475 else {
476 if( couldBeGrid )
477 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=LOOP,init=UNIFORM,loopdepth=4,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
478 else {
479 if( task == "PR" )
480 ia = newInfAlgFromString( "HAK[doubleloop=1,clusters=MIN,init=UNIFORM,tol=1e-9" + maxtimestr + ",maxiter=100000,verbose=" + toString(ia_verbose) + "]", fgs[ev] );
481 else
482 ia = newInfAlgFromString( "GIBBS[maxiter=1000000000,burnin=1000,restart=10000000" + maxtimestr + ",verbose=" + toString(ia_verbose) + "]", fgs[ev] );
483 }
484 }
485 }
486
487 // initialize
488 if( verbose )
489 cout << " Initializing inference algorithm..." << endl;
490 ia->init();
491 // run
492 if( verbose )
493 cout << " Running inference algorithm..." << endl;
494 ia->run();
495 if( verbose )
496 cout << " Inference algorithm finished..." << endl;
497 } catch( Exception &e ) {
498 failed = true;
499 if( verbose ) {
500 cout << " Inference algorithm failed...!" << endl;
501 cout << " Exception: " << e.what() << endl;
502 }
503 }
504
505 if( verbose )
506 cout << " Used time: " << toc() - tic << endl;
507 if( !failed && verbose ) {
508 try {
509 cout << " Number of iterations: " << ia->Iterations() << endl;
510 } catch( Exception &e ) {
511 cout << " Number of iterations: N/A" << endl;
512 }
513 try {
514 cout << " Final maxdiff: " << ia->maxDiff() << endl;
515 } catch( Exception &e ) {
516 cout << " Final maxdiff: N/A" << endl;
517 }
518 }
519
520 // update results for inference task
521 if( !failed ) {
522 if( task == "PR" ) {
523 PRbest cur;
524
525 // calculate PR value
526 cur.value = ia->logZ() / dai::log((Real)10.0);
527
528 // get maxdiff
529 try {
530 cur.maxdiff = ia->maxDiff();
531 } catch( Exception &e ) {
532 cur.maxdiff = 1e-9;
533 }
534
535 // only update if this run has converged
536 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestPR[ev].maxdiff)) && !dai::isnan(cur.value) ) {
537 // if this was exact inference, we are ready
538 if( solver == 1 ) {
539 ev2go.erase( ev2go.begin() + _ev );
540 _ev--;
541 cur.ready = true;
542 }
543
544 if( verbose )
545 cout << " Replacing best PR value so far (" << bestPR[ev].value << ") with new value " << cur.value << endl;
546 bestPR[ev] = cur;
547 improved = true;
548 } else {
549 if( verbose )
550 cout << " Discarding PR value " << cur.value << endl;
551 }
552 } else if( task == "MAR" ) {
553 MARbest cur;
554
555 // get variable beliefs
556 bool hasnans = false;
557 cur.beliefs.reserve( fgs[ev].nrVars() );
558 for( size_t i = 0; i < fgs[ev].nrVars(); i++ ) {
559 cur.beliefs.push_back( ia->beliefV(i) );
560 if( cur.beliefs.back().hasNaNs() )
561 hasnans = true;
562 }
563
564 // get maxdiff
565 try {
566 cur.maxdiff = ia->maxDiff();
567 } catch( Exception &e ) {
568 cur.maxdiff = 1e-9;
569 }
570
571 // only update if this run has converged
572 if( ((cur.maxdiff <= 1e-9) || (cur.maxdiff <= bestMAR[ev].maxdiff)) && !hasnans ) {
573 // if this was exact inference, we are ready
574 if( solver == 1 ) {
575 ev2go.erase( ev2go.begin() + _ev );
576 _ev--;
577 cur.ready = true;
578 }
579
580 if( verbose )
581 cout << " Replacing best beliefs so far with new beliefs" << endl;
582 bestMAR[ev] = cur;
583 improved = true;
584 } else {
585 if( verbose )
586 cout << " Discarding beliefs" << endl;
587 }
588 } else if( task == "MPE" ) {
589 MPEbest cur;
590
591 // calculate MPE state
592 cur.state = ia->findMaximum();
593
594 // calculate MPE value
595 cur.value = fgs[ev].logScore( cur.state );
596
597 // update best MPE state and value
598 if( cur.value > bestMPE[ev].value && !dai::isnan(cur.value) ) {
599 // if this was exact inference, we are ready
600 if( solver == 1 ) {
601 ev2go.erase( ev2go.begin() + _ev );
602 _ev--;
603 cur.ready = true;
604 }
605
606 if( verbose )
607 cout << " Replacing best MPE value so far (" << bestMPE[ev].value << ") with new value " << cur.value << endl;
608 bestMPE[ev] = cur;
609 improved = true;
610 } else {
611 if( verbose )
612 cout << " New MPE value " << cur.value << " not better than best one so far " << bestMPE[ev].value << endl;
613 }
614 }
615 }
616
617 // remove inference algorithm
618 if( verbose )
619 cout << " Cleaning up..." << endl;
620 if( !failed )
621 delete ia;
622
623 // write current best output to stream
624 if( improved ) {
625 if( verbose )
626 cout << " Writing output..." << endl;
627 if( first )
628 first = false;
629 else
630 os << "-BEGIN-" << endl;
631 os << evid.size() << endl;
632 for( size_t ev = 0; ev < evid.size(); ev++ ) {
633 if( task == "PR" ) {
634 // output probability of evidence
635 os << bestPR[ev].value << endl;
636 } else if( task == "MAR" ) {
637 // output variable marginals
638 os << bestMAR[ev].beliefs.size() << " ";
639 for( size_t i = 0; i < bestMAR[ev].beliefs.size(); i++ ) {
640 os << bestMAR[ev].beliefs[i].nrStates() << " ";
641 for( size_t s = 0; s < bestMAR[ev].beliefs[i].nrStates(); s++ )
642 os << bestMAR[ev].beliefs[i][s] << " ";
643 }
644 os << endl;
645 } else if( task == "MPE" ) {
646 // output MPE state
647 os << fgs[ev].nrVars() << " ";
648 for( size_t i = 0; i < fgs[ev].nrVars(); i++ )
649 os << bestMPE[ev].state[i] << " ";
650 os << endl;
651 }
652 }
653 os.flush();
654 }
655
656 if( verbose )
657 cout << " Done..." << endl;
658
659 if( !improved )
660 subsolver++;
661 if( improved || subsolver >= nrsubsolvers[solver] || couldBeGrid ) {
662 subsolver = 0;
663 _ev++;
664 }
665 }
666
667 if( task == "MPE" && solver == 2 && (toc() - starttic) < UAI_time && MPEdamping != 0.0 ) {
668 MPEdamping /= 2.0;
669 solver--; // repeat this one
670 }
671 if( ev2go.size() == 0 )
672 break;
673 }
674
675 // close output file
676 if( verbose )
677 cout << "Closing output file..." << endl;
678 os.close();
679
680 if( verbose )
681 cout << "Done!" << endl;
682 }
683
684 return 0;
685 }