1 /* This file is part of libDAI - http://www.libdai.org/
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
9 #include <dai/dai_config.h>
18 #include <dai/gibbs.h>
20 #include <dai/properties.h>
29 void Gibbs::setProperties( const PropertySet
&opts
) {
30 DAI_ASSERT( opts
.hasKey("maxiter") );
31 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
33 if( opts
.hasKey("restart") )
34 props
.restart
= opts
.getStringAs
<size_t>("restart");
36 props
.restart
= props
.maxiter
;
37 if( opts
.hasKey("burnin") )
38 props
.burnin
= opts
.getStringAs
<size_t>("burnin");
41 if( opts
.hasKey("maxtime") )
42 props
.maxtime
= opts
.getStringAs
<Real
>("maxtime");
44 props
.maxtime
= INFINITY
;
45 if( opts
.hasKey("verbose") )
46 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
52 PropertySet
Gibbs::getProperties() const {
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
);
63 string
Gibbs::printProperties() const {
64 stringstream
s( stringstream::out
);
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
<< "]";
75 void Gibbs::construct() {
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 ) );
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 ) );
91 _state
.resize( nrVars(), 0 );
94 _max_state
.resize( nrVars(), 0 );
96 _max_score
= logScore( _max_state
);
100 void Gibbs::updateCounts() {
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
) {
114 size_t Gibbs::getFactorEntry( size_t I
) {
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
];
127 size_t Gibbs::getFactorEntryDiff( size_t I
, size_t i
) {
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
];
136 skip
*= var(j
).states();
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 );
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
] );
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
);
163 i_given_MB
.normalize();
168 void Gibbs::resampleVar( size_t i
) {
169 _state
[i
] = getVarDist(i
).draw();
173 void Gibbs::randomizeState() {
174 for( size_t i
= 0; i
< nrVars(); i
++ )
175 _state
[i
] = rnd( var(i
).states() );
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 );
190 if( props
.verbose
>= 1 )
191 cerr
<< "Starting " << identify() << "...";
192 if( props
.verbose
>= 3 )
197 for( ; _iters
< props
.maxiter
&& (toc() - tic
) < props
.maxtime
; _iters
++ ) {
198 if( (_iters
% props
.restart
) == 0 )
200 for( size_t i
= 0; i
< nrVars(); i
++ )
202 if( (_iters
% props
.restart
) > props
.burnin
)
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
;
213 if( props
.verbose
>= 3 )
214 cerr
<< name() << "::run: ran " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
219 return std::pow( _iters
, -0.5 );
223 Factor
Gibbs::beliefV( size_t i
) const {
224 if( _sample_count
== 0 )
225 return Factor( var(i
) );
227 return Factor( var(i
), _var_counts
[i
] ).normalized();
231 Factor
Gibbs::beliefF( size_t I
) const {
232 if( _sample_count
== 0 )
233 return Factor( factor(I
).vars() );
235 return Factor( factor(I
).vars(), _factor_counts
[I
] ).normalized();
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
) );
249 Factor
Gibbs::belief( const VarSet
&ns
) const {
252 else if( ns
.size() == 1 )
253 return beliefV( findVar( *(ns
.begin()) ) );
256 for( I
= 0; I
< nrFactors(); I
++ )
257 if( factor(I
).vars() >> ns
)
259 if( I
== nrFactors() )
260 DAI_THROW(BELIEF_NOT_AVAILABLE
);
261 return beliefF(I
).marginal(ns
);
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
);
273 return gibbs
.state();
277 } // end of namespace dai