Multiple changes: changes in build system, one workaround and one bug fix
[libdai.git] / 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