include/dai/gibbs.h

   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 /// \file
  10 /// \brief Defines class Gibbs, which implements Gibbs sampling
  11
  12
  13 #ifndef __defined_libdai_gibbs_h
  14 #define __defined_libdai_gibbs_h
  15
  16
  17 #include <dai/dai_config.h>
  18 #ifdef DAI_WITH_GIBBS
  19
  20
  21 #include <dai/daialg.h>
  22 #include <dai/factorgraph.h>
  23 #include <dai/properties.h>
  24
  25
  26 namespace dai {
  27
  28
  29 /// Approximate inference algorithm "Gibbs sampling"
  30 /** \author Frederik Eaton
  31  */
  32 class Gibbs : public DAIAlgFG {
  33     private:
  34         /// Type used to store the counts of various states
  35         typedef std::vector<size_t> _count_t;
  36         /// Type used to store the joint state of all variables
  37         typedef std::vector<size_t> _state_t;
  38         /// Number of samples counted so far (excluding burn-in periods)
  39         size_t _sample_count;
  40         /// State counts for each variable
  41         std::vector<_count_t> _var_counts;
  42         /// State counts for each factor
  43         std::vector<_count_t> _factor_counts;
  44         /// Number of iterations done (including burn-in periods)
  45         size_t _iters;
  46         /// Current joint state of all variables
  47         _state_t _state;
  48         /// Joint state with maximum probability seen so far
  49         _state_t _max_state;
  50         /// Highest score so far
  51         Real _max_score;
  52
  53     public:
  54         /// Parameters for Gibbs
  55         struct Properties {
  56             /// Maximum number of iterations
  57             size_t maxiter;
  58
  59             /// Maximum time (in seconds)
  60             double maxtime;
  61
  62             /// Number of iterations after which a random restart is made
  63             size_t restart;
  64
  65             /// Number of "burn-in" iterations after each (re)start (for which no statistics are gathered)
  66             size_t burnin;
  67
  68             /// Verbosity (amount of output sent to stderr)
  69             size_t verbose;
  70         } props;
  71
  72     public:
  73         /// Default constructor
  74         Gibbs() : DAIAlgFG(), _sample_count(0), _var_counts(), _factor_counts(), _iters(0), _state(), _max_state(), _max_score(-INFINITY) {}
  75
  76         /// Construct from FactorGraph \a fg and PropertySet \a opts
  77         /** \param fg Factor graph.
  78          *  \param opts Parameters @see Properties
  79          */
  80         Gibbs( const FactorGraph &fg, const PropertySet &opts ) : DAIAlgFG(fg), _sample_count(0), _var_counts(), _factor_counts(), _iters(0), _state(), _max_state(), _max_score(-INFINITY) {
  81             setProperties( opts );
  82             construct();
  83         }
  84
  85
  86     /// \name General InfAlg interface
  87     //@{
  88         virtual Gibbs* clone() const { return new Gibbs(*this); }
  89         virtual Gibbs* construct( const FactorGraph &fg, const PropertySet &opts ) const { return new Gibbs( fg, opts ); }
  90         virtual std::string name() const { return "GIBBS"; }
  91         virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
  92         virtual Factor belief( const VarSet &vs ) const;
  93         virtual Factor beliefV( size_t i ) const;
  94         virtual Factor beliefF( size_t I ) const;
  95         virtual std::vector<Factor> beliefs() const;
  96         virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
  97         std::vector<size_t> findMaximum() const { return _max_state; }
  98         virtual void init();
  99         virtual void init( const VarSet &/*ns*/ ) { init(); }
 100         virtual Real run();
 101         virtual Real maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
 102         virtual size_t Iterations() const { return _iters; }
 103         virtual void setMaxIter( size_t maxiter ) { props.maxiter = maxiter; }
 104         virtual void setProperties( const PropertySet &opts );
 105         virtual PropertySet getProperties() const;
 106         virtual std::string printProperties() const;
 107     //@}
 108
 109
 110     /// \name Additional interface specific for Gibbs
 111     //@{
 112         /// Draw the current joint state of all variables from a uniform random distribution
 113         void randomizeState();
 114         /// Return reference to current state of all variables
 115         std::vector<size_t>& state() { return _state; }
 116         /// Return constant reference to current state of all variables
 117         const std::vector<size_t>& state() const { return _state; }
 118     //@}
 119
 120     private:
 121         /// Helper function for constructors
 122         void construct();
 123         /// Updates all counts (_sample_count, _var_counts, _factor_counts) based on current state
 124         void updateCounts();
 125         /// Calculate conditional distribution of variable \a i, given the current state
 126         Prob getVarDist( size_t i );
 127         /// Draw state of variable \a i randomly from its conditional distribution and update the current state
 128         void resampleVar( size_t i );
 129         /// Calculates linear index into factor \a I corresponding to the current state
 130         size_t getFactorEntry( size_t I );
 131         /// Calculates the differences between linear indices into factor \a I corresponding with a state change of variable \a i
 132         size_t getFactorEntryDiff( size_t I, size_t i );
 133 };
 134
 135
 136 /// Runs Gibbs sampling for \a maxiter iterations (of which \a burnin for burn-in) on FactorGraph \a fg, and returns the resulting state
 137 /** \relates Gibbs
 138  */
 139 std::vector<size_t> getGibbsState( const FactorGraph &fg, size_t maxiter );
 140
 141
 142 } // end of namespace dai
 143
 144
 145 /** \example example_sprinkler_gibbs.cpp
 146  *  This example shows how to use the Gibbs class.
 147  */
 148
 149
 150 #endif
 151
 152
 153 #endif