Changed license from GPL v2+ to FreeBSD (aka BSD 2-clause) license
[libdai.git] / include / dai / mr.h
index a6050b6..62bb955 100644 (file)
@@ -1,18 +1,13 @@
 /*  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) 2006-2011, The libDAI authors. All rights reserved.
  *
- *  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
+ *  Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
  */
 
 
 /// \file
-/// \brief Defines class MR
-/// \todo Improve documentation
+/// \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 <dai/graph.h>
 #include <boost/dynamic_bitset.hpp>
 
 
 namespace dai {
 
 
-/// Approximate inference algorithm by Montanari and Rizzo
+/// Approximate inference algorithm by Montanari and Rizzo [\ref MoR05]
+/** \author Bastian Wemmenhove wrote the original implementation before it was merged into libDAI
+ */
 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<Real> >                 tJ;        // tJ[i][_j] is the tanh of the interaction between spin i and its neighbour nb[i][_j]
-        std::vector<Real>                               theta;     // theta[i] is the local field on spin i
-        std::vector<std::vector<Real> >                 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<Real> > >   cors;
+        /// Is the underlying factor graph supported?
+        bool supported;
 
-        static const size_t kmax = 31;
-        typedef boost::dynamic_bitset<> sub_nb;
+        /// The interaction graph (Markov graph)
+        GraphAL G;
 
-        size_t N;
+        /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour G.nb(i,_j)
+        std::vector<std::vector<Real> >                 tJ;
+        /// theta[i] is the local field on spin \a i
+        std::vector<Real>                               theta;
 
-        std::vector<Real> Mag;
+        /// M[i][_j] is \f$ M^{(i)}_j \f$
+        std::vector<std::vector<Real> >                 M;
+        /// Cavity correlations
+        std::vector<std::vector<std::vector<Real> > >   cors;
 
+        /// Type used for managing a subset of neighbors
+        typedef boost::dynamic_bitset<> sub_nb;
+        
+        /// Magnetizations
+        std::vector<Real> Mag;
+        
+        /// Maximum difference encountered so far
         Real _maxdiff;
+        
+        /// Number of iterations needed
         size_t _iters;
 
     public:
-        /// Parameters of this inference algorithm
+        /// 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
+            /// Verbosity (amount of output sent to stderr)
             size_t verbose;
 
-            /// Tolerance
+            /// Tolerance for convergence test
             Real tol;
 
             /// Update equations
@@ -77,23 +92,26 @@ class MR : public DAIAlgFG {
             InitType inits;
         } props;
 
-        /// Name of this inference method
-        static const char *Name;
-
     public:
         /// Default constructor
-        MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
+        MR() : DAIAlgFG(), supported(), G(), tJ(), theta(), M(), cors(), Mag(), _maxdiff(), _iters(), props() {}
 
-        /// Construct from FactorGraph fg and PropertySet opts
+        /// Construct from FactorGraph \a fg and PropertySet \a opts
+        /** \param fg Factor graph.
+         *  \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 &n ) const;
-        virtual Factor belief( const VarSet &/*ns*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
+        virtual std::string name() const { return "MR"; }
+        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 std::vector<Factor> beliefs() const;
         virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
         virtual void init() {}
@@ -101,41 +119,63 @@ class MR : public DAIAlgFG {
         virtual Real run();
         virtual Real maxDiff() const { return _maxdiff; }
         virtual size_t Iterations() const { return _iters; }
-    //@}
-
-    /// \name Managing parameters (which are stored in MR::props)
-    //@{
-        /// Set parameters of this inference algorithm.
-        /** The parameters are set according to \a opts. 
-         *  The values can be stored either as std::string or as the type of the corresponding MR::props member.
-         */
-        void setProperties( const PropertySet &opts );
-        /// Returns parameters of this inference algorithm converted into a PropertySet.
-        PropertySet getProperties() const;
-        /// Returns parameters of this inference algorithm formatted as a string in the format "[key1=val1,key2=val2,...,keyn=valn]".
-        std::string printProperties() const;
+        virtual void setProperties( const PropertySet &opts );
+        virtual PropertySet getProperties() const;
+        virtual std::string printProperties() const;
     //@}
 
     private:
-        void init(size_t Nin, Real *_w, Real *_th);
-        void makekindex();
-        void init_cor();
-        Real init_cor_resp();
-        void solvemcav();
-        void solveM();
-
+        /// Initialize cors
+        Real calcCavityCorrelations();
+        
+        /// Iterate update equations for cavity fields
+        void propagateCavityFields();
+        
+        /// Calculate magnetizations
+        void calcMagnetizations();
+
+        /// 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);
-
+        
+        /// 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);
-
-        Real sign(Real a) { return (a >= 0) ? 1.0 : -1.0; }
 };