Added Gibbs code provided by Frederik Eaton
authorJoris Mooij <joris@jorismooij.nl>
Fri, 14 Nov 2008 15:55:13 +0000 (16:55 +0100)
committerJoris Mooij <joris@jorismooij.nl>
Fri, 14 Nov 2008 15:55:13 +0000 (16:55 +0100)
Makefile
Makefile.shared
Makefile.win
include/dai/alldai.h
include/dai/gibbs.h [new file with mode: 0644]
src/alldai.cpp
src/gibbs.cpp [new file with mode: 0644]

index e0fa522..79addfd 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -27,6 +27,7 @@ WITH_LC=true
 WITH_TREEEP=true
 WITH_JTREE=true
 WITH_MR=true
+WITH_GIBBS=true
 # Build with debug info?
 DEBUG=true
 # Build matlab interface?
index e2790d4..c137f74 100644 (file)
@@ -67,6 +67,10 @@ ifdef WITH_MR
        CCFLAGS:=$(CCFLAGS) -DDAI_WITH_MR
        OBJECTS:=$(OBJECTS) mr$(OE)
 endif
+ifdef WITH_GIBBS
+       CCFLAGS:=$(CCFLAGS) -DDAI_WITH_GIBBS
+       OBJECTS:=$(OBJECTS) gibbs$(OE)
+endif
 
 
 HEADERS=$(INC)/bipgraph.h $(INC)/index.h $(INC)/var.h $(INC)/factor.h $(INC)/varset.h $(INC)/smallset.h $(INC)/prob.h $(INC)/daialg.h $(INC)/properties.h $(INC)/alldai.h $(INC)/enum.h $(INC)/exceptions.h
@@ -128,6 +132,9 @@ weightedgraph$(OE) : $(SRC)/weightedgraph.cpp $(INC)/weightedgraph.h $(HEADERS)
 mr$(OE) : $(SRC)/mr.cpp $(INC)/mr.h $(HEADERS)
        $(CC) $(CCFLAGS) -c $(SRC)/mr.cpp
 
+gibbs$(OE) : $(SRC)/gibbs.cpp $(INC)/gibbs.h $(HEADERS)
+       $(CC) $(CCFLAGS) -c $(SRC)/gibbs.cpp
+
 properties$(OE) : $(SRC)/properties.cpp $(HEADERS)
        $(CC) $(CCFLAGS) -c $(SRC)/properties.cpp
 
index d0281d0..459e41a 100755 (executable)
@@ -27,6 +27,7 @@ WITH_LC=true
 WITH_TREEEP=true\r
 WITH_JTREE=true\r
 WITH_MR=true\r
+WITH_GIBBS=true\r
 # Build with debug info?\r
 DEBUG=true\r
 # Build matlab interface?\r
index d13114b..430e8eb 100644 (file)
@@ -54,6 +54,9 @@
 #ifdef DAI_WITH_MR
     #include <dai/mr.h>
 #endif
+#ifdef DAI_WITH_GIBBS
+    #include <dai/gibbs.h>
+#endif
 
 
 /// Namespace for libDAI
@@ -92,6 +95,9 @@ static const char* DAINames[] = {
 #endif
 #ifdef DAI_WITH_MR
     MR::Name,
+#endif
+#ifdef DAI_WITH_GIBBS
+    Gibbs::Name,
 #endif
     ""
 };
diff --git a/include/dai/gibbs.h b/include/dai/gibbs.h
new file mode 100644 (file)
index 0000000..fa9e296
--- /dev/null
@@ -0,0 +1,108 @@
+/*  Copyright (C) 2008  Frederik Eaton [frederik at ofb dot net]
+
+    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
+*/
+
+
+#ifndef __defined_libdai_gibbs_h
+#define __defined_libdai_gibbs_h
+
+
+#include <dai/daialg.h>
+#include <dai/factorgraph.h>
+#include <dai/properties.h>
+
+
+namespace dai {
+
+
+class Gibbs : public DAIAlgFG {
+    public:
+        /// Parameters of this inference algorithm
+        struct Properties {
+            /// Number of iterations
+            size_t iters;
+
+            /// Verbosity
+            size_t verbose;
+        } props;
+
+        /// Name of this inference algorithm
+        static const char *Name;
+
+    protected:
+        typedef std::vector<size_t> _count_t;
+        size_t _sample_count;
+        std::vector<_count_t> _var_counts;
+        std::vector<_count_t> _factor_counts;
+
+        typedef std::vector<size_t> _state_t;
+        void update_counts(_state_t &st);
+        void randomize_state(_state_t &st);
+        Prob get_var_dist(_state_t &st, size_t i);
+        void resample_var(_state_t &st, size_t i);
+        size_t get_factor_entry(const _state_t &st, int factor);
+
+    public:
+        // default constructor
+        Gibbs() : DAIAlgFG() {}
+        // copy constructor
+        Gibbs(const Gibbs & x) : DAIAlgFG(x), _sample_count(x._sample_count), _var_counts(x._var_counts), _factor_counts(x._factor_counts) {}
+        // construct Gibbs object from FactorGraph
+        Gibbs( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg) {
+            setProperties( opts );
+            construct();
+        }
+        // assignment operator
+        Gibbs & operator=(const Gibbs & x) {
+            if(this!=&x) {
+                DAIAlgFG::operator=(x);
+                _sample_count = x._sample_count;
+                _var_counts = x._var_counts;
+                _factor_counts = x._factor_counts;
+            }
+            return *this;
+        }
+        
+        virtual Gibbs* clone() const { return new Gibbs(*this); }
+        virtual Gibbs* create() const { return new Gibbs(); }
+        virtual std::string identify() const { return std::string(Name) + printProperties(); }
+        virtual Factor belief( const Var &n ) const;
+        virtual Factor belief( const VarSet &ns ) const;
+        virtual std::vector<Factor> beliefs() const;
+        virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
+        virtual void init();
+        virtual void init( const VarSet &/*ns*/ ) { init(); }
+        virtual double run();
+        virtual double maxDiff() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
+        virtual size_t Iterations() const { return props.iters; }
+
+        Factor beliefV( size_t i ) const;
+        Factor beliefF( size_t I ) const;
+
+        void construct();
+        /// Set Props according to the PropertySet opts, where the values can be stored as std::strings or as the type of the corresponding Props member
+        void setProperties( const PropertySet &opts );
+        PropertySet getProperties() const;
+        std::string printProperties() const;
+};
+
+
+} // end of namespace dai
+
+
+#endif
index 1f8d898..643c917 100644 (file)
@@ -62,6 +62,10 @@ InfAlg *newInfAlg( const string &name, const FactorGraph &fg, const PropertySet
 #ifdef DAI_WITH_MR
     if( name == MR::Name )
         return new MR (fg, opts);
+#endif
+#ifdef DAI_WITH_GIBBS
+    if( name == Gibbs::Name )
+        return new Gibbs (fg, opts);
 #endif
     DAI_THROW(UNKNOWN_DAI_ALGORITHM);
 }
diff --git a/src/gibbs.cpp b/src/gibbs.cpp
new file mode 100644 (file)
index 0000000..ca430d1
--- /dev/null
@@ -0,0 +1,260 @@
+/*  Copyright (C) 2008  Frederik Eaton [frederik at ofb dot net]
+
+    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
+*/
+
+
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <set>
+#include <algorithm>
+#include <dai/gibbs.h>
+#include <dai/util.h>
+#include <dai/properties.h>
+
+
+namespace dai {
+
+
+using namespace std;
+
+
+const char *Gibbs::Name = "GIBBS";
+
+
+void Gibbs::setProperties( const PropertySet &opts ) {
+    assert( opts.hasKey("iters") );
+    props.iters = opts.getStringAs<size_t>("iters");
+
+    if( opts.hasKey("verbose") )
+        props.verbose = opts.getStringAs<size_t>("verbose");
+    else
+        props.verbose = 0;
+}
+
+
+PropertySet Gibbs::getProperties() const {
+    PropertySet opts;
+    opts.Set( "iters", props.iters );
+    opts.Set( "verbose", props.verbose );
+    return opts;
+}
+
+
+string Gibbs::printProperties() const {
+    stringstream s( stringstream::out );
+    s << "[";
+    s << "iters=" << props.iters << ",";
+    s << "verbose=" << props.verbose << "]";
+    return s.str();
+}
+
+
+void Gibbs::construct() {
+    _var_counts.clear();
+    _var_counts.reserve( nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        _var_counts.push_back( _count_t( var(i).states(), 0 ) );
+    _factor_counts.clear();
+    _factor_counts.reserve( nrFactors() );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        _factor_counts.push_back( _count_t( factor(I).states(), 0 ) );
+    _sample_count = 0;
+}
+
+
+void Gibbs::update_counts( _state_t &st ) {
+    for( size_t i = 0; i < nrVars(); i++ )
+        _var_counts[i][st[i]]++;
+    for( size_t I = 0; I < nrFactors(); I++ ) {
+        if( 0 ) {
+/*            multind mi( factor(I).vars() );
+            _state_t f_st( factor(I).vars().size() );
+            int k = 0;
+            foreach( size_t j, nbF(I) )
+                f_st[k++] = st[j];
+            _factor_counts[I][mi.li(f_st)]++;*/
+        } else {
+            size_t ent = get_factor_entry(st, I);
+            _factor_counts[I][ent]++;
+        }
+    }
+    _sample_count++;
+}
+
+
+inline
+size_t Gibbs::get_factor_entry(const _state_t &st, int factor) {
+  size_t f_entry=0;
+  int rank = nbF(factor).size();
+  for(int j=rank-1; j>=0; j--) {
+      int jn = nbF(factor)[j];
+      f_entry *= var(jn).states();
+      f_entry += st[jn];
+  }
+  return f_entry;
+}
+
+
+Prob Gibbs::get_var_dist( _state_t &st, size_t i ) {
+    assert( st.size() == vars().size() );
+    assert( i < nrVars() );
+    if( 1 ) {
+        // use markov blanket of n to calculate distribution
+        size_t dim = var(i).states();
+        Neighbors &facts = nbV(i);
+
+        Prob values( dim, 1.0 );
+
+        for( size_t I = 0; I < facts.size(); I++ ) {
+            size_t fa = facts[I];
+            const Factor &f = factor(fa);
+            int save_ind = st[i];
+            for( size_t k = 0; k < dim; k++ ) {
+                st[i] = k;
+                int f_entry = get_factor_entry(st, fa);
+                values[k] *= f[f_entry];
+            }
+            st[i] = save_ind;
+        }
+
+        return values.normalized();
+    } else {
+/*        Var vi = var(i);
+        Factor d(vi);
+        assert(vi.states()>0);
+        assert(vi.label()>=0);
+        // loop over factors containing i (nbV(i)):
+        foreach(size_t I, nbV(i)) {
+            // use multind to find linear state for variables != i in factor
+            assert(I<nrFactors());
+            assert(factor(I).vars().size() > 0);
+            VarSet vs (factor(I).vars() / vi);
+            multind mi(vs);
+            _state_t I_st(vs.size());
+            int k=0;
+            foreach(size_t l, nbF(I)) {
+                if(l!=i) I_st[k++] = st[l];
+            }
+            // use slice(ns,ns_state) to get beliefs for variable i
+            // multiply all these beliefs together
+            d *= factor(I).slice(vs, mi.li(I_st));
+        }
+        d.p().normalize();
+        return d.p();*/
+    }
+}
+
+
+void Gibbs::resample_var( _state_t &st, size_t i ) {
+    // draw randomly from conditional distribution and update 'st'
+    st[i] = get_var_dist( st, i ).draw();
+}
+
+
+void Gibbs::randomize_state( _state_t &st ) {
+    assert( st.size() == nrVars() );
+    for( size_t i = 0; i < nrVars(); i++ )
+        st[i] = rnd_int( 0, var(i).states() - 1 );
+}
+
+
+void Gibbs::init() {
+    for( size_t i = 0; i < nrVars(); i++ )
+        fill( _var_counts[i].begin(), _var_counts[i].end(), 0 );
+    for( size_t I = 0; I < nrFactors(); I++ )
+        fill( _factor_counts[I].begin(), _factor_counts[I].end(), 0 );
+    _sample_count = 0;
+}
+
+
+double Gibbs::run() {
+    if( props.verbose >= 1 )
+        cout << "Starting " << identify() << "...";
+    if( props.verbose >= 3 )
+        cout << endl;
+
+    double tic = toc();
+    
+    vector<size_t> state( nrVars() );
+    randomize_state( state );
+
+    for( size_t iter = 0; iter < props.iters; iter++ ) {
+        for( size_t j = 0; j < nrVars(); j++ )
+            resample_var( state, j );
+        update_counts( state );
+    }
+
+    if( props.verbose >= 3 ) {
+        for( size_t i = 0; i < nrVars(); i++ ) {
+            cerr << "belief for variable " << var(i) << ": " << beliefV(i) << endl;
+            cerr << "counts for variable " << var(i) << ": " << Prob( _var_counts[i].begin(), _var_counts[i].end() ) << endl;
+        }
+    }
+    
+    if( props.verbose >= 3 )
+        cout << "Gibbs::run:  ran " << props.iters << " passes (" << toc() - tic << " clocks)." << endl;
+
+    return 0.0;
+}
+
+
+Factor Gibbs::beliefV( size_t i ) const {
+    Prob p( _var_counts[i].begin(), _var_counts[i].end() );
+    p.normalize();
+    return( Factor( var(i), p ) );
+}
+
+
+Factor Gibbs::beliefF( size_t I ) const {
+    Prob p( _factor_counts[I].begin(), _factor_counts[I].end() );
+    p.normalize();
+    return( Factor( factor(I).vars(), p ) );
+}
+
+
+Factor Gibbs::belief( const Var &n ) const {
+    return( beliefV( findVar( n ) ) );
+}
+
+
+vector<Factor> Gibbs::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;
+}
+
+
+Factor Gibbs::belief( const VarSet &ns ) const {
+    if( ns.size() == 1 )
+        return belief( *(ns.begin()) );
+    else {
+        size_t I;
+        for( I = 0; I < nrFactors(); I++ )
+            if( factor(I).vars() >> ns )
+                break;
+        assert( I != nrFactors() );
+        return beliefF(I).marginal(ns);
+    }
+}
+
+
+} // end of namespace dai