New git HEAD version
[libdai.git] / src / factorgraph.cpp
index 416ecf8..3cff7eb 100644 (file)
@@ -1,11 +1,8 @@
 /*  This file is part of libDAI - http://www.libdai.org/
  *
 /*  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) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
- *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
  */
 
 
  */
 
 
@@ -112,7 +109,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
     is >> nr_Factors;
     if( is.fail() )
         DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
     is >> nr_Factors;
     if( is.fail() )
         DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Cannot read number of factors");
-    if( verbose >= 2 )
+    if( verbose >= 1 )
         cerr << "Reading " << nr_Factors << " factors..." << endl;
 
     getline (is,line);
         cerr << "Reading " << nr_Factors << " factors..." << endl;
 
     getline (is,line);
@@ -121,13 +118,13 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
 
     map<long,size_t> vardims;
     for( size_t I = 0; I < nr_Factors; I++ ) {
 
     map<long,size_t> vardims;
     for( size_t I = 0; I < nr_Factors; I++ ) {
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "Reading factor " << I << "..." << endl;
         size_t nr_members;
         while( (is.peek()) == '#' )
             getline(is,line);
         is >> nr_members;
             cerr << "Reading factor " << I << "..." << endl;
         size_t nr_members;
         while( (is.peek()) == '#' )
             getline(is,line);
         is >> nr_members;
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  nr_members: " << nr_members << endl;
 
         vector<long> labels;
             cerr << "  nr_members: " << nr_members << endl;
 
         vector<long> labels;
@@ -138,7 +135,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             is >> mi_label;
             labels.push_back(mi_label);
         }
             is >> mi_label;
             labels.push_back(mi_label);
         }
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  labels: " << labels << endl;
 
         vector<size_t> dims;
             cerr << "  labels: " << labels << endl;
 
         vector<size_t> dims;
@@ -149,7 +146,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             is >> mi_dim;
             dims.push_back(mi_dim);
         }
             is >> mi_dim;
             dims.push_back(mi_dim);
         }
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  dimensions: " << dims << endl;
 
         // add the Factor
             cerr << "  dimensions: " << dims << endl;
 
         // add the Factor
@@ -166,6 +163,8 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             Ivars.push_back( Var(labels[mi], dims[mi]) );
         }
         facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
             Ivars.push_back( Var(labels[mi], dims[mi]) );
         }
         facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
+        if( verbose >= 2 )
+            cerr << "  vardims: " << vardims << endl;
 
         // calculate permutation object
         Permute permindex( Ivars );
 
         // calculate permutation object
         Permute permindex( Ivars );
@@ -175,7 +174,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
         while( (is.peek()) == '#' )
             getline(is,line);
         is >> nr_nonzeros;
         while( (is.peek()) == '#' )
             getline(is,line);
         is >> nr_nonzeros;
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  nonzeroes: " << nr_nonzeros << endl;
         for( size_t k = 0; k < nr_nonzeros; k++ ) {
             size_t li;
             cerr << "  nonzeroes: " << nr_nonzeros << endl;
         for( size_t k = 0; k < nr_nonzeros; k++ ) {
             size_t li;
@@ -204,8 +203,8 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
 VarSet FactorGraph::Delta( size_t i ) const {
     // calculate Markov Blanket
     VarSet Del;
 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
-        foreach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        bforeach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
             Del |= var(j);
 
     return Del;
             Del |= var(j);
 
     return Del;
@@ -220,15 +219,34 @@ VarSet FactorGraph::Delta( const VarSet &ns ) const {
 }
 
 
 }
 
 
+SmallSet<size_t> FactorGraph::Deltai( size_t i ) const {
+    // calculate Markov Blanket
+    SmallSet<size_t> Del;
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        bforeach( const Neighbor &j, nbF(I) ) // for all neighboring variables j of I
+            Del |= j;
+
+    return Del;
+}
+
+
 void FactorGraph::makeCavity( size_t i, bool backup ) {
     // fills all Factors that include var(i) with ones
     map<size_t,Factor> newFacs;
 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
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
         newFacs[I] = Factor( factor(I).vars(), (Real)1 );
     setFactors( newFacs, backup );
 }
 
 
         newFacs[I] = Factor( factor(I).vars(), (Real)1 );
     setFactors( newFacs, backup );
 }
 
 
+void FactorGraph::makeRegionCavity( std::vector<size_t> facInds, bool backup ) {
+    map<size_t,Factor> newFacs;
+    for( size_t I = 0; I < facInds.size(); I++ )
+       newFacs[facInds[I]] = Factor(factor(facInds[I]).vars(), (Real)1);
+    setFactors( newFacs, backup );
+}
+
+
 void FactorGraph::ReadFromFile( const char *filename ) {
     ifstream infile;
     infile.open( filename );
 void FactorGraph::ReadFromFile( const char *filename ) {
     ifstream infile;
     infile.open( filename );
@@ -261,7 +279,7 @@ void FactorGraph::printDot( std::ostream &os ) const {
     for( size_t I = 0; I < nrFactors(); I++ )
         os << "\tf" << I << ";" << endl;
     for( size_t i = 0; i < nrVars(); i++ )
     for( size_t I = 0; I < nrFactors(); I++ )
         os << "\tf" << I << ";" << endl;
     for( size_t i = 0; i < nrVars(); i++ )
-        foreach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
+        bforeach( const Neighbor &I, nbV(i) )  // for all neighboring factors I of i
             os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
     os << "}" << endl;
 }
             os << "\tv" << var(i).label() << " -- f" << I << ";" << endl;
     os << "}" << endl;
 }
@@ -270,8 +288,8 @@ void FactorGraph::printDot( std::ostream &os ) const {
 GraphAL FactorGraph::MarkovGraph() const {
     GraphAL G( nrVars() );
     for( size_t i = 0; i < nrVars(); i++ )
 GraphAL FactorGraph::MarkovGraph() const {
     GraphAL G( nrVars() );
     for( size_t i = 0; i < nrVars(); i++ )
-        foreach( const Neighbor &I, nbV(i) )
-            foreach( const Neighbor &j, nbF(I) )
+        bforeach( const Neighbor &I, nbV(i) )
+            bforeach( const Neighbor &j, nbF(I) )
                 if( i < j )
                     G.addEdge( i, j, true );
     return G;
                 if( i < j )
                     G.addEdge( i, j, true );
     return G;
@@ -289,8 +307,8 @@ bool FactorGraph::isMaximal( size_t I ) const {
                     return false;
         return true;
     } else {
                     return false;
         return true;
     } else {
-        foreach( const Neighbor& i, nbF(I) ) {
-            foreach( const Neighbor& J, nbV(i) ) {
+        bforeach( const Neighbor& i, nbF(I) ) {
+            bforeach( const Neighbor& J, nbV(i) ) {
                 if( J != I )
                     if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
                         return false;
                 if( J != I )
                     if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
                         return false;
@@ -312,8 +330,8 @@ size_t FactorGraph::maximalFactor( size_t I ) const {
                     return maximalFactor( J );
         return I;
     } else {
                     return maximalFactor( J );
         return I;
     } else {
-        foreach( const Neighbor& i, nbF(I) ) {
-            foreach( const Neighbor& J, nbV(i) ) {
+        bforeach( const Neighbor& i, nbF(I) ) {
+            bforeach( const Neighbor& J, nbV(i) ) {
                 if( J != I )
                     if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
                         return maximalFactor( J );
                 if( J != I )
                     if( (factor(J).vars() >> I_vars) && (factor(J).vars().size() != I_size) )
                         return maximalFactor( J );
@@ -337,7 +355,7 @@ vector<VarSet> FactorGraph::maximalFactorDomains() const {
 }
 
 
 }
 
 
-Real FactorGraph::logScore( const std::vector<size_t>& statevec ) {
+Real FactorGraph::logScore( const std::vector<size_t>& statevec ) const {
     // Construct a State object that represents statevec
     // This decouples the representation of the joint state in statevec from the factor graph
     map<Var, size_t> statemap;
     // Construct a State object that represents statevec
     // This decouples the representation of the joint state in statevec from the factor graph
     map<Var, size_t> statemap;
@@ -349,7 +367,7 @@ Real FactorGraph::logScore( const std::vector<size_t>& statevec ) {
     // by summing the log factor entries of the factors that correspond to this joint configuration
     Real lS = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ )
     // by summing the log factor entries of the factors that correspond to this joint configuration
     Real lS = 0.0;
     for( size_t I = 0; I < nrFactors(); I++ )
-        lS += dai::log( factor(I)[S(factor(I).vars())] );
+        lS += dai::log( factor(I)[BigInt_size_t(S(factor(I).vars()))] );
     return lS;
 }
 
     return lS;
 }
 
@@ -360,7 +378,7 @@ void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
     mask.set( x, (Real)1 );
 
     map<size_t, Factor> newFacs;
     mask.set( x, (Real)1 );
 
     map<size_t, Factor> newFacs;
-    foreach( const Neighbor &I, nbV(i) )
+    bforeach( const Neighbor &I, nbV(i) )
         newFacs[I] = factor(I) * mask;
     setFactors( newFacs, backup );
 
         newFacs[I] = factor(I) * mask;
     setFactors( newFacs, backup );
 
@@ -372,13 +390,13 @@ void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
     Var n = var(i);
     Factor mask_n( n, (Real)0 );
 
     Var n = var(i);
     Factor mask_n( n, (Real)0 );
 
-    foreach( size_t i, is ) {
+    bforeach( size_t i, is ) {
         DAI_ASSERT( i <= n.states() );
         mask_n.set( i, (Real)1 );
     }
 
     map<size_t, Factor> newFacs;
         DAI_ASSERT( i <= n.states() );
         mask_n.set( i, (Real)1 );
     }
 
     map<size_t, Factor> newFacs;
-    foreach( const Neighbor &I, nbV(i) )
+    bforeach( const Neighbor &I, nbV(i) )
         newFacs[I] = factor(I) * mask_n;
     setFactors( newFacs, backup );
 }
         newFacs[I] = factor(I) * mask_n;
     setFactors( newFacs, backup );
 }
@@ -388,7 +406,7 @@ void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup )
     size_t st = factor(I).nrStates();
     Factor newF( factor(I).vars(), (Real)0 );
 
     size_t st = factor(I).nrStates();
     Factor newF( factor(I).vars(), (Real)0 );
 
-    foreach( size_t i, is ) {
+    bforeach( size_t i, is ) {
         DAI_ASSERT( i <= st );
         newF.set( i, factor(I)[i] );
     }
         DAI_ASSERT( i <= st );
         newF.set( i, factor(I)[i] );
     }