Merge branch 'joris'
[libdai.git] / 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 algoritm
58 jt.init();
59 // Run junction tree algorithm
60 jt.run();
61
62 // Construct a BP (belief propagation) object from the FactorGraph fg
63 // using the parameters specified by opts and two additional properties,
64 // specifying the type of updates the BP algorithm should perform and
65 // whether they should be done in the real or in the logdomain
66 BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false));
67 // Initialize belief propagation algorithm
68 bp.init();
69 // Run belief propagation algorithm
70 bp.run();
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 //
77 // Note that inference is set to MAXPROD, which means that the object
78 // will perform the max-product algorithm instead of the sum-product algorithm
79 BP mp(fg, opts("updates",string("SEQFIX"))("logdomain",false)("inference",string("MAXPROD")));
80 // Initialize max-product algorithm
81 mp.init();
82 // Run max-product algorithm
83 mp.run();
84 // Calculate joint state of all variables that has maximum probability
85 // based on the max-product result
86 vector<size_t> mpstate = mp.findMaximum();
87
88 // Report variable marginals for fg, calculated by the junction tree algorithm
89 cout << "Exact variable marginals:" << endl;
90 for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
91 cout << jt.belief(fg.var(i)) << endl; // display the "belief" of jt for that variable
92
93 // Report variable marginals for fg, calculated by the belief propagation algorithm
94 cout << "Approximate (loopy belief propagation) variable marginals:" << endl;
95 for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
96 cout << bp.belief(fg.var(i)) << endl; // display the belief of bp for that variable
97
98 // Report factor marginals for fg, calculated by the junction tree algorithm
99 cout << "Exact factor marginals:" << endl;
100 for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
101 cout << jt.belief(fg.factor(I).vars()) << endl; // display the "belief" of jt for the variables in that factor
102
103 // Report factor marginals for fg, calculated by the belief propagation algorithm
104 cout << "Approximate (loopy belief propagation) factor marginals:" << endl;
105 for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
106 cout << bp.belief(fg.factor(I).vars()) << endl; // display the belief of bp for the variables in that factor
107
108 // Report log partition sum (normalizing constant) of fg, calculated by the junction tree algorithm
109 cout << "Exact log partition sum: " << jt.logZ() << endl;
110
111 // Report log partition sum of fg, approximated by the belief propagation algorithm
112 cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
113
114 // Report max-product variable marginals
115 cout << "Max-product variable marginals:" << endl;
116 for( size_t i = 0; i < fg.nrVars(); i++ )
117 cout << mp.belief(fg.var(i)) << endl;
118
119 // Report max-product factor marginals
120 cout << "Max-product factor marginals:" << endl;
121 for( size_t I = 0; I < fg.nrFactors(); I++ )
122 cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
123
124 // Report max-product joint state
125 cout << "Max-product state:" << endl;
126 for( size_t i = 0; i < mpstate.size(); i++ )
127 cout << fg.var(i) << ": " << mpstate[i] << endl;
128 }
129
130 return 0;
131 }