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 }