Added utils/uai2fg, ExactInf::findMaximum(), and fixed two bugs
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 20 Apr 2010 11:01:07 +0000 (13:01 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 20 Apr 2010 11:01:07 +0000 (13:01 +0200)
ChangeLog
Makefile
examples/example.cpp
include/dai/exactinf.h
src/exactinf.cpp
src/jtree.cpp
utils/uai2fg.cpp [new file with mode: 0644]

index 9fbb776..52effc4 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -6,12 +6,16 @@ git master HEAD
 * [Stefano Pellegrini] Fixed bug in JTree::findMaximum()
 * Fixed some bugs in the MatLab interface build system
 * Fixed a bug in utils/fginfo.cpp
+* Added utils/uai2fg, a utility to convert the UAI inference competition file
+  formats to the libDAI factor graph file format
 * utils/createfg now has a new factor type ('ISINGUNIFORM') and the old
   'ISING' factor type has been renamed to 'ISINGGAUSS'
 * tests/testdai option "marginals" now has five possible values: NONE
   outputs no marginals, VAR only variable marginals, FAC only factor
   marginals, VARFAC both types of marginals, and ALL outputs all
   marginals calculated by the algorithm.
+* Improved exactinf.h/cpp:
+  - Added ExactInf::findMaximum()
 * Improved treeep.h/cpp:
   - changed TreeEP::construct( const RootedTree& ) into
     TreeEP::construct( const FactorGraph&, const RootedTree& )
@@ -22,6 +26,7 @@ git master HEAD
   - changed JTree::GenerateJT( const std::vector<VarSet> & )
     into JTree::GenerateJT( const FactorGraph &, const std::vector<VarSet> & )
   - Now also supports disconnected factor graphs
+  - Fixed a bug in JTree::findMaximum() and simplified the code
 * Improved regiongraph.h/cpp:
   - Made (previously public) members RegionGraph::G, RegionGraph::ORs,
     RegionGraph::IRs and RegionGraph::fac2OR protected.
index 98096d6..7bd6ed4 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -144,7 +144,7 @@ unittests : tests/unit/var_test$(EE) tests/unit/smallset_test$(EE) tests/unit/va
 
 tests : tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE) $(unittests)
 
-utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE)
+utils : utils/createfg$(EE) utils/fg2dot$(EE) utils/fginfo$(EE) utils/uai2fg$(EE)
 
 lib: $(LIB)/libdai$(LE)
 
@@ -233,6 +233,9 @@ utils/fg2dot$(EE) : utils/fg2dot.cpp $(HEADERS) $(LIB)/libdai$(LE)
 utils/fginfo$(EE) : utils/fginfo.cpp $(HEADERS) $(LIB)/libdai$(LE)
        $(CC) $(CCO)$@ $< $(LIBS)
 
+utils/uai2fg$(EE) : utils/uai2fg.cpp $(HEADERS) $(LIB)/libdai$(LE)
+       $(CC) $(CCO)$@ $< $(LIBS)
+
 
 # LIBRARY
 ##########
@@ -294,7 +297,7 @@ clean :
        -rm tests/testdai$(EE) tests/testem/testem$(EE) tests/testbbp$(EE)
        -rm tests/unit/var_test$(EE) tests/unit/smallset_test$(EE) tests/unit/varset_test$(EE) tests/unit/graph_test$(EE) tests/unit/bipgraph_test$(EE) tests/unit/weightedgraph_test$(EE) tests/unit/enum_test$(EE) tests/unit/util_test$(EE) tests/unit/exceptions_test$(EE) tests/unit/properties_test$(EE) tests/unit/index_test$(EE) tests/unit/prob_test$(EE) tests/unit/factor_test$(EE) tests/unit/factorgraph_test$(EE) tests/unit/clustergraph_test$(EE) tests/unit/regiongraph_test$(EE) tests/unit/daialg_test$(EE) tests/unit/alldai_test$(EE)
        -rm factorgraph_test.fg alldai_test.aliases
-       -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE)
+       -rm utils/fg2dot$(EE) utils/createfg$(EE) utils/fginfo$(EE) utils/uai2fg$(EE)
        -rm -R doc
        -rm -R lib
 else
@@ -307,24 +310,22 @@ clean :
        -del examples\*$(EE).manifest
        -del examples\*.ilk
        -del examples\*.pdb
-       -del tests\testdai$(EE)
-       -del tests\testbbp$(EE)
-       -del tests\testdai$(EE).manifest
-       -del tests\testbbp$(EE).manifest
-       -del tests\testem\testem$(EE)
-       -del tests\testem\testem$(EE).manifest
+       -del tests\*$(EE)
+       -del tests\*$(EE).manifest
        -del tests\*.pdb
        -del tests\*.ilk
+       -del tests\testem\*$(EE)
+       -del tests\testem\*$(EE).manifest
        -del tests\testem\*.pdb
        -del tests\testem\*.ilk
        -del utils\*$(EE)
        -del utils\*$(EE).manifest
        -del utils\*.pdb
        -del utils\*.ilk
-       -del tests\unit\*_test.ilk
-       -del tests\unit\*_test.pdb
        -del tests\unit\*_test$(EE)
        -del tests\unit\*_test$(EE).manifest
+       -del tests\unit\*_test.pdb
+       -del tests\unit\*_test.ilk
        -del factorgraph_test.fg
        -del alldai_test.aliases
        -del $(LIB)\libdai$(LE)
index 68e509f..390c4da 100644 (file)
@@ -10,6 +10,7 @@
 
 
 #include <iostream>
+#include <map>
 #include <dai/alldai.h>  // Include main libDAI header file
 
 
@@ -17,6 +18,25 @@ using namespace dai;
 using namespace std;
 
 
+// Evaluates the log probability of the joint configuration vstate of the variables in factor graph fg
+Real evalState( const FactorGraph& fg, const std::vector<size_t> vstate ) {
+    // First, construct a map from Var objects to their states
+    // This decouples the representation of the joint state in vstate
+    // from the factor graph
+    map<Var, size_t> state;
+    for( size_t i = 0; i < vstate.size(); i++ )
+        state[fg.var(i)] = vstate[i];
+
+    // Evaluate the log probability of the joint configuration in state
+    // by summing the log factor entries of the factors in factor graph fg
+    // that correspond to the joint configuration
+    Real logZ = 0.0;
+    for( size_t I = 0; I < fg.nrFactors(); I++ )
+        logZ += dai::log( fg.factor(I)[calcLinearState( fg.factor(I).vars(), state )] );
+    return logZ;
+}
+
+
 int main( int argc, char *argv[] ) {
     if ( argc != 2 ) {
         cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
@@ -131,12 +151,12 @@ int main( int argc, char *argv[] ) {
             cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
 
         // Report exact MAP joint state
-        cout << "Exact MAP state:" << endl;
+        cout << "Exact MAP state (log probability = " << evalState( fg, jtmapstate ) << "):" << endl;
         for( size_t i = 0; i < jtmapstate.size(); i++ )
             cout << fg.var(i) << ": " << jtmapstate[i] << endl;
 
         // Report max-product MAP joint state
-        cout << "Approximate (max-product) MAP state:" << endl;
+        cout << "Approximate (max-product) MAP state (log probability = " << evalState( fg, mpstate ) << "):" << endl;
         for( size_t i = 0; i < mpstate.size(); i++ )
             cout << fg.var(i) << ": " << mpstate[i] << endl;
     }
index 86de319..cf9c411 100644 (file)
@@ -87,14 +87,18 @@ class ExactInf : public DAIAlgFG {
         virtual std::string printProperties() const;
     //@}
 
-    /// \name Additional interface specific for JTree
+    /// \name Additional interface specific for ExactInf
     //@{
         /// Calculates marginal probability distribution for variables \a vs
         /** \note The complexity of this calculation is exponential in the number of variables.
          */
         Factor calcMarginal( const VarSet &vs ) const;
-    //@}
 
+        /// Calculates the joint state of all variables that has maximum probability
+        /** \note The complexity of this calculation is exponential in the number of variables.
+         */
+        std::vector<std::size_t> findMaximum() const;
+    //@}
 
     private:
         /// Helper function for constructors
index 980c9c5..97a1a9c 100644 (file)
@@ -97,6 +97,25 @@ Factor ExactInf::calcMarginal( const VarSet &vs ) const {
     return P.marginal( vs, true );
 }
 
+        
+std::vector<std::size_t> ExactInf::findMaximum() const {
+    Factor P;
+    for( size_t I = 0; I < nrFactors(); I++ )
+        P *= factor(I);
+    size_t linearState = P.p().argmax().first;
+
+    // convert to state
+    map<Var, size_t> state = calcState( P.vars(), linearState );
+
+    // convert to desired output data structure
+    vector<size_t> mapState;
+    mapState.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        mapState.push_back( state[var(i)] );
+
+    return mapState;
+}
+
 
 vector<Factor> ExactInf::beliefs() const {
     vector<Factor> result = _beliefsV;
index f106e90..44bdb61 100644 (file)
@@ -117,25 +117,25 @@ void JTree::construct( const FactorGraph &fg, const std::vector<VarSet> &cl, boo
     // Construct a weighted graph (each edge is weighted with the cardinality
     // of the intersection of the nodes, where the nodes are the elements of cl).
     WeightedGraph<int> JuncGraph;
-    std::vector<bool> connected( cl.size(), false );
+    // Start by connecting all clusters with cluster zero, and weight zero,
+    // in order to get a connected weighted graph
+    for( size_t i = 1; i < cl.size(); i++ )
+        JuncGraph[UEdge(i,0)] = 0;
     for( size_t i = 0; i < cl.size(); i++ ) {
         for( size_t j = i + 1; j < cl.size(); j++ ) {
             size_t w = (cl[i] & cl[j]).size();
-            if( w ) {
+            if( w )
                 JuncGraph[UEdge(i,j)] = w;
-                connected[i] = true;
-                connected[j] = true;
-            }
         }
     }
-    // for clusters that have no overlap with other clusters,
-    // connect them with the zeroth cluster
-    for( size_t i = 1; i < cl.size(); i++ )
-        if( !connected[i] )
-            JuncGraph[UEdge(i,0)] = 0;
+    if( props.verbose >= 3 )
+        cerr << "Weightedgraph: " << JuncGraph << endl;
 
     // Construct maximal spanning tree using Prim's algorithm
     RTree = MaxSpanningTree( JuncGraph, true );
+    if( props.verbose >= 3 )
+        cerr << "Spanning tree: " << RTree << endl;
+    DAI_DEBASSERT( RTree.size() == cl.size() - 1 );
 
     // Construct corresponding region graph
 
@@ -223,9 +223,12 @@ Factor JTree::belief( const VarSet &vs ) const {
     for( beta = Qb.begin(); beta != Qb.end(); beta++ )
         if( beta->vars() >> vs )
             break;
-    if( beta != Qb.end() )
-        return( beta->marginal(vs) );
-    else {
+    if( beta != Qb.end() ) {
+        if( props.inference == Properties::InfType::SUMPROD )
+            return( beta->marginal(vs) );
+        else
+            return( beta->maxMarginal(vs) );
+    } else {
         vector<Factor>::const_iterator alpha;
         for( alpha = Qa.begin(); alpha != Qa.end(); alpha++ )
             if( alpha->vars() >> vs )
@@ -233,8 +236,12 @@ Factor JTree::belief( const VarSet &vs ) const {
         if( alpha == Qa.end() ) {
             DAI_THROW(BELIEF_NOT_AVAILABLE);
             return Factor();
-        } else
-            return( alpha->marginal(vs) );
+        } else {
+            if( props.inference == Properties::InfType::SUMPROD )
+                return( alpha->marginal(vs) );
+            else
+                return( alpha->maxMarginal(vs) );
+        }
     }
 }
 
@@ -454,16 +461,22 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
     for( beta = Qb.begin(); beta != Qb.end(); beta++ )
         if( beta->vars() >> vs )
             break;
-    if( beta != Qb.end() )
-        return( beta->marginal(vs) );
-    else {
+    if( beta != Qb.end() ) {
+        if( props.inference == Properties::InfType::SUMPROD )
+            return( beta->marginal(vs) );
+        else
+            return( beta->maxMarginal(vs) );
+    } else {
         vector<Factor>::const_iterator alpha;
         for( alpha = Qa.begin(); alpha != Qa.end(); alpha++ )
             if( alpha->vars() >> vs )
                 break;
-        if( alpha != Qa.end() )
-            return( alpha->marginal(vs) );
-        else {
+        if( alpha != Qa.end() ) {
+            if( props.inference == Properties::InfType::SUMPROD )
+                return( alpha->marginal(vs) );
+            else
+                return( alpha->maxMarginal(vs) );
+        } else {
             // Find subtree to do efficient inference
             RootedTree T;
             size_t Tsize = findEfficientTree( vs, T );
@@ -509,7 +522,11 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
                             Qa[T[i].second] *= piet;
                         }
 
-                    Factor new_Qb = Qa[T[i].second].marginal( IR( b[i] ), false );
+                    Factor new_Qb;
+                    if( props.inference == Properties::InfType::SUMPROD )
+                        new_Qb = Qa[T[i].second].marginal( IR( b[i] ), false );
+                    else
+                        new_Qb = Qa[T[i].second].maxMarginal( IR( b[i] ), false );
                     logZ += log(new_Qb.normalize());
                     Qa[T[i].first] *= new_Qb / Qb[b[i]];
                     Qb[b[i]] = new_Qb;
@@ -518,7 +535,10 @@ Factor JTree::calcMarginal( const VarSet& vs ) {
 
                 Factor piet( vsrem, 0.0 );
                 piet.set( s, exp(logZ) );
-                Pvs += piet * Qa[T[0].first].marginal( vs / vsrem, false );      // OPTIMIZE ME
+                if( props.inference == Properties::InfType::SUMPROD )
+                    Pvs += piet * Qa[T[0].first].marginal( vs / vsrem, false );      // OPTIMIZE ME
+                else
+                    Pvs += piet * Qa[T[0].first].maxMarginal( vs / vsrem, false );      // OPTIMIZE ME
 
                 // Restore clamped beliefs
                 for( map<size_t,Factor>::const_iterator alpha = Qa_old.begin(); alpha != Qa_old.end(); alpha++ )
@@ -562,80 +582,38 @@ std::pair<size_t,double> boundTreewidth( const FactorGraph &fg, greedyVariableEl
 
 
 std::vector<size_t> JTree::findMaximum() const {
-    vector<size_t> maximum( nrVars() );
-    vector<bool> visitedVars( nrVars(), false );
-    vector<bool> visitedFactors( nrFactors(), false );
-    stack<size_t> scheduledFactors;
-    for( size_t i = 0; i < nrVars(); ++i ) {
-        if( visitedVars[i] )
-            continue;
-        visitedVars[i] = true;
-
-        // Maximise with respect to variable i
-        Prob prod = beliefV(i).p();
-        maximum[i] = prod.argmax().first;
-
-        foreach( const Neighbor &I, nbV(i) )
-            if( !visitedFactors[I] )
-                scheduledFactors.push(I);
-
-        while( !scheduledFactors.empty() ){
-            size_t I = scheduledFactors.top();
-            scheduledFactors.pop();
-            if( visitedFactors[I] )
-                continue;
-            visitedFactors[I] = true;
-
-            // Evaluate if some neighboring variables still need to be fixed; if not, we're done
-            bool allDetermined = true;
-            foreach( const Neighbor &j, nbF(I) )
-                if( !visitedVars[j.node] ) {
-                    allDetermined = false;
-                    break;
-                }
-            if( allDetermined )
-                continue;
-
-            // Calculate product of incoming messages on factor I
-            Prob prod2 = beliefF(I).p();
-
-            // The allowed configuration is restrained according to the variables assigned so far:
-            // pick the argmax amongst the allowed states
-            Real maxProb = -numeric_limits<Real>::max();
-            State maxState( factor(I).vars() );
-            for( State s( factor(I).vars() ); s.valid(); ++s ){
-                // First, calculate whether this state is consistent with variables that
-                // have been assigned already
-                bool allowedState = true;
-                foreach( const Neighbor &j, nbF(I) )
-                    if( visitedVars[j.node] && maximum[j.node] != s(var(j.node)) ) {
-                        allowedState = false;
-                        break;
-                    }
-                // If it is consistent, check if its probability is larger than what we have seen so far
-                if( allowedState && prod2[s] > maxProb ) {
-                    maxState = s;
-                    maxProb = prod2[s];
-                }
+#ifdef DAI_DEBUG
+    // check consistency of variables and factors
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        size_t linearState = beliefF(I).p().argmax().first;
+        map<Var,size_t> state = calcState( factor(I).vars(), linearState );
+        foreach( const Neighbor& i, nbF(I) ) {
+            if( state[var(i)] != beliefV(i).p().argmax().first ) {
+                cerr << "ERROR: inconsistency between MAP beliefs: factor belief " << beliefF(I) << ": " << state << " is inconsistent with " << beliefV(i) << endl;
+                DAI_THROW(RUNTIME_ERROR);
             }
-
-            // Decode the argmax
-            foreach( const Neighbor &j, nbF(I) ) {
-                if( visitedVars[j.node] ) {
-                    // We have already visited j earlier - hopefully our state is consistent
-                    if( maximum[j.node] != maxState(var(j.node)) && props.verbose >= 1 )
-                        cerr << "JTree::findMaximum - warning: maximum not consistent due to loops." << endl;
-                } else {
-                    // We found a consistent state for variable j
-                    visitedVars[j.node] = true;
-                    maximum[j.node] = maxState( var(j.node) );
-                    foreach( const Neighbor &J, nbV(j) )
-                        if( !visitedFactors[J] )
-                            scheduledFactors.push(J);
+        }
+    }
+    // check consistency of 
+    for( size_t alpha = 0; alpha < nrORs(); alpha++ ) {
+        size_t linearStateQa = Qa[alpha].p().argmax().first;
+        map<Var,size_t> stateQa = calcState( Qa[alpha].vars(), linearStateQa );
+        foreach( const Neighbor& beta, nbOR(alpha) ) {
+            size_t linearStateQb = Qb[beta].p().argmax().first;
+            map<Var,size_t> stateQb = calcState( Qb[beta].vars(), linearStateQb );
+            for( map<Var,size_t>::const_iterator it = stateQb.begin(); it != stateQb.end(); it++ )
+                if( stateQa[it->first] != it->second ) {
+                    cerr << "ERROR: inconsistency between MAP beliefs: OR belief " << Qa[alpha] << " (with MAP " << stateQa << ") is inconsistent with IR belief " << Qb[beta] << " (with MAP " << stateQb << ")" << endl;
+                    DAI_THROW(RUNTIME_ERROR);
                 }
-            }
         }
     }
+#endif
+
+    vector<size_t> maximum;
+    maximum.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        maximum.push_back( beliefV(i).p().argmax().first );
     return maximum;
 }
 
diff --git a/utils/uai2fg.cpp b/utils/uai2fg.cpp
new file mode 100644 (file)
index 0000000..6c2347b
--- /dev/null
@@ -0,0 +1,271 @@
+/*  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-2010  Joris Mooij  [joris dot mooij at libdai dot org]
+ */
+
+
+#include <iostream>
+#include <fstream>
+#include <dai/alldai.h>
+#include <dai/util.h>
+#include <dai/index.h>
+#include <dai/jtree.h>
+
+
+using namespace std;
+using namespace dai;
+
+
+map<size_t, size_t> ReadUAIEvidenceFile( char *filename ) {
+    map<size_t, size_t> evid;
+
+    ifstream is;
+    is.open( filename );
+    if( is.is_open() ) {
+        size_t nr_evid;
+        is >> nr_evid;
+        if( is.fail() )
+            DAI_THROWE(INVALID_EVIDENCE_FILE,"Cannot read number of observed variables");
+
+        for( size_t i = 0; i < nr_evid; i++ ) {
+            size_t label, val;
+            is >> label;
+            is >> val;
+            evid[label] = val;
+        }
+
+        is.close();
+    } else
+        DAI_THROWE(CANNOT_READ_FILE,"Cannot read from file " + std::string(filename));
+
+    return evid;
+}
+
+
+pair<vector<Var>, vector<Factor> > ReadUAIFGFile( const char *filename, long verbose ) {
+    vector<Factor> factors;
+    vector<Var> vars;
+
+    ifstream is;
+    is.open( filename );
+    if( is.is_open() ) {
+        try { // FIXME: no exceptions are generated by iostreams by default
+            size_t nr_f, nr_v;
+            string line;
+            
+            getline(is,line);
+            if( line != "BAYES" && line != "MARKOV" && line != "BAYES\r" && line != "MARKOV\r" )
+                DAI_THROW(INVALID_FACTORGRAPH_FILE);
+            if( verbose >= 2 )
+                cout << "Reading " << line << " network..." << endl;
+
+            is >> nr_v;
+            if( is.fail() )
+                DAI_THROW(INVALID_FACTORGRAPH_FILE);
+            if( verbose >= 2 )
+                cout << "Reading " << nr_v << " variables..." << endl;
+
+            for( size_t i = 0; i < nr_v; i++ ) {
+                size_t dim;
+                is >> dim;
+                if( is.fail() )
+                    DAI_THROW(INVALID_FACTORGRAPH_FILE);
+                vars.push_back( Var( i, dim ) );
+            }
+
+            is >> nr_f;
+            if( is.fail() )
+                DAI_THROW(INVALID_FACTORGRAPH_FILE);
+            if( verbose >= 2 )
+                cout << "Reading " << nr_f << " factors..." << endl;
+
+            vector<vector<long> > labels;
+            for( size_t I = 0; I < nr_f; I++ ) {
+                if( verbose >= 3 )
+                    cout << "Reading factor " << I << "..." << endl;
+                size_t I_nrvars;
+                is >> I_nrvars;
+                if( verbose >= 3 )
+                    cout << "  I_nrvars: " << I_nrvars << endl;
+
+                vector<long> I_labels;
+                vector<size_t> I_dims;
+                VarSet I_vars;
+                for( size_t _i = 0; _i < I_nrvars; _i++ ) {
+                    long label;
+                    is >> label;
+                    I_labels.push_back(label);
+                    I_dims.push_back(vars[label].states());
+                    I_vars |= vars[label];
+                }
+                if( verbose >= 3 )
+                    cout << "  labels: " << I_labels << ", dimensions " << I_dims << endl;
+
+                // add the Factor and the labels
+                factors.push_back(Factor(I_vars,0.0));
+                labels.push_back(I_labels);
+                
+            }
+
+            for( size_t I = 0; I < nr_f; I++ ) {
+                if( verbose >= 3 )
+                    cout << "Reading factor " << I << "..." << endl;
+
+                reverse(labels[I].begin(), labels[I].end());    // last label is least significant
+
+                size_t I_nrvars = factors[I].vars().size();
+                vector<size_t> I_dims;
+                for( size_t _i = 0; _i < I_nrvars; _i++ )
+                    I_dims.push_back(vars[labels[I][_i]].states());
+                if( verbose >= 3 )
+                    cout << "  labels: " << labels[I] << ", dimensions " << I_dims << endl;
+
+                // calculate permutation sigma (internally, members are sorted)
+                vector<size_t> sigma(I_nrvars,0);
+                VarSet::const_iterator j = factors[I].vars().begin();
+                for( size_t mi = 0; mi < I_nrvars; mi++,j++ ) {
+                    long search_for = j->label();
+                    vector<long>::iterator j_loc = find(labels[I].begin(),labels[I].end(),search_for);
+                    sigma[mi] = j_loc - labels[I].begin();
+                }
+                if( verbose >= 3 )
+                    cout << "  sigma: " << sigma << endl;
+
+                // calculate multindices
+                Permute permindex( I_dims, sigma );
+                
+                // read values
+                size_t nr_nonzeros;
+                is >> nr_nonzeros;
+                if( verbose >= 3 ) 
+                    cout << "  nonzeroes: " << nr_nonzeros << endl;
+                assert( nr_nonzeros == factors[I].states() );
+                for( size_t li = 0; li < nr_nonzeros; li++ ) {
+                    double val;
+                    is >> val;
+                    factors[I][permindex.convertLinearIndex( li )] = val;
+                }
+            }
+
+            if( verbose >= 3 )
+                cout << "factors:" << factors << endl;
+
+
+            is.close();
+        } catch (...) {
+            is.close();
+            throw;
+        }
+    } else
+        DAI_THROW(CANNOT_READ_FILE);
+
+    return pair<vector<Var>, vector<Factor> >(vars,factors);
+}
+
+
+int main( int argc, char *argv[] ) {
+    if ( argc != 6 ) {
+        cout << "Usage: " << argv[0] << " <filename.uai> <filename.uai.evid> <filename.fg> <type> <run_jtree>" << endl << endl;
+        cout << "Converts input files from UAI approximate inference evaluation to" << endl;
+        cout << "libDAI factor graph format." << endl;
+        cout << "Reads factor graph <filename.uai> and evidence <filename.uai.evid>" << endl;
+        cout << "and writes the clamped factor graph to <filename.fg>." << endl;
+        cout << "If type==0, uses surgery, otherwise, uses safe method." << endl;
+        cout << "If run_jtree!=0, runs a junction tree and reports the results in the UAI 2008 results file format." << endl;
+        return 1;
+    } else {
+        long verbose = 0;
+        long type = atoi( argv[4] );
+        bool run_jtree = atoi( argv[5] );
+
+        pair<vector<Var>, vector<Factor> > input;
+        input = ReadUAIFGFile( argv[1], verbose );
+        FactorGraph fg0( input.second );
+        map<size_t,size_t> evid = ReadUAIEvidenceFile( argv[2] );
+        double logZcorr = 0.0;
+        if( type == 0 ) {
+            for( size_t I = 0; I < input.second.size(); I++ ) {
+                for( map<size_t,size_t>::const_iterator e = evid.begin(); e != evid.end(); e++ ) {
+                    if( input.second[I].vars() >> input.first[e->first] ) {
+                        if( verbose >= 2 )
+                            cout << "clamping " << e->first << " to value " << e->second << " in factor " << I << " = " << input.second[I].vars() << endl;
+                        input.second[I] = input.second[I].slice( input.first[e->first], e->second );
+                        if( verbose >= 2 )
+                            cout << "...remaining vars: " << input.second[I].vars() << endl;
+                    }
+                }
+            }
+            // remove empty factors
+            for( vector<Factor>::iterator I = input.second.begin(); I != input.second.end(); )
+                if( I->vars().size() == 0 ) {
+                    logZcorr += std::log( (Real)(*I)[0] );
+                    I = input.second.erase( I );
+                } else
+                    I++;
+            // multiply with logZcorr constant
+            if( input.second.size() == 0 )
+                input.second.push_back( Factor( VarSet(), std::exp(logZcorr) ) );
+            else
+                input.second.front() *= std::exp(logZcorr);
+        }
+        for( map<size_t,size_t>::const_iterator e = evid.begin(); e != evid.end(); e++ ) {
+            // add delta factors
+            input.second.push_back( Factor( input.first[e->first], 0.0 ) );
+            input.second.back()[e->second] = 1.0;
+        }
+        FactorGraph fg(input.second);
+
+        fg.WriteToFile( argv[3] );
+
+        if( run_jtree ) {
+            JTree jt0( fg0, PropertySet()("updates",string("HUGIN")) );
+            jt0.init();
+            jt0.run();
+
+            JTree jt( fg, PropertySet()("updates",string("HUGIN")) );
+            jt.init();
+            jt.run();
+
+            cout.precision( 8 );
+            if( evid.size() )
+                cout << "z " << (jt.logZ() - jt0.logZ()) / dai::log((Real)10.0) << endl;
+            else
+                cout << "z " << jt.logZ() / dai::log((Real)10.0) << endl;
+
+            cout << "m " << jt.nrVars() << " ";
+            for( size_t i = 0; i < jt.nrVars(); i++ ) {
+                cout << jt.var(i).states() << " ";
+                for( size_t s = 0; s < jt.var(i).states(); s++ )
+                    cout << jt.beliefV(i)[s] << " ";
+            }
+            cout << endl;
+
+            cout << "s ";
+            jt.props.inference = JTree::Properties::InfType::MAXPROD;
+            jt.init();
+            jt.run();
+            vector<size_t> MAP = jt.findMaximum();
+            map<Var, size_t> state;
+            for( size_t i = 0; i < MAP.size(); i++ )
+                state[jt.var(i)] = MAP[i];
+            double log_MAP_prob = 0.0;
+            for( size_t I = 0; I < jt.nrFactors(); I++ )
+                log_MAP_prob += dai::log( jt.factor(I)[calcLinearState( jt.factor(I).vars(), state )] );
+            if( evid.size() )
+                cout << (log_MAP_prob - jt0.logZ()) / dai::log((Real)10.0) << " ";
+            else
+                cout << log_MAP_prob / dai::log((Real)10.0) << " ";
+
+            cout << jt.nrVars() << " ";
+            for( size_t i = 0; i < jt.nrVars(); i++ )
+                cout << MAP[i] << " ";
+            cout << endl;
+        }
+    }
+
+    return 0;
+}