src/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/dai_config.h>
  10 #ifdef DAI_WITH_GIBBS
  11
  12
  13 #include <iostream>
  14 #include <sstream>
  15 #include <map>
  16 #include <set>
  17 #include <algorithm>
  18 #include <dai/gibbs.h>
  19 #include <dai/util.h>
  20 #include <dai/properties.h>
  21
  22
  23 namespace dai {
  24
  25
  26 using namespace std;
  27
  28
  29 void Gibbs::setProperties( const PropertySet &opts ) {
  30     DAI_ASSERT( opts.hasKey("maxiter") );
  31     props.maxiter = opts.getStringAs<size_t>("maxiter");
  32
  33     if( opts.hasKey("restart") )
  34         props.restart = opts.getStringAs<size_t>("restart");
  35     else
  36         props.restart = props.maxiter;
  37     if( opts.hasKey("burnin") )
  38         props.burnin = opts.getStringAs<size_t>("burnin");
  39     else
  40         props.burnin = 0;
  41     if( opts.hasKey("maxtime") )
  42         props.maxtime = opts.getStringAs<Real>("maxtime");
  43     else
  44         props.maxtime = INFINITY;
  45     if( opts.hasKey("verbose") )
  46         props.verbose = opts.getStringAs<size_t>("verbose");
  47     else
  48         props.verbose = 0;
  49 }
  50
  51
  52 PropertySet Gibbs::getProperties() const {
  53     PropertySet opts;
  54     opts.set( "maxiter", props.maxiter );
  55     opts.set( "maxtime", props.maxtime );
  56     opts.set( "restart", props.restart );
  57     opts.set( "burnin", props.burnin );
  58     opts.set( "verbose", props.verbose );
  59     return opts;
  60 }
  61
  62
  63 string Gibbs::printProperties() const {
  64     stringstream s( stringstream::out );
  65     s << "[";
  66     s << "maxiter=" << props.maxiter << ",";
  67     s << "maxtime=" << props.maxtime << ",";
  68     s << "restart=" << props.restart << ",";
  69     s << "burnin=" << props.burnin << ",";
  70     s << "verbose=" << props.verbose << "]";
  71     return s.str();
  72 }
  73
  74
  75 void Gibbs::construct() {
  76     _sample_count = 0;
  77
  78     _var_counts.clear();
  79     _var_counts.reserve( nrVars() );
  80     for( size_t i = 0; i < nrVars(); i++ )
  81         _var_counts.push_back( _count_t( var(i).states(), 0 ) );
  82
  83     _factor_counts.clear();
  84     _factor_counts.reserve( nrFactors() );
  85     for( size_t I = 0; I < nrFactors(); I++ )
  86         _factor_counts.push_back( _count_t( factor(I).nrStates(), 0 ) );
  87
  88     _iters = 0;
  89
  90     _state.clear();
  91     _state.resize( nrVars(), 0 );
  92
  93     _max_state.clear();
  94     _max_state.resize( nrVars(), 0 );
  95
  96     _max_score = logScore( _max_state );
  97 }
  98
  99
 100 void Gibbs::updateCounts() {
 101     _sample_count++;
 102     for( size_t i = 0; i < nrVars(); i++ )
 103         _var_counts[i][_state[i]]++;
 104     for( size_t I = 0; I < nrFactors(); I++ )
 105         _factor_counts[I][getFactorEntry(I)]++;
 106     Real score = logScore( _state );
 107     if( score > _max_score ) {
 108         _max_state = _state;
 109         _max_score = score;
 110     }
 111 }
 112
 113
 114 size_t Gibbs::getFactorEntry( size_t I ) {
 115     size_t f_entry = 0;
 116     for( int _j = nbF(I).size() - 1; _j >= 0; _j-- ) {
 117         // note that iterating over nbF(I) yields the same ordering
 118         // of variables as iterating over factor(I).vars()
 119         size_t j = nbF(I)[_j];
 120         f_entry *= var(j).states();
 121         f_entry += _state[j];
 122     }
 123     return f_entry;
 124 }
 125
 126
 127 size_t Gibbs::getFactorEntryDiff( size_t I, size_t i ) {
 128     size_t skip = 1;
 129     for( size_t _j = 0; _j < nbF(I).size(); _j++ ) {
 130         // note that iterating over nbF(I) yields the same ordering
 131         // of variables as iterating over factor(I).vars()
 132         size_t j = nbF(I)[_j];
 133         if( i == j )
 134             break;
 135         else
 136             skip *= var(j).states();
 137     }
 138     return skip;
 139 }
 140
 141
 142 Prob Gibbs::getVarDist( size_t i ) {
 143     DAI_ASSERT( i < nrVars() );
 144     size_t i_states = var(i).states();
 145     Prob i_given_MB( i_states, 1.0 );
 146
 147     // use Markov blanket of var(i) to calculate distribution
 148     bforeach( const Neighbor &I, nbV(i) ) {
 149         const Factor &f_I = factor(I);
 150         size_t I_skip = getFactorEntryDiff( I, i );
 151         size_t I_entry = getFactorEntry(I) - (_state[i] * I_skip);
 152         for( size_t st_i = 0; st_i < i_states; st_i++ ) {
 153             i_given_MB.set( st_i, i_given_MB[st_i] * f_I[I_entry] );
 154             I_entry += I_skip;
 155         }
 156     }
 157
 158     if( i_given_MB.sum() == 0.0 )
 159         // If no state of i is allowed, use uniform distribution
 160         // FIXME is that indeed the right thing to do?
 161         i_given_MB = Prob( i_states );
 162     else
 163         i_given_MB.normalize();
 164     return i_given_MB;
 165 }
 166
 167
 168 void Gibbs::resampleVar( size_t i ) {
 169     _state[i] = getVarDist(i).draw();
 170 }
 171
 172
 173 void Gibbs::randomizeState() {
 174     for( size_t i = 0; i < nrVars(); i++ )
 175         _state[i] = rnd( var(i).states() );
 176 }
 177
 178
 179 void Gibbs::init() {
 180     _sample_count = 0;
 181     for( size_t i = 0; i < nrVars(); i++ )
 182         fill( _var_counts[i].begin(), _var_counts[i].end(), 0 );
 183     for( size_t I = 0; I < nrFactors(); I++ )
 184         fill( _factor_counts[I].begin(), _factor_counts[I].end(), 0 );
 185     _iters = 0;
 186 }
 187
 188
 189 Real Gibbs::run() {
 190     if( props.verbose >= 1 )
 191         cerr << "Starting " << identify() << "...";
 192     if( props.verbose >= 3 )
 193         cerr << endl;
 194
 195     double tic = toc();
 196
 197     for( ; _iters < props.maxiter && (toc() - tic) < props.maxtime; _iters++ ) {
 198         if( (_iters % props.restart) == 0 )
 199             randomizeState();
 200         for( size_t i = 0; i < nrVars(); i++ )
 201             resampleVar( i );
 202         if( (_iters % props.restart) > props.burnin )
 203             updateCounts();
 204     }
 205
 206     if( props.verbose >= 3 ) {
 207         for( size_t i = 0; i < nrVars(); i++ ) {
 208             cerr << "Belief for variable " << var(i) << ": " << beliefV(i) << endl;
 209             cerr << "Counts for variable " << var(i) << ": " << Prob( _var_counts[i] ) << endl;
 210         }
 211     }
 212
 213     if( props.verbose >= 3 )
 214         cerr << name() << "::run:  ran " << _iters << " passes (" << toc() - tic << " seconds)." << endl;
 215
 216     if( _iters == 0 )
 217         return INFINITY;
 218     else
 219         return std::pow( _iters, -0.5 );
 220 }
 221
 222
 223 Factor Gibbs::beliefV( size_t i ) const {
 224     if( _sample_count == 0 )
 225         return Factor( var(i) );
 226     else
 227         return Factor( var(i), _var_counts[i] ).normalized();
 228 }
 229
 230
 231 Factor Gibbs::beliefF( size_t I ) const {
 232     if( _sample_count == 0 )
 233         return Factor( factor(I).vars() );
 234     else
 235         return Factor( factor(I).vars(), _factor_counts[I] ).normalized();
 236 }
 237
 238
 239 vector<Factor> Gibbs::beliefs() const {
 240     vector<Factor> result;
 241     for( size_t i = 0; i < nrVars(); ++i )
 242         result.push_back( beliefV(i) );
 243     for( size_t I = 0; I < nrFactors(); ++I )
 244         result.push_back( beliefF(I) );
 245     return result;
 246 }
 247
 248
 249 Factor Gibbs::belief( const VarSet &ns ) const {
 250     if( ns.size() == 0 )
 251         return Factor();
 252     else if( ns.size() == 1 )
 253         return beliefV( findVar( *(ns.begin()) ) );
 254     else {
 255         size_t I;
 256         for( I = 0; I < nrFactors(); I++ )
 257             if( factor(I).vars() >> ns )
 258                 break;
 259         if( I == nrFactors() )
 260             DAI_THROW(BELIEF_NOT_AVAILABLE);
 261         return beliefF(I).marginal(ns);
 262     }
 263 }
 264
 265
 266 std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t maxiter ) {
 267     PropertySet gibbsProps;
 268     gibbsProps.set( "maxiter", maxiter );
 269     gibbsProps.set( "burnin", size_t(0) );
 270     gibbsProps.set( "verbose", size_t(0) );
 271     Gibbs gibbs( fg, gibbsProps );
 272     gibbs.run();
 273     return gibbs.state();
 274 }
 275
 276
 277 } // end of namespace dai
 278
 279
 280 #endif