Replaced doubles by Reals, fixed two bugs
[libdai.git] / examples / example.cpp
index 0705bde..210a293 100644 (file)
@@ -1,23 +1,12 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
-    Radboud University Nijmegen, The Netherlands /
-    Max Planck Institute for Biological Cybernetics, Germany
-
-    This file is part of libDAI.
-
-    libDAI is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    libDAI is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  libDAI is licensed under the terms of the GNU General Public License version
+ *  2, or (at your option) any later version. libDAI is distributed without any
+ *  warranty. See the file COPYING for more details.
+ *
+ *  Copyright (C) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ */
 
 
 #include <iostream>
@@ -40,9 +29,9 @@ int main( int argc, char *argv[] ) {
         fg.ReadFromFile(argv[1]);
 
         // Set some constants
-        size_t  maxiter = 10000;
-        double  tol = 1e-9;
-        size_t  verb = 1;
+        size_t maxiter = 10000;
+        Real   tol = 1e-9;
+        size_t verb = 1;
 
         // Store the constants in a PropertySet object
         PropertySet opts;
@@ -54,11 +43,21 @@ int main( int argc, char *argv[] ) {
         // using the parameters specified by opts and an additional property
         // that specifies the type of updates the JTree algorithm should perform
         JTree jt( fg, opts("updates",string("HUGIN")) );
-        // Initialize junction tree algoritm
+        // Initialize junction tree algorithm
         jt.init();
         // Run junction tree algorithm
         jt.run();
 
+        // Construct another JTree (junction tree) object that is used to calculate
+        // the joint configuration of variables that has maximum probability (MAP state)
+        JTree jtmap( fg, opts("updates",string("HUGIN"))("inference",string("MAXPROD")) );
+        // Initialize junction tree algorithm
+        jtmap.init();
+        // Run junction tree algorithm
+        jtmap.run();
+        // Calculate joint state of all variables that has maximum probability
+        vector<size_t> jtmapstate = jtmap.findMaximum();
+
         // Construct a BP (belief propagation) object from the FactorGraph fg
         // using the parameters specified by opts and two additional properties,
         // specifying the type of updates the BP algorithm should perform and
@@ -111,18 +110,33 @@ int main( int argc, char *argv[] ) {
         // Report log partition sum of fg, approximated by the belief propagation algorithm
         cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
 
+        // Report exact MAP variable marginals
+        cout << "Exact MAP variable marginals:" << endl;
+        for( size_t i = 0; i < fg.nrVars(); i++ )
+            cout << jtmap.belief(fg.var(i)) << endl;
+
         // Report max-product variable marginals
-        cout << "Max-product variable marginals:" << endl;
+        cout << "Approximate (max-product) MAP variable marginals:" << endl;
         for( size_t i = 0; i < fg.nrVars(); i++ )
             cout << mp.belief(fg.var(i)) << endl;
 
+        // Report exact MAP factor marginals
+        cout << "Exact MAP factor marginals:" << endl;
+        for( size_t I = 0; I < fg.nrFactors(); I++ )
+            cout << jtmap.belief(fg.factor(I).vars()) << "=" << jtmap.beliefF(I) << endl;
+
         // Report max-product factor marginals
-        cout << "Max-product factor marginals:" << endl;
+        cout << "Approximate (max-product) MAP factor marginals:" << endl;
         for( size_t I = 0; I < fg.nrFactors(); I++ )
             cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
 
-        // Report max-product joint state
-        cout << "Max-product state:" << endl;
+        // Report exact MAP joint state
+        cout << "Exact MAP state:" << endl;
+        for( size_t i = 0; i < jtmapstate.size(); i++ )
+            cout << fg.var(i) << ": " << jtmapstate[i] << endl;
+
+        // Report max-product MAP joint state
+        cout << "Approximate (max-product) MAP state:" << endl;
         for( size_t i = 0; i < mpstate.size(); i++ )
             cout << fg.var(i) << ": " << mpstate[i] << endl;
     }