Merge branch 'pletscher'
[libdai.git] / src / gibbs.cpp
index 96a8a17..c1bff37 100644 (file)
@@ -1,21 +1,11 @@
-/*  Copyright (C) 2008  Frederik Eaton [frederik at ofb dot net]
-
-    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/
+ *
+ *  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) 2008  Frederik Eaton  [frederik at ofb dot net]
+ */
 
 
 #include <iostream>
@@ -38,7 +28,7 @@ const char *Gibbs::Name = "GIBBS";
 
 
 void Gibbs::setProperties( const PropertySet &opts ) {
-    assert( opts.hasKey("iters") );
+    DAI_ASSERT( opts.hasKey("iters") );
     props.iters = opts.getStringAs<size_t>("iters");
 
     if( opts.hasKey("verbose") )
@@ -70,7 +60,7 @@ void Gibbs::construct() {
     _var_counts.reserve( nrVars() );
     for( size_t i = 0; i < nrVars(); i++ )
         _var_counts.push_back( _count_t( var(i).states(), 0 ) );
-    
+
     _factor_counts.clear();
     _factor_counts.reserve( nrFactors() );
     for( size_t I = 0; I < nrFactors(); I++ )
@@ -121,11 +111,11 @@ inline size_t Gibbs::getFactorEntryDiff( size_t I, size_t i ) {
 
 
 Prob Gibbs::getVarDist( size_t i ) {
-    assert( i < nrVars() );
+    DAI_ASSERT( i < nrVars() );
     size_t i_states = var(i).states();
     Prob i_given_MB( i_states, 1.0 );
 
-    // use markov blanket of var(i) to calculate distribution
+    // use Markov blanket of var(i) to calculate distribution
     foreach( const Neighbor &I, nbV(i) ) {
         const Factor &f_I = factor(I);
         size_t I_skip = getFactorEntryDiff( I, i );
@@ -136,7 +126,13 @@ Prob Gibbs::getVarDist( size_t i ) {
         }
     }
 
-    return i_given_MB.normalized();
+    if( i_given_MB.sum() == 0.0 )
+        // If no state of i is allowed, use uniform distribution
+        // FIXME is that indeed the right thing to do?
+        i_given_MB = Prob( i_states );
+    else
+        i_given_MB.normalize();
+    return i_given_MB;
 }
 
 
@@ -163,12 +159,12 @@ void Gibbs::init() {
 
 double Gibbs::run() {
     if( props.verbose >= 1 )
-        cout << "Starting " << identify() << "...";
+        cerr << "Starting " << identify() << "...";
     if( props.verbose >= 3 )
-        cout << endl;
+        cerr << endl;
 
     double tic = toc();
-    
+
     randomizeState();
 
     for( size_t iter = 0; iter < props.iters; iter++ ) {
@@ -183,20 +179,20 @@ double Gibbs::run() {
             cerr << "counts for variable " << var(i) << ": " << Prob( _var_counts[i].begin(), _var_counts[i].end() ) << endl;
         }
     }
-    
+
     if( props.verbose >= 3 )
-        cout << "Gibbs::run:  ran " << props.iters << " passes (" << toc() - tic << " clocks)." << endl;
+        cerr << Name << "::run:  ran " << props.iters << " passes (" << toc() - tic << " clocks)." << endl;
 
     return 0.0;
 }
 
 
-inline Factor Gibbs::beliefV( size_t i ) const {
+Factor Gibbs::beliefV( size_t i ) const {
     return Factor( var(i), _var_counts[i].begin() ).normalized();
 }
 
 
-inline Factor Gibbs::beliefF( size_t I ) const {
+Factor Gibbs::beliefF( size_t I ) const {
     return Factor( factor(I).vars(), _factor_counts[I].begin() ).normalized();
 }
 
@@ -208,9 +204,9 @@ Factor Gibbs::belief( const Var &n ) const {
 
 vector<Factor> Gibbs::beliefs() const {
     vector<Factor> result;
-    for( size_t i = 0; i < nrVars(); i++ )
+    for( size_t i = 0; i < nrVars(); ++i )
         result.push_back( beliefV(i) );
-    for( size_t I = 0; I < nrFactors(); I++ )
+    for( size_t I = 0; I < nrFactors(); ++I )
         result.push_back( beliefF(I) );
     return result;
 }
@@ -224,7 +220,7 @@ Factor Gibbs::belief( const VarSet &ns ) const {
         for( I = 0; I < nrFactors(); I++ )
             if( factor(I).vars() >> ns )
                 break;
-        assert( I != nrFactors() );
+        DAI_ASSERT( I != nrFactors() );
         return beliefF(I).marginal(ns);
     }
 }