1 /* Copyright (C) 2008 Frederik Eaton [frederik at ofb dot net],
2 Joris Mooij [joris dot mooij at tuebingen dot mpg dot de]
4 This file is part of libDAI.
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 #include <dai/gibbs.h>
29 #include <dai/properties.h>
38 const char *Gibbs::Name
= "GIBBS";
41 void Gibbs::setProperties( const PropertySet
&opts
) {
42 assert( opts
.hasKey("iters") );
43 props
.iters
= opts
.getStringAs
<size_t>("iters");
45 if( opts
.hasKey("verbose") )
46 props
.verbose
= opts
.getStringAs
<size_t>("verbose");
52 PropertySet
Gibbs::getProperties() const {
54 opts
.Set( "iters", props
.iters
);
55 opts
.Set( "verbose", props
.verbose
);
60 string
Gibbs::printProperties() const {
61 stringstream
s( stringstream::out
);
63 s
<< "iters=" << props
.iters
<< ",";
64 s
<< "verbose=" << props
.verbose
<< "]";
69 void Gibbs::construct() {
71 _var_counts
.reserve( nrVars() );
72 for( size_t i
= 0; i
< nrVars(); i
++ )
73 _var_counts
.push_back( _count_t( var(i
).states(), 0 ) );
75 _factor_counts
.clear();
76 _factor_counts
.reserve( nrFactors() );
77 for( size_t I
= 0; I
< nrFactors(); I
++ )
78 _factor_counts
.push_back( _count_t( factor(I
).states(), 0 ) );
83 _state
.resize( nrVars(), 0 );
87 void Gibbs::updateCounts() {
88 for( size_t i
= 0; i
< nrVars(); i
++ )
89 _var_counts
[i
][_state
[i
]]++;
90 for( size_t I
= 0; I
< nrFactors(); I
++ )
91 _factor_counts
[I
][getFactorEntry(I
)]++;
96 inline size_t Gibbs::getFactorEntry( size_t I
) {
98 for( int _j
= nbF(I
).size() - 1; _j
>= 0; _j
-- ) {
99 // note that iterating over nbF(I) yields the same ordering
100 // of variables as iterating over factor(I).vars()
101 size_t j
= nbF(I
)[_j
];
102 f_entry
*= var(j
).states();
103 f_entry
+= _state
[j
];
109 inline size_t Gibbs::getFactorEntryDiff( size_t I
, size_t i
) {
111 for( size_t _j
= 0; _j
< nbF(I
).size(); _j
++ ) {
112 // note that iterating over nbF(I) yields the same ordering
113 // of variables as iterating over factor(I).vars()
114 size_t j
= nbF(I
)[_j
];
118 skip
*= var(j
).states();
124 Prob
Gibbs::getVarDist( size_t i
) {
125 assert( i
< nrVars() );
126 size_t i_states
= var(i
).states();
127 Prob
i_given_MB( i_states
, 1.0 );
129 // use Markov blanket of var(i) to calculate distribution
130 foreach( const Neighbor
&I
, nbV(i
) ) {
131 const Factor
&f_I
= factor(I
);
132 size_t I_skip
= getFactorEntryDiff( I
, i
);
133 size_t I_entry
= getFactorEntry(I
) - (_state
[i
] * I_skip
);
134 for( size_t st_i
= 0; st_i
< i_states
; st_i
++ ) {
135 i_given_MB
[st_i
] *= f_I
[I_entry
];
140 if( i_given_MB
.sum() == 0.0 )
141 // If no state of i is allowed, use uniform distribution
142 // FIXME is that indeed the right thing to do?
143 i_given_MB
= Prob( i_states
);
145 i_given_MB
.normalize();
150 inline void Gibbs::resampleVar( size_t i
) {
151 // draw randomly from conditional distribution and update _state
152 _state
[i
] = getVarDist(i
).draw();
156 void Gibbs::randomizeState() {
157 for( size_t i
= 0; i
< nrVars(); i
++ )
158 _state
[i
] = rnd_int( 0, var(i
).states() - 1 );
163 for( size_t i
= 0; i
< nrVars(); i
++ )
164 fill( _var_counts
[i
].begin(), _var_counts
[i
].end(), 0 );
165 for( size_t I
= 0; I
< nrFactors(); I
++ )
166 fill( _factor_counts
[I
].begin(), _factor_counts
[I
].end(), 0 );
171 double Gibbs::run() {
172 if( props
.verbose
>= 1 )
173 cerr
<< "Starting " << identify() << "...";
174 if( props
.verbose
>= 3 )
181 for( size_t iter
= 0; iter
< props
.iters
; iter
++ ) {
182 for( size_t j
= 0; j
< nrVars(); j
++ )
187 if( props
.verbose
>= 3 ) {
188 for( size_t i
= 0; i
< nrVars(); i
++ ) {
189 cerr
<< "belief for variable " << var(i
) << ": " << beliefV(i
) << endl
;
190 cerr
<< "counts for variable " << var(i
) << ": " << Prob( _var_counts
[i
].begin(), _var_counts
[i
].end() ) << endl
;
194 if( props
.verbose
>= 3 )
195 cerr
<< Name
<< "::run: ran " << props
.iters
<< " passes (" << toc() - tic
<< " clocks)." << endl
;
201 Factor
Gibbs::beliefV( size_t i
) const {
202 return Factor( var(i
), _var_counts
[i
].begin() ).normalized();
206 Factor
Gibbs::beliefF( size_t I
) const {
207 return Factor( factor(I
).vars(), _factor_counts
[I
].begin() ).normalized();
211 Factor
Gibbs::belief( const Var
&n
) const {
212 return( beliefV( findVar( n
) ) );
216 vector
<Factor
> Gibbs::beliefs() const {
217 vector
<Factor
> result
;
218 for( size_t i
= 0; i
< nrVars(); ++i
)
219 result
.push_back( beliefV(i
) );
220 for( size_t I
= 0; I
< nrFactors(); ++I
)
221 result
.push_back( beliefF(I
) );
226 Factor
Gibbs::belief( const VarSet
&ns
) const {
228 return belief( *(ns
.begin()) );
231 for( I
= 0; I
< nrFactors(); I
++ )
232 if( factor(I
).vars() >> ns
)
234 assert( I
!= nrFactors() );
235 return beliefF(I
).marginal(ns
);
240 } // end of namespace dai