3dca7e7fe4d8b8be49f21cfcf0e129f9352b122c
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.
14 #include <dai/gibbs.h>
16 #include <dai/properties.h>
25 void Gibbs::setProperties( const PropertySet
&opts
) {
26 DAI_ASSERT( opts
.hasKey("maxiter") );
27 props
.maxiter
= opts
.getStringAs
<size_t>("maxiter");
29 if( opts
.hasKey("restart") )
30 props
.restart
= opts
.getStringAs
<size_t>("restart");
32 props
.restart
= props
.maxiter
;
33 if( opts
.hasKey("burnin") )
34 props
.burnin
= opts
.getStringAs
<size_t>("burnin");
37 if( opts
.hasKey("maxtime") )
38 props
.maxtime
= opts
.getStringAs
<Real
>("maxtime");
40 props
.maxtime
= INFINITY
;
41 if( opts
.hasKey("verbose") )
42 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
48 PropertySet
Gibbs::getProperties() const {
50 opts
.set( "maxiter", props
.maxiter
);
51 opts
.set( "maxtime", props
.maxtime
);
52 opts
.set( "restart", props
.restart
);
53 opts
.set( "burnin", props
.burnin
);
54 opts
.set( "verbose", props
.verbose
);
59 string
Gibbs::printProperties() const {
60 stringstream
s( stringstream::out
);
62 s
<< "maxiter=" << props
.maxiter
<< ",";
63 s
<< "maxtime=" << props
.maxtime
<< ",";
64 s
<< "restart=" << props
.restart
<< ",";
65 s
<< "burnin=" << props
.burnin
<< ",";
66 s
<< "verbose=" << props
.verbose
<< "]";
71 void Gibbs::construct() {
75 _var_counts
.reserve( nrVars() );
76 for( size_t i
= 0; i
< nrVars(); i
++ )
77 _var_counts
.push_back( _count_t( var(i
).states(), 0 ) );
79 _factor_counts
.clear();
80 _factor_counts
.reserve( nrFactors() );
81 for( size_t I
= 0; I
< nrFactors(); I
++ )
82 _factor_counts
.push_back( _count_t( factor(I
).nrStates(), 0 ) );
87 _state
.resize( nrVars(), 0 );
90 _max_state
.resize( nrVars(), 0 );
92 _max_score
= logScore( _max_state
);
96 void Gibbs::updateCounts() {
98 for( size_t i
= 0; i
< nrVars(); i
++ )
99 _var_counts
[i
][_state
[i
]]++;
100 for( size_t I
= 0; I
< nrFactors(); I
++ )
101 _factor_counts
[I
][getFactorEntry(I
)]++;
102 Real score
= logScore( _state
);
103 if( score
> _max_score
) {
110 size_t Gibbs::getFactorEntry( size_t I
) {
112 for( int _j
= nbF(I
).size() - 1; _j
>= 0; _j
-- ) {
113 // note that iterating over nbF(I) yields the same ordering
114 // of variables as iterating over factor(I).vars()
115 size_t j
= nbF(I
)[_j
];
116 f_entry
*= var(j
).states();
117 f_entry
+= _state
[j
];
123 size_t Gibbs::getFactorEntryDiff( size_t I
, size_t i
) {
125 for( size_t _j
= 0; _j
< nbF(I
).size(); _j
++ ) {
126 // note that iterating over nbF(I) yields the same ordering
127 // of variables as iterating over factor(I).vars()
128 size_t j
= nbF(I
)[_j
];
132 skip
*= var(j
).states();
138 Prob
Gibbs::getVarDist( size_t i
) {
139 DAI_ASSERT( i
< nrVars() );
140 size_t i_states
= var(i
).states();
141 Prob
i_given_MB( i_states
, 1.0 );
143 // use Markov blanket of var(i) to calculate distribution
144 bforeach( const Neighbor
&I
, nbV(i
) ) {
145 const Factor
&f_I
= factor(I
);
146 size_t I_skip
= getFactorEntryDiff( I
, i
);
147 size_t I_entry
= getFactorEntry(I
) - (_state
[i
] * I_skip
);
148 for( size_t st_i
= 0; st_i
< i_states
; st_i
++ ) {
149 i_given_MB
.set( st_i
, i_given_MB
[st_i
] * f_I
[I_entry
] );
154 if( i_given_MB
.sum() == 0.0 )
155 // If no state of i is allowed, use uniform distribution
156 // FIXME is that indeed the right thing to do?
157 i_given_MB
= Prob( i_states
);
159 i_given_MB
.normalize();
164 void Gibbs::resampleVar( size_t i
) {
165 _state
[i
] = getVarDist(i
).draw();
169 void Gibbs::randomizeState() {
170 for( size_t i
= 0; i
< nrVars(); i
++ )
171 _state
[i
] = rnd( var(i
).states() );
177 for( size_t i
= 0; i
< nrVars(); i
++ )
178 fill( _var_counts
[i
].begin(), _var_counts
[i
].end(), 0 );
179 for( size_t I
= 0; I
< nrFactors(); I
++ )
180 fill( _factor_counts
[I
].begin(), _factor_counts
[I
].end(), 0 );
186 if( props
.verbose
>= 1 )
187 cerr
<< "Starting " << identify() << "...";
188 if( props
.verbose
>= 3 )
193 for( ; _iters
< props
.maxiter
&& (toc() - tic
) < props
.maxtime
; _iters
++ ) {
194 if( (_iters
% props
.restart
) == 0 )
196 for( size_t i
= 0; i
< nrVars(); i
++ )
198 if( (_iters
% props
.restart
) > props
.burnin
)
202 if( props
.verbose
>= 3 ) {
203 for( size_t i
= 0; i
< nrVars(); i
++ ) {
204 cerr
<< "Belief for variable " << var(i
) << ": " << beliefV(i
) << endl
;
205 cerr
<< "Counts for variable " << var(i
) << ": " << Prob( _var_counts
[i
] ) << endl
;
209 if( props
.verbose
>= 3 )
210 cerr
<< name() << "::run: ran " << _iters
<< " passes (" << toc() - tic
<< " seconds)." << endl
;
215 return std::pow( _iters
, -0.5 );
219 Factor
Gibbs::beliefV( size_t i
) const {
220 if( _sample_count
== 0 )
221 return Factor( var(i
) );
223 return Factor( var(i
), _var_counts
[i
] ).normalized();
227 Factor
Gibbs::beliefF( size_t I
) const {
228 if( _sample_count
== 0 )
229 return Factor( factor(I
).vars() );
231 return Factor( factor(I
).vars(), _factor_counts
[I
] ).normalized();
235 vector
<Factor
> Gibbs::beliefs() const {
236 vector
<Factor
> result
;
237 for( size_t i
= 0; i
< nrVars(); ++i
)
238 result
.push_back( beliefV(i
) );
239 for( size_t I
= 0; I
< nrFactors(); ++I
)
240 result
.push_back( beliefF(I
) );
245 Factor
Gibbs::belief( const VarSet
&ns
) const {
248 else if( ns
.size() == 1 )
249 return beliefV( findVar( *(ns
.begin()) ) );
252 for( I
= 0; I
< nrFactors(); I
++ )
253 if( factor(I
).vars() >> ns
)
255 if( I
== nrFactors() )
256 DAI_THROW(BELIEF_NOT_AVAILABLE
);
257 return beliefF(I
).marginal(ns
);
262 std::vector
<size_t> getGibbsState( const FactorGraph
&fg
, size_t maxiter
) {
263 PropertySet gibbsProps
;
264 gibbsProps
.set( "maxiter", maxiter
);
265 gibbsProps
.set( "burnin", size_t(0) );
266 gibbsProps
.set( "verbose", size_t(0) );
267 Gibbs
gibbs( fg
, gibbsProps
);
269 return gibbs
.state();
273 } // end of namespace dai