Cleaned up variable elimination code in ClusterGraph
[libdai.git] / include / dai / mr.h
index 2b55e96..ce8a7c6 100644 (file)
@@ -1,22 +1,17 @@
-/*  Copyright (C) 2006-2008  Joris Mooij  [j dot mooij at science dot ru dot nl]
-    Radboud University Nijmegen, The Netherlands
-    
-    This file is part of libDAI.
+/*  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) 2007       Bastian Wemmenhove
+ *  Copyright (C) 2007-2009  Joris Mooij         [joris dot mooij at libdai dot org]
+ *  Copyright (C) 2007       Radboud University Nijmegen, The Netherlands
+ */
 
-    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
-*/
+/// \file
+/// \brief Defines class MR, which implements loop corrections as proposed by Montanari and Rizzo
 
 
 #ifndef __defined_libdai_mr_h
 #include <dai/enum.h>
 #include <dai/properties.h>
 #include <dai/exceptions.h>
+#include <boost/dynamic_bitset.hpp>
 
 
 namespace dai {
 
 
-class sub_nb;
-
-
+/// Approximate inference algorithm by Montanari and Rizzo [\ref MoR05]
+/** \author Bastian Wemmenhove wrote the original implementation before it was merged into libDAI
+ *  \todo Clean up code (use a BipartiteGraph-like implementation for the graph structure)
+ */
 class MR : public DAIAlgFG {
     private:
-        bool supported;                                            // is the underlying factor graph supported?
-
-        std::vector<size_t>                             con;       // con[i] = connectivity of spin i
-        std::vector<std::vector<size_t> >               nb;        // nb[i] are the neighbours of spin i
-        std::vector<std::vector<double> >               tJ;        // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
-        std::vector<double>                             theta;     // theta[i] is the local field on spin i
-        std::vector<std::vector<double> >               M;         // M[i][_j] is M^{(i)}_j
-        std::vector<std::vector<size_t> >               kindex;    // the _j'th neighbour of spin i has spin i as its kindex[i][_j]'th neighbour
-        std::vector<std::vector<std::vector<double> > > cors;
-    
+        /// Is the underlying factor graph supported?
+        bool supported;
+        /// con[i] = connectivity of spin \a i
+        std::vector<size_t>                             con;
+        /// nb[i] are the neighbours of spin \a i
+        std::vector<std::vector<size_t> >               nb;
+        /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour nb[i][_j]
+        std::vector<std::vector<Real> >                 tJ;
+        /// theta[i] is the local field on spin \a i
+        std::vector<Real>                               theta;
+        /// M[i][_j] is \f$ M^{(i)}_j \f$
+        std::vector<std::vector<Real> >                 M;
+        /// The \a _j 'th neighbour of spin \a i has spin \a i as its kindex[i][_j]'th neighbour
+        std::vector<std::vector<size_t> >               kindex;
+        /// Cavity correlations
+        std::vector<std::vector<std::vector<Real> > >   cors;
+        /// Maximum connectivity
         static const size_t kmax = 31;
-        
+        /// Type used for managing a subset of neighbors
+        typedef boost::dynamic_bitset<> sub_nb;
+        /// Number of variables (spins)
         size_t N;
-
-        std::vector<double> Mag;
+        /// Magnetizations
+        std::vector<Real> Mag;
+        /// Maximum difference encountered so far
+        Real _maxdiff;
+        /// Number of iterations needed
+        size_t _iters;
 
     public:
+        /// Parameters for MR
         struct Properties {
+            /// Enumeration of different types of update equations
+            /** The possible update equations are:
+             *  - FULL full updates, slow but accurate
+             *  - LINEAR linearized updates, faster but less accurate
+             */
+            DAI_ENUM(UpdateType,FULL,LINEAR);
+
+            /// Enumeration of different ways of initializing the cavity correlations
+            /** The possible cavity initializations are:
+             *  - RESPPROP using response propagation ("linear response")
+             *  - CLAMPING using clamping and BP
+             *  - EXACT using JunctionTree
+             */
+            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
+
+            /// Verbosity (amount of output sent to stderr)
             size_t verbose;
-            double tol;
-            DAI_ENUM(UpdateType,FULL,LINEAR)
-            DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT)
+
+            /// Tolerance for convergence test
+            Real tol;
+
+            /// Update equations
             UpdateType updates;
+
+            /// How to initialize the cavity correlations
             InitType inits;
         } props;
-        double maxdiff;
+
+        /// Name of this inference method
+        static const char *Name;
 
     public:
-        MR() {}
-        MR( const FactorGraph & fg, const PropertySet &opts );
-        void init(size_t Nin, double *_w, double *_th);
+        /// Default constructor
+        MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
+
+        /// Construct from FactorGraph \a fg and PropertySet \a opts
+        /** \param opts Parameters @see Properties
+         *  \note This implementation only deals with binary variables and pairwise interactions.
+         *  \throw NOT_IMPLEMENTED if \a fg has factors depending on three or more variables or has variables with more than two possible states.
+         */
+        MR( const FactorGraph &fg, const PropertySet &opts );
+
+
+    /// \name General InfAlg interface
+    //@{
+        virtual MR* clone() const { return new MR(*this); }
+        virtual std::string identify() const;
+        virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
+        virtual Factor belief( const VarSet &/*vs*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
+        virtual Factor beliefV( size_t i ) 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*/ ) { DAI_THROW(NOT_IMPLEMENTED); }
+        virtual Real run();
+        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;
+    //@}
+
+    private:
+        /// Returns the signum of \a a
+        Real sign(Real a) { return (a >= 0) ? 1.0 : -1.0; }
+
+        /// Initialize N, con, nb, tJ, theta
+        void init(size_t Nin, Real *_w, Real *_th);
+        
+        /// Initialize kindex
         void makekindex();
-        void read_files();
+        
+        /// Initialize cors
         void init_cor();
-        double init_cor_resp();
+        
+        /// Calculate cors using response propagation
+        Real init_cor_resp();
+        
+        /// Iterate update equations for cavity fields
         void solvemcav();
+        
+        /// Calculate magnetizations
         void solveM();
-        double run();
-        Factor belief( const Var &n ) const;
-        Factor belief( const VarSet &/*ns*/ ) const { 
-            DAI_THROW(NOT_IMPLEMENTED);
-            return Factor(); 
-        }
-        std::vector<Factor> beliefs() const;
-        Real logZ() const { 
-            DAI_THROW(NOT_IMPLEMENTED);
-            return 0.0; 
-        }
-        void init() {}
-        /// Clear messages and beliefs corresponding to the nodes in ns
-        virtual void init( const VarSet &/*ns*/ ) {
-            DAI_THROW(NOT_IMPLEMENTED);
-        }
-        static const char *Name;
-        std::string identify() const;
-        double _tJ(size_t i, sub_nb A);
-
-        double Omega(size_t i, size_t _j, size_t _l);
-        double T(size_t i, sub_nb A);
-        double T(size_t i, size_t _j);
-        double Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
-        double Gamma(size_t i, size_t _l1, size_t _l2);
-
-        double appM(size_t i, sub_nb A);
-        void sum_subs(size_t j, sub_nb A, double *sum_even, double *sum_odd);
-
-        double sign(double a) { return (a >= 0) ? 1.0 : -1.0; }
-        MR* clone() const { return new MR(*this); }
-        /// Create (virtual constructor)
-        virtual MR* create() const { return new MR(); }
 
-        void setProperties( const PropertySet &opts );
-        PropertySet getProperties() const;
-        std::string printProperties() const;
-        double maxDiff() const { return maxdiff; }
-}; 
-
-
-// represents a subset of nb[i] as a binary number
-// the elements of a subset should be thought of as indices in nb[i]
-class sub_nb {
-    private:
-        size_t s;
-        size_t bits;
-    
-    public:
-        // construct full subset containing nr_elmt elements
-        sub_nb(size_t nr_elmt) {
-#ifdef DAI_DEBUG
-            assert( nr_elmt < sizeof(size_t) / sizeof(char) * 8 );
-#endif
-            bits = nr_elmt;
-            s = (1U << bits) - 1;
-        }
-
-        // copy constructor
-        sub_nb( const sub_nb & x ) : s(x.s), bits(x.bits) {}
-
-        // assignment operator 
-        sub_nb & operator=( const sub_nb &x ) {
-            if( this != &x ) {
-                s = x.s;
-                bits = x.bits;
-            }
-            return *this;
-        }
-
-        // returns number of elements
-        size_t size() {
-            size_t size = 0;
-            for(size_t bit = 0; bit < bits; bit++)
-                if( s & (1U << bit) )
-                    size++;
-            return size;
-        }
-
-        // increases s by one (for enumeration in lexicographical order)
-        sub_nb operator++() { 
-            s++; 
-            if( s >= (1U << bits) )
-                s = 0;
-            return *this; 
-        }
+        /// Calculate the product of all tJ[i][_j] for _j in A
+        /** \param i variable index
+         *  \param A subset of neighbors of variable \a i
+         */
+        Real _tJ(size_t i, sub_nb A);
         
-        // return i'th element of this subset
-        size_t operator[](size_t i) { 
-            size_t bit;
-            for(bit = 0; bit < bits; bit++ )
-                if( s & (1U << bit) ) {
-                    if( i == 0 )
-                        break;
-                    else
-                        i--;
-                }
-#ifdef DAI_DEBUG
-            assert( bit < bits );
-#endif
-            return bit;
-        }
-
-        // add index _j to this subset
-        sub_nb &operator +=(size_t _j) {
-            s |= (1U << _j); 
-            return *this;
-        }
-
-        // return copy with index _j
-        sub_nb operator+(size_t _j) {
-            sub_nb x = *this;
-            x += _j;
-            return x;
-        }
-
-        // delete index _j from this subset
-        sub_nb &operator -=(size_t _j) {
-            s &= ~(1U << _j); 
-            return *this;
-        }
-
-        // return copy without index _j
-        sub_nb operator-(size_t _j) {
-            sub_nb x = *this;
-            x -= _j;
-            return x;
-        }
-
-        // empty this subset
-        sub_nb & clear() {
-            s = 0;
-            return *this;
-        }
-
-        // returns true if subset is empty
-        bool empty() { return (s == 0); }
-
-        // return 1 if _j is contained, 0 otherwise ("is element of")
-        size_t operator&(size_t _j) { return s & (1U << _j); }
-
-        friend std::ostream& operator<< (std::ostream& os, const sub_nb x) {
-            if( x.bits == 0 )
-                os << "empty";
-            else {
-                for(size_t bit = x.bits; bit > 0; bit-- )
-                    if( x.s & (1U << (bit-1)) )
-                        os << "1";
-                    else
-                        os << "0";
-            }
-            return os;
-        }
+        /// Calculate \f$ \Omega^{(i)}_{j,l} \f$ as defined in [\ref MoR05] eqn. (2.15)
+        Real Omega(size_t i, size_t _j, size_t _l);
+        
+        /// Calculate \f$ T^{(i)}_A \f$ as defined in [\ref MoR05] eqn. (2.17) with \f$ A = \{l_1,l_2,\dots\} \f$
+        /** \param i variable index
+         *  \param A subset of neighbors of variable \a i
+         */
+        Real T(size_t i, sub_nb A);
+        
+        /// Calculates \f$ T^{(i)}_j \f$ where \a j is the \a _j 'th neighbor of \a i
+        Real T(size_t i, size_t _j);
+        
+        /// Calculates \f$ \Gamma^{(i)}_{j,l_1l_2} \f$ as defined in [\ref MoR05] eqn. (2.16)
+        Real Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
+        
+        /// Calculates \f$ \Gamma^{(i)}_{l_1l_2} \f$ as defined in [\ref MoK07] on page 1141
+        Real Gamma(size_t i, size_t _l1, size_t _l2);
+        
+        /// Approximates moments of variables in \a A
+        /** Calculate the moment of variables in \a A from M and cors, neglecting higher order cumulants,
+         *  defined as the sum over all partitions of A into subsets of cardinality two at most of the
+         *  product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
+         *  entries of the partitions.
+         *
+         *  \param i variable index
+         *  \param A subset of neighbors of variable \a i
+         */
+        Real appM(size_t i, sub_nb A);
+        
+        /// Calculate sum over all even/odd subsets B of \a A of _tJ(j,B) appM(j,B)
+        /** \param j variable index
+         *  \param A subset of neighbors of variable \a j
+         *  \param sum_even on return, will contain the sum over all even subsets
+         *  \param sum_odd on return, will contain the sum over all odd subsets
+         */
+        void sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd);
 };