Added decimation algorithm DECMAP
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Aug 2010 12:38:15 +0000 (14:38 +0200)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Tue, 3 Aug 2010 12:38:15 +0000 (14:38 +0200)
13 files changed:
ChangeLog
Makefile
Makefile.ALL
README
examples/example.cpp
include/dai/alldai.h
include/dai/daialg.h
include/dai/decmap.h [new file with mode: 0644]
include/dai/doc.h
include/dai/factor.h
src/alldai.cpp
src/decmap.cpp [new file with mode: 0644]
src/factor.cpp

index b39f939..8f59471 100644 (file)
--- 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).
index 8b84af2..d76c56f 100644 (file)
--- 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
 ###########
index a13c11c..df10601 100644 (file)
@@ -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 (file)
--- 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
index 390c4da..1a51931 100644 (file)
@@ -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<size_t> 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<size_t> 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;
index bdf1bd3..98b6d69 100644 (file)
@@ -61,6 +61,9 @@
 #ifdef DAI_WITH_CBP
     #include <dai/cbp.h>
 #endif
+#ifdef DAI_WITH_DECMAP
+    #include <dai/decmap.h>
+#endif
 
 
 /// Namespace for libDAI
@@ -152,6 +155,9 @@ static const char* DAINames[] = {
 #ifdef DAI_WITH_CBP
     CBP::Name,
 #endif
+#ifdef DAI_WITH_DECMAP
+    DecMAP::Name,
+#endif
     ""
 };
 
index 0242c79..249172f 100644 (file)
@@ -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<std::size_t> 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 (file)
index 0000000..4229aa5
--- /dev/null
@@ -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 <dai/daialg.h>
+
+
+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<size_t> _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<Factor> 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<size_t> 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
index b08a337..c44be0b 100644 (file)
@@ -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
  *  - 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
index 24b6d89..e354f50 100644 (file)
@@ -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
 
 
index 46a5a39..266227f 100644 (file)
@@ -69,6 +69,10 @@ InfAlg *newInfAlg( const std::string &name, const FactorGraph &fg, const Propert
     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 (file)
index 0000000..dc9300f
--- /dev/null
@@ -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 <dai/alldai.h>
+
+
+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<string>("ianame");
+    props.iaopts = opts.getStringAs<PropertySet>("iaopts");
+    if( opts.hasKey("verbose") )
+        props.verbose = opts.getStringAs<size_t>("verbose");
+    else
+        props.verbose = 0;
+    if( opts.hasKey("reinit") )
+        props.reinit = opts.getStringAs<bool>("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<size_t>( 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<Var, size_t> 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<Factor> DecMAP::beliefs() const {
+    vector<Factor> 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<size_t> 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<size_t> varsToClamp;
+
+        // schedule clamping for the free variables with zero entropy
+        for( SmallSet<size_t>::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<size_t> Istate(1,0);
+        Istate[0] = clamped->beliefF(bestI).p().argmax().first;
+        map<Var, size_t> 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<Var, size_t> 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
index bcc4905..624363d 100644 (file)
@@ -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