Some small documentation updates
[libdai.git] / 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 <map>
14 #include <dai/alldai.h> // Include main libDAI header file
15
16
17 using namespace dai;
18 using namespace std;
19
20
21 int main( int argc, char *argv[] ) {
22 if ( argc != 2 ) {
23 cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
24 cout << "Reads factor graph <filename.fg> and runs" << endl;
25 cout << "Belief Propagation and JunctionTree on it." << endl << endl;
26 return 1;
27 } else {
28 // Read FactorGraph from the file specified by the first command line argument
29 FactorGraph fg;
30 fg.ReadFromFile(argv[1]);
31
32 // Set some constants
33 size_t maxiter = 10000;
34 Real tol = 1e-9;
35 size_t verb = 1;
36
37 // Store the constants in a PropertySet object
38 PropertySet opts;
39 opts.set("maxiter",maxiter); // Maximum number of iterations
40 opts.set("tol",tol); // Tolerance for convergence
41 opts.set("verbose",verb); // Verbosity (amount of output generated)
42
43 // Construct a JTree (junction tree) object from the FactorGraph fg
44 // using the parameters specified by opts and an additional property
45 // that specifies the type of updates the JTree algorithm should perform
46 JTree jt( fg, opts("updates",string("HUGIN")) );
47 // Initialize junction tree algorithm
48 jt.init();
49 // Run junction tree algorithm
50 jt.run();
51
52 // Construct another JTree (junction tree) object that is used to calculate
53 // the joint configuration of variables that has maximum probability (MAP state)
54 JTree jtmap( fg, opts("updates",string("HUGIN"))("inference",string("MAXPROD")) );
55 // Initialize junction tree algorithm
56 jtmap.init();
57 // Run junction tree algorithm
58 jtmap.run();
59 // Calculate joint state of all variables that has maximum probability
60 vector<size_t> jtmapstate = jtmap.findMaximum();
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("SEQRND"))("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("SEQRND"))("logdomain",false)("inference",string("MAXPROD"))("damping",string("0.1")));
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 // Construct a decimation algorithm object from the FactorGraph fg
89 // using the parameters specified by opts and three additional properties,
90 // specifying that the decimation algorithm should use the max-product
91 // algorithm and should completely reinitalize its state at every step
92 DecMAP decmap(fg, opts("reinit",true)("ianame",string("BP"))("iaopts",string("[damping=0.1,inference=MAXPROD,logdomain=0,maxiter=1000,tol=1e-9,updates=SEQRND,verbose=1]")) );
93 decmap.init();
94 decmap.run();
95 vector<size_t> decmapstate = decmap.findMaximum();
96
97 // Report variable marginals for fg, calculated by the junction tree algorithm
98 cout << "Exact variable marginals:" << endl;
99 for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
100 cout << jt.belief(fg.var(i)) << endl; // display the "belief" of jt for that variable
101
102 // Report variable marginals for fg, calculated by the belief propagation algorithm
103 cout << "Approximate (loopy belief propagation) variable marginals:" << endl;
104 for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
105 cout << bp.belief(fg.var(i)) << endl; // display the belief of bp for that variable
106
107 // Report factor marginals for fg, calculated by the junction tree algorithm
108 cout << "Exact factor marginals:" << endl;
109 for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
110 cout << jt.belief(fg.factor(I).vars()) << endl; // display the "belief" of jt for the variables in that factor
111
112 // Report factor marginals for fg, calculated by the belief propagation algorithm
113 cout << "Approximate (loopy belief propagation) factor marginals:" << endl;
114 for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
115 cout << bp.belief(fg.factor(I).vars()) << endl; // display the belief of bp for the variables in that factor
116
117 // Report log partition sum (normalizing constant) of fg, calculated by the junction tree algorithm
118 cout << "Exact log partition sum: " << jt.logZ() << endl;
119
120 // Report log partition sum of fg, approximated by the belief propagation algorithm
121 cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
122
123 // Report exact MAP variable marginals
124 cout << "Exact MAP variable marginals:" << endl;
125 for( size_t i = 0; i < fg.nrVars(); i++ )
126 cout << jtmap.belief(fg.var(i)) << endl;
127
128 // Report max-product variable marginals
129 cout << "Approximate (max-product) MAP variable marginals:" << endl;
130 for( size_t i = 0; i < fg.nrVars(); i++ )
131 cout << mp.belief(fg.var(i)) << endl;
132
133 // Report exact MAP factor marginals
134 cout << "Exact MAP factor marginals:" << endl;
135 for( size_t I = 0; I < fg.nrFactors(); I++ )
136 cout << jtmap.belief(fg.factor(I).vars()) << "=" << jtmap.beliefF(I) << endl;
137
138 // Report max-product factor marginals
139 cout << "Approximate (max-product) MAP factor marginals:" << endl;
140 for( size_t I = 0; I < fg.nrFactors(); I++ )
141 cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
142
143 // Report exact MAP joint state
144 cout << "Exact MAP state (log score = " << fg.logScore( jtmapstate ) << "):" << endl;
145 for( size_t i = 0; i < jtmapstate.size(); i++ )
146 cout << fg.var(i) << ": " << jtmapstate[i] << endl;
147
148 // Report max-product MAP joint state
149 cout << "Approximate (max-product) MAP state (log score = " << fg.logScore( mpstate ) << "):" << endl;
150 for( size_t i = 0; i < mpstate.size(); i++ )
151 cout << fg.var(i) << ": " << mpstate[i] << endl;
152
153 // Report DecMAP joint state
154 cout << "Approximate DecMAP state (log score = " << fg.logScore( decmapstate ) << "):" << endl;
155 for( size_t i = 0; i < decmapstate.size(); i++ )
156 cout << fg.var(i) << ": " << decmapstate[i] << endl;
157 }
158
159 return 0;
160 }