include/dai/bp.h

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