examples/example.cpp

   1 /*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
   2     Radboud University Nijmegen, The Netherlands /
   3     Max Planck Institute for Biological Cybernetics, Germany
   4
   5     This file is part of libDAI.
   6
   7     libDAI is free software; you can redistribute it and/or modify
   8     it under the terms of the GNU General Public License as published by
   9     the Free Software Foundation; either version 2 of the License, or
  10     (at your option) any later version.
  11
  12     libDAI is distributed in the hope that it will be useful,
  13     but WITHOUT ANY WARRANTY; without even the implied warranty of
  14     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15     GNU General Public License for more details.
  16
  17     You should have received a copy of the GNU General Public License
  18     along with libDAI; if not, write to the Free Software
  19     Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 */
  21
  22
  23 #include <iostream>
  24 #include <dai/alldai.h>  // Include main libDAI header file
  25
  26
  27 using namespace dai;
  28 using namespace std;
  29
  30
  31 int main( int argc, char *argv[] ) {
  32     if ( argc != 2 ) {
  33         cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
  34         cout << "Reads factor graph <filename.fg> and runs" << endl;
  35         cout << "Belief Propagation and JunctionTree on it." << endl << endl;
  36         return 1;
  37     } else {
  38         // Read FactorGraph from the file specified by the first command line argument
  39         FactorGraph fg;
  40         fg.ReadFromFile(argv[1]);
  41
  42         // Set some constants
  43         size_t  maxiter = 10000;
  44         double  tol = 1e-9;
  45         size_t  verb = 1;
  46
  47         // Store the constants in a PropertySet object
  48         PropertySet opts;
  49         opts.Set("maxiter",maxiter);  // Maximum number of iterations
  50         opts.Set("tol",tol);          // Tolerance for convergence
  51         opts.Set("verbose",verb);     // Verbosity (amount of output generated)
  52
  53         // Construct a JTree (junction tree) object from the FactorGraph fg
  54         // using the parameters specified by opts and an additional property
  55         // that specifies the type of updates the JTree algorithm should perform
  56         JTree jt( fg, opts("updates",string("HUGIN")) );
  57         // Initialize junction tree algorithm
  58         jt.init();
  59         // Run junction tree algorithm
  60         jt.run();
  61
  62         // Construct another JTree (junction tree) object that is used to calculate
  63         // the joint configuration of variables that has maximum probability (MAP state)
  64         JTree jtmap( fg, opts("updates",string("HUGIN"))("inference",string("MAXPROD")) );
  65         // Initialize junction tree algorithm
  66         jtmap.init();
  67         // Run junction tree algorithm
  68         jtmap.run();
  69         // Calculate joint state of all variables that has maximum probability
  70         vector<size_t> jtmapstate = jtmap.findMaximum();
  71
  72         // Construct a BP (belief propagation) object from the FactorGraph fg
  73         // using the parameters specified by opts and two additional properties,
  74         // specifying the type of updates the BP algorithm should perform and
  75         // whether they should be done in the real or in the logdomain
  76         BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false));
  77         // Initialize belief propagation algorithm
  78         bp.init();
  79         // Run belief propagation algorithm
  80         bp.run();
  81
  82         // Construct a BP (belief propagation) object from the FactorGraph fg
  83         // using the parameters specified by opts and two additional properties,
  84         // specifying the type of updates the BP algorithm should perform and
  85         // whether they should be done in the real or in the logdomain
  86         //
  87         // Note that inference is set to MAXPROD, which means that the object
  88         // will perform the max-product algorithm instead of the sum-product algorithm
  89         BP mp(fg, opts("updates",string("SEQFIX"))("logdomain",false)("inference",string("MAXPROD")));
  90         // Initialize max-product algorithm
  91         mp.init();
  92         // Run max-product algorithm
  93         mp.run();
  94         // Calculate joint state of all variables that has maximum probability
  95         // based on the max-product result
  96         vector<size_t> mpstate = mp.findMaximum();
  97
  98         // Report variable marginals for fg, calculated by the junction tree algorithm
  99         cout << "Exact variable marginals:" << endl;
 100         for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
 101             cout << jt.belief(fg.var(i)) << endl; // display the "belief" of jt for that variable
 102
 103         // Report variable marginals for fg, calculated by the belief propagation algorithm
 104         cout << "Approximate (loopy belief propagation) variable marginals:" << endl;
 105         for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
 106             cout << bp.belief(fg.var(i)) << endl; // display the belief of bp for that variable
 107
 108         // Report factor marginals for fg, calculated by the junction tree algorithm
 109         cout << "Exact factor marginals:" << endl;
 110         for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 111             cout << jt.belief(fg.factor(I).vars()) << endl;  // display the "belief" of jt for the variables in that factor
 112
 113         // Report factor marginals for fg, calculated by the belief propagation algorithm
 114         cout << "Approximate (loopy belief propagation) factor marginals:" << endl;
 115         for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 116             cout << bp.belief(fg.factor(I).vars()) << endl; // display the belief of bp for the variables in that factor
 117
 118         // Report log partition sum (normalizing constant) of fg, calculated by the junction tree algorithm
 119         cout << "Exact log partition sum: " << jt.logZ() << endl;
 120
 121         // Report log partition sum of fg, approximated by the belief propagation algorithm
 122         cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
 123
 124         // Report exact MAP variable marginals
 125         cout << "Exact MAP variable marginals:" << endl;
 126         for( size_t i = 0; i < fg.nrVars(); i++ )
 127             cout << jtmap.belief(fg.var(i)) << endl;
 128
 129         // Report max-product variable marginals
 130         cout << "Approximate (max-product) MAP variable marginals:" << endl;
 131         for( size_t i = 0; i < fg.nrVars(); i++ )
 132             cout << mp.belief(fg.var(i)) << endl;
 133
 134         // Report exact MAP factor marginals
 135         cout << "Exact MAP factor marginals:" << endl;
 136         for( size_t I = 0; I < fg.nrFactors(); I++ )
 137             cout << jtmap.belief(fg.factor(I).vars()) << "=" << jtmap.beliefF(I) << endl;
 138
 139         // Report max-product factor marginals
 140         cout << "Approximate (max-product) MAP factor marginals:" << endl;
 141         for( size_t I = 0; I < fg.nrFactors(); I++ )
 142             cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
 143
 144         // Report exact MAP joint state
 145         cout << "Exact MAP state:" << endl;
 146         for( size_t i = 0; i < jtmapstate.size(); i++ )
 147             cout << fg.var(i) << ": " << jtmapstate[i] << endl;
 148
 149         // Report max-product MAP joint state
 150         cout << "Approximate (max-product) MAP state:" << endl;
 151         for( size_t i = 0; i < mpstate.size(); i++ )
 152             cout << fg.var(i) << ": " << mpstate[i] << endl;
 153     }
 154
 155     return 0;
 156 }