examples/example.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-2009  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 <dai/alldai.h>  // Include main libDAI header file
  14
  15
  16 using namespace dai;
  17 using namespace std;
  18
  19
  20 int main( int argc, char *argv[] ) {
  21     if ( argc != 2 ) {
  22         cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
  23         cout << "Reads factor graph <filename.fg> and runs" << endl;
  24         cout << "Belief Propagation and JunctionTree on it." << endl << endl;
  25         return 1;
  26     } else {
  27         // Read FactorGraph from the file specified by the first command line argument
  28         FactorGraph fg;
  29         fg.ReadFromFile(argv[1]);
  30
  31         // Set some constants
  32         size_t maxiter = 10000;
  33         Real   tol = 1e-9;
  34         size_t verb = 1;
  35
  36         // Store the constants in a PropertySet object
  37         PropertySet opts;
  38         opts.Set("maxiter",maxiter);  // Maximum number of iterations
  39         opts.Set("tol",tol);          // Tolerance for convergence
  40         opts.Set("verbose",verb);     // Verbosity (amount of output generated)
  41
  42         // Construct a JTree (junction tree) object from the FactorGraph fg
  43         // using the parameters specified by opts and an additional property
  44         // that specifies the type of updates the JTree algorithm should perform
  45         JTree jt( fg, opts("updates",string("HUGIN")) );
  46         // Initialize junction tree algorithm
  47         jt.init();
  48         // Run junction tree algorithm
  49         jt.run();
  50
  51         // Construct another JTree (junction tree) object that is used to calculate
  52         // the joint configuration of variables that has maximum probability (MAP state)
  53         JTree jtmap( fg, opts("updates",string("HUGIN"))("inference",string("MAXPROD")) );
  54         // Initialize junction tree algorithm
  55         jtmap.init();
  56         // Run junction tree algorithm
  57         jtmap.run();
  58         // Calculate joint state of all variables that has maximum probability
  59         vector<size_t> jtmapstate = jtmap.findMaximum();
  60
  61         // Construct a BP (belief propagation) object from the FactorGraph fg
  62         // using the parameters specified by opts and two additional properties,
  63         // specifying the type of updates the BP algorithm should perform and
  64         // whether they should be done in the real or in the logdomain
  65         BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false));
  66         // Initialize belief propagation algorithm
  67         bp.init();
  68         // Run belief propagation algorithm
  69         bp.run();
  70
  71         // Construct a BP (belief propagation) object from the FactorGraph fg
  72         // using the parameters specified by opts and two additional properties,
  73         // specifying the type of updates the BP algorithm should perform and
  74         // whether they should be done in the real or in the logdomain
  75         //
  76         // Note that inference is set to MAXPROD, which means that the object
  77         // will perform the max-product algorithm instead of the sum-product algorithm
  78         BP mp(fg, opts("updates",string("SEQFIX"))("logdomain",false)("inference",string("MAXPROD")));
  79         // Initialize max-product algorithm
  80         mp.init();
  81         // Run max-product algorithm
  82         mp.run();
  83         // Calculate joint state of all variables that has maximum probability
  84         // based on the max-product result
  85         vector<size_t> mpstate = mp.findMaximum();
  86
  87         // Report variable marginals for fg, calculated by the junction tree algorithm
  88         cout << "Exact variable marginals:" << endl;
  89         for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
  90             cout << jt.belief(fg.var(i)) << endl; // display the "belief" of jt for that variable
  91
  92         // Report variable marginals for fg, calculated by the belief propagation algorithm
  93         cout << "Approximate (loopy belief propagation) variable marginals:" << endl;
  94         for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
  95             cout << bp.belief(fg.var(i)) << endl; // display the belief of bp for that variable
  96
  97         // Report factor marginals for fg, calculated by the junction tree algorithm
  98         cout << "Exact factor marginals:" << endl;
  99         for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 100             cout << jt.belief(fg.factor(I).vars()) << endl;  // display the "belief" of jt for the variables in that factor
 101
 102         // Report factor marginals for fg, calculated by the belief propagation algorithm
 103         cout << "Approximate (loopy belief propagation) factor marginals:" << endl;
 104         for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 105             cout << bp.belief(fg.factor(I).vars()) << endl; // display the belief of bp for the variables in that factor
 106
 107         // Report log partition sum (normalizing constant) of fg, calculated by the junction tree algorithm
 108         cout << "Exact log partition sum: " << jt.logZ() << endl;
 109
 110         // Report log partition sum of fg, approximated by the belief propagation algorithm
 111         cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
 112
 113         // Report exact MAP variable marginals
 114         cout << "Exact MAP variable marginals:" << endl;
 115         for( size_t i = 0; i < fg.nrVars(); i++ )
 116             cout << jtmap.belief(fg.var(i)) << endl;
 117
 118         // Report max-product variable marginals
 119         cout << "Approximate (max-product) MAP variable marginals:" << endl;
 120         for( size_t i = 0; i < fg.nrVars(); i++ )
 121             cout << mp.belief(fg.var(i)) << endl;
 122
 123         // Report exact MAP factor marginals
 124         cout << "Exact MAP factor marginals:" << endl;
 125         for( size_t I = 0; I < fg.nrFactors(); I++ )
 126             cout << jtmap.belief(fg.factor(I).vars()) << "=" << jtmap.beliefF(I) << endl;
 127
 128         // Report max-product factor marginals
 129         cout << "Approximate (max-product) MAP factor marginals:" << endl;
 130         for( size_t I = 0; I < fg.nrFactors(); I++ )
 131             cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
 132
 133         // Report exact MAP joint state
 134         cout << "Exact MAP state:" << endl;
 135         for( size_t i = 0; i < jtmapstate.size(); i++ )
 136             cout << fg.var(i) << ": " << jtmapstate[i] << endl;
 137
 138         // Report max-product MAP joint state
 139         cout << "Approximate (max-product) MAP state:" << endl;
 140         for( size_t i = 0; i < mpstate.size(); i++ )
 141             cout << fg.var(i) << ": " << mpstate[i] << endl;
 142     }
 143
 144     return 0;
 145 }