Improved documentation of include/dai/index.h
[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 <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 double 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 }