include/dai/mr.h

   1 /*  This file is part of libDAI - http://www.libdai.org/
   2  *
   3  *  libDAI is licensed under the terms of the GNU General Public License version
   4  *  2, or (at your option) any later version. libDAI is distributed without any
   5  *  warranty. See the file COPYING for more details.
   6  *
   7  *  Copyright (C) 2007       Bastian Wemmenhove
   8  *  Copyright (C) 2007-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   9  *  Copyright (C) 2007       Radboud University Nijmegen, The Netherlands
  10  */
  11
  12
  13 /// \file
  14 /// \brief Defines class MR, which implements loop corrections as proposed by Montanari and Rizzo
  15
  16
  17 #ifndef __defined_libdai_mr_h
  18 #define __defined_libdai_mr_h
  19
  20
  21 #include <vector>
  22 #include <string>
  23 #include <dai/factorgraph.h>
  24 #include <dai/daialg.h>
  25 #include <dai/enum.h>
  26 #include <dai/properties.h>
  27 #include <dai/exceptions.h>
  28 #include <boost/dynamic_bitset.hpp>
  29
  30
  31 namespace dai {
  32
  33
  34 /// Approximate inference algorithm by Montanari and Rizzo [\ref MoR05]
  35 /** \author Bastian Wemmenhove wrote the original implementation before it was merged into libDAI
  36  *  \todo Clean up code (use a BipartiteGraph-like implementation for the graph structure, and BBP for response propagation).
  37  */
  38 class MR : public DAIAlgFG {
  39     private:
  40         /// Is the underlying factor graph supported?
  41         bool supported;
  42         /// con[i] = connectivity of spin \a i
  43         std::vector<size_t>                             con;
  44         /// nb[i] are the neighbours of spin \a i
  45         std::vector<std::vector<size_t> >               nb;
  46         /// tJ[i][_j] is the hyperbolic tangent of the interaction between spin \a i and its neighbour nb[i][_j]
  47         std::vector<std::vector<Real> >                 tJ;
  48         /// theta[i] is the local field on spin \a i
  49         std::vector<Real>                               theta;
  50         /// M[i][_j] is \f$ M^{(i)}_j \f$
  51         std::vector<std::vector<Real> >                 M;
  52         /// The \a _j 'th neighbour of spin \a i has spin \a i as its kindex[i][_j]'th neighbour
  53         std::vector<std::vector<size_t> >               kindex;
  54         /// Cavity correlations
  55         std::vector<std::vector<std::vector<Real> > >   cors;
  56         /// Maximum connectivity
  57         static const size_t kmax = 31;
  58         /// Type used for managing a subset of neighbors
  59         typedef boost::dynamic_bitset<> sub_nb;
  60         /// Number of variables (spins)
  61         size_t N;
  62         /// Magnetizations
  63         std::vector<Real> Mag;
  64         /// Maximum difference encountered so far
  65         Real _maxdiff;
  66         /// Number of iterations needed
  67         size_t _iters;
  68
  69     public:
  70         /// Parameters for MR
  71         struct Properties {
  72             /// Enumeration of different types of update equations
  73             /** The possible update equations are:
  74              *  - FULL full updates, slow but accurate
  75              *  - LINEAR linearized updates, faster but less accurate
  76              */
  77             DAI_ENUM(UpdateType,FULL,LINEAR);
  78
  79             /// Enumeration of different ways of initializing the cavity correlations
  80             /** The possible cavity initializations are:
  81              *  - RESPPROP using response propagation ("linear response")
  82              *  - CLAMPING using clamping and BP
  83              *  - EXACT using JunctionTree
  84              */
  85             DAI_ENUM(InitType,RESPPROP,CLAMPING,EXACT);
  86
  87             /// Verbosity (amount of output sent to stderr)
  88             size_t verbose;
  89
  90             /// Tolerance for convergence test
  91             Real tol;
  92
  93             /// Update equations
  94             UpdateType updates;
  95
  96             /// How to initialize the cavity correlations
  97             InitType inits;
  98         } props;
  99
 100         /// Name of this inference method
 101         static const char *Name;
 102
 103     public:
 104         /// Default constructor
 105         MR() : DAIAlgFG(), supported(), con(), nb(), tJ(), theta(), M(), kindex(), cors(), N(), Mag(), _maxdiff(), _iters(), props() {}
 106
 107         /// Construct from FactorGraph \a fg and PropertySet \a opts
 108         /** \param opts Parameters @see Properties
 109          *  \note This implementation only deals with binary variables and pairwise interactions.
 110          *  \throw NOT_IMPLEMENTED if \a fg has factors depending on three or more variables or has variables with more than two possible states.
 111          */
 112         MR( const FactorGraph &fg, const PropertySet &opts );
 113
 114
 115     /// \name General InfAlg interface
 116     //@{
 117         virtual MR* clone() const { return new MR(*this); }
 118         virtual std::string identify() const;
 119         virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
 120         virtual Factor belief( const VarSet &/*vs*/ ) const { DAI_THROW(NOT_IMPLEMENTED); return Factor(); }
 121         virtual Factor beliefV( size_t i ) const;
 122         virtual std::vector<Factor> beliefs() const;
 123         virtual Real logZ() const { DAI_THROW(NOT_IMPLEMENTED); return 0.0; }
 124         virtual void init() {}
 125         virtual void init( const VarSet &/*ns*/ ) { DAI_THROW(NOT_IMPLEMENTED); }
 126         virtual Real run();
 127         virtual Real maxDiff() const { return _maxdiff; }
 128         virtual size_t Iterations() const { return _iters; }
 129         virtual void setProperties( const PropertySet &opts );
 130         virtual PropertySet getProperties() const;
 131         virtual std::string printProperties() const;
 132     //@}
 133
 134     private:
 135         /// Returns the signum of \a a
 136         Real sign(Real a) { return (a >= 0) ? 1.0 : -1.0; }
 137
 138         /// Initialize N, con, nb, tJ, theta
 139         void init(size_t Nin, Real *_w, Real *_th);
 140
 141         /// Initialize kindex
 142         void makekindex();
 143
 144         /// Initialize cors
 145         Real init_cor();
 146
 147         /// Calculate cors using response propagation
 148         Real init_cor_resp();
 149
 150         /// Iterate update equations for cavity fields
 151         void solvemcav();
 152
 153         /// Calculate magnetizations
 154         void solveM();
 155
 156         /// Calculate the product of all tJ[i][_j] for _j in A
 157         /** \param i variable index
 158          *  \param A subset of neighbors of variable \a i
 159          */
 160         Real _tJ(size_t i, sub_nb A);
 161
 162         /// Calculate \f$ \Omega^{(i)}_{j,l} \f$ as defined in [\ref MoR05] eqn. (2.15)
 163         Real Omega(size_t i, size_t _j, size_t _l);
 164
 165         /// Calculate \f$ T^{(i)}_A \f$ as defined in [\ref MoR05] eqn. (2.17) with \f$ A = \{l_1,l_2,\dots\} \f$
 166         /** \param i variable index
 167          *  \param A subset of neighbors of variable \a i
 168          */
 169         Real T(size_t i, sub_nb A);
 170
 171         /// Calculates \f$ T^{(i)}_j \f$ where \a j is the \a _j 'th neighbor of \a i
 172         Real T(size_t i, size_t _j);
 173
 174         /// Calculates \f$ \Gamma^{(i)}_{j,l_1l_2} \f$ as defined in [\ref MoR05] eqn. (2.16)
 175         Real Gamma(size_t i, size_t _j, size_t _l1, size_t _l2);
 176
 177         /// Calculates \f$ \Gamma^{(i)}_{l_1l_2} \f$ as defined in [\ref MoK07] on page 1141
 178         Real Gamma(size_t i, size_t _l1, size_t _l2);
 179
 180         /// Approximates moments of variables in \a A
 181         /** Calculate the moment of variables in \a A from M and cors, neglecting higher order cumulants,
 182          *  defined as the sum over all partitions of A into subsets of cardinality two at most of the
 183          *  product of the cumulants (either first order, i.e. M, or second order, i.e. cors) of the
 184          *  entries of the partitions.
 185          *
 186          *  \param i variable index
 187          *  \param A subset of neighbors of variable \a i
 188          */
 189         Real appM(size_t i, sub_nb A);
 190
 191         /// Calculate sum over all even/odd subsets B of \a A of _tJ(j,B) appM(j,B)
 192         /** \param j variable index
 193          *  \param A subset of neighbors of variable \a j
 194          *  \param sum_even on return, will contain the sum over all even subsets
 195          *  \param sum_odd on return, will contain the sum over all odd subsets
 196          */
 197         void sum_subs(size_t j, sub_nb A, Real *sum_even, Real *sum_odd);
 198 };
 199
 200
 201 } // end of namespace dai
 202
 203
 204 #endif