New git HEAD version
[libdai.git] / src / emalg.cpp
index 3e50ec2..68969f8 100644 (file)
@@ -1,11 +1,8 @@
 /*  This file is part of libDAI - http://www.libdai.org/
  *
- *  libDAI is licensed under the terms of the GNU General Public License version
- *  2, or (at your option) any later version. libDAI is distributed without any
- *  warranty. See the file COPYING for more details.
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
  *
- *  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
- *  Copyright (C) 2009  University of California, Santa Cruz
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
  */
 
 
@@ -37,6 +34,8 @@ ParameterEstimation* ParameterEstimation::construct( const std::string &method,
 }
 
 
+std::string CondProbEstimation::_name = "CondProbEstimation";
+
 ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
     size_t target_dimension = p.getStringAs<size_t>("target_dim");
     size_t total_dimension = p.getStringAs<size_t>("total_dim");
@@ -46,35 +45,39 @@ ParameterEstimation* CondProbEstimation::factory( const PropertySet &p ) {
     return new CondProbEstimation( target_dimension, Prob( total_dimension, pseudo_count ) );
 }
 
-
 CondProbEstimation::CondProbEstimation( size_t target_dimension, const Prob &pseudocounts )
-  : _target_dim(target_dimension), _stats(pseudocounts), _initial_stats(pseudocounts)
+  : _target_dim(target_dimension), _initial_stats(pseudocounts), _props()
 {
-    DAI_ASSERT( !(_stats.size() % _target_dim) );
+    DAI_ASSERT( !(_initial_stats.size() % _target_dim) );
+    _props.setAsString<size_t>("target_dim", _target_dim);
+    _props.setAsString<size_t>("total_dim", _initial_stats.size());
+    _props.setAsString<Real>("pseudo_count", _initial_stats.get(0));
 }
 
 
-void CondProbEstimation::addSufficientStatistics( const Prob &p ) {
-    _stats += p;
+Prob CondProbEstimation::parametersToFactor(const Prob& p) {
+    Prob result(p);
+    return result;
 }
 
-
-Prob CondProbEstimation::estimate() {
+// In the case of a conditional probability table, the 
+// parameters are identical to the estimated factor
+Prob CondProbEstimation::parameters(const Prob& p) {
+    Prob stats = p + _initial_stats;
     // normalize pseudocounts
-    for( size_t parent = 0; parent < _stats.size(); parent += _target_dim ) {
+    for( size_t parent = 0; parent < stats.size(); parent += _target_dim ) {
         // calculate norm
         size_t top = parent + _target_dim;
-        Real norm = std::accumulate( &(_stats[parent]), &(_stats[parent]) + _target_dim, 0.0 );
+        Real norm = 0.0;
+        for( size_t i = parent; i < top; ++i )
+            norm += stats[i];
         if( norm != 0.0 )
             norm = 1.0 / norm;
         // normalize
         for( size_t i = parent; i < top; ++i )
-            _stats[i] *= norm;
+            stats.set( i, stats[i] * norm );
     }
-    // reset _stats to _initial_stats
-    Prob result = _stats;
-    _stats = _initial_stats;
-    return result;
+    return stats;
 }
 
 
@@ -88,19 +91,20 @@ void SharedParameters::setPermsAndVarSetsFromVarOrders() {
     if( _varorders.size() == 0 )
         return;
     DAI_ASSERT( _estimation != NULL );
+    _expectations = new Prob(_estimation->probSize(), 0);
 
     // 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, vs );
         _varsets[foi->first] = vs;
-        DAI_ASSERT( _estimation->probSize() == vs.nrStates() );
+        DAI_ASSERT( (BigInt)_estimation->probSize() == vs.nrStates() );
     }
 }
 
 
 SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool ownPE )
-  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE)
+  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE), _expectations(NULL)
 {
     // Calculate the necessary permutations and varsets
     setPermsAndVarSetsFromVarOrders();
@@ -108,7 +112,7 @@ SharedParameters::SharedParameters( const FactorOrientations &varorders, Paramet
 
 
 SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
-  : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true)
+  : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true), _expectations(NULL)
 {
     // Read the desired parameter estimation method from the stream
     std::string est_method;
@@ -127,8 +131,7 @@ SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
         while( line.size() == 0 && getline(is, line) )
             ;
 
-        std::vector<std::string> fields;
-        tokenizeString(line, fields, " \t");
+        std::vector<std::string> fields = tokenizeString( line, true, " \t" );
 
         // Lookup the factor in the factorgraph
         if( fields.size() < 1 )
@@ -165,29 +168,29 @@ SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
 }
 
 
-void SharedParameters::collectSufficientStatistics( InfAlg &alg ) {
+void SharedParameters::collectExpectations( InfAlg &alg ) {
     for( std::map< FactorIndex, Permute >::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
         Permute &perm = i->second;
         VarSet &vs = _varsets[i->first];
 
         Factor b = alg.belief(vs);
-        Prob p( b.states(), 0.0 );
-        for( size_t entry = 0; entry < b.states(); ++entry )
-            p[entry] = b[perm.convertLinearIndex(entry)]; // apply inverse permutation
-        _estimation->addSufficientStatistics( p );
+        Prob p( b.nrStates(), 0.0 );
+        for( size_t entry = 0; entry < b.nrStates(); ++entry )
+            p.set( entry, b[perm.convertLinearIndex(entry)] ); // apply inverse permutation
+        (*_expectations) += p;
     }
 }
 
 
 void SharedParameters::setParameters( FactorGraph &fg ) {
-    Prob p = _estimation->estimate();
+    Prob p = _estimation->estimate(this->currentExpectations());
     for( std::map<FactorIndex, Permute>::iterator i = _perms.begin(); i != _perms.end(); ++i ) {
         Permute &perm = i->second;
         VarSet &vs = _varsets[i->first];
 
         Factor f( vs, 0.0 );
-        for( size_t entry = 0; entry < f.states(); ++entry )
-            f[perm.convertLinearIndex(entry)] = p[entry];
+        for( size_t entry = 0; entry < f.nrStates(); ++entry )
+            f.set( perm.convertLinearIndex(entry), p[entry] );
 
         fg.setFactor( i->first, f );
     }
@@ -205,7 +208,7 @@ MaximizationStep::MaximizationStep( std::istream &is, const FactorGraph &fg_varl
 
 void MaximizationStep::addExpectations( InfAlg &alg ) {
     for( size_t i = 0; i < _params.size(); ++i )
-        _params[i].collectSufficientStatistics( alg );
+        _params[i].collectExpectations( alg );
 }
 
 
@@ -215,6 +218,11 @@ void MaximizationStep::maximize( FactorGraph &fg ) {
 }
 
 
+void MaximizationStep::clear( ) {
+    for( size_t i = 0; i < _params.size(); ++i )
+        _params[i].clear( );
+}
+
 const std::string EMAlg::MAX_ITERS_KEY("max_iters");
 const std::string EMAlg::LOG_Z_TOL_KEY("log_z_tol");
 const size_t EMAlg::MAX_ITERS_DEFAULT = 30;
@@ -268,6 +276,8 @@ Real EMAlg::iterate( MaximizationStep &mstep ) {
     Real logZ = 0;
     Real likelihood = 0;
 
+    mstep.clear();
+
     _estep.run();
     logZ = _estep.logZ();