Fixed two problems related to g++ 4.0.0 on Darwin 9.8.0
[libdai.git] / tests / testdai.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 <fstream>
14 #include <map>
15 #include <numeric>
16 #include <cmath>
17 #include <cstdlib>
18 #include <cstring>
19 #include <boost/program_options.hpp>
20 #include <dai/util.h>
21 #include <dai/alldai.h>
22
23
24 using namespace std;
25 using namespace dai;
26 namespace po = boost::program_options;
27
28
29 std::vector<Real> calcDists( const vector<Factor> &x, const vector<Factor> &y ) {
30 vector<Real> errs;
31 errs.reserve( x.size() );
32 DAI_ASSERT( x.size() == y.size() );
33 for( size_t i = 0; i < x.size(); i++ )
34 errs.push_back( dist( x[i], y[i], DISTTV ) );
35 return errs;
36 }
37
38
39 /// Wrapper class for DAI approximate inference algorithms
40 class TestDAI {
41 protected:
42 /// Stores a pointer to an InfAlg object, managed by this class
43 InfAlg *obj;
44 /// Stores the name of the InfAlg algorithm
45 string name;
46 /// Stores the total variation distances of the variable marginals
47 vector<Real> varErr;
48 /// Stores the total variation distances of the factor marginals
49 vector<Real> facErr;
50
51 public:
52 /// Stores the variable marginals
53 vector<Factor> varMarginals;
54 /// Stores the factor marginals
55 vector<Factor> facMarginals;
56 /// Stores all marginals
57 vector<Factor> allMarginals;
58 /// Stores the logarithm of the partition sum
59 Real logZ;
60 /// Stores the maximum difference in the last iteration
61 Real maxdiff;
62 /// Stores the computation time (in seconds)
63 double time;
64 /// Stores the number of iterations needed
65 size_t iters;
66 /// Does the InfAlg support logZ()?
67 bool has_logZ;
68 /// Does the InfAlg support maxDiff()?
69 bool has_maxdiff;
70 /// Does the InfAlg support Iterations()?
71 bool has_iters;
72
73 /// Construct from factor graph \a fg, name \a _name, and set of properties \a opts
74 TestDAI( const FactorGraph &fg, const string &_name, const PropertySet &opts ) : obj(NULL), name(_name), varErr(), facErr(), varMarginals(), facMarginals(), allMarginals(), logZ(0.0), maxdiff(0.0), time(0), iters(0U), has_logZ(false), has_maxdiff(false), has_iters(false) {
75 double tic = toc();
76
77 if( name == "LDPC" ) {
78 // special case: simulating a Low Density Parity Check code
79 Real zero[2] = {1.0, 0.0};
80 for( size_t i = 0; i < fg.nrVars(); i++ )
81 varMarginals.push_back( Factor(fg.var(i), zero) );
82 allMarginals = varMarginals;
83 logZ = 0.0;
84 maxdiff = 0.0;
85 iters = 1;
86 has_logZ = false;
87 has_maxdiff = false;
88 has_iters = false;
89 } else
90 // create a corresponding InfAlg object
91 obj = newInfAlg( name, fg, opts );
92
93 // Add the time needed to create the object
94 time += toc() - tic;
95 }
96
97 /// Destructor
98 ~TestDAI() {
99 if( obj != NULL )
100 delete obj;
101 }
102
103 /// Identify
104 string identify() const {
105 if( obj != NULL )
106 return obj->identify();
107 else
108 return "NULL";
109 }
110
111 /// Run the algorithm and store its results
112 void doDAI() {
113 double tic = toc();
114 if( obj != NULL ) {
115 // Initialize
116 obj->init();
117 // Run
118 obj->run();
119 // Record the time
120 time += toc() - tic;
121
122 // Store logarithm of the partition sum (if supported)
123 try {
124 logZ = obj->logZ();
125 has_logZ = true;
126 } catch( Exception &e ) {
127 if( e.code() == Exception::NOT_IMPLEMENTED )
128 has_logZ = false;
129 else
130 throw;
131 }
132
133 // Store maximum difference encountered in last iteration (if supported)
134 try {
135 maxdiff = obj->maxDiff();
136 has_maxdiff = true;
137 } catch( Exception &e ) {
138 if( e.code() == Exception::NOT_IMPLEMENTED )
139 has_maxdiff = false;
140 else
141 throw;
142 }
143
144 // Store number of iterations needed (if supported)
145 try {
146 iters = obj->Iterations();
147 has_iters = true;
148 } catch( Exception &e ) {
149 if( e.code() == Exception::NOT_IMPLEMENTED )
150 has_iters = false;
151 else
152 throw;
153 }
154
155 // Store variable marginals
156 varMarginals.clear();
157 for( size_t i = 0; i < obj->fg().nrVars(); i++ )
158 varMarginals.push_back( obj->beliefV( i ) );
159
160 // Store factor marginals
161 facMarginals.clear();
162 for( size_t I = 0; I < obj->fg().nrFactors(); I++ )
163 try {
164 facMarginals.push_back( obj->beliefF( I ) );
165 } catch( Exception &e ) {
166 if( e.code() == Exception::BELIEF_NOT_AVAILABLE )
167 facMarginals.push_back( Factor( obj->fg().factor(I).vars(), INFINITY ) );
168 else
169 throw;
170 }
171
172 // Store all marginals calculated by the method
173 allMarginals = obj->beliefs();
174 };
175 }
176
177 /// Calculate total variation distance of variable and factor marginals with respect to those in \a varMargs and \a facMargs
178 void calcErrors( const vector<Factor>& varMargs, const vector<Factor>& facMargs ) {
179 varErr = calcDists( varMarginals, varMargs );
180 facErr = calcDists( facMarginals, facMargs );
181 }
182
183 /// Return maximum variable error
184 Real maxVarErr() {
185 return( *max_element( varErr.begin(), varErr.end() ) );
186 }
187
188 /// Return average variable error
189 Real avgVarErr() {
190 return( accumulate( varErr.begin(), varErr.end(), 0.0 ) / varErr.size() );
191 }
192
193 /// Return maximum factor error
194 Real maxFacErr() {
195 return( *max_element( facErr.begin(), facErr.end() ) );
196 }
197
198 /// Return average factor error
199 Real avgFacErr() {
200 return( accumulate( facErr.begin(), facErr.end(), 0.0 ) / facErr.size() );
201 }
202 };
203
204
205 /// Clips a real number: if the absolute value of \a x is less than \a minabs, return \a minabs, else return \a x
206 Real clipReal( Real x, Real minabs ) {
207 if( abs(x) < minabs )
208 return minabs;
209 else
210 return x;
211 }
212
213
214 /// Which marginals to outpu (none, only variable, only factor, variable and factor, all)
215 DAI_ENUM(MarginalsOutputType,NONE,VAR,FAC,VARFAC,ALL);
216
217
218 /// Main function
219 int main( int argc, char *argv[] ) {
220 // Variables to store command line options
221 // Filename of factor graph
222 string filename;
223 // Filename for aliases
224 string aliases;
225 // Approximate Inference methods to use
226 vector<string> methods;
227 // Which marginals to output
228 MarginalsOutputType marginals;
229 // Output number of iterations?
230 bool report_iters = true;
231 // Output calculation time?
232 bool report_time = true;
233
234 // Define required command line options
235 po::options_description opts_required("Required options");
236 opts_required.add_options()
237 ("filename", po::value< string >(&filename), "Filename of factor graph")
238 ("methods", po::value< vector<string> >(&methods)->multitoken(), "DAI methods to perform")
239 ;
240
241 // Define allowed command line options
242 po::options_description opts_optional("Allowed options");
243 opts_optional.add_options()
244 ("help", "Produce help message")
245 ("aliases", po::value< string >(&aliases), "Filename for aliases")
246 ("marginals", po::value< MarginalsOutputType >(&marginals), "Output marginals? (NONE/VAR/FAC/VARFAC/ALL, default=NONE)")
247 ("report-time", po::value< bool >(&report_time), "Output calculation time (default==1)?")
248 ("report-iters", po::value< bool >(&report_iters), "Output iterations needed (default==1)?")
249 ;
250
251 // Define all command line options
252 po::options_description cmdline_options;
253 cmdline_options.add(opts_required).add(opts_optional);
254
255 // Parse command line
256 po::variables_map vm;
257 po::store(po::parse_command_line(argc, argv, cmdline_options), vm);
258 po::notify(vm);
259
260 // Display help message if necessary
261 if( vm.count("help") || !(vm.count("filename") && vm.count("methods")) ) {
262 cout << "This program is part of libDAI - http://www.libdai.org/" << endl << endl;
263 cout << "Usage: ./testdai --filename <filename.fg> --methods <method1> [<method2> <method3> ...]" << endl << endl;
264 cout << "Reads factor graph <filename.fg> and performs the approximate inference algorithms" << endl;
265 cout << "<method*>, reporting for each method:" << endl;
266 cout << " o the calculation time needed, in seconds (if report-time == 1);" << endl;
267 cout << " o the number of iterations needed (if report-iters == 1);" << endl;
268 cout << " o the maximum (over all variables) total variation error in the variable marginals;" << endl;
269 cout << " o the average (over all variables) total variation error in the variable marginals;" << endl;
270 cout << " o the maximum (over all factors) total variation error in the factor marginals;" << endl;
271 cout << " o the average (over all factors) total variation error in the factor marginals;" << endl;
272 cout << " o the error (difference) of the logarithm of the partition sums;" << endl << endl;
273 cout << "All errors are calculated by comparing the results of the current method with" << endl;
274 cout << "the results of the first method (the base method). If marginals==VAR, additional" << endl;
275 cout << "output consists of the variable marginals, if marginals==FAC, the factor marginals" << endl;
276 cout << "if marginals==VARFAC, both variable and factor marginals, and if marginals==ALL, all" << endl;
277 cout << "marginals calculated by the method are reported." << endl << endl;
278 cout << "<method*> should be a list of one or more methods, seperated by spaces, in the format:" << endl << endl;
279 cout << " name[key1=val1,key2=val2,key3=val3,...,keyn=valn]" << endl << endl;
280 cout << "where name should be the name of an algorithm in libDAI (or an alias, if an alias" << endl;
281 cout << "filename is provided), followed by a list of properties (surrounded by rectangular" << endl;
282 cout << "brackets), where each property consists of a key=value pair and the properties are" << endl;
283 cout << "seperated by commas. If an alias file is specified, alias substitution is performed." << endl;
284 cout << "This is done by looking up the name in the alias file and substituting the alias" << endl;
285 cout << "by its corresponding method as defined in the alias file. Properties are parsed from" << endl;
286 cout << "left to right, so if a property occurs repeatedly, the right-most value is used." << endl << endl;
287 cout << opts_required << opts_optional << endl;
288 #ifdef DAI_DEBUG
289 cout << "Note: this is a debugging build of libDAI." << endl << endl;
290 #endif
291 cout << "Example: ./testdai --filename testfast.fg --aliases aliases.conf --methods JTREE_HUGIN BP_SEQFIX BP_PARALL[maxiter=5]" << endl;
292 return 1;
293 }
294
295 try {
296 // Read aliases
297 map<string,string> Aliases;
298 if( !aliases.empty() )
299 Aliases = readAliasesFile( aliases );
300
301 // Read factor graph
302 FactorGraph fg;
303 fg.ReadFromFile( filename.c_str() );
304
305 // Declare variables used for storing variable factor marginals and log partition sum of base method
306 vector<Factor> varMarginals0;
307 vector<Factor> facMarginals0;
308 Real logZ0 = 0.0;
309
310 // Output header
311 cout.setf( ios_base::scientific );
312 cout.precision( 3 );
313 cout << "# " << filename << endl;
314 cout.width( 39 );
315 cout << left << "# METHOD" << "\t";
316 if( report_time )
317 cout << right << "SECONDS " << "\t";
318 if( report_iters )
319 cout << "ITERS" << "\t";
320 cout << "MAX VAR ERR" << "\t";
321 cout << "AVG VAR ERR" << "\t";
322 cout << "MAX FAC ERR" << "\t";
323 cout << "AVG FAC ERR" << "\t";
324 cout << "LOGZ ERROR" << "\t";
325 cout << "MAXDIFF" << "\t";
326 cout << endl;
327
328 // For each method...
329 for( size_t m = 0; m < methods.size(); m++ ) {
330 // Parse method
331 pair<string, PropertySet> meth = parseNameProperties( methods[m], Aliases );
332
333 // Check whether name is valid
334 size_t n = 0;
335 for( ; strlen( DAINames[n] ) != 0; n++ )
336 if( meth.first == DAINames[n] )
337 break;
338 if( strlen( DAINames[n] ) == 0 )
339 DAI_THROWE(UNKNOWN_DAI_ALGORITHM,string("Unknown DAI algorithm \"") + meth.first + string("\" in \"") + methods[m] + string("\""));
340
341 // Construct object for running the method
342 TestDAI testdai(fg, meth.first, meth.second );
343
344 // Run the method
345 testdai.doDAI();
346
347 // For the base method, store its variable marginals and logarithm of the partition sum
348 if( m == 0 ) {
349 varMarginals0 = testdai.varMarginals;
350 facMarginals0 = testdai.facMarginals;
351 logZ0 = testdai.logZ;
352 }
353
354 // Calculate errors relative to base method
355 testdai.calcErrors( varMarginals0, facMarginals0 );
356
357 // Output method name
358 cout.width( 39 );
359 cout << left << methods[m] << "\t";
360 // Output calculation time, if requested
361 if( report_time )
362 cout << right << testdai.time << "\t";
363 // Output number of iterations, if requested
364 if( report_iters ) {
365 if( testdai.has_iters ) {
366 cout << testdai.iters << "\t";
367 } else {
368 cout << "N/A \t";
369 }
370 }
371
372 // If this is not the base method
373 if( m > 0 ) {
374 cout.setf( ios_base::scientific );
375 cout.precision( 3 );
376
377 // Output maximum error in variable marginals
378 Real mev = clipReal( testdai.maxVarErr(), 1e-9 );
379 cout << mev << "\t";
380
381 // Output average error in variable marginals
382 Real aev = clipReal( testdai.avgVarErr(), 1e-9 );
383 cout << aev << "\t";
384
385 // Output maximum error in factor marginals
386 Real mef = clipReal( testdai.maxFacErr(), 1e-9 );
387 if( mef == INFINITY )
388 cout << "N/A \t";
389 else
390 cout << mef << "\t";
391
392 // Output average error in factor marginals
393 Real aef = clipReal( testdai.avgFacErr(), 1e-9 );
394 if( aef == INFINITY )
395 cout << "N/A \t";
396 else
397 cout << aef << "\t";
398
399 // Output error in log partition sum
400 if( testdai.has_logZ ) {
401 cout.setf( ios::showpos );
402 Real le = clipReal( testdai.logZ - logZ0, 1e-9 );
403 cout << le << "\t";
404 cout.unsetf( ios::showpos );
405 } else
406 cout << "N/A \t";
407
408 // Output maximum difference in last iteration
409 if( testdai.has_maxdiff ) {
410 Real md = clipReal( testdai.maxdiff, 1e-9 );
411 if( dai::isnan( mev ) )
412 md = mev;
413 if( dai::isnan( aev ) )
414 md = aev;
415 if( md == INFINITY )
416 md = 1.0;
417 cout << md << "\t";
418 } else
419 cout << "N/A \t";
420 }
421 cout << endl;
422
423 // Output marginals, if requested
424 if( marginals == MarginalsOutputType::VAR || marginals == MarginalsOutputType::VARFAC )
425 for( size_t i = 0; i < testdai.varMarginals.size(); i++ )
426 cout << "# " << testdai.varMarginals[i] << endl;
427 if( marginals == MarginalsOutputType::FAC || marginals == MarginalsOutputType::VARFAC )
428 for( size_t I = 0; I < testdai.facMarginals.size(); I++ )
429 cout << "# " << testdai.facMarginals[I] << endl;
430 if( marginals == MarginalsOutputType::ALL )
431 for( size_t I = 0; I < testdai.allMarginals.size(); I++ )
432 cout << "# " << testdai.allMarginals[I] << endl;
433 }
434
435 return 0;
436 } catch( string &s ) {
437 // Abort with error message
438 cerr << "Exception: " << s << endl;
439 return 2;
440 }
441 }