Updated copyright headers
[libdai.git] / include / dai / emalg.h
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 #ifndef __defined_libdai_emalg_h
13 #define __defined_libdai_emalg_h
14
15
16 #include <vector>
17 #include <map>
18
19 #include <dai/factor.h>
20 #include <dai/daialg.h>
21 #include <dai/evidence.h>
22 #include <dai/index.h>
23 #include <dai/properties.h>
24
25
26 /// \file
27 /// \brief Defines classes related to Expectation Maximization: EMAlg, ParameterEstimation, CondProbEstimation and SharedParameters
28 /// \todo Describe EM file format
29
30
31 namespace dai {
32
33
34 /// Interface for a parameter estimation method.
35 /** This parameter estimation interface is based on sufficient statistics.
36 * Implementations are responsible for collecting data from a probability
37 * vector passed to it from a SharedParameters container object.
38 *
39 * Implementations of this interface should register a factory function
40 * via the static ParameterEstimation::registerMethod function.
41 * The default registry only contains CondProbEstimation, named
42 * "ConditionalProbEstimation".
43 */
44 class ParameterEstimation {
45 public:
46 /// A pointer to a factory function.
47 typedef ParameterEstimation* (*ParamEstFactory)( const PropertySet& );
48
49 /// Virtual destructor for deleting pointers to derived classes.
50 virtual ~ParameterEstimation() {}
51 /// Virtual copy constructor.
52 virtual ParameterEstimation* clone() const = 0;
53
54 /// General factory method for construction of ParameterEstimation subclasses.
55 static ParameterEstimation* construct( const std::string &method, const PropertySet &p );
56
57 /// Register a subclass so that it can be used with ParameterEstimation::construct.
58 static void registerMethod( const std::string &method, const ParamEstFactory &f ) {
59 if( _registry == NULL )
60 loadDefaultRegistry();
61 (*_registry)[method] = f;
62 }
63
64 /// Estimate the factor using the accumulated sufficient statistics and reset.
65 virtual Prob estimate() = 0;
66
67 /// Accumulate the sufficient statistics for p.
68 virtual void addSufficientStatistics( const Prob &p ) = 0;
69
70 /// Returns the size of the Prob that should be passed to addSufficientStatistics.
71 virtual size_t probSize() const = 0;
72
73 private:
74 /// A static registry containing all methods registered so far.
75 static std::map<std::string, ParamEstFactory> *_registry;
76
77 /// Registers default ParameterEstimation subclasses (currently, only CondProbEstimation).
78 static void loadDefaultRegistry();
79 };
80
81
82 /// Estimates the parameters of a conditional probability table, using pseudocounts.
83 class CondProbEstimation : private ParameterEstimation {
84 private:
85 /// Number of states of the variable of interest
86 size_t _target_dim;
87 /// Current pseudocounts
88 Prob _stats;
89 /// Initial pseudocounts
90 Prob _initial_stats;
91
92 public:
93 /// Constructor
94 /** For a conditional probability \f$ P( X | Y ) \f$,
95 * \param target_dimension should equal \f$ | X | \f$
96 * \param pseudocounts has length \f$ |X| \cdot |Y| \f$
97 */
98 CondProbEstimation( size_t target_dimension, const Prob &pseudocounts );
99
100 /// Virtual constructor, using a PropertySet.
101 /** Some keys in the PropertySet are required.
102 * For a conditional probability \f$ P( X | Y ) \f$,
103 * - target_dimension should be equal to \f$ | X | \f$
104 * - total_dimension should be equal to \f$ |X| \cdot |Y| \f$
105 *
106 * An optional key is:
107 * - pseudo_count, which specifies the initial counts (defaults to 1)
108 */
109 static ParameterEstimation* factory( const PropertySet &p );
110
111 /// Virtual copy constructor
112 virtual ParameterEstimation* clone() const { return new CondProbEstimation( _target_dim, _initial_stats ); }
113
114 /// Virtual destructor
115 virtual ~CondProbEstimation() {}
116
117 /// Returns an estimate of the conditional probability distribution.
118 /** The format of the resulting Prob keeps all the values for
119 * \f$ P(X | Y=y) \f$ in sequential order in the array.
120 */
121 virtual Prob estimate();
122
123 /// Accumulate sufficient statistics from the expectations in p.
124 virtual void addSufficientStatistics( const Prob &p );
125
126 /// Returns the required size for arguments to addSufficientStatistics.
127 virtual size_t probSize() const { return _stats.size(); }
128 };
129
130
131 /// A single factor or set of factors whose parameters should be estimated.
132 /** To ensure that parameters can be shared between different factors during
133 * EM learning, each factor's values are reordered to match a desired variable
134 * ordering. The ordering of the variables in a factor may therefore differ
135 * from the canonical ordering used in libDAI. The SharedParameters
136 * class couples one or more factors (together with the specified orderings
137 * of the variables) with a ParameterEstimation object, taking care of the
138 * necessary permutations of the factor entries / parameters.
139 */
140 class SharedParameters {
141 public:
142 /// Convenience label for an index into a factor in a FactorGraph.
143 typedef size_t FactorIndex;
144 /// Convenience label for a grouping of factor orientations.
145 typedef std::map<FactorIndex, std::vector<Var> > FactorOrientations;
146
147 private:
148 /// Maps factor indices to the corresponding VarSets
149 std::map<FactorIndex, VarSet> _varsets;
150 /// Maps factor indices to the corresponding Permute objects that permute the desired ordering into the canonical ordering
151 std::map<FactorIndex, Permute> _perms;
152 /// Maps factor indices to the corresponding desired variable orderings
153 FactorOrientations _varorders;
154 /// Parameter estimation method to be used
155 ParameterEstimation *_estimation;
156 /// Indicates whether the object pointed to by _estimation should be deleted upon destruction
157 bool _deleteEstimation;
158
159 /// Calculates the permutation that permutes the variables in varorder into the canonical ordering
160 static Permute calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS );
161
162 /// Initializes _varsets and _perms from _varorders
163 void setPermsAndVarSetsFromVarOrders();
164
165 public:
166 /// Copy constructor
167 SharedParameters( const SharedParameters &sp );
168
169 /// Constructor
170 /** \param varorders all the factor orientations for this parameter
171 * \param estimation a pointer to the parameter estimation method
172 * \param deletePE whether the parameter estimation object should be deleted in the destructor
173 */
174 SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool deletePE=false );
175
176 /// Constructor for making an object from a stream and a factor graph
177 SharedParameters( std::istream &is, const FactorGraph &fg_varlookup );
178
179 /// Destructor
180 ~SharedParameters() {
181 if( _deleteEstimation )
182 delete _estimation;
183 }
184
185 /// Collect the necessary statistics from expected values
186 void collectSufficientStatistics( InfAlg &alg );
187
188 /// Estimate and set the shared parameters
189 void setParameters( FactorGraph &fg );
190
191 /// Returns the parameters
192 void collectParameters( const FactorGraph &fg, std::vector<Real> &outVals, std::vector<Var> &outVarOrder );
193 };
194
195
196 /// A MaximizationStep groups together several parameter estimation tasks into a single unit.
197 class MaximizationStep {
198 private:
199 std::vector<SharedParameters> _params;
200
201 public:
202 /// Default constructor
203 MaximizationStep() : _params() {}
204
205 /// Constructor from a vector of SharedParameters objects
206 MaximizationStep( std::vector<SharedParameters> &maximizations ) : _params(maximizations) {}
207
208 /// Constructor from an input stream and a corresponding factor graph
209 MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup );
210
211 /// Collect the beliefs from this InfAlg as expectations for the next Maximization step.
212 void addExpectations( InfAlg &alg );
213
214 /// Using all of the currently added expectations, make new factors with maximized parameters and set them in the FactorGraph.
215 void maximize( FactorGraph &fg );
216
217 /// @name Iterator interface
218 //@{
219 typedef std::vector<SharedParameters>::iterator iterator;
220 typedef std::vector<SharedParameters>::const_iterator const_iterator;
221 iterator begin() { return _params.begin(); }
222 const_iterator begin() const { return _params.begin(); }
223 iterator end() { return _params.end(); }
224 const_iterator end() const { return _params.end(); }
225 //@}
226 };
227
228
229 /// EMAlg performs Expectation Maximization to learn factor parameters.
230 /** This requires specifying:
231 * - Evidence (instances of observations from the graphical model),
232 * - InfAlg for performing the E-step, which includes the factor graph,
233 * - a vector of MaximizationSteps steps to be performed.
234 *
235 * This implementation can perform incremental EM by using multiple
236 * MaximizationSteps. An expectation step is performed between execution
237 * of each MaximizationStep. A call to iterate() will cycle through all
238 * MaximizationSteps.
239 *
240 * Having multiple and separate maximization steps allows for maximizing some
241 * parameters, performing another E step, and then maximizing separate
242 * parameters, which may result in faster convergence in some cases.
243 */
244 class EMAlg {
245 private:
246 /// All the data samples used during learning
247 const Evidence &_evidence;
248
249 /// How to do the expectation step
250 InfAlg &_estep;
251
252 /// The maximization steps to take
253 std::vector<MaximizationStep> _msteps;
254
255 /// Number of iterations done
256 size_t _iters;
257
258 /// History of likelihoods
259 std::vector<Real> _lastLogZ;
260
261 /// Maximum number of iterations
262 size_t _max_iters;
263
264 /// Convergence tolerance
265 Real _log_z_tol;
266
267 public:
268 /// Key for setting maximum iterations @see setTermConditions
269 static const std::string MAX_ITERS_KEY;
270 /// Default maximum iterations @see setTermConditions
271 static const size_t MAX_ITERS_DEFAULT;
272 /// Key for setting likelihood termination condition @see setTermConditions
273 static const std::string LOG_Z_TOL_KEY;
274 /// Default likelihood tolerance @see setTermConditions
275 static const Real LOG_Z_TOL_DEFAULT;
276
277 /// Construct an EMAlg from all these objects
278 EMAlg( const Evidence &evidence, InfAlg &estep, std::vector<MaximizationStep> &msteps, const PropertySet &termconditions )
279 : _evidence(evidence), _estep(estep), _msteps(msteps), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT)
280 {
281 setTermConditions( termconditions );
282 }
283
284 /// Construct an EMAlg from an Evidence object, an InfAlg object, and an input stream
285 EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &mstep_file );
286
287 /// Change the conditions for termination
288 /** There are two possible parameters in the PropertySet
289 * - max_iters maximum number of iterations
290 * - log_z_tol proportion of increase in logZ
291 *
292 * \see hasSatisifiedTermConditions()
293 */
294 void setTermConditions( const PropertySet &p );
295
296 /// Determine if the termination conditions have been met.
297 /** There are two sufficient termination conditions:
298 * -# the maximum number of iterations has been performed
299 * -# the ratio of logZ increase over previous logZ is less than the
300 * tolerance, i.e.,
301 * \f$ \frac{\log(Z_t) - \log(Z_{t-1})}{| \log(Z_{t-1}) | } < \mathrm{tol} \f$.
302 */
303 bool hasSatisfiedTermConditions() const;
304
305 /// Return the last calculated log likelihood
306 Real getLogZ() const { return _lastLogZ.back(); }
307
308 /// Returns number of iterations done so far
309 size_t getCurrentIters() const { return _iters; }
310
311 /// Get the iteration method used
312 const InfAlg& eStep() const { return _estep; }
313
314 /// Perform an iteration over all maximization steps
315 Real iterate();
316
317 /// Perform an iteration over a single MaximizationStep
318 Real iterate( MaximizationStep &mstep );
319
320 /// Iterate until termination conditions are satisfied
321 void run();
322
323 /// @name Iterator interface
324 //@{
325 typedef std::vector<MaximizationStep>::iterator s_iterator;
326 typedef std::vector<MaximizationStep>::const_iterator const_s_iterator;
327 s_iterator s_begin() { return _msteps.begin(); }
328 const_s_iterator s_begin() const { return _msteps.begin(); }
329 s_iterator s_end() { return _msteps.end(); }
330 const_s_iterator s_end() const { return _msteps.end(); }
331 //@}
332 };
333
334
335 } // end of namespace dai
336
337
338 #endif