tests/test.cpp

   1 /*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
   2     Radboud University Nijmegen, The Netherlands
   3
   4     This file is part of libDAI.
   5
   6     libDAI is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     libDAI is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with libDAI; if not, write to the Free Software
  18     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19 */
  20
  21
  22 #include <iostream>
  23 #include <fstream>
  24 #include <map>
  25 #include <numeric>
  26 #include <cmath>
  27 #include <cstdlib>
  28 #include <cstring>
  29 #include <boost/program_options.hpp>
  30 #include <dai/util.h>
  31 #include <dai/alldai.h>
  32
  33
  34 using namespace std;
  35 using namespace dai;
  36 namespace po = boost::program_options;
  37
  38
  39 class TestDAI {
  40     protected:
  41         InfAlg          *obj;
  42         string          name;
  43         vector<double>  err;
  44
  45     public:
  46         vector<Factor>  q;
  47         double          logZ;
  48         double          maxdiff;
  49         double          time;
  50         size_t          iters;
  51         bool            has_logZ;
  52         bool            has_maxdiff;
  53         bool            has_iters;
  54
  55         TestDAI( const FactorGraph &fg, const string &_name, const PropertySet &opts ) : obj(NULL), name(_name), err(), q(), logZ(0.0), maxdiff(0.0), time(0), iters(0U), has_logZ(false), has_maxdiff(false), has_iters(false) {
  56             double tic = toc();
  57             if( name == "LDPC" ) {
  58                 double zero[2] = {1.0, 0.0};
  59                 q.clear();
  60                 for( size_t i = 0; i < fg.nrVars(); i++ )
  61                     q.push_back( Factor(Var(i,2), zero) );
  62                 logZ = 0.0;
  63                 maxdiff = 0.0;
  64                 iters = 1;
  65                 has_logZ = false;
  66                 has_maxdiff = false;
  67                 has_iters = false;
  68             } else
  69                 obj = newInfAlg( name, fg, opts );
  70             time += toc() - tic;
  71         }
  72
  73         ~TestDAI() {
  74             if( obj != NULL )
  75                 delete obj;
  76         }
  77
  78         string identify() {
  79             if( obj != NULL )
  80                 return obj->identify();
  81             else
  82                 return "NULL";
  83         }
  84
  85         vector<Factor> allBeliefs() {
  86             vector<Factor> result;
  87             for( size_t i = 0; i < obj->fg().nrVars(); i++ )
  88                 result.push_back( obj->belief( obj->fg().var(i) ) );
  89             return result;
  90         }
  91
  92         void doDAI() {
  93             double tic = toc();
  94             if( obj != NULL ) {
  95                 obj->init();
  96                 obj->run();
  97                 time += toc() - tic;
  98                 try {
  99                     logZ = obj->logZ();
 100                     has_logZ = true;
 101                 } catch( Exception &e ) {
 102                     has_logZ = false;
 103                 }
 104                 try {
 105                     maxdiff = obj->maxDiff();
 106                     has_maxdiff = true;
 107                 } catch( Exception &e ) {
 108                     has_maxdiff = false;
 109                 }
 110                 has_iters = false;
 111                 iters = 0;
 112                 q = allBeliefs();
 113             };
 114         }
 115
 116         void calcErrs( const TestDAI &x ) {
 117             err.clear();
 118             err.reserve( q.size() );
 119             for( size_t i = 0; i < q.size(); i++ )
 120                 err.push_back( dist( q[i], x.q[i], Prob::DISTTV ) );
 121         }
 122
 123         void calcErrs( const vector<Factor> &x ) {
 124             err.clear();
 125             err.reserve( q.size() );
 126             for( size_t i = 0; i < q.size(); i++ )
 127                 err.push_back( dist( q[i], x[i], Prob::DISTTV ) );
 128         }
 129
 130         double maxErr() {
 131             return( *max_element( err.begin(), err.end() ) );
 132         }
 133
 134         double avgErr() {
 135             return( accumulate( err.begin(), err.end(), 0.0 ) / err.size() );
 136         }
 137 };
 138
 139
 140 pair<string, PropertySet> parseMethod( const string &_s, const map<string,string> & aliases ) {
 141     // s = first part of _s, until '['
 142     string::size_type pos = _s.find_first_of('[');
 143     string s;
 144     if( pos == string::npos )
 145         s = _s;
 146     else
 147         s = _s.substr(0,pos);
 148
 149     // if the first part is an alias, substitute
 150     if( aliases.find(s) != aliases.end() )
 151         s = aliases.find(s)->second;
 152
 153     // attach second part, merging properties if necessary
 154     if( pos != string::npos ) {
 155         if( s.at(s.length()-1) == ']' ) {
 156             s = s.erase(s.length()-1,1) + ',' + _s.substr(pos+1);
 157         } else
 158             s = s + _s.substr(pos);
 159     }
 160
 161     pair<string, PropertySet> result;
 162     string & name = result.first;
 163     PropertySet & opts = result.second;
 164
 165     pos = s.find_first_of('[');
 166     if( pos == string::npos )
 167         throw "Malformed method";
 168     name = s.substr( 0, pos );
 169     size_t n = 0;
 170     for( ; strlen( DAINames[n] ) != 0; n++ )
 171         if( name == DAINames[n] )
 172             break;
 173     if( strlen( DAINames[n] ) == 0 && (name != "LDPC") )
 174         DAI_THROW(UNKNOWN_DAI_ALGORITHM);
 175
 176     stringstream ss;
 177     ss << s.substr(pos,s.length());
 178     ss >> opts;
 179
 180     return result;
 181 }
 182
 183
 184 double clipdouble( double x, double minabs ) {
 185     if( fabs(x) < minabs )
 186         return minabs;
 187     else
 188         return x;
 189 }
 190
 191
 192 int main( int argc, char *argv[] ) {
 193     try {
 194         string filename;
 195         string aliases;
 196         vector<string> methods;
 197         double tol;
 198         size_t maxiter;
 199         size_t verbose;
 200         bool marginals = false;
 201         bool report_time = true;
 202
 203         po::options_description opts_required("Required options");
 204         opts_required.add_options()
 205             ("filename", po::value< string >(&filename), "Filename of FactorGraph")
 206             ("methods", po::value< vector<string> >(&methods)->multitoken(), "DAI methods to test")
 207         ;
 208
 209         po::options_description opts_optional("Allowed options");
 210         opts_optional.add_options()
 211             ("help", "produce help message")
 212             ("aliases", po::value< string >(&aliases), "Filename for aliases")
 213             ("tol", po::value< double >(&tol), "Override tolerance")
 214             ("maxiter", po::value< size_t >(&maxiter), "Override maximum number of iterations")
 215             ("verbose", po::value< size_t >(&verbose), "Override verbosity")
 216             ("marginals", po::value< bool >(&marginals), "Output single node marginals?")
 217             ("report-time", po::value< bool >(&report_time), "Report calculation time")
 218         ;
 219
 220         po::options_description cmdline_options;
 221         cmdline_options.add(opts_required).add(opts_optional);
 222
 223         po::variables_map vm;
 224         po::store(po::parse_command_line(argc, argv, cmdline_options), vm);
 225         po::notify(vm);
 226
 227         if( vm.count("help") || !(vm.count("filename") && vm.count("methods")) ) {
 228             cout << "Reads factorgraph <filename.fg> and performs the approximate" << endl;
 229             cout << "inference algorithms <method*>, reporting calculation time, max and average" << endl;
 230             cout << "error and relative logZ error (comparing with the results of" << endl;
 231             cout << "<method0>, the base method, for which one can use JTREE_HUGIN)." << endl << endl;
 232             cout << opts_required << opts_optional << endl;
 233             return 1;
 234         }
 235
 236         // Read aliases
 237         map<string,string> Aliases;
 238         if( !aliases.empty() ) {
 239             ifstream infile;
 240             infile.open (aliases.c_str());
 241             if (infile.is_open()) {
 242                 while( true ) {
 243                     string line;
 244                     getline(infile,line);
 245                     if( infile.fail() )
 246                         break;
 247                     if( (!line.empty()) && (line[0] != '#') ) {
 248                         string::size_type pos = line.find(':',0);
 249                         if( pos == string::npos )
 250                             throw "Invalid alias";
 251                         else {
 252                             string::size_type posl = line.substr(0, pos).find_last_not_of(" \t");
 253                             string key = line.substr(0, posl + 1);
 254                             string::size_type posr = line.substr(pos + 1, line.length()).find_first_not_of(" \t");
 255                             string val = line.substr(pos + 1 + posr, line.length());
 256                             Aliases[key] = val;
 257                         }
 258                     }
 259                 }
 260                 infile.close();
 261             } else
 262                 throw "Error opening aliases file";
 263         }
 264
 265         FactorGraph fg;
 266         fg.ReadFromFile( filename.c_str() );
 267
 268         vector<Factor> q0;
 269         double logZ0 = 0.0;
 270
 271         cout << "# " << filename << endl;
 272         cout.width( 40 );
 273         cout << left << "# METHOD" << "  ";
 274         if( report_time ) {
 275             cout.width( 10 );
 276             cout << right << "SECONDS" << "   ";
 277         }
 278         cout.width( 10 );
 279         cout << "MAX ERROR" << "  ";
 280         cout.width( 10 );
 281         cout << "AVG ERROR" << "  ";
 282         cout.width( 10 );
 283         cout << "LOGZ ERROR" << "  ";
 284         cout.width( 10 );
 285         cout << "MAXDIFF" << "  ";
 286         cout.width( 10 );
 287         cout << "ITERS" << endl;
 288
 289         for( size_t m = 0; m < methods.size(); m++ ) {
 290             pair<string, PropertySet> meth = parseMethod( methods[m], Aliases );
 291
 292             if( vm.count("tol") )
 293                 meth.second.Set("tol",tol);
 294             if( vm.count("maxiter") )
 295                 meth.second.Set("maxiter",maxiter);
 296             if( vm.count("verbose") )
 297                 meth.second.Set("verbose",verbose);
 298             TestDAI piet(fg, meth.first, meth.second );
 299             piet.doDAI();
 300             if( m == 0 ) {
 301                 q0 = piet.q;
 302                 logZ0 = piet.logZ;
 303             }
 304             piet.calcErrs(q0);
 305
 306             cout.width( 40 );
 307 //                cout << left << piet.identify() << "  ";
 308             cout << left << methods[m] << "  ";
 309             if( report_time ) {
 310                 cout.width( 10 );
 311                 cout << right << piet.time << "    ";
 312             }
 313
 314             if( m > 0 ) {
 315                 cout.setf( ios_base::scientific );
 316                 cout.precision( 3 );
 317                 cout.width( 10 );
 318                 double me = clipdouble( piet.maxErr(), 1e-9 );
 319                 cout << me << "  ";
 320                 cout.width( 10 );
 321                 double ae = clipdouble( piet.avgErr(), 1e-9 );
 322                 cout << ae << "  ";
 323                 cout.width( 10 );
 324                 if( piet.has_logZ ) {
 325                     double le = clipdouble( piet.logZ / logZ0 - 1.0, 1e-9 );
 326                     cout << le << "  ";
 327                 } else {
 328                     cout << "N/A         ";
 329                 }
 330                 cout.width( 10 );
 331                 if( piet.has_maxdiff ) {
 332                     double md = clipdouble( piet.maxdiff, 1e-9 );
 333                     if( isnan( me ) )
 334                         md = me;
 335                     if( isnan( ae ) )
 336                         md = ae;
 337                     cout << md << "  ";
 338                 } else {
 339                     cout << "N/A         ";
 340                 }
 341                 cout.width( 10 );
 342                 if( piet.has_iters ) {
 343                     cout << piet.iters << "  " << endl;
 344                 } else {
 345                     cout << "N/A         ";
 346                 }
 347             }
 348             cout << endl;
 349
 350             if( marginals ) {
 351                 for( size_t i = 0; i < piet.q.size(); i++ )
 352                     cout << "# " << piet.q[i] << endl;
 353             }
 354         }
 355     } catch(const char *e) {
 356         cerr << "Exception: " << e << endl;
 357         return 1;
 358     } catch(exception& e) {
 359         cerr << "Exception: " << e.what() << endl;
 360         return 1;
 361     }
 362     catch(...) {
 363         cerr << "Exception of unknown type!" << endl;
 364     }
 365
 366     return 0;
 367 }