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