Polished EM code
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 28 Jul 2009 14:20:52 +0000 (16:20 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 28 Jul 2009 14:20:52 +0000 (16:20 +0200)
Makefile
include/dai/emalg.h
include/dai/evidence.h
src/emalg.cpp
tests/testem/runtests.bat [new file with mode: 0755]
tests/testem/testem.cpp
tests/testregression

index 37732cb..12eff2f 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -252,19 +252,21 @@ endif
 # REGRESSION TESTS
 ###################
 
-ifneq ($(OS),WINDOWS)
 testregression : tests/testdai$(EE)
        @echo Starting regression test...this can take a minute or so!
+ifneq ($(OS),WINDOWS)
        cd tests && ./testregression && cd ..
 else
-testregression : tests/testdai$(EE)
-       @echo Starting regression test...this can take a minute or so!
        cd tests && testregression.bat && cd ..
 endif
 
 testem : tests/testem/testem$(EE)
        @echo Starting EM tests
+ifneq ($(OS),WINDOWS)
        cd tests/testem && ./runtests && cd ../..
+else
+       cd tests\testem && runtests && cd ..\..
+endif
 
 
 # DOCUMENTATION
@@ -288,11 +290,12 @@ clean :
        -rm matlab/*$(ME)
        -rm examples/example$(EE) examples/example_bipgraph$(EE) examples/example_varset$(EE) examples/example_sprinkler$(EE)
        -rm tests/testdai$(EE)
+       -rm tests/testem/testem$(EE)
        -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE)
        -rm -R doc
        -rm -R lib
 else
 .PHONY : clean
 clean :
-       -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
+       -del *$(OE) *.ilk *.pdb *$(EE) matlab\*$(ME) examples\*$(EE) examples\*.ilk examples\*.pdb tests\testdai$(EE) tests\testem\testem$(EE) tests\*.pdb tests\*.ilk utils\*$(EE) utils\*.pdb utils\*.ilk $(LIB)\libdai$(LE)
 endif
index 1e353ff..b483612 100644 (file)
@@ -35,7 +35,7 @@
 
 /// \file 
 /** \brief Defines classes related to Expectation Maximization:
- *  EMAlg, ParameterEstimate, and FactorOrientations
+ *  EMAlg, ParameterEstimation, CondProbEstimation and SharedParameters
  */
 
 
@@ -87,38 +87,38 @@ class CondProbEstimation : private ParameterEstimation {
         Prob _stats;
         Prob _initial_stats;
     public:
-    /** For a conditional probability \f$ Pr( X | Y ) \f$, 
-     *  \param target_dimension should equal \f$ | X | \f$
-     *  \param pseudocounts has length \f$ |X| \cdot |Y| \f$
-     */
-    CondProbEstimation( size_t target_dimension, Prob pseudocounts );
-
-    /// Virtual constructor, using a PropertySet.
-    /** Some keys in the PropertySet are required:
-     *     - target_dimension, which should be equal to \f$ | X | \f$
-     *     - total_dimension, which sholud be equal to \f$ |X| \cdot |Y| \f$
-     *  
-     *  An optional key is:
-     *     - pseudo_count which specifies the initial counts (defaults to 1)
-     */
-    static ParameterEstimation* factory( const PropertySet& p );
-    /// Virtual destructor
-    virtual ~CondProbEstimation() {}
-    /// Returns an estimate of the conditional probability distribution.
-    /** The format of the resulting Prob keeps all the values for 
-     *  \f$ P(X | Y=a) \f$ sequential in the array.
-     */
-    virtual Prob estimate();
-    /// Accumulate sufficient statistics from the expectations in p.
-    virtual void addSufficientStatistics( Prob& p );
-    /// Returns the required size for arguments to addSufficientStatistics
-    virtual size_t probSize() const { 
-        return _stats.size(); 
-    }
-    /// Virtual copy constructor.
-    virtual ParameterEstimation* clone() const {
-        return new CondProbEstimation( _target_dim, _initial_stats );
-    }
+        /** For a conditional probability \f$ Pr( X | Y ) \f$, 
+         *  \param target_dimension should equal \f$ | X | \f$
+         *  \param pseudocounts has length \f$ |X| \cdot |Y| \f$
+         */
+        CondProbEstimation( size_t target_dimension, Prob pseudocounts );
+
+        /// Virtual constructor, using a PropertySet.
+        /** Some keys in the PropertySet are required:
+         *     - target_dimension, which should be equal to \f$ | X | \f$
+         *     - total_dimension, which should be equal to \f$ |X| \cdot |Y| \f$
+         *  
+         *  An optional key is:
+         *     - pseudo_count, which specifies the initial counts (defaults to 1)
+         */
+        static ParameterEstimation* factory( const PropertySet& p );
+        /// Virtual destructor
+        virtual ~CondProbEstimation() {}
+        /// Returns an estimate of the conditional probability distribution.
+        /** The format of the resulting Prob keeps all the values for 
+         *  \f$ P(X | Y=a) \f$ sequential in the array.
+         */
+        virtual Prob estimate();
+        /// Accumulate sufficient statistics from the expectations in p.
+        virtual void addSufficientStatistics( Prob& p );
+        /// Returns the required size for arguments to addSufficientStatistics
+        virtual size_t probSize() const { 
+            return _stats.size(); 
+        }
+        /// Virtual copy constructor.
+        virtual ParameterEstimation* clone() const {
+            return new CondProbEstimation( _target_dim, _initial_stats );
+        }
 };
 
 
@@ -127,7 +127,7 @@ class CondProbEstimation : private ParameterEstimation {
  *  canonical variable ordering.  This canonical variable ordering
  *  will likely not be the order of variables required to make two
  *  factors parameters isomorphic.  Therefore, this ordering of the
- *  variables must be specified for ever factory to ensure that
+ *  variables must be specified for each factory to ensure that
  *  parameters can be shared between different factors during EM.
  */
 class SharedParameters {
@@ -143,7 +143,9 @@ class SharedParameters {
         ParameterEstimation* _estimation;
         bool _deleteEstimation;
 
-        static Permute calculatePermutation( const std::vector<Var>& varorder, const std::vector<size_t>& dims, VarSet& outVS );
+        /// Calculates the permutation that permutes the variables in varorder into the canonical ordering
+        static Permute calculatePermutation( const std::vector<Var>& varorder, VarSet& outVS );
+        /// Initializes _varsets and _perms from _varorders
         void setPermsAndVarSetsFromVarOrders();
 
     public:
@@ -179,9 +181,10 @@ class MaximizationStep {
     private:
         std::vector<SharedParameters> _params;
     public:
+        /// Default constructor
         MaximizationStep() : _params() {}
 
-        /// Construct an step object that contains all these estimation problems
+        /// Construct a step object that contains all these estimation problems
         MaximizationStep( std::vector<SharedParameters>& maximizations ) : _params(maximizations) {}  
 
         /// Construct a step from an input stream
@@ -209,6 +212,10 @@ class MaximizationStep {
  *  MaximizationSteps.  An expectation step is performed between execution
  *  of each MaximizationStep.  A call to iterate() will cycle through all
  *  MaximizationSteps.
+ *
+ *  Having multiple and separate maximization steps allows for maximizing some
+ *  parameters, performing another E step, and then maximizing separate
+ *  parameters, which may result in faster convergence in some cases.
  */  
 class EMAlg {
     private:
@@ -220,10 +227,17 @@ class EMAlg {
 
         /// The maximization steps to take
         std::vector<MaximizationStep> _msteps;
+
+        /// Number of iterations done
         size_t _iters;
+
+        /// History of likelihoods
         std::vector<Real> _lastLogZ;
 
+        /// Maximum number of iterations
         size_t _max_iters;
+
+        /// Convergence tolerance
         Real _log_z_tol;
 
     public:
@@ -237,7 +251,7 @@ class EMAlg {
         static const Real LOG_Z_TOL_DEFAULT;
 
         /// Construct an EMAlg from all these objects
-        EMAlg( const Evidence& evidence, InfAlg& estep, std::vector<MaximizationStep>& msteps, PropertySet* termconditions = NULL) :
+        EMAlg( const Evidence& evidence, InfAlg& estep, std::vector<MaximizationStep>& msteps, const PropertySet &termconditions ) :
           _evidence(evidence), _estep(estep), _msteps(msteps), _iters(0), _lastLogZ(), _max_iters(MAX_ITERS_DEFAULT), _log_z_tol(LOG_Z_TOL_DEFAULT) { 
               setTermConditions( termconditions );
         }
@@ -245,14 +259,14 @@ class EMAlg {
         /// Construct an EMAlg from an input stream
         EMAlg( const Evidence& evidence, InfAlg& estep, std::istream& mstep_file );
 
-        /// Change the coditions for termination
-        /** There are two possible parameters in the PropertySety
+        /// Change the conditions for termination
+        /** There are two possible parameters in the PropertySet
          *    - max_iters  maximum number of iterations (default 30)
          *    - log_z_tol proportion of increase in logZ (default 0.01)
          *
          *  \see hasSatisifiedTermConditions()
          */
-        void setTermConditions( const PropertySetp );
+        void setTermConditions( const PropertySet &p );
 
         /// Determine if the termination conditions have been met.
         /** There are two sufficient termination conditions:
@@ -264,6 +278,7 @@ class EMAlg {
          */
         bool hasSatisfiedTermConditions() const;
 
+        /// Returns number of iterations done so far
         size_t getCurrentIters() const { 
             return _iters; 
         }
index 25ae67d..44a7984 100644 (file)
@@ -53,6 +53,9 @@ class Observation {
 
 
 /// Stores multiple observations of variables
+/** The Evidence class contains multiple samples, where a sample is the joint
+ *  observation of the states of some variables.
+ */
 class Evidence {
     private:
         std::vector<Observation> _samples;
index 0831599..3ef4be4 100644 (file)
@@ -87,28 +87,24 @@ Prob CondProbEstimation::estimate() {
 }
 
 
-Permute SharedParameters::calculatePermutation(const std::vector< Var >& varorder, const std::vector< size_t >& dims, VarSet& outVS) {
-    std::vector<long> labels(dims.size());
+Permute SharedParameters::calculatePermutation( const std::vector< Var >& varorder, VarSet& outVS ) {
+    std::vector<size_t> dims;
+    std::vector<long> labels;
   
-    // Check that the variable set is compatible
-    if (varorder.size() != dims.size())
-        DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
-  
-    // Collect all labels, and order them in vs
-    for (size_t di = 0; di < dims.size(); ++di) {
-        if (dims[di] != varorder[di].states())
-            DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
-        outVS |= varorder[di];
-        labels[di] = varorder[di].label();
+    // Collect all labels and dimensions, and order them in vs
+    dims.reserve( varorder.size() );
+    labels.reserve( varorder.size() );
+    for( size_t i = 0; i < varorder.size(); i++ ) {
+        dims.push_back( varorder[i].states() );
+        labels.push_back( varorder[i].label() );
+        outVS |= varorder[i];
     }
   
     // Construct the sigma array for the permutation object
-    std::vector<size_t> sigma(dims.size(), 0);
-    VarSet::iterator set_iterator = outVS.begin();
-    for (size_t vs_i = 0; vs_i < dims.size(); ++vs_i, ++set_iterator) {
-        std::vector< long >::iterator location = find(labels.begin(), labels.end(), set_iterator->label());
-        sigma[vs_i] = location - labels.begin();
-    }
+    std::vector<size_t> sigma;
+    sigma.reserve( dims.size() );
+    for (VarSet::iterator set_iterator = outVS.begin(); sigma.size() < dims.size(); ++set_iterator)
+        sigma.push_back( find(labels.begin(), labels.end(), set_iterator->label()) - labels.begin() );
   
     return Permute(dims, sigma);
 }
@@ -117,28 +113,22 @@ Permute SharedParameters::calculatePermutation(const std::vector< Var >& varorde
 void SharedParameters::setPermsAndVarSetsFromVarOrders() {
     if (_varorders.size() == 0)
         return;
-    FactorOrientations::const_iterator foi = _varorders.begin();
-    std::vector< size_t > dims(foi->second.size());
-    size_t total_dim = 1;
-    for (size_t i = 0; i < dims.size(); ++i) {
-        dims[i] = foi->second[i].states();
-        total_dim *= dims[i];
-    }
-  
-    // Construct the permutation objects
-    for ( ; foi != _varorders.end(); ++foi) {
+    if( _estimation == NULL )
+        DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
+
+    // Construct the permutation objects and the varsets
+    for( FactorOrientations::const_iterator foi = _varorders.begin(); foi != _varorders.end(); ++foi ) {
         VarSet vs;
-        _perms[foi->first] = calculatePermutation(foi->second, dims, vs);
+        _perms[foi->first] = calculatePermutation(foi->second, vs);
         _varsets[foi->first] = vs;
+        if( _estimation->probSize() != vs.nrStates() )
+            DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
     }
-  
-    if (_estimation == NULL || _estimation->probSize() != total_dim)
-        DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
 }
 
 
 SharedParameters::SharedParameters(std::istream& is, const FactorGraph& fg_varlookup)
-  : _varsets(), _perms(), _varorders(), _estimation(NULL), _deleteEstimation(1
+  : _varsets(), _perms(), _varorders(), _estimation(NULL), _deleteEstimation(true
 {
     std::string est_method;
     PropertySet props;
@@ -149,27 +139,28 @@ SharedParameters::SharedParameters(std::istream& is, const FactorGraph& fg_varlo
 
     size_t num_factors;
     is >> num_factors;
-    for (size_t sp_i = 0; sp_i < num_factors; ++sp_i) {
+    for( size_t sp_i = 0; sp_i < num_factors; ++sp_i ) {
         std::string line;
-        std::vector< std::string > fields;
-        size_t factor;
-        std::vector< Var > var_order;
-        std::istringstream iss;
-
         while(line.size() == 0 && getline(is, line))
             ;
+
+        std::vector< std::string > fields;
         tokenizeString(line, fields, " \t");
 
         // Lookup the factor in the factorgraph
         if (fields.size() < 1)
             DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+        std::istringstream iss;
         iss.str(fields[0]);
+        size_t factor;
         iss >> factor;
         const VarSet& vs = fg_varlookup.factor(factor).vars();
         if (fields.size() != vs.size() + 1)
             DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
 
         // Construct the vector of Vars
+        std::vector< Var > var_order;
+        var_order.reserve( vs.size() );
         for (size_t fi = 1; fi < fields.size(); ++fi) {
             // Lookup a single variable by label
             long label;
@@ -192,21 +183,20 @@ SharedParameters::SharedParameters(std::istream& is, const FactorGraph& fg_varlo
 SharedParameters::SharedParameters( const SharedParameters& sp ) 
   : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _deleteEstimation(sp._deleteEstimation)
 {
-  if (_deleteEstimation)
-      _estimation = _estimation->clone();
+    if (_deleteEstimation)
+        _estimation = _estimation->clone();
 }
 
 
 SharedParameters::SharedParameters( const FactorOrientations& varorders, ParameterEstimation* estimation )
-  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _deleteEstimation(0
+  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _deleteEstimation(false
 {
-  setPermsAndVarSetsFromVarOrders();
+    setPermsAndVarSetsFromVarOrders();
 }
 
 
 void SharedParameters::collectSufficientStatistics(InfAlg& alg) {
-    std::map< FactorIndex, Permute >::iterator i = _perms.begin();
-    for ( ; i != _perms.end(); ++i) {
+    for ( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i) {
         Permute& perm = i->second;
         VarSet& vs = _varsets[i->first];
         
@@ -221,8 +211,7 @@ void SharedParameters::collectSufficientStatistics(InfAlg& alg) {
 
 void SharedParameters::setParameters(FactorGraph& fg) {
     Prob p = _estimation->estimate();
-    std::map< FactorIndex, Permute >::iterator i = _perms.begin();
-    for ( ; i != _perms.end(); ++i) {
+    for ( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i) {
         Permute& perm = i->second;
         VarSet& vs = _varsets[i->first];
         
@@ -280,46 +269,45 @@ EMAlg::EMAlg(const Evidence& evidence, InfAlg& estep, std::istream& msteps_file)
 }      
 
 
-void EMAlg::setTermConditions(const PropertySet* p) {
-    if (NULL == p)
-        return;
-    if (p->hasKey(MAX_ITERS_KEY))
-        _max_iters = p->getStringAs<size_t>(MAX_ITERS_KEY);
-    if (p->hasKey(LOG_Z_TOL_KEY))
-        _log_z_tol = p->getStringAs<Real>(LOG_Z_TOL_KEY);
+void EMAlg::setTermConditions(const PropertySet &p) {
+    if (p.hasKey(MAX_ITERS_KEY))
+        _max_iters = p.getStringAs<size_t>(MAX_ITERS_KEY);
+    if (p.hasKey(LOG_Z_TOL_KEY))
+        _log_z_tol = p.getStringAs<Real>(LOG_Z_TOL_KEY);
 }
 
 
 bool EMAlg::hasSatisfiedTermConditions() const {
     if (_iters >= _max_iters) {
-        return 1;
+        return true;
     } else if (_lastLogZ.size() < 3) { 
         // need at least 2 to calculate ratio
         // Also, throw away first iteration, as the parameters may not
         // have been normalized according to the estimation method
-        return 0;
+        return false;
     } else {
         Real current = _lastLogZ[_lastLogZ.size() - 1];
         Real previous = _lastLogZ[_lastLogZ.size() - 2];
-        if (previous == 0) return 0;
+        if (previous == 0) 
+            return false;
         Real diff = current - previous;
         if (diff < 0) {
             std::cerr << "Error: in EM log-likehood decreased from " << previous << " to " << current << std::endl;
-            return 1;
+            return true;
         }
         return diff / abs(previous) <= _log_z_tol;
     }
 }
 
 
-Real EMAlg::iterate(MaximizationStep& mstep) {
-    Evidence::const_iterator e = _evidence.begin();
+Real EMAlg::iterate( MaximizationStep& mstep ) {
     Real logZ = 0;
 
     // Expectation calculation
-    for ( ; e != _evidence.end(); ++e) {
+    for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
         InfAlg* clamped = _estep.clone();
         e->applyEvidence(*clamped);
+        clamped->init();
         clamped->run();
       
         logZ += clamped->logZ();
diff --git a/tests/testem/runtests.bat b/tests/testem/runtests.bat
new file mode 100755 (executable)
index 0000000..44bdf2e
--- /dev/null
@@ -0,0 +1,7 @@
+testem 2var.fg 2var_data.tab 2var.em > testem.out.tmp
+testem 3var.fg 2var_data.tab 3var.em >> testem.out.tmp
+testem ..\hoi1.fg hoi1_data.tab hoi1_share_f0_f2.em >> testem.out.tmp
+testem ..\hoi1.fg hoi1_data.tab hoi1_share_f0_f1_f2.em >> testem.out.tmp
+diff -s testem.out.tmp testem.out
+
+del testem.out.tmp
index 74e5c99..a6a62f2 100644 (file)
@@ -42,7 +42,7 @@ void usage( const string& msg ) {
 
 int main( int argc, char** argv ) {
     if( argc != 4 )
-      usage("Incorrect number of arguments.");
+        usage("Incorrect number of arguments.");
     
     FactorGraph fg;
     fg.ReadFromFile( argv[1] );
@@ -58,8 +58,7 @@ int main( int argc, char** argv ) {
     e.addEvidenceTabFile( estream, fg );
 
     cout << "Number of samples: " << e.nrSamples() << endl;
-    Evidence::iterator ps = e.begin();
-    for( ; ps != e.end(); ps++ )
+    for( Evidence::iterator ps = e.begin(); ps != e.end(); ps++ )
         cout << "Sample #" << (ps - e.begin()) << " has " << ps->observations().size() << " observations." << endl;
 
     ifstream emstream( argv[3] );
index 55630a1..cdae018 100755 (executable)
@@ -3,6 +3,6 @@ TMPFILE1=`mktemp /var/tmp/testfast.XXXXXX`
 trap 'rm -f $TMP_FILE' 0 1 15
 
 ./testall testfast.fg > $TMPFILE1
-diff -s $TMPFILE1 testfast.out
+diff -s $TMPFILE1 testfast.out || exit 1
 
 rm -f $TMPFILE1