From 3b28074fdffcca0570d8c9424a5558b6e18c82f2 Mon Sep 17 00:00:00 2001 From: Joris Mooij Date: Tue, 3 Aug 2010 14:38:15 +0200 Subject: [PATCH] Added decimation algorithm DECMAP --- ChangeLog | 2 + Makefile | 7 ++ Makefile.ALL | 1 + README | 3 +- examples/example.cpp | 18 ++++- include/dai/alldai.h | 6 ++ include/dai/daialg.h | 10 ++- include/dai/decmap.h | 97 ++++++++++++++++++++++ include/dai/doc.h | 6 +- include/dai/factor.h | 7 ++ src/alldai.cpp | 4 + src/decmap.cpp | 189 +++++++++++++++++++++++++++++++++++++++++++ src/factor.cpp | 8 ++ 13 files changed, 351 insertions(+), 7 deletions(-) create mode 100644 include/dai/decmap.h create mode 100644 src/decmap.cpp diff --git a/ChangeLog b/ChangeLog index b39f939..8f59471 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,5 +1,7 @@ git HEAD -------- +* Added DECMAP algorithm. +* Added createFactorDelta( const VarSet& vs, size_t state ) * [Peter Rockett] Improved Makefiles for image segmentation example. * Removed deprecated interfaces. * Added FactorGraph::isMaximal(size_t) and FactorGraph::maximalFactor(size_t). diff --git a/Makefile b/Makefile index 8b84af2..d76c56f 100644 --- a/Makefile +++ b/Makefile @@ -89,6 +89,10 @@ ifdef WITH_CBP WITHFLAGS:=$(WITHFLAGS) -DDAI_WITH_CBP NAMES:=$(NAMES) bbp cbp bp_dual endif +ifdef WITH_DECMAP + WITHFLAGS:=$(WITHFLAGS) -DDAI_WITH_DECMAP + NAMES:=$(NAMES) decmap +endif # Define standard libDAI header dependencies, source file names and object file names HEADERS=$(foreach name,graph dag bipgraph index var factor varset smallset prob daialg properties alldai enum exceptions util,$(INC)/$(name).h) @@ -175,6 +179,9 @@ treeep$(OE) : $(SRC)/treeep.cpp $(INC)/treeep.h $(HEADERS) $(INC)/weightedgraph. emalg$(OE) : $(SRC)/emalg.cpp $(INC)/emalg.h $(INC)/evidence.h $(HEADERS) $(CC) -c $< +decmap$(OE) : $(SRC)/decmap.cpp $(INC)/decmap.h $(HEADERS) + $(CC) -c $< + # EXAMPLES ########### diff --git a/Makefile.ALL b/Makefile.ALL index a13c11c..df10601 100644 --- a/Makefile.ALL +++ b/Makefile.ALL @@ -31,6 +31,7 @@ WITH_JTREE=true WITH_MR=true WITH_GIBBS=true WITH_CBP=true +WITH_DECMAP=true # Build doxygen documentation? WITH_DOC= diff --git a/README b/README index 1fa7815..068c130 100644 --- a/README +++ b/README @@ -87,7 +87,8 @@ Currently, libDAI supports the following (approximate) inference methods: * Double-loop GBP [HAK03]; * Various variants of Loop Corrected Belief Propagation [MoK07, MoR05]; * Gibbs sampler; - * Conditioned Belief Propagation [EaG09]. + * Conditioned Belief Propagation [EaG09]; + * Decimation algorithm. These inference methods can be used to calculate partition sums, marginals over subsets of variables, and MAP states (the joint state of variables that has diff --git a/examples/example.cpp b/examples/example.cpp index 390c4da..1a51931 100644 --- a/examples/example.cpp +++ b/examples/example.cpp @@ -82,7 +82,7 @@ int main( int argc, char *argv[] ) { // using the parameters specified by opts and two additional properties, // specifying the type of updates the BP algorithm should perform and // whether they should be done in the real or in the logdomain - BP bp(fg, opts("updates",string("SEQFIX"))("logdomain",false)); + BP bp(fg, opts("updates",string("SEQRND"))("logdomain",false)); // Initialize belief propagation algorithm bp.init(); // Run belief propagation algorithm @@ -95,7 +95,7 @@ int main( int argc, char *argv[] ) { // // Note that inference is set to MAXPROD, which means that the object // will perform the max-product algorithm instead of the sum-product algorithm - BP mp(fg, opts("updates",string("SEQFIX"))("logdomain",false)("inference",string("MAXPROD"))); + BP mp(fg, opts("updates",string("SEQRND"))("logdomain",false)("inference",string("MAXPROD"))("damping",string("0.1"))); // Initialize max-product algorithm mp.init(); // Run max-product algorithm @@ -104,6 +104,15 @@ int main( int argc, char *argv[] ) { // based on the max-product result vector mpstate = mp.findMaximum(); + // Construct a decimation algorithm object from the FactorGraph fg + // using the parameters specified by opts and three additional properties, + // specifying that the decimation algorithm should use the max-product + // algorithm and should completely reinitalize its state at every step + DecMAP decmap(fg, opts("reinit",true)("ianame",string("BP"))("iaopts",string("[damping=0.1,inference=MAXPROD,logdomain=0,maxiter=1000,tol=1e-9,updates=SEQRND,verbose=1]")) ); + decmap.init(); + decmap.run(); + vector decmapstate = decmap.findMaximum(); + // Report variable marginals for fg, calculated by the junction tree algorithm cout << "Exact variable marginals:" << endl; for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg @@ -159,6 +168,11 @@ int main( int argc, char *argv[] ) { 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; + + // Report DecMAP joint state + cout << "Approximate DecMAP state (log probability = " << evalState( fg, decmapstate ) << "):" << endl; + for( size_t i = 0; i < decmapstate.size(); i++ ) + cout << fg.var(i) << ": " << decmapstate[i] << endl; } return 0; diff --git a/include/dai/alldai.h b/include/dai/alldai.h index bdf1bd3..98b6d69 100644 --- a/include/dai/alldai.h +++ b/include/dai/alldai.h @@ -61,6 +61,9 @@ #ifdef DAI_WITH_CBP #include #endif +#ifdef DAI_WITH_DECMAP + #include +#endif /// Namespace for libDAI @@ -151,6 +154,9 @@ static const char* DAINames[] = { #endif #ifdef DAI_WITH_CBP CBP::Name, +#endif +#ifdef DAI_WITH_DECMAP + DecMAP::Name, #endif "" }; diff --git a/include/dai/daialg.h b/include/dai/daialg.h index 0242c79..249172f 100644 --- a/include/dai/daialg.h +++ b/include/dai/daialg.h @@ -111,15 +111,21 @@ class InfAlg { */ virtual Real logZ() const = 0; + /// Calculates the joint state of all variables that has maximum probability + /** \note Before this method is called, run() should have been called. + * \throw NOT_IMPLEMENTED if not implemented/supported + */ + virtual std::vector findMaximum() const { DAI_THROW(NOT_IMPLEMENTED); } + /// Returns maximum difference between single variable beliefs in the last iteration. /** \throw NOT_IMPLEMENTED if not implemented/supported */ - virtual Real maxDiff() const = 0; + virtual Real maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); }; /// Returns number of iterations done (one iteration passes over the complete factorgraph). /** \throw NOT_IMPLEMENTED if not implemented/supported */ - virtual size_t Iterations() const = 0; + virtual size_t Iterations() const { DAI_THROW(NOT_IMPLEMENTED); }; //@} /// \name Changing the factor graph diff --git a/include/dai/decmap.h b/include/dai/decmap.h new file mode 100644 index 0000000..4229aa5 --- /dev/null +++ b/include/dai/decmap.h @@ -0,0 +1,97 @@ +/* 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) 2010 Joris Mooij [joris dot mooij at libdai dot org] + */ + + +/// \file +/// \brief Defines class DecMAP, which constructs a MAP state by decimation + + +#ifndef __defined_libdai_decmap_h +#define __defined_libdai_decmap_h + + +#include + + +namespace dai { + + +/// Approximate inference algorithm DecMAP, which constructs a MAP state by decimation +/** Decimation involves repeating the following two steps until no free variables remain: + * - run an approximate inference algorithm, + * - clamp the factor with the lowest entropy to its most probable state + */ +class DecMAP : public DAIAlgFG { + private: + /// Stores the final MAP state + std::vector _state; + /// Stores the log probability of the MAP state + Real _logp; + /// Maximum difference encountered so far + Real _maxdiff; + /// Number of iterations needed + size_t _iters; + + public: + /// Parameters for DecMAP + struct Properties { + /// Verbosity (amount of output sent to stderr) + size_t verbose; + + /// Complete or partial reinitialization of clamped subgraphs? + bool reinit; + + /// Name of the algorithm used to calculate the beliefs on clamped subgraphs + std::string ianame; + + /// Parameters for the algorithm used to calculate the beliefs on clamped subgraphs + PropertySet iaopts; + } props; + + /// Name of this inference algorithm + static const char *Name; + + public: + /// Default constructor + DecMAP() : DAIAlgFG(), _state(), _logp(), _maxdiff(), _iters(), props() {} + + /// Construct from FactorGraph \a fg and PropertySet \a opts + /** \param fg Factor graph. + * \param opts Parameters @see Properties + */ + DecMAP( const FactorGraph &fg, const PropertySet &opts ); + + + /// \name General InfAlg interface + //@{ + virtual DecMAP* clone() const { return new DecMAP(*this); } + virtual std::string identify() const; + virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); } + virtual Factor belief( const VarSet &/*vs*/ ) const; + virtual Factor beliefV( size_t i ) const; + virtual Factor beliefF( size_t I ) const { return belief( factor(I).vars() ); } + virtual std::vector beliefs() const; + virtual Real logZ() const { return _logp; } + virtual void init() { _maxdiff = 0.0; _iters = 0; } + virtual void init( const VarSet &/*ns*/ ) { init(); } + virtual Real run(); + virtual std::vector findMaximum() const { return _state; } + virtual Real maxDiff() const { return _maxdiff; } + virtual size_t Iterations() const { return _iters; } + virtual void setProperties( const PropertySet &opts ); + virtual PropertySet getProperties() const; + virtual std::string printProperties() const; + //@} +}; + + +} // end of namespace dai + + +#endif diff --git a/include/dai/doc.h b/include/dai/doc.h index b08a337..c44be0b 100644 --- a/include/dai/doc.h +++ b/include/dai/doc.h @@ -77,7 +77,8 @@ * - Various variants of Loop Corrected Belief Propagation * [\ref MoK07, \ref MoR05]; * - Gibbs sampler; - * - Conditioned Belief Propagation [\ref EaG09]. + * - Conditioned Belief Propagation [\ref EaG09]; + * - Decimation algorithm. * * These inference methods can be used to calculate partition sums, marginals * over subsets of variables, and MAP states (the joint state of variables that @@ -479,9 +480,10 @@ * - Loop Corrected Belief Propagation: dai::MR [\ref MoR05] and dai::LC [\ref MoK07] * - Gibbs sampling: dai::Gibbs * - Conditioned Belief Propagation: dai::CBP [\ref EaG09] + * - Decimation algorithm: dai::DECMAP * * Not all inference tasks are implemented by each method: calculating MAP states - * is only possible with dai::JTree and dai::BP, calculating partition sums is + * is only possible with dai::JTree, dai::BP and dai::DECMAP; calculating partition sums is * not possible with dai::MR, dai::LC and dai::Gibbs. * * \section terminology-learning Parameter learning diff --git a/include/dai/factor.h b/include/dai/factor.h index 24b6d89..e354f50 100644 --- a/include/dai/factor.h +++ b/include/dai/factor.h @@ -669,6 +669,13 @@ Factor createFactorPotts( const Var &x1, const Var &x2, Real J ); Factor createFactorDelta( const Var &v, size_t state ); +/// Returns a Kronecker delta point mass +/** \param vs Set of variables + * \param state The state of \a vs that should get value 1 + */ +Factor createFactorDelta( const VarSet& vs, size_t state ); + + } // end of namespace dai diff --git a/src/alldai.cpp b/src/alldai.cpp index 46a5a39..266227f 100644 --- a/src/alldai.cpp +++ b/src/alldai.cpp @@ -68,6 +68,10 @@ InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const Propert #ifdef DAI_WITH_CBP if( name == CBP::Name ) return new CBP (fg, opts); +#endif +#ifdef DAI_WITH_DECMAP + if( name == DecMAP::Name ) + return new DecMAP (fg, opts); #endif DAI_THROWE(UNKNOWN_DAI_ALGORITHM,"Unknown libDAI algorithm: " + name); } diff --git a/src/decmap.cpp b/src/decmap.cpp new file mode 100644 index 0000000..dc9300f --- /dev/null +++ b/src/decmap.cpp @@ -0,0 +1,189 @@ +/* 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) 2010 Joris Mooij [joris dot mooij at libdai dot org] + */ + + +#include + + +namespace dai { + + +using namespace std; + + +const char *DecMAP::Name = "DECMAP"; + + +void DecMAP::setProperties( const PropertySet &opts ) { + DAI_ASSERT( opts.hasKey("ianame") ); + DAI_ASSERT( opts.hasKey("iaopts") ); + + props.ianame = opts.getStringAs("ianame"); + props.iaopts = opts.getStringAs("iaopts"); + if( opts.hasKey("verbose") ) + props.verbose = opts.getStringAs("verbose"); + else + props.verbose = 0; + if( opts.hasKey("reinit") ) + props.reinit = opts.getStringAs("reinit"); + else + props.reinit = true; +} + + +PropertySet DecMAP::getProperties() const { + PropertySet opts; + opts.set( "verbose", props.verbose ); + opts.set( "reinit", props.reinit ); + opts.set( "ianame", props.ianame ); + opts.set( "iaopts", props.iaopts ); + return opts; +} + + +string DecMAP::printProperties() const { + stringstream s( stringstream::out ); + s << "["; + s << "verbose=" << props.verbose << ","; + s << "reinit=" << props.reinit << ","; + s << "ianame=" << props.ianame << ","; + s << "iaopts=" << props.iaopts << "]"; + return s.str(); +} + + +DecMAP::DecMAP( const FactorGraph& fg, const PropertySet& opts ) : DAIAlgFG(fg), _state(), _logp(), _maxdiff(), _iters(), props() { + setProperties( opts ); + + _state = vector( nrVars(), 0 ); + _logp = -INFINITY; +} + + +string DecMAP::identify() const { + return string(Name) + printProperties(); +} + + +Factor DecMAP::belief( const VarSet& vs ) const { + if( vs.size() == 0 ) + return Factor(); + else { + map state; + for( VarSet::const_iterator v = vs.begin(); v != vs.end(); v++ ) + state[*v] = _state[findVar(*v)]; + return createFactorDelta( vs, calcLinearState( vs, state ) ); + } +} + + +Factor DecMAP::beliefV( size_t i ) const { + return createFactorDelta( var(i), _state[i] ); +} + + +vector DecMAP::beliefs() const { + vector result; + for( size_t i = 0; i < nrVars(); ++i ) + result.push_back( beliefV(i) ); + for( size_t I = 0; I < nrFactors(); ++I ) + result.push_back( beliefF(I) ); + return result; +} + + +Real DecMAP::run() { + if( props.verbose >= 1 ) + cerr << "Starting " << identify() << "..."; + if( props.verbose >= 2 ) + cerr << endl; + + // the variables which have not been clamped yet + SmallSet freeVars; + for( size_t i = 0; i < nrVars(); i++ ) + freeVars |= i; + + // prepare the inference algorithm object + InfAlg *clamped = newInfAlg( props.ianame, fg(), props.iaopts ); + + // decimate until no free variables remain + while( freeVars.size() ) { + Real md = clamped->run(); + if( md > _maxdiff ) + _maxdiff = md; + _iters += clamped->Iterations(); + + // store the variables that need initialization + VarSet varsToInit; + SmallSet varsToClamp; + + // schedule clamping for the free variables with zero entropy + for( SmallSet::const_iterator it = freeVars.begin(); it != freeVars.end(); ) { + if( clamped->beliefV( *it ).entropy() == 0.0 ) { + // this variable should be clamped + varsToInit |= var( *it ); + varsToClamp |= *it; + _state[*it] = clamped->beliefV( *it ).p().argmax().first; + freeVars.erase( *it ); + } else + it++; + } + + // find the free factor with lowest entropy + size_t bestI = 0; + Real bestEnt = INFINITY; + for( size_t I = 0; I < nrFactors(); I++ ) { + // check if the factor is still free + if( freeVars.intersects( bipGraph().nb2Set(I) ) ) { + Real EntI = clamped->beliefF(I).entropy(); + if( EntI < bestEnt ) { + bestI = I; + bestEnt = EntI; + } + } + } + + // schedule clamping for the factor with lowest entropy + vector Istate(1,0); + Istate[0] = clamped->beliefF(bestI).p().argmax().first; + map Istatemap = calcState( factor(bestI).vars(), Istate[0] ); + foreach( size_t i, bipGraph().nb2Set(bestI) & freeVars ) { + varsToInit |= var(i); + varsToClamp |= i; + _state[i] = Istatemap[var(i)]; + freeVars.erase(i); + } + + // clamp all variables scheduled for clamping + foreach( size_t i, varsToClamp ) + clamped->clamp( i, _state[i], false ); + + // initialize clamped for the next run + if( props.reinit ) + clamped->init(); + else + clamped->init( varsToInit ); + } + + // calculate MAP state + map state; + for( size_t i = 0; i < nrVars(); i++ ) + state[var(i)] = _state[i]; + _logp = 0.0; + for( size_t I = 0; I < nrFactors(); I++ ) + _logp += dai::log( factor(I)[calcLinearState( factor(I).vars(), state )] ); + + // clean up + delete clamped; + + return _maxdiff; +} + + +} // end of namespace dai diff --git a/src/factor.cpp b/src/factor.cpp index bcc4905..624363d 100644 --- a/src/factor.cpp +++ b/src/factor.cpp @@ -63,4 +63,12 @@ Factor createFactorDelta( const Var &v, size_t state ) { } +Factor createFactorDelta( const VarSet& vs, size_t state ) { + Factor fac( vs, 0.0 ); + DAI_ASSERT( state < vs.nrStates() ); + fac.set( state, 1.0 ); + return fac; +} + + } // end of namespace dai -- 2.20.1