New git HEAD version
[libdai.git] / examples / example_sprinkler_gibbs.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <dai/factorgraph.h>
10 #include <dai/gibbs.h>
11 #include <dai/properties.h>
12 #include <iostream>
13 #include <fstream>
14
15 using namespace std;
16 using namespace dai;
17
18 int main() {
19 #ifdef DAI_WITH_GIBBS
20 // This example program illustrates how to use Gibbs sampling
21 // to sample from a joint probability distribution described
22 // by a factor graph, using the sprinkler network example discussed at
23 // http://www.cs.ubc.ca/~murphyk/Bayes/bnintro.html
24 //
25 // The file sprinkler.fg has to be generated by first running
26 // example_sprinkler
27
28 // Read the factorgraph from the file
29 FactorGraph SprinklerNetwork;
30 SprinklerNetwork.ReadFromFile( "sprinkler.fg" );
31 cout << "Sprinkler network read from sprinkler.fg" << endl;
32
33 // Output some information about the factorgraph
34 cout << SprinklerNetwork.nrVars() << " variables" << endl;
35 cout << SprinklerNetwork.nrFactors() << " factors" << endl;
36
37 // Prepare a Gibbs sampler
38 PropertySet gibbsProps;
39 gibbsProps.set("maxiter", size_t(100)); // number of Gibbs sampler iterations
40 gibbsProps.set("burnin", size_t(0));
41 gibbsProps.set("verbose", size_t(0));
42 Gibbs gibbsSampler( SprinklerNetwork, gibbsProps );
43
44 // Open a .tab file for writing
45 ofstream outfile;
46 outfile.open( "sprinkler.tab" );
47 if( !outfile.is_open() )
48 throw "Cannot write to file!";
49
50 // Write header (consisting of variable labels)
51 for( size_t i = 0; i < SprinklerNetwork.nrVars(); i++ )
52 outfile << (i == 0 ? "" : "\t") << SprinklerNetwork.var(i).label();
53 outfile << endl << endl;
54
55 // Draw samples from joint distribution using Gibbs sampling
56 // and write them to the .tab file
57 size_t nrSamples = 1000;
58 std::vector<size_t> state;
59 for( size_t t = 0; t < nrSamples; t++ ) {
60 gibbsSampler.init();
61 gibbsSampler.run();
62 state = gibbsSampler.state();
63 for( size_t i = 0; i < state.size(); i++ )
64 outfile << (i == 0 ? "" : "\t") << state[i];
65 outfile << endl;
66 }
67 cout << nrSamples << " samples written to sprinkler.tab" << endl;
68
69 // Close the .tab file
70 outfile.close();
71
72 return 0;
73 #else
74 cout << "libDAI was configured without Gibbs (this can be changed in include/dai/dai_config.h)." << endl;
75 #endif
76 }