* [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& )
- 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.
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)
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
##########
-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
-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)
#include <iostream>
+#include <map>
#include <dai/alldai.h> // Include main libDAI header file
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;
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;
}
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
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;
// 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
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 )
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) );
+ }
}
}
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 );
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;
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++ )
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;
}
--- /dev/null
+/* 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;
+}
projects / libdai.git / commitdiff