New git HEAD version
[libdai.git] / src / emalg.cpp
index e0fcab3..68969f8 100644 (file)
@@ -1,22 +1,9 @@
-/*  Copyright (C) 2009  Charles Vaske  [cvaske at soe dot ucsc dot edu]
-    University of California Santa Cruz
-
-    This file is part of libDAI.
-
-    libDAI is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    libDAI is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with libDAI; if not, write to the Free Software
-    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
-*/
+/*  This file is part of libDAI - http://www.libdai.org/
+ *
+ *  Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
+ *
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
+ */
 
 
 #include <dai/util.h>
 namespace dai {
 
 
+// Initialize static private member of ParameterEstimation
 std::map<std::string, ParameterEstimation::ParamEstFactory> *ParameterEstimation::_registry = NULL;
 
 
 void ParameterEstimation::loadDefaultRegistry() {
     _registry = new std::map<std::string, ParamEstFactory>();
-    (*_registry)["ConditionalProbEstimation"] = CondProbEstimation::factory;
+    (*_registry)["CondProbEstimation"] = CondProbEstimation::factory;
 }
 
 
@@ -40,12 +28,14 @@ ParameterEstimation* ParameterEstimation::construct( const std::string &method,
         loadDefaultRegistry();
     std::map<std::string, ParamEstFactory>::iterator i = _registry->find(method);
     if( i == _registry->end() )
-        DAI_THROW(UNKNOWN_PARAMETER_ESTIMATION_METHOD);
+        DAI_THROWE(UNKNOWN_PARAMETER_ESTIMATION_METHOD, "Unknown parameter estimation method: " + method);
     ParamEstFactory factory = i->second;
     return factory(p);
 }
 
 
+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");
@@ -55,82 +45,74 @@ 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) 
+CondProbEstimation::CondProbEstimation( size_t target_dimension, const Prob &pseudocounts )
+  : _target_dim(target_dimension), _initial_stats(pseudocounts), _props()
 {
-    if( _stats.size() % _target_dim )
-        DAI_THROW(MALFORMED_PROPERTY);
+    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
-        Real norm = 0.0;
         size_t top = parent + _target_dim;
+        Real norm = 0.0;
         for( size_t i = parent; i < top; ++i )
-            norm += _stats[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;
 }
 
 
-Permute SharedParameters::calculatePermutation( const std::vector<Var> &varorder, VarSet &outVS ) {
-    // Collect all labels and dimensions, and order them in vs
-    std::vector<size_t> dims;
-    dims.reserve( varorder.size() );
-    std::vector<long> labels;
-    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;
-    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 );
+Permute SharedParameters::calculatePermutation( const std::vector<Var> &varOrder, VarSet &outVS ) {
+    outVS = VarSet( varOrder.begin(), varOrder.end(), varOrder.size() );
+    return Permute( varOrder );
 }
 
 
 void SharedParameters::setPermsAndVarSetsFromVarOrders() {
     if( _varorders.size() == 0 )
         return;
-    if( _estimation == NULL )
-        DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
+    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;
-        if( _estimation->probSize() != vs.nrStates() )
-            DAI_THROW(INVALID_SHARED_PARAMETERS_ORDER);
+        DAI_ASSERT( (BigInt)_estimation->probSize() == vs.nrStates() );
     }
 }
 
 
-SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg_varlookup )
-  : _varsets(), _perms(), _varorders(), _estimation(NULL), _deleteEstimation(true) 
+SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool ownPE )
+  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _ownEstimation(ownPE), _expectations(NULL)
+{
+    // Calculate the necessary permutations and varsets
+    setPermsAndVarSetsFromVarOrders();
+}
+
+
+SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg )
+  : _varsets(), _perms(), _varorders(), _estimation(NULL), _ownEstimation(true), _expectations(NULL)
 {
     // Read the desired parameter estimation method from the stream
     std::string est_method;
@@ -149,34 +131,33 @@ SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg_varl
         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 )
-            DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+            DAI_THROWE(INVALID_EMALG_FILE,"Empty line unexpected");
         std::istringstream iss;
         iss.str( fields[0] );
         size_t factor;
         iss >> factor;
-        const VarSet &vs = fg_varlookup.factor(factor).vars();
+        const VarSet &vs = fg.factor(factor).vars();
         if( fields.size() != vs.size() + 1 )
-            DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+            DAI_THROWE(INVALID_EMALG_FILE,"Number of fields does not match factor size");
 
         // 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;
+            size_t label;
             std::istringstream labelparse( fields[fi] );
             labelparse >> label;
             VarSet::const_iterator vsi = vs.begin();
             for( ; vsi != vs.end(); ++vsi )
-                if( vsi->label() == label ) 
+                if( vsi->label() == label )
                     break;
             if( vsi == vs.end() )
-                DAI_THROW(INVALID_SHARED_PARAMETERS_INPUT_LINE);
+                DAI_THROWE(INVALID_EMALG_FILE,"Specified variables do not match the factor variables");
             var_order.push_back( *vsi );
         }
         _varorders[factor] = var_order;
@@ -187,72 +168,35 @@ SharedParameters::SharedParameters( std::istream &is, const FactorGraph &fg_varl
 }
 
 
-SharedParameters::SharedParameters( const SharedParameters &sp ) 
-  : _varsets(sp._varsets), _perms(sp._perms), _varorders(sp._varorders), _estimation(sp._estimation), _deleteEstimation(sp._deleteEstimation)
-{
-    // If sp owns its _estimation object, we should clone it instead
-    if( _deleteEstimation )
-        _estimation = _estimation->clone();
-}
-
-
-SharedParameters::SharedParameters( const FactorOrientations &varorders, ParameterEstimation *estimation, bool deleteParameterEstimationInDestructor )
-  : _varsets(), _perms(), _varorders(varorders), _estimation(estimation), _deleteEstimation(deleteParameterEstimationInDestructor) 
-{
-    // Calculate the necessary permutations
-    setPermsAndVarSetsFromVarOrders();
-}
-
-
-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.convert_linear_index(entry)];
-        _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.convert_linear_index(entry)] = p[entry];
+        for( size_t entry = 0; entry < f.nrStates(); ++entry )
+            f.set( perm.convertLinearIndex(entry), p[entry] );
 
         fg.setFactor( i->first, f );
     }
 }
 
 
-void SharedParameters::collectParameters( const FactorGraph& fg, std::vector< Real >& outVals, std::vector< Var >& outVarOrder ) {
-    FactorOrientations::iterator it = _varorders.begin();
-    if (it == _varorders.end()) {
-       return;
-    }
-    FactorIndex i = it->first;
-    std::vector< Var >::iterator var_it = _varorders[i].begin();
-    std::vector< Var >::iterator var_stop = _varorders[i].end();
-    for ( ; var_it != var_stop; ++var_it) {
-       outVarOrder.push_back(*var_it);
-    }
-    const Factor& f = fg.factor(i);
-    assert(f.vars() == _varsets[i]);
-    const Permute& perm = _perms[i];
-    for (size_t val_index = 0; val_index < f.states(); ++val_index) {
-       outVals.push_back(f[perm.convert_linear_index(val_index)]);
-    }
-}
-
-
 MaximizationStep::MaximizationStep( std::istream &is, const FactorGraph &fg_varlookup ) : _params() {
     size_t num_params = -1;
     is >> num_params;
@@ -264,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 );
 }
 
 
@@ -274,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;
@@ -289,7 +238,7 @@ EMAlg::EMAlg( const Evidence &evidence, InfAlg &estep, std::istream &msteps_file
     _msteps.reserve(num_msteps);
     for( size_t i = 0; i < num_msteps; ++i )
         _msteps.push_back( MaximizationStep( msteps_file, estep.fg() ) );
-}      
+}
 
 
 void EMAlg::setTermConditions( const PropertySet &p ) {
@@ -325,29 +274,33 @@ bool EMAlg::hasSatisfiedTermConditions() const {
 
 Real EMAlg::iterate( MaximizationStep &mstep ) {
     Real logZ = 0;
-       Real likelihood = 0;
-       
-       _estep.run();
-       logZ = _estep.logZ();
-       
+    Real likelihood = 0;
+
+    mstep.clear();
+
+    _estep.run();
+    logZ = _estep.logZ();
+
     // Expectation calculation
     for( Evidence::const_iterator e = _evidence.begin(); e != _evidence.end(); ++e ) {
         InfAlg* clamped = _estep.clone();
-        e->applyEvidence( *clamped );
+        // Apply evidence
+        for( Evidence::Observation::const_iterator i = e->begin(); i != e->end(); ++i )
+            clamped->clamp( clamped->fg().findVar(i->first), i->second );
         clamped->init();
         clamped->run();
-               
-               likelihood += clamped->logZ() - logZ;
-               
+
+        likelihood += clamped->logZ() - logZ;
+
         mstep.addExpectations( *clamped );
-               
+
         delete clamped;
     }
-    
+
     // Maximization of parameters
     mstep.maximize( _estep.fg() );
-       
-       return likelihood;
+
+    return likelihood;
 }