example_sprinkler_gibbs.cpp

/*  This file is part of libDAI - http://www.libdai.org/
 *
 *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
 *
 *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
 */


#include <dai/factorgraph.h>
#include <dai/gibbs.h>
#include <dai/properties.h>
#include <iostream>
#include <fstream>

using namespace std;
using namespace dai;

int main() {
#ifdef DAI_WITH_GIBBS
    // This example program illustrates how to use Gibbs sampling
    // to sample from a joint probability distribution described
    // by a factor graph, using the sprinkler network example discussed at
    // http://www.cs.ubc.ca/~murphyk/Bayes/bnintro.html
    //
    // The file sprinkler.fg has to be generated by first running
    // example_sprinkler

    // Read the factorgraph from the file
    FactorGraph SprinklerNetwork;
    SprinklerNetwork.ReadFromFile( "sprinkler.fg" );
    cout << "Sprinkler network read from sprinkler.fg" << endl;

    // Output some information about the factorgraph
    cout << SprinklerNetwork.nrVars() << " variables" << endl;
    cout << SprinklerNetwork.nrFactors() << " factors" << endl;

    // Prepare a Gibbs sampler
    PropertySet gibbsProps;
    gibbsProps.set("maxiter", size_t(100));   // number of Gibbs sampler iterations
    gibbsProps.set("burnin", size_t(0));
    gibbsProps.set("verbose", size_t(0));
    Gibbs gibbsSampler( SprinklerNetwork, gibbsProps );

    // Open a .tab file for writing
    ofstream outfile;
    outfile.open( "sprinkler.tab" );
    if( !outfile.is_open() )
        throw "Cannot write to file!";

    // Write header (consisting of variable labels)
    for( size_t i = 0; i < SprinklerNetwork.nrVars(); i++ )
        outfile << (i == 0 ? "" : "\t") << SprinklerNetwork.var(i).label();
    outfile << endl << endl;

    // Draw samples from joint distribution using Gibbs sampling
    // and write them to the .tab file
    size_t nrSamples = 1000;
    std::vector<size_t> state;
    for( size_t t = 0; t < nrSamples; t++ ) {
        gibbsSampler.init();
        gibbsSampler.run();
        state = gibbsSampler.state();
        for( size_t i = 0; i < state.size(); i++ )
            outfile << (i == 0 ? "" : "\t") << state[i];
        outfile << endl;
    }
    cout << nrSamples << " samples written to sprinkler.tab" << endl;

    // Close the .tab file
    outfile.close();

    return 0;
#else
    cout << "libDAI was configured without Gibbs (this can be changed in include/dai/dai_config.h)." << endl;
#endif
}