New git HEAD version
[libdai.git] / src / factorgraph.cpp
index 8813e6d..3cff7eb 100644 (file)
@@ -1,23 +1,9 @@
-/*  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/
+ *
+ *  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 <iostream>
@@ -41,7 +27,7 @@ namespace dai {
 using namespace std;
 
 
-FactorGraph::FactorGraph( const std::vector<Factor> &P ) : G(), _backup() {
+FactorGraph::FactorGraph( const std::vector<Factor> &P ) : _G(), _backup() {
     // add factors, obtain variables
     set<Var> varset;
     _factors.reserve( P.size() );
@@ -79,7 +65,7 @@ void FactorGraph::constructGraph( size_t nrEdges ) {
     }
 
     // create bipartite graph
-    G.construct( nrVars(), nrFactors(), edges.begin(), edges.end() );
+    _G.construct( nrVars(), nrFactors(), edges.begin(), edges.end() );
 }
 
 
@@ -97,12 +83,12 @@ std::ostream& operator<< ( std::ostream &os, const FactorGraph &fg ) {
             os << i->states() << " ";
         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 )
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
+            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 )
+        for( size_t k = 0; k < fg.factor(I).nrStates(); k++ )
+            if( fg.factor(I)[k] != (Real)0 )
                 os << k << " " << setw(os.precision()+4) << fg.factor(I)[k] << endl;
     }
 
@@ -123,22 +109,22 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
     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);
-    if( is.fail() )
-        DAI_THROW(INVALID_FACTORGRAPH_FILE);
+    if( is.fail() || line.size() > 0 )
+        DAI_THROWE(INVALID_FACTORGRAPH_FILE,"Expecting empty line");
 
     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;
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  nr_members: " << nr_members << endl;
 
         vector<long> labels;
@@ -149,7 +135,7 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             is >> mi_label;
             labels.push_back(mi_label);
         }
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             cerr << "  labels: " << labels << endl;
 
         vector<size_t> dims;
@@ -160,11 +146,12 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
             is >> mi_dim;
             dims.push_back(mi_dim);
         }
-        if( verbose >= 3 )
+        if( verbose >= 2 )
             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() ) {
@@ -173,34 +160,25 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
                     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]);
-        }
-        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();
+            Ivars.push_back( Var(labels[mi], dims[mi]) );
         }
-        if( verbose >= 3 )
-            cerr << "  sigma: " << sigma << endl;
+        facs.push_back( Factor( VarSet( Ivars.begin(), Ivars.end(), Ivars.size() ), (Real)0 ) );
+        if( verbose >= 2 )
+            cerr << "  vardims: " << vardims << endl;
 
-        // calculate multindices
-        Permute permindex( dims, sigma );
+        // calculate permutation object
+        Permute permindex( Ivars );
 
         // read values
         size_t 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;
-            double val;
+            Real val;
             while( (is.peek()) == '#' )
                 getline(is,line);
             is >> li;
@@ -208,9 +186,8 @@ std::istream& operator>> ( std::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().set( permindex.convertLinearIndex( li ), val );
         }
     }
 
@@ -223,16 +200,11 @@ std::istream& operator>> ( std::istream& is, FactorGraph &fg ) {
 }
 
 
-VarSet FactorGraph::delta( unsigned 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
-        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;
@@ -242,16 +214,35 @@ 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++ )
-        result |= Delta(findVar(*n));
+        result |= Delta( findVar(*n) );
     return result;
 }
 
 
-void FactorGraph::makeCavity( unsigned i, bool backup ) {
+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;
-    foreach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
-        newFacs[I] = Factor(factor(I).vars(), 1.0);
+    bforeach( const Neighbor &I, nbV(i) ) // for all neighboring factors I of i
+        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 );
 }
 
@@ -280,7 +271,7 @@ void FactorGraph::WriteToFile( const char *filename, size_t precision ) const {
 
 
 void FactorGraph::printDot( std::ostream &os ) const {
-    os << "graph G {" << endl;
+    os << "graph FactorGraph {" << endl;
     os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
     for( size_t i = 0; i < nrVars(); i++ )
         os << "\tv" << var(i).label() << ";" << endl;
@@ -288,36 +279,106 @@ 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++ )
-        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;
 }
 
 
-vector<VarSet> FactorGraph::Cliques() const {
-    vector<VarSet> result;
+GraphAL FactorGraph::MarkovGraph() const {
+    GraphAL G( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        bforeach( const Neighbor &I, nbV(i) )
+            bforeach( const Neighbor &j, nbF(I) )
+                if( i < j )
+                    G.addEdge( i, j, true );
+    return G;
+}
 
-    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() );
+bool FactorGraph::isMaximal( size_t I ) const {
+    const VarSet& I_vars = factor(I).vars();
+    size_t I_size = I_vars.size();
+
+    if( I_size == 0 ) {
+        for( size_t J = 0; J < nrFactors(); J++ ) 
+            if( J != I )
+                if( factor(J).vars().size() > 0 )
+                    return false;
+        return true;
+    } else {
+        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;
+            }
+        }
+        return true;
+    }
+}
+
+
+size_t FactorGraph::maximalFactor( size_t I ) const {
+    const VarSet& I_vars = factor(I).vars();
+    size_t I_size = I_vars.size();
+
+    if( I_size == 0 ) {
+        for( size_t J = 0; J < nrFactors(); J++ )
+            if( J != I )
+                if( factor(J).vars().size() > 0 )
+                    return maximalFactor( J );
+        return I;
+    } else {
+        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 );
+            }
+        }
+        return I;
     }
+}
+
 
+vector<VarSet> FactorGraph::maximalFactorDomains() const {
+    vector<VarSet> result;
+
+    for( size_t I = 0; I < nrFactors(); I++ )
+        if( isMaximal( I ) )
+            result.push_back( factor(I).vars() );
+
+    if( result.size() == 0 )
+        result.push_back( VarSet() );
     return result;
 }
 
 
+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;
+    for( size_t i = 0; i < statevec.size(); i++ )
+        statemap[var(i)] = statevec[i];
+    State S(statemap);
+
+    // Evaluate the log probability of the joint configuration in 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++ )
+        lS += dai::log( factor(I)[BigInt_size_t(S(factor(I).vars()))] );
+    return lS;
+}
+
+
 void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
-    assert( x <= var(i).states() );
-    Factor mask( var(i), 0.0 );
-    mask[x] = 1.0;
+    DAI_ASSERT( x <= var(i).states() );
+    Factor mask( var(i), (Real)0 );
+    mask.set( x, (Real)1 );
 
     map<size_t, Factor> newFacs;
-    foreach( const BipartiteGraph::Neighbor &I, nbV(i) )
+    bforeach( const Neighbor &I, nbV(i) )
         newFacs[I] = factor(I) * mask;
     setFactors( newFacs, backup );
 
@@ -327,27 +388,27 @@ void FactorGraph::clamp( size_t i, size_t x, bool backup ) {
 
 void FactorGraph::clampVar( size_t i, const vector<size_t> &is, bool backup ) {
     Var n = var(i);
-    Factor mask_n( n, 0.0 );
+    Factor mask_n( n, (Real)0 );
 
-    foreach( size_t i, is ) {
-        assert( i <= n.states() );
-        mask_n[i] = 1.0;
+    bforeach( size_t i, is ) {
+        DAI_ASSERT( i <= n.states() );
+        mask_n.set( i, (Real)1 );
     }
 
     map<size_t, Factor> newFacs;
-    foreach( const BipartiteGraph::Neighbor &I, nbV(i) )
+    bforeach( const Neighbor &I, nbV(i) )
         newFacs[I] = factor(I) * mask_n;
     setFactors( newFacs, backup );
 }
 
 
 void FactorGraph::clampFactor( size_t I, const vector<size_t> &is, bool backup ) {
-    size_t st = factor(I).states();
-    Factor newF( factor(I).vars(), 0.0 );
+    size_t st = factor(I).nrStates();
+    Factor newF( factor(I).vars(), (Real)0 );
 
-    foreach( size_t i, is ) {
-        assert( i <= st );
-        newF[i] = factor(I)[i];
+    bforeach( size_t i, is ) {
+        DAI_ASSERT( i <= st );
+        newF.set( i, factor(I)[i] );
     }
 
     setFactor( I, newF, backup );
@@ -367,7 +428,8 @@ void FactorGraph::restoreFactor( size_t I ) {
     if( it != _backup.end() ) {
         setFactor(I, it->second);
         _backup.erase(it);
-    }
+    } else
+        DAI_THROW(OBJECT_NOT_FOUND);
 }
 
 
@@ -421,9 +483,11 @@ bool FactorGraph::isBinary() const {
 }
 
 
-FactorGraph FactorGraph::clamped( const Var &v, 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;
+    clamped_facs.push_back( createFactorDelta( v, state ) );
     for( size_t I = 0; I < nrFactors(); I++ ) {
         VarSet v_I = factor(I).vars();
         Factor new_factor;