Working on documentation of include/dai/emalg.h
[libdai.git] / src / emalg.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
6 *
7 * Copyright (C) 2009 Charles Vaske [cvaske at soe dot ucsc dot edu]
8 * Copyright (C) 2009 University of California, Santa Cruz
9 */
10
11
12 #include <dai/util.h>
13 #include <dai/emalg.h>
14
15
16 namespace dai {
17
18
19 // Initialize static private member of ParameterEstimation
20 std::map<std::string, ParameterEstimation::ParamEstFactory> *ParameterEstimation::_registry = NULL;
21
22
23 void ParameterEstimation::loadDefaultRegistry() {
24 _registry = new std::map<std::string, ParamEstFactory>();
25 (*_registry)["ConditionalProbEstimation"] = CondProbEstimation::factory;
26 }
27
28
29 ParameterEstimation* ParameterEstimation::construct( const std::string &method, const PropertySet &p ) {
30 if( _registry == NULL )
31 loadDefaultRegistry();
32 std::map<std::string, ParamEstFactory>::iterator i = _registry->find(method);
33 if( i == _registry->end() )
34 DAI_THROWE(UNKNOWN_PARAMETER_ESTIMATION_METHOD, "Unknown parameter estimation method: " + method);
35 ParamEstFactory factory = i->second;
36 return factory(p);
37 }
38
39
40 ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
41 size_t target_dimension = p.getStringAs<size_t>("target_dim");
42 size_t total_dimension = p.getStringAs<size_t>("total_dim");
43 Real pseudo_count = 1;
44 if( p.hasKey("pseudo_count") )
45 pseudo_count = p.getStringAs<Real>("pseudo_count");
46 return new CondProbEstimation( target_dimension, Prob( total_dimension, pseudo_count ) );
47 }
48
49
50 CondProbEstimation::CondProbEstimation( size_t target_dimension, const Prob &pseudocounts )
51 : _target_dim(target_dimension), _stats(pseudocounts), _initial_stats(pseudocounts)
52 {
53 DAI_ASSERT( !(_stats.size() % _target_dim) );
54 }
55
56
57 void CondProbEstimation::addSufficientStatistics( const Prob &p ) {
58 _stats += p;
59 }
60
61
62 Prob CondProbEstimation::estimate() {
63 // normalize pseudocounts
64 for( size_t parent = 0; parent < _stats.size(); parent += _target_dim ) {
65 // calculate norm
66 size_t top = parent + _target_dim;
67 Real norm = std::accumulate( &(_stats[parent]), &(_stats[top]), 0.0 );
68 if( norm != 0.0 )
69 norm = 1.0 / norm;
70 // normalize
71 for( size_t i = parent; i < top; ++i )
72 _stats[i] *= norm;
73 }
74 // reset _stats to _initial_stats
75 Prob result = _stats;
76 _stats = _initial_stats;
77 return result;
78 }
79
80
81 Permute SharedParameters::calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS ) {
82 // Collect all labels and dimensions, and order them in vs
83 std::vector<size_t> dims;
84 dims.reserve( varorder.size() );
85 std::vector<long> labels;
86 labels.reserve( varorder.size() );
87 for( size_t i = 0; i < varorder.size(); i++ ) {
88 dims.push_back( varorder[i].states() );
89 labels.push_back( varorder[i].label() );
90 outVS |= varorder[i];
91 }
92
93 // Construct the sigma array for the permutation object
94 std::vector<size_t> sigma;
95 sigma.reserve( dims.size() );
96 for( VarSet::iterator set_iterator = outVS.begin(); sigma.size() < dims.size(); ++set_iterator )
97 sigma.push_back( find(labels.begin(), labels.end(), set_iterator->label()) - labels.begin() );
98
99 return Permute( dims, sigma );
100 }
101
102
103 void SharedParameters::setPermsAndVarSetsFromVarOrders() {
104 if( _varorders.size() == 0 )
105 return;
106 DAI_ASSERT( _estimation != NULL );
107
108 // Construct the permutation objects and the varsets
109 for( FactorOrientations::const_iterator foi = _varorders.begin(); foi != _varorders.end(); ++foi ) {
110 VarSet vs;
111 _perms[foi->first] = calculatePermutation( foi->second, vs );
112 _varsets[foi->first] = vs;
113 DAI_ASSERT( _estimation->probSize() == vs.nrStates() );
114 }
115 }
116
117
118 SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool ownPE )
119 : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE)
120 {
121 // Calculate the necessary permutations and varsets
122 setPermsAndVarSetsFromVarOrders();
123 }
124
125
126 SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
127 : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true)
128 {
129 // Read the desired parameter estimation method from the stream
130 std::string est_method;
131 PropertySet props;
132 is >> est_method;
133 is >> props;
134
135 // Construct a corresponding object
136 _estimation = ParameterEstimation::construct( est_method, props );
137
138 // Read in the factors that are to be estimated
139 size_t num_factors;
140 is >> num_factors;
141 for( size_t sp_i = 0; sp_i < num_factors; ++sp_i ) {
142 std::string line;
143 while( line.size() == 0 && getline(is, line) )
144 ;
145
146 std::vector<std::string> fields;
147 tokenizeString(line, fields, " \t");
148
149 // Lookup the factor in the factorgraph
150 if( fields.size() < 1 )
151 DAI_THROW(INVALID_EMALG_FILE);
152 std::istringstream iss;
153 iss.str( fields[0] );
154 size_t factor;
155 iss >> factor;
156 const VarSet &vs = fg.factor(factor).vars();
157 if( fields.size() != vs.size() + 1 )
158 DAI_THROW(INVALID_EMALG_FILE);
159
160 // Construct the vector of Vars
161 std::vector<Var> var_order;
162 var_order.reserve( vs.size() );
163 for( size_t fi = 1; fi < fields.size(); ++fi ) {
164 // Lookup a single variable by label
165 long label;
166 std::istringstream labelparse( fields[fi] );
167 labelparse >> label;
168 VarSet::const_iterator vsi = vs.begin();
169 for( ; vsi != vs.end(); ++vsi )
170 if( vsi->label() == label )
171 break;
172 if( vsi == vs.end() )
173 DAI_THROW(INVALID_EMALG_FILE);
174 var_order.push_back( *vsi );
175 }
176 _varorders[factor] = var_order;
177 }
178
179 // Calculate the necessary permutations
180 setPermsAndVarSetsFromVarOrders();
181 }
182
183
184 void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
185 for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
186 Permute &perm = i->second;
187 VarSet &vs = _varsets[i->first];
188
189 Factor b = alg.belief(vs);
190 Prob p( b.states(), 0.0 );
191 for( size_t entry = 0; entry < b.states(); ++entry )
192 p[entry] = b[perm.convertLinearIndex(entry)];
193 _estimation->addSufficientStatistics( p );
194 }
195 }
196
197
198 void SharedParameters::setParameters( FactorGraph &fg ) {
199 Prob p = _estimation->estimate();
200 for( std::map<FactorIndex, Permute>::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
201 Permute &perm = i->second;
202 VarSet &vs = _varsets[i->first];
203
204 Factor f( vs, 0.0 );
205 for( size_t entry = 0; entry < f.states(); ++entry )
206 f[perm.convertLinearIndex(entry)] = p[entry];
207
208 fg.setFactor( i->first, f );
209 }
210 }
211
212
213 void SharedParameters::collectParameters( const FactorGraph &fg, std::vector<Real> &outVals, std::vector<Var> &outVarOrder ) {
214 FactorOrientations::iterator it = _varorders.begin();
215 if( it == _varorders.end() )
216 return;
217 FactorIndex I = it->first;
218 for( std::vector<Var>::const_iterator var_it = _varorders[I].begin(); var_it != _varorders[I].end(); ++var_it )
219 outVarOrder.push_back( *var_it );
220
221 const Factor &f = fg.factor(I);
222 DAI_ASSERT( f.vars() == _varsets[I] );
223 const Permute &perm = _perms[I];
224 for( size_t val_index = 0; val_index < f.states(); ++val_index )
225 outVals.push_back( f[perm.convertLinearIndex(val_index)] );
226 }
227
228
229 MaximizationStep::MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup ) : _params() {
230 size_t num_params = -1;
231 is >> num_params;
232 _params.reserve( num_params );
233 for( size_t i = 0; i < num_params; ++i )
234 _params.push_back( SharedParameters( is, fg_varlookup ) );
235 }
236
237
238 void MaximizationStep::addExpectations( InfAlg &alg ) {
239 for( size_t i = 0; i < _params.size(); ++i )
240 _params[i].collectSufficientStatistics( alg );
241 }
242
243
244 void MaximizationStep::maximize( FactorGraph &fg ) {
245 for( size_t i = 0; i < _params.size(); ++i )
246 _params[i].setParameters( fg );
247 }
248
249
250 const std::string EMAlg::MAX_ITERS_KEY("max_iters");
251 const std::string EMAlg::LOG_Z_TOL_KEY("log_z_tol");
252 const size_t EMAlg::MAX_ITERS_DEFAULT = 30;
253 const Real EMAlg::LOG_Z_TOL_DEFAULT = 0.01;
254
255
256 EMAlg::EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &msteps_file )
257 : _evidence(evidence), _estep(estep), _msteps(), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
258 {
259 msteps_file.exceptions( std::istream::eofbit | std::istream::failbit | std::istream::badbit );
260 size_t num_msteps = -1;
261 msteps_file >> num_msteps;
262 _msteps.reserve(num_msteps);
263 for( size_t i = 0; i < num_msteps; ++i )
264 _msteps.push_back( MaximizationStep( msteps_file, estep.fg() ) );
265 }
266
267
268 void EMAlg::setTermConditions( const PropertySet &p ) {
269 if( p.hasKey(MAX_ITERS_KEY) )
270 _max_iters = p.getStringAs<size_t>(MAX_ITERS_KEY);
271 if( p.hasKey(LOG_Z_TOL_KEY) )
272 _log_z_tol = p.getStringAs<Real>(LOG_Z_TOL_KEY);
273 }
274
275
276 bool EMAlg::hasSatisfiedTermConditions() const {
277 if( _iters >= _max_iters )
278 return true;
279 else if( _lastLogZ.size() < 3 )
280 // need at least 2 to calculate ratio
281 // Also, throw away first iteration, as the parameters may not
282 // have been normalized according to the estimation method
283 return false;
284 else {
285 Real current = _lastLogZ[_lastLogZ.size() - 1];
286 Real previous = _lastLogZ[_lastLogZ.size() - 2];
287 if( previous == 0 )
288 return false;
289 Real diff = current - previous;
290 if( diff < 0 ) {
291 std::cerr << "Error: in EM log-likehood decreased from " << previous << " to " << current << std::endl;
292 return true;
293 }
294 return (diff / fabs(previous)) <= _log_z_tol;
295 }
296 }
297
298
299 Real EMAlg::iterate( MaximizationStep &mstep ) {
300 Real logZ = 0;
301 Real likelihood = 0;
302
303 _estep.run();
304 logZ = _estep.logZ();
305
306 // Expectation calculation
307 for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
308 InfAlg* clamped = _estep.clone();
309 // Apply evidence
310 for( Observation::const_iterator i = e->begin(); i != e->end(); ++i )
311 clamped->clamp( clamped->fg().findVar(i->first), i->second );
312 clamped->init();
313 clamped->run();
314
315 likelihood += clamped->logZ() - logZ;
316
317 mstep.addExpectations( *clamped );
318
319 delete clamped;
320 }
321
322 // Maximization of parameters
323 mstep.maximize( _estep.fg() );
324
325 return likelihood;
326 }
327
328
329 Real EMAlg::iterate() {
330 Real likelihood;
331 for( size_t i = 0; i < _msteps.size(); ++i )
332 likelihood = iterate( _msteps[i] );
333 _lastLogZ.push_back( likelihood );
334 ++_iters;
335 return likelihood;
336 }
337
338
339 void EMAlg::run() {
340 while( !hasSatisfiedTermConditions() )
341 iterate();
342 }
343
344
345 } // end of namespace dai