Documented all exceptions and did some general cleanups
[libdai.git] / src / factorgraph.cpp
index 57614c2..41088f5 100644 (file)
@@ -1,26 +1,16 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
-    Radboud University Nijmegen, The Netherlands /
-    Max Planck Institute for Biological Cybernetics, Germany
-
-    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) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ */
 
 
 #include <iostream>
+#include <iomanip>
 #include <iterator>
 #include <map>
 #include <set>
@@ -31,6 +21,7 @@
 #include <dai/factorgraph.h>
 #include <dai/util.h>
 #include <dai/exceptions.h>
+#include <boost/lexical_cast.hpp>
 
 
 namespace dai {
@@ -63,10 +54,10 @@ FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _backup() {
 void FactorGraph::constructGraph( size_t nrEdges ) {
     // create a mapping for indices
     hash_map<size_t, size_t> hashmap;
-    
+
     for( size_t i = 0; i < vars().size(); i++ )
         hashmap[var(i).label()] = i;
-    
+
     // create edge list
     vector<Edge> edges;
     edges.reserve( nrEdges );
@@ -82,7 +73,7 @@ void FactorGraph::constructGraph( size_t nrEdges ) {
 
 
 /// Writes a FactorGraph to an output stream
-ostream& operator << (ostream& os, const FactorGraph& fg) {
+std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) {
     os << fg.nrFactors() << endl;
 
     for( size_t I = 0; I < fg.nrFactors(); I++ ) {
@@ -96,15 +87,12 @@ ostream& operator << (ostream& os, const FactorGraph& fg) {
         os << endl;
         size_t nr_nonzeros = 0;
         for( size_t k = 0; k < fg.factor(I).states(); k++ )
-            if( fg.factor(I)[k] != 0.0 )
+            if( fg.factor(I)[k] != (Real)0 )
                 nr_nonzeros++;
         os << nr_nonzeros << endl;
         for( size_t k = 0; k < fg.factor(I).states(); k++ )
-            if( fg.factor(I)[k] != 0.0 ) {
-                char buf[20];
-                sprintf(buf,"%18.14g", fg.factor(I)[k]);
-                os << k << " " << buf << endl;
-            }
+            if( fg.factor(I)[k] != (Real)0 )
+                os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
     }
 
     return(os);
@@ -112,24 +100,24 @@ ostream& operator << (ostream& os, const FactorGraph& fg) {
 
 
 /// Reads a FactorGraph from an input stream
-istream& operator >> (istream& is, FactorGraph& fg) {
+std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
     long verbose = 0;
 
     vector<Factor> facs;
     size_t nr_Factors;
     string line;
-    
+
     while( (is.peek()) == '#' )
         getline(is,line);
     is >> nr_Factors;
     if( is.fail() )
-        DAI_THROW(INVALID_FACTORGRAPH_FILE);
+        DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
     if( verbose >= 2 )
         cerr << "Reading " << nr_Factors << " factors..." << endl;
 
     getline (is,line);
     if( is.fail() )
-        DAI_THROW(INVALID_FACTORGRAPH_FILE);
+        DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Expecting empty line");
 
     map<long,size_t> vardims;
     for( size_t I = 0; I < nr_Factors; I++ ) {
@@ -165,43 +153,33 @@ istream& operator >> (istream& is, FactorGraph& fg) {
             cerr << "  dimensions: " << dims << endl;
 
         // add the Factor
-        VarSet I_vars;
+        vector<Var> Ivars;
+        Ivars.reserve( nr_members );
         for( size_t mi = 0; mi < nr_members; mi++ ) {
             map<long,size_t>::iterator vdi = vardims.find( labels[mi] );
             if( vdi != vardims.end() ) {
                 // check whether dimensions are consistent
                 if( vdi->second != dims[mi] )
-                    DAI_THROW(INVALID_FACTORGRAPH_FILE);
+                    DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Variable with label " + boost::lexical_cast<string>(labels[mi]) + " has inconsistent dimensions.");
             } else
                 vardims[labels[mi]] = dims[mi];
-            I_vars |= Var(labels[mi], dims[mi]);
+            Ivars.push_back( Var(labels[mi], dims[mi]) );
         }
-        facs.push_back( Factor( I_vars, 0.0 ) );
-        
-        // calculate permutation sigma (internally, members are sorted)
-        vector<size_t> sigma(nr_members,0);
-        VarSet::iterator j = I_vars.begin();
-        for( size_t mi = 0; mi < nr_members; mi++,j++ ) {
-            long search_for = j->label();
-            vector<long>::iterator j_loc = find(labels.begin(),labels.end(),search_for);
-            sigma[mi] = j_loc - labels.begin();
-        }
-        if( verbose >= 3 )
-            cerr << "  sigma: " << sigma << endl;
+        facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
+
+        // calculate permutation object
+        Permute permindex( Ivars );
 
-        // calculate multindices
-        Permute permindex( dims, sigma );
-        
         // read values
         size_t nr_nonzeros;
         while( (is.peek()) == '#' )
             getline(is,line);
         is >> nr_nonzeros;
-        if( verbose >= 3 ) 
+        if( verbose >= 3 )
             cerr << "  nonzeroes: " << nr_nonzeros << endl;
         for( size_t k = 0; k < nr_nonzeros; k++ ) {
             size_t li;
-            double val;
+            Real val;
             while( (is.peek()) == '#' )
                 getline(is,line);
             is >> li;
@@ -209,9 +187,8 @@ istream& operator >> (istream& is, FactorGraph& fg) {
                 getline(is,line);
             is >> val;
 
-            // store value, but permute indices first according
-            // to internal representation
-            facs.back()[permindex.convert_linear_index( li  )] = val;
+            // store value, but permute indices first according to internal representation
+            facs.back()[permindex.convertLinearIndex( li )] = val;
         }
     }
 
@@ -224,12 +201,12 @@ istream& operator >> (istream& is, FactorGraph& fg) {
 }
 
 
-VarSet FactorGraph::delta( unsigned i ) const {
+VarSet FactorGraph::delta( size_t i ) const {
     return( Delta(i) / var(i) );
 }
 
 
-VarSet FactorGraph::Delta( unsigned i ) const {
+VarSet FactorGraph::Delta( size_t i ) const {
     // calculate Markov Blanket
     VarSet Del;
     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
@@ -242,17 +219,17 @@ VarSet FactorGraph::Delta( unsigned i ) const {
 
 VarSet FactorGraph::Delta( const VarSet &ns ) const {
     VarSet result;
-    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ ) 
+    for( VarSet::const_iterator n = ns.begin(); n != ns.end(); n++ )
         result |= Delta(findVar(*n));
     return result;
 }
 
 
-void FactorGraph::makeCavity( unsigned i, bool backup ) {
+void FactorGraph::makeCavity( size_t i, bool backup ) {
     // fills all Factors that include var(i) with ones
     map<size_t,Factor> newFacs;
     foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
-        newFacs[I] = Factor(factor(I).vars(), 1.0);
+        newFacs[I] = Factor( factor(I).vars(), (Real)1 );
     setFactors( newFacs, backup );
 }
 
@@ -264,18 +241,19 @@ void FactorGraph::ReadFromFile( const char *filename ) {
         infile >> *this;
         infile.close();
     } else
-        DAI_THROW(CANNOT_READ_FILE);
+        DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
 }
 
 
-void FactorGraph::WriteToFile( const char *filename ) const {
+void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
     ofstream outfile;
     outfile.open( filename );
     if( outfile.is_open() ) {
+        outfile.precision( precision );
         outfile << *this;
         outfile.close();
     } else
-        DAI_THROW(CANNOT_WRITE_FILE);
+        DAI_THROWE(CANNOT_WRITE_FILE,"Cannot write to file " + std::string(filename));
 }
 
 
@@ -296,13 +274,13 @@ void FactorGraph::printDot( std::ostream &os ) const {
 
 vector<VarSet> FactorGraph::Cliques() const {
     vector<VarSet> result;
-    
+
     for( size_t I = 0; I < nrFactors(); I++ ) {
         bool maximal = true;
         for( size_t J = 0; (J < nrFactors()) && maximal; J++ )
             if( (factor(J).vars() >> factor(I).vars()) && (factor(J).vars() != factor(I).vars()) )
                 maximal = false;
-        
+
         if( maximal )
             result.push_back( factor(I).vars() );
     }
@@ -311,30 +289,53 @@ vector<VarSet> FactorGraph::Cliques() const {
 }
 
 
-void FactorGraph::clamp( const Var & n, size_t i, bool backup ) {
-    assert( i <= n.states() );
+void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
+    DAI_ASSERT( x <= var(i).states() );
+    Factor mask( var(i), (Real)0 );
+    mask[x] = (Real)1;
 
-    // Multiply each factor that contains the variable with a delta function
+    map<size_t, Factor> newFacs;
+    foreach( const Neighbor &I, nbV(i) )
+        newFacs[I] = factor(I) * mask;
+    setFactors( newFacs, backup );
+
+    return;
+}
 
-    Factor delta_n_i(n,0.0);
-    delta_n_i[i] = 1.0;
+
+void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
+    Var n = var(i);
+    Factor mask_n( n, (Real)0 );
+
+    foreach( size_t i, is ) {
+        DAI_ASSERT( i <= n.states() );
+        mask_n[i] = (Real)1;
+    }
 
     map<size_t, Factor> newFacs;
-    // For all factors that contain n
-    for( size_t I = 0; I < nrFactors(); I++ ) 
-        if( factor(I).vars().contains( n ) )
-            // Multiply it with a delta function
-            newFacs[I] = factor(I) * delta_n_i;
+    foreach( const Neighbor &I, nbV(i) )
+        newFacs[I] = factor(I) * mask_n;
     setFactors( newFacs, backup );
+}
 
-    return;
+
+void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
+    size_t st = factor(I).states();
+    Factor newF( factor(I).vars(), (Real)0 );
+
+    foreach( size_t i, is ) {
+        DAI_ASSERT( i <= st );
+        newF[i] = factor(I)[i];
+    }
+
+    setFactor( I, newF, backup );
 }
 
 
 void FactorGraph::backupFactor( size_t I ) {
     map<size_t,Factor>::iterator it = _backup.find( I );
     if( it != _backup.end() )
-        DAI_THROW( MULTIPLE_UNDO );
+        DAI_THROW(MULTIPLE_UNDO);
     _backup[I] = factor(I);
 }
 
@@ -373,6 +374,7 @@ void FactorGraph::restoreFactors() {
     _backup.clear();
 }
 
+
 void FactorGraph::backupFactors( const std::set<size_t> & facs ) {
     for( std::set<size_t>::const_iterator fac = facs.begin(); fac != facs.end(); fac++ )
         backupFactor( *fac );
@@ -397,14 +399,15 @@ bool FactorGraph::isBinary() const {
 }
 
 
-FactorGraph FactorGraph::clamped( const Var & v_i, size_t state ) const {
-    Real zeroth_order = 1.0;
+FactorGraph FactorGraph::clamped( size_t i, size_t state ) const {
+    Var v = var( i );
+    Real zeroth_order = (Real)1;
     vector<Factor> clamped_facs;
     for( size_t I = 0; I < nrFactors(); I++ ) {
         VarSet v_I = factor(I).vars();
         Factor new_factor;
-        if( v_I.intersects( v_i ) )
-            new_factor = factor(I).slice( v_i, state );
+        if( v_I.intersects( v ) )
+            new_factor = factor(I).slice( v, state );
         else
             new_factor = factor(I);