include/dai/bp.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) 2006-2009  Joris Mooij  [joris dot mooij at libdai dot org]
   8  *  Copyright (C) 2006-2007  Radboud University Nijmegen, The Netherlands
   9  */
  10
  11
  12 /// \file
  13 /// \brief Defines class BP, which implements (Loopy) Belief Propagation
  14
  15
  16 #ifndef __defined_libdai_bp_h
  17 #define __defined_libdai_bp_h
  18
  19
  20 #include <string>
  21 #include <dai/daialg.h>
  22 #include <dai/factorgraph.h>
  23 #include <dai/properties.h>
  24 #include <dai/enum.h>
  25
  26
  27 namespace dai {
  28
  29
  30 /// Approximate inference algorithm "(Loopy) Belief Propagation"
  31 /** The Loopy Belief Propagation algorithm uses message passing
  32  *  to approximate marginal probability distributions ("beliefs") for variables
  33  *  and factors (more precisely, for the subset of variables depending on the factor).
  34  *  There are two variants, the sum-product algorithm (corresponding to
  35  *  finite temperature) and the max-product algorithm (corresponding to
  36  *  zero temperature).
  37  *
  38  *  The messages \f$m_{I\to i}(x_i)\f$ are passed from factors \f$I\f$ to variables \f$i\f$.
  39  *  In case of the sum-product algorith, the update equation is:
  40  *    \f[ m_{I\to i}(x_i) \propto \sum_{x_{N_I\setminus\{i\}}} f_I(x_I) \prod_{j\in N_I\setminus\{i\}} \prod_{J\in N_j\setminus\{I\}} m_{J\to j}\f]
  41  *  and in case of the max-product algorithm:
  42  *    \f[ m_{I\to i}(x_i) \propto \max_{x_{N_I\setminus\{i\}}} f_I(x_I) \prod_{j\in N_I\setminus\{i\}} \prod_{J\in N_j\setminus\{I\}} m_{J\to j}\f]
  43  *  In order to improve convergence, the updates can be damped. For improved numerical stability,
  44  *  the updates can be done in the log-domain alternatively.
  45  *
  46  *  After convergence, the variable beliefs are calculated by:
  47  *    \f[ b_i(x_i) \propto \prod_{I\in N_i} m_{I\to i}(x_i)\f]
  48  *  and the factor beliefs are calculated by:
  49  *    \f[ b_I(x_I) \propto f_I(x_I) \prod_{j\in N_I} \prod_{J\in N_j\setminus\{I\}} m_{J\to j}(x_j) \f]
  50  *  The logarithm of the partition sum is calculated by:
  51  *    \f[ \log Z = \sum_i (1 - |N_i|) \sum_{x_i} b_i(x_i) \log b_i(x_i) - \sum_I \sum_{x_I} b_I(x_I) \log \frac{b_I(x_I)}{f_I(x_I)} \f]
  52  *
  53  *  For the max-product algorithm, a heuristic way of finding the MAP state (the
  54  *  joint configuration of all variables which has maximum probability) is provided
  55  *  by the findMaximum() method, which can be called after convergence.
  56  *
  57  *  \note There are two implementations, an optimized one (the default) which caches IndexFor objects,
  58  *  and a slower, less complicated one which is easier to maintain/understand. The slower one can be
  59  *  enabled by defining DAI_BP_FAST as false in the source file.
  60  */
  61 class BP : public DAIAlgFG {
  62     protected:
  63         /// Type used for index cache
  64         typedef std::vector<size_t> ind_t;
  65         /// Type used for storing edge properties
  66         struct EdgeProp {
  67             /// Index cached for this edge
  68             ind_t  index;
  69             /// Old message living on this edge
  70             Prob   message;
  71             /// New message living on this edge
  72             Prob   newMessage;
  73             /// Residual for this edge
  74             Real   residual;
  75         };
  76         /// Stores all edge properties
  77         std::vector<std::vector<EdgeProp> > _edges;
  78         /// Type of lookup table (only used for maximum-residual BP)
  79         typedef std::multimap<Real, std::pair<std::size_t, std::size_t> > LutType;
  80         /// Lookup table (only used for maximum-residual BP)
  81         std::vector<std::vector<LutType::iterator> > _edge2lut;
  82         /// Lookup table (only used for maximum-residual BP)
  83         LutType _lut;
  84         /// Maximum difference between variable beliefs encountered so far
  85         Real _maxdiff;
  86         /// Number of iterations needed
  87         size_t _iters;
  88         /// The history of message updates (only recorded if \a recordSentMessages is \c true)
  89         std::vector<std::pair<std::size_t, std::size_t> > _sentMessages;
  90
  91     public:
  92         /// Parameters for BP
  93         struct Properties {
  94             /// Enumeration of possible update schedules
  95             /** The following update schedules have been defined:
  96              *  - PARALL parallel updates
  97              *  - SEQFIX sequential updates using a fixed sequence
  98              *  - SEQRND sequential updates using a random sequence
  99              *  - SEQMAX maximum-residual updates [\ref EMK06]
 100              */
 101             DAI_ENUM(UpdateType,SEQFIX,SEQRND,SEQMAX,PARALL);
 102
 103             /// Enumeration of inference variants
 104             /** There are two inference variants:
 105              *  - SUMPROD Sum-Product
 106              *  - MAXPROD Max-Product (equivalent to Min-Sum)
 107              */
 108             DAI_ENUM(InfType,SUMPROD,MAXPROD);
 109
 110             /// Verbosity (amount of output sent to stderr)
 111             size_t verbose;
 112
 113             /// Maximum number of iterations
 114             size_t maxiter;
 115
 116             /// Tolerance for convergence test
 117             Real tol;
 118
 119             /// Whether updates should be done in logarithmic domain or not
 120             bool logdomain;
 121
 122             /// Damping constant (0.0 means no damping, 1.0 is maximum damping)
 123             Real damping;
 124
 125             /// Message update schedule
 126             UpdateType updates;
 127
 128             /// Inference variant
 129             InfType inference;
 130         } props;
 131
 132         /// Name of this inference algorithm
 133         static const char *Name;
 134
 135         /// Specifies whether the history of message updates should be recorded
 136         bool recordSentMessages;
 137
 138     public:
 139     /// \name Constructors/destructors
 140     //@{
 141         /// Default constructor
 142         BP() : DAIAlgFG(), _edges(), _edge2lut(), _lut(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {}
 143
 144         /// Construct from FactorGraph \a fg and PropertySet \a opts
 145         /** \param opts Parameters @see Properties
 146          */
 147         BP( const FactorGraph & fg, const PropertySet &opts ) : DAIAlgFG(fg), _edges(), _maxdiff(0.0), _iters(0U), _sentMessages(), props(), recordSentMessages(false) {
 148             setProperties( opts );
 149             construct();
 150         }
 151
 152         /// Copy constructor
 153         BP( const BP &x ) : DAIAlgFG(x), _edges(x._edges), _edge2lut(x._edge2lut),
 154             _lut(x._lut), _maxdiff(x._maxdiff), _iters(x._iters), _sentMessages(x._sentMessages),
 155             props(x.props), recordSentMessages(x.recordSentMessages)
 156         {
 157             for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
 158                 _edge2lut[l->second.first][l->second.second] = l;
 159         }
 160
 161         /// Assignment operator
 162         BP& operator=( const BP &x ) {
 163             if( this != &x ) {
 164                 DAIAlgFG::operator=( x );
 165                 _edges = x._edges;
 166                 _lut = x._lut;
 167                 for( LutType::iterator l = _lut.begin(); l != _lut.end(); ++l )
 168                     _edge2lut[l->second.first][l->second.second] = l;
 169                 _maxdiff = x._maxdiff;
 170                 _iters = x._iters;
 171                 _sentMessages = x._sentMessages;
 172                 props = x.props;
 173                 recordSentMessages = x.recordSentMessages;
 174             }
 175             return *this;
 176         }
 177     //@}
 178
 179     /// \name General InfAlg interface
 180     //@{
 181         virtual BP* clone() const { return new BP(*this); }
 182         virtual std::string identify() const;
 183         virtual Factor belief( const Var &v ) const { return beliefV( findVar( v ) ); }
 184         virtual Factor belief( const VarSet &vs ) const;
 185         virtual Factor beliefV( size_t i ) const;
 186         virtual Factor beliefF( size_t I ) const;
 187         virtual std::vector<Factor> beliefs() const;
 188         virtual Real logZ() const;
 189         virtual void init();
 190         virtual void init( const VarSet &ns );
 191         virtual Real run();
 192         virtual Real maxDiff() const { return _maxdiff; }
 193         virtual size_t Iterations() const { return _iters; }
 194         virtual void setProperties( const PropertySet &opts );
 195         virtual PropertySet getProperties() const;
 196         virtual std::string printProperties() const;
 197     //@}
 198
 199     /// \name Additional interface specific for BP
 200     //@{
 201         /// Calculates the joint state of all variables that has maximum probability
 202         /** \pre Assumes that run() has been called and that \a props.inference == \c MAXPROD
 203          */
 204         std::vector<std::size_t> findMaximum() const;
 205
 206         /// Returns history of which messages have been updated
 207         const std::vector<std::pair<std::size_t, std::size_t> >& getSentMessages() const {
 208             return _sentMessages;
 209         }
 210
 211         /// Clears history of which messages have been updated
 212         void clearSentMessages() { _sentMessages.clear(); }
 213     //@}
 214
 215     protected:
 216         /// Returns constant reference to message from the \a _I 'th neighbor of variable \a i to variable \a i
 217         const Prob & message(size_t i, size_t _I) const { return _edges[i][_I].message; }
 218         /// Returns reference to message from the \a _I 'th neighbor of variable \a i to variable \a i
 219         Prob & message(size_t i, size_t _I) { return _edges[i][_I].message; }
 220         /// Returns constant reference to updated message from the \a _I 'th neighbor of variable \a i to variable \a i
 221         const Prob & newMessage(size_t i, size_t _I) const { return _edges[i][_I].newMessage; }
 222         /// Returns reference to updated message from the \a _I 'th neighbor of variable \a i to variable \a i
 223         Prob & newMessage(size_t i, size_t _I) { return _edges[i][_I].newMessage; }
 224         /// Returns constant reference to cached index for the edge between variable \a i and its \a _I 'th neighbor
 225         const ind_t & index(size_t i, size_t _I) const { return _edges[i][_I].index; }
 226         /// Returns reference to cached index for the edge between variable \a i and its \a _I 'th neighbor
 227         ind_t & index(size_t i, size_t _I) { return _edges[i][_I].index; }
 228         /// Returns constant reference to residual for the edge between variable \a i and its \a _I 'th neighbor
 229         const Real & residual(size_t i, size_t _I) const { return _edges[i][_I].residual; }
 230         /// Returns reference to residual for the edge between variable \a i and its \a _I 'th neighbor
 231         Real & residual(size_t i, size_t _I) { return _edges[i][_I].residual; }
 232
 233         /// Calculate the updated message from the \a _I 'th neighbor of variable \a i to variable \a i
 234         virtual void calcNewMessage( size_t i, size_t _I );
 235         /// Replace the "old" message from the \a _I 'th neighbor of variable \a i to variable \a i by the "new" (updated) message
 236         void updateMessage( size_t i, size_t _I );
 237         /// Set the residual (difference between new and old message) for the edge between variable \a i and its \a _I 'th neighbor to \a r
 238         void updateResidual( size_t i, size_t _I, Real r );
 239         /// Finds the edge which has the maximum residual (difference between new and old message)
 240         void findMaxResidual( size_t &i, size_t &_I );
 241         /// Calculates unnormalized belief of variable \a i
 242         virtual void calcBeliefV( size_t i, Prob &p ) const;
 243         /// Calculates unnormalized belief of factor \a I
 244         virtual void calcBeliefF( size_t I, Prob &p ) const;
 245
 246         /// Helper function for constructors
 247         virtual void construct();
 248 };
 249
 250
 251 } // end of namespace dai
 252
 253
 254 #endif