Merge branch 'joris'
authorJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 2 Mar 2009 19:51:34 +0000 (20:51 +0100)
committerJoris Mooij <joris.mooij@tuebingen.mpg.de>
Mon, 2 Mar 2009 19:51:34 +0000 (20:51 +0100)
Conflicts:
example.cpp
include/dai/factor.h
include/dai/var.h
src/bp.cpp

1  2 
examples/example.cpp
include/dai/bp.h
include/dai/factor.h
include/dai/factorgraph.h
include/dai/index.h
include/dai/prob.h
include/dai/var.h
src/bp.cpp

index eba3630,0000000..0705bde
mode 100644,000000..100644
--- /dev/null
@@@ -1,100 -1,0 +1,131 @@@
-         // Report single-variable marginals for fg, calculated by the junction tree algorithm
-         cout << "Exact single-variable marginals:" << endl;
 +/*  Copyright (C) 2006-2008  Joris Mooij  [joris dot mooij at tuebingen dot mpg dot de]
 +    Radboud University Nijmegen, The Netherlands /
 +    Max Planck Institute for Biological Cybernetics, Germany
 +
 +    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 <dai/alldai.h>  // Include main libDAI header file
 +
 +
 +using namespace dai;
 +using namespace std;
 +
 +
 +int main( int argc, char *argv[] ) {
 +    if ( argc != 2 ) {
 +        cout << "Usage: " << argv[0] << " <filename.fg>" << endl << endl;
 +        cout << "Reads factor graph <filename.fg> and runs" << endl;
 +        cout << "Belief Propagation and JunctionTree on it." << endl << endl;
 +        return 1;
 +    } else {
 +        // Read FactorGraph from the file specified by the first command line argument
 +        FactorGraph fg;
 +        fg.ReadFromFile(argv[1]);
 +
 +        // Set some constants
 +        size_t  maxiter = 10000;
 +        double  tol = 1e-9;
 +        size_t  verb = 1;
 +
 +        // Store the constants in a PropertySet object
 +        PropertySet opts;
 +        opts.Set("maxiter",maxiter);  // Maximum number of iterations
 +        opts.Set("tol",tol);          // Tolerance for convergence
 +        opts.Set("verbose",verb);     // Verbosity (amount of output generated)
 +
 +        // Construct a JTree (junction tree) object from the FactorGraph fg
 +        // using the parameters specified by opts and an additional property
 +        // that specifies the type of updates the JTree algorithm should perform
 +        JTree jt( fg, opts("updates",string("HUGIN")) );
 +        // Initialize junction tree algoritm
 +        jt.init();
 +        // Run junction tree algorithm
 +        jt.run();
 +
 +        // Construct a BP (belief propagation) object from the FactorGraph fg
 +        // 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));
 +        // Initialize belief propagation algorithm
 +        bp.init();
 +        // Run belief propagation algorithm
 +        bp.run();
 +
-         // Report single-variable marginals for fg, calculated by the belief propagation algorithm
-         cout << "Approximate (loopy belief propagation) single-variable marginals:" << endl;
++        // Construct a BP (belief propagation) object from the FactorGraph fg
++        // 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
++        //
++        // 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")));
++        // Initialize max-product algorithm
++        mp.init();
++        // Run max-product algorithm
++        mp.run();
++        // Calculate joint state of all variables that has maximum probability
++        // based on the max-product result
++        vector<size_t> mpstate = mp.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
 +            cout << jt.belief(fg.var(i)) << endl; // display the "belief" of jt for that variable
 +
++        // Report variable marginals for fg, calculated by the belief propagation algorithm
++        cout << "Approximate (loopy belief propagation) variable marginals:" << endl;
 +        for( size_t i = 0; i < fg.nrVars(); i++ ) // iterate over all variables in fg
 +            cout << bp.belief(fg.var(i)) << endl; // display the belief of bp for that variable
 +
 +        // Report factor marginals for fg, calculated by the junction tree algorithm
 +        cout << "Exact factor marginals:" << endl;
 +        for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 +            cout << jt.belief(fg.factor(I).vars()) << endl;  // display the "belief" of jt for the variables in that factor
 +
 +        // Report factor marginals for fg, calculated by the belief propagation algorithm
 +        cout << "Approximate (loopy belief propagation) factor marginals:" << endl;
 +        for( size_t I = 0; I < fg.nrFactors(); I++ ) // iterate over all factors in fg
 +            cout << bp.belief(fg.factor(I).vars()) << endl; // display the belief of bp for the variables in that factor
 +
 +        // Report log partition sum (normalizing constant) of fg, calculated by the junction tree algorithm
 +        cout << "Exact log partition sum: " << jt.logZ() << endl;
 +
 +        // Report log partition sum of fg, approximated by the belief propagation algorithm
 +        cout << "Approximate (loopy belief propagation) log partition sum: " << bp.logZ() << endl;
++
++        // Report max-product variable marginals
++        cout << "Max-product variable marginals:" << endl;
++        for( size_t i = 0; i < fg.nrVars(); i++ )
++            cout << mp.belief(fg.var(i)) << endl;
++
++        // Report max-product factor marginals
++        cout << "Max-product factor marginals:" << endl;
++        for( size_t I = 0; I < fg.nrFactors(); I++ )
++            cout << mp.belief(fg.factor(I).vars()) << "=" << mp.beliefF(I) << endl;
++
++        // Report max-product joint state
++        cout << "Max-product state:" << endl;
++        for( size_t i = 0; i < mpstate.size(); i++ )
++            cout << fg.var(i) << ": " << mpstate[i] << endl;
 +    }
 +
 +    return 0;
 +}
Simple merge
@@@ -72,24 -68,25 +72,30 @@@ template <typename T> class TFactor 
          TProb<T>    _p;
  
      public:
 -        /// Construct Factor with empty VarSet
+         /// Iterator over factor entries
+               typedef typename TProb<T>::iterator iterator;
+         /// Const iterator over factor entries
+               typedef typename TProb<T>::const_iterator const_iterator;
 +        /// Constructs TFactor depending on no variables, with value p
          TFactor ( Real p = 1.0 ) : _vs(), _p(1,p) {}
  
 -        /// Construct Factor from VarSet
 +        /// Constructs TFactor depending on variables in ns, with uniform distribution
          TFactor( const VarSet& ns ) : _vs(ns), _p(_vs.nrStates()) {}
          
 -        /// Construct Factor from VarSet and initial value
 +        /// Constructs TFactor depending on variables in ns, with all values set to p
          TFactor( const VarSet& ns, Real p ) : _vs(ns), _p(_vs.nrStates(),p) {}
          
 -        /// Construct Factor from VarSet and initial array
 -        TFactor( const VarSet& ns, const Real *p ) : _vs(ns), _p(_vs.nrStates(),p) {}
 -
 -        /// Construct Factor from VarSet and TProb<T>
 +        /// Constructs TFactor depending on variables in ns, copying the values from the range starting at begin
 +        /** \param ns contains the variables that the new TFactor should depend on.
 +         *  \tparam Iterator Iterates over instances of type T; should support addition of size_t.
 +         *  \param begin Points to first element to be added.
 +         */
 +        template<typename TIterator>
 +        TFactor( const VarSet& ns, TIterator begin ) : _vs(ns), _p(begin, begin + _vs.nrStates(), _vs.nrStates()) {}
 +
 +        /// Constructs TFactor depending on variables in ns, with values set to the TProb p
          TFactor( const VarSet& ns, const TProb<T>& p ) : _vs(ns), _p(p) {
  #ifdef DAI_DEBUG
              assert( _vs.nrStates() == _p.size() );
          const VarSet & vars() const { return _vs; }
  
          /// Returns the number of possible joint states of the variables
 +        /** \note This is equal to the length of the value vector.
 +         */
          size_t states() const { return _p.size(); }
  
 -        /// Returns a copy of the i'th probability value
 +        /// Returns a copy of the i'th entry of the value vector
          T operator[] (size_t i) const { return _p[i]; }
  
 -        /// Returns a reference to the i'th probability value
 +        /// Returns a reference to the i'th entry of the value vector
          T& operator[] (size_t i) { return _p[i]; }
+         
+         /// Returns iterator pointing to first entry
+         iterator begin() { return _p.begin(); }
+         /// Returns const iterator pointing to first entry
+               const_iterator begin() const { return _p.begin(); }
+               /// Returns iterator pointing beyond last entry
+               iterator end() { return _p.end(); }
+               /// Returns const iterator pointing beyond last entry
+               const_iterator end() const { return _p.end(); }
  
 -        /// Sets all probability entries to p
 +        /// Sets all values to p
          TFactor<T> & fill (T p) { _p.fill( p ); return(*this); }
  
 -        /// Fills all probability entries with random values
 +        /// Draws all values i.i.d. from a uniform distribution on [0,1)
          TFactor<T> & randomize () { _p.randomize(); return(*this); }
  
 -        /// Returns product of *this with x
 -        TFactor<T> operator* (T x) const {
 -            Factor result = *this;
 -            result.p() *= x;
 -            return result;
 +
 +        /// Multiplies *this with scalar t
 +        TFactor<T>& operator*= (T t) {
 +            _p *= t;
 +            return *this;
          }
  
 -        /// Multiplies each probability entry with x
 -        TFactor<T>& operator*= (T x) {
 -            _p *= x;
 +        /// Divides *this by scalar t
 +        TFactor<T>& operator/= (T t) {
 +            _p /= t;
              return *this;
          }
  
Simple merge
Simple merge
@@@ -106,10 -108,22 +111,22 @@@ template <typename T> class TProb 
  #endif
          }
          
 -        /// Returns a reference to the i'th probability entry
 +        /// Returns reference to the i'th entry
          T& operator[]( size_t i ) { return _p[i]; }
+         
+         /// Returns iterator pointing to first entry
+         iterator begin() { return _p.begin(); }
+         /// Returns const iterator pointing to first entry
+         const_iterator begin() const { return _p.begin(); }
+         /// Returns iterator pointing beyond last entry
+         iterator end() { return _p.end(); }
+         /// Returns const iterator pointing beyond last entry
+         const_iterator end() const { return _p.end(); }
  
 -        /// Sets all elements to x
 +        /// Sets all entries to x
          TProb<T> & fill(T x) { 
              std::fill( _p.begin(), _p.end(), x );
              return *this;
@@@ -73,18 -66,42 +74,42 @@@ class Var 
          /// Returns reference to number of states
          size_t& states () { return _states; }
  
 -        /// Smaller-than operator (only compares labels)
 +        /// Smaller-than operator (compares only labels)
          bool operator < ( const Var& n ) const { return( _label <  n._label ); }
 -        /// Larger-than operator (only compares labels)
 +        /// Larger-than operator (compares only labels)
          bool operator > ( const Var& n ) const { return( _label >  n._label ); }
-         /// Smaller-than-or-equal-to operator (compares only labels)
-         bool operator <= ( const Var& n ) const { return( _label <= n._label ); }
-         /// Larger-than-or-equal-to operator (compares only labels)
-         bool operator >= ( const Var& n ) const { return( _label >= n._label ); }
-         /// Not-equal-to operator (compares only labels)
-         bool operator != ( const Var& n ) const { return( _label != n._label ); }
-         /// Equal-to operator (compares only labels)
-         bool operator == ( const Var& n ) const { return( _label == n._label ); }
+         /// Smaller-than-or-equal-to operator (only compares labels)
+         bool operator <= ( const Var& n ) const { 
+ #ifdef DAI_DEBUG
+             if( _label == n._label )
+                 assert( _states == n._states );
+ #endif
+             return( _label <= n._label ); 
+         }
+         /// Larger-than-or-equal-to operator (only compares labels)
+         bool operator >= ( const Var& n ) const { 
+ #ifdef DAI_DEBUG
+             if( _label == n._label )
+                 assert( _states == n._states );
+ #endif
+             return( _label >= n._label ); 
+         }
+         /// Not-equal-to operator (only compares labels)
+         bool operator != ( const Var& n ) const { 
+ #ifdef DAI_DEBUG
+             if( _label == n._label )
+                 assert( _states == n._states );
+ #endif
+             return( _label != n._label ); 
+         }
+         /// Equal-to operator (only compares labels)
+         bool operator == ( const Var& n ) const { 
+ #ifdef DAI_DEBUG
+             if( _label == n._label )
+                 assert( _states == n._states );
+ #endif            
+             return( _label == n._label ); 
+         }
  
          /// Writes a Var to an output stream
          friend std::ostream& operator << ( std::ostream& os, const Var& n ) {
diff --cc src/bp.cpp
@@@ -375,8 -410,8 +416,8 @@@ Factor BP::belief( const VarSet &ns ) c
  }
  
  
- Factor BP::beliefF (size_t I) const {
+ Factor BP::beliefF( size_t I ) const {
 -    if( 0 == 1 ) {
 +    if( !DAI_BP_FAST ) {
          /*  UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION */
  
          Factor prod( factor(I) );