include/dai/bbp.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) 2009  Frederik Eaton [frederik at ofb dot net]
   8  */
   9
  10
  11 /// \file
  12 /// \brief Defines class BBP [\ref EaG09]
  13 /// \todo Fit more closely into libDAI framework
  14 /// \todo Improve documentation
  15 /// \author Frederik Eaton
  16
  17
  18 #ifndef ___defined_libdai_bbp_h
  19 #define ___defined_libdai_bbp_h
  20
  21
  22 #include <vector>
  23 #include <utility>
  24
  25 #include <dai/prob.h>
  26 #include <dai/daialg.h>
  27 #include <dai/factorgraph.h>
  28 #include <dai/enum.h>
  29 #include <dai/bp_dual.h>
  30
  31
  32 namespace dai {
  33
  34
  35 /// Computes the adjoint of the unnormed probability vector from the normalizer and the adjoint of the normalized probability vector @see eqn. (13) in [\ref EaG09]
  36 Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w );
  37
  38 /// Runs Gibbs sampling for \a iters iterations on ia.fg(), and returns state
  39 std::vector<size_t> getGibbsState( const InfAlg &ia, size_t iters );
  40
  41
  42 /// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
  43 /** \author Frederik Eaton
  44  */
  45 class BBP {
  46     protected:
  47     /// \name Inputs
  48     //@{
  49         BP_dual _bp_dual;
  50         const FactorGraph *_fg;
  51         const InfAlg *_ia;
  52     //@}
  53
  54         /// Number of iterations done
  55         size_t _iters;
  56
  57     /// \name Outputs
  58     //@{
  59         /// Variable factor adjoints
  60         std::vector<Prob> _adj_psi_V;
  61         /// Factor adjoints
  62         std::vector<Prob> _adj_psi_F;
  63         /// Variable->factor message adjoints (indexed [i][_I])
  64         std::vector<std::vector<Prob> > _adj_n;
  65         /// Factor->variable message adjoints (indexed [i][_I])
  66         std::vector<std::vector<Prob> > _adj_m;
  67         /// Normalized variable belief adjoints
  68         std::vector<Prob> _adj_b_V;
  69         /// Normalized factor belief adjoints
  70         std::vector<Prob> _adj_b_F;
  71     //@}
  72
  73     /// \name Helper quantities computed from the BP messages
  74     //@{
  75         /// _T[i][_I] (see eqn. (41) in [\ref EaG09])
  76         std::vector<std::vector<Prob > > _T;
  77         /// _U[I][_i] (see eqn. (42) in [\ref EaG09])
  78         std::vector<std::vector<Prob > > _U;
  79         /// _S[i][_I][_j] (see eqn. (43) in [\ref EaG09])
  80         std::vector<std::vector<std::vector<Prob > > > _S;
  81         /// _R[I][_i][_J] (see eqn. (44) in [\ref EaG09])
  82         std::vector<std::vector<std::vector<Prob > > > _R;
  83     //@}
  84
  85         /// Unnormalized variable belief adjoints
  86         std::vector<Prob> _adj_b_V_unnorm;
  87         /// Unnormalized factor belief adjoints
  88         std::vector<Prob> _adj_b_F_unnorm;
  89
  90         /// Initial variable factor adjoints
  91         std::vector<Prob> _init_adj_psi_V;
  92         /// Initial factor adjoints
  93         std::vector<Prob> _init_adj_psi_F;
  94
  95         /// Unnormalized variable->factor message adjoint (indexed [i][_I])
  96         std::vector<std::vector<Prob> > _adj_n_unnorm;
  97         /// Unnormalized factor->variable message adjoint (indexed [i][_I])
  98         std::vector<std::vector<Prob> > _adj_m_unnorm;
  99         /// Updated normalized variable->factor message adjoint (indexed [i][_I])
 100         std::vector<std::vector<Prob> > _new_adj_n;
 101         /// Updated normalized factor->variable message adjoint (indexed [i][_I])
 102         std::vector<std::vector<Prob> > _new_adj_m;
 103
 104     /// \name Optimized indexing (for performance)
 105     //@{
 106         /// Calculates _indices, which is a cache of IndexFor @see bp.cpp
 107         void RegenerateInds();
 108
 109         /// Index type
 110         typedef std::vector<size_t>  _ind_t;
 111         /// Cached indices (indexed [i][_I])
 112         std::vector<std::vector<_ind_t> >  _indices;
 113         /// Returns an index from the cache
 114         const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
 115     //@}
 116
 117     /// \name Initialization
 118     //@{
 119         /// Calculate T values; see eqn. (41) in [\ref EaG09]
 120         void RegenerateT();
 121         /// Calculate U values; see eqn. (42) in [\ref EaG09]
 122         void RegenerateU();
 123         /// Calculate S values; see eqn. (43) in [\ref EaG09]
 124         void RegenerateS();
 125         /// Calculate R values; see eqn. (44) in [\ref EaG09]
 126         void RegenerateR();
 127         /// Calculate _adj_b_V_unnorm and _adj_b_F_unnorm from _adj_b_V and _adj_b_F
 128         void RegenerateInputs();
 129         /// Initialise members for factor adjoints (call after RegenerateInputs)
 130         void RegeneratePsiAdjoints();
 131         /// Initialise members for message adjoints (call after RegenerateInputs) for parallel algorithm
 132         void RegenerateParMessageAdjoints();
 133         /// Initialise members for message adjoints (call after RegenerateInputs) for sequential algorithm
 134         /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
 135          *  than updating _adj_n (and friends) which are unused in the sequential algorithm.
 136          */
 137         void RegenerateSeqMessageAdjoints();
 138     //@}
 139
 140         /// Returns reference to T value; see eqn. (41) in [\ref EaG09]
 141         Prob & T(size_t i, size_t _I) { return _T[i][_I]; }
 142         /// Returns constant reference to T value; see eqn. (41) in [\ref EaG09]
 143         const Prob & T(size_t i, size_t _I) const { return _T[i][_I]; }
 144         /// Returns reference to U value; see eqn. (42) in [\ref EaG09]
 145         Prob & U(size_t I, size_t _i) { return _U[I][_i]; }
 146         /// Returns constant reference to U value; see eqn. (42) in [\ref EaG09]
 147         const Prob & U(size_t I, size_t _i) const { return _U[I][_i]; }
 148         /// Returns reference to S value; see eqn. (43) in [\ref EaG09]
 149         Prob & S(size_t i, size_t _I, size_t _j) { return _S[i][_I][_j]; }
 150         /// Returns constant reference to S value; see eqn. (43) in [\ref EaG09]
 151         const Prob & S(size_t i, size_t _I, size_t _j) const { return _S[i][_I][_j]; }
 152         /// Returns reference to R value; see eqn. (44) in [\ref EaG09]
 153         Prob & R(size_t I, size_t _i, size_t _J) { return _R[I][_i][_J]; }
 154         /// Returns constant reference to R value; see eqn. (44) in [\ref EaG09]
 155         const Prob & R(size_t I, size_t _i, size_t _J) const { return _R[I][_i][_J]; }
 156
 157     /// \name Parallel algorithm
 158     //@{
 159         /// Calculates new variable->factor message adjoint
 160         /** Increases variable factor adjoint according to eqn. (27) in [\ref EaG09] and
 161          *  calculates the new variable->factor message adjoint according to eqn. (29) in [\ref EaG09].
 162          */
 163         void calcNewN( size_t i, size_t _I );
 164         /// Calculates new factor->variable message adjoint
 165         /** Increases factor adjoint according to eqn. (28) in [\ref EaG09] and
 166          *  calculates the new factor->variable message adjoint according to the r.h.s. of eqn. (30) in [\ref EaG09].
 167          */
 168         void calcNewM( size_t i, size_t _I );
 169         /// Calculates unnormalized variable->factor message adjoint from the normalized one
 170         void calcUnnormMsgN( size_t i, size_t _I );
 171         /// Calculates unnormalized factor->variable message adjoint from the normalized one
 172         void calcUnnormMsgM( size_t i, size_t _I );
 173         /// Updates (un)normalized variable->factor message adjoints
 174         void upMsgN( size_t i, size_t _I );
 175         /// Updates (un)normalized factor->variable message adjoints
 176         void upMsgM( size_t i, size_t _I );
 177         /// Do one parallel update of all message adjoints
 178         void doParUpdate();
 179     //@}
 180
 181     /// \name Sequential algorithm
 182     //@{
 183         /// Helper function for sendSeqMsgM: increases factor->variable message adjoint by p and calculates the corresponding unnormalized adjoint
 184         void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
 185         //  DISABLED BECAUSE IT IS BUGGY:
 186         //  void updateSeqMsgM( size_t i, size_t _I );
 187         /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint
 188         void setSeqMsgM( size_t i, size_t _I, const Prob &p );
 189         /// Implements routine Send-n in Figure 5 in [\ref EaG09]
 190         void sendSeqMsgN( size_t i, size_t _I, const Prob &f );
 191         /// Implements routine Send-m in Figure 5 in [\ref EaG09]
 192         void sendSeqMsgM( size_t i, size_t _I );
 193     //@}
 194
 195         /// Calculates averaged L-1 norm of unnormalized message adjoints
 196         Real getUnMsgMag();
 197         /// Calculates averaged L-1 norms of current and new normalized message adjoints
 198         void getMsgMags( Real &s, Real &new_s );
 199
 200         /// Sets all vectors _adj_b_F to zero
 201         void zero_adj_b_F() {
 202             _adj_b_F.clear();
 203             _adj_b_F.reserve( _fg->nrFactors() );
 204             for( size_t I = 0; I < _fg->nrFactors(); I++ )
 205                 _adj_b_F.push_back( Prob( _fg->factor(I).states(), Real( 0.0 ) ) );
 206         }
 207
 208         /// Returns indices and magnitude of the largest normalized factor->variable message adjoint
 209         void getArgmaxMsgM( size_t &i, size_t &_I, Real &mag );
 210         /// Returns magnitude of the largest (in L1-norm) normalized factor->variable message adjoint
 211         Real getMaxMsgM();
 212         /// Calculates sum of L1 norms of all normalized factor->variable message adjoints
 213         Real getTotalMsgM();
 214         /// Calculates sum of L1 norms of all updated normalized factor->variable message adjoints
 215         Real getTotalNewMsgM();
 216         /// Calculates sum of L1 norms of all normalized variable->factor message adjoints
 217         Real getTotalMsgN();
 218
 219     public:
 220         /// Called by \a init, recalculates intermediate values
 221         void Regenerate();
 222
 223         /// Constructor
 224         BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
 225             props.set(opts);
 226         }
 227
 228         /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph fg
 229         std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
 230         /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph fg
 231         std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
 232
 233         /// Initializes belief adjoints and initial factor adjoints and regenerates
 234         void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F, const std::vector<Prob> &adj_psi_V, const std::vector<Prob> &adj_psi_F ) {
 235             _adj_b_V = adj_b_V;
 236             _adj_b_F = adj_b_F;
 237             _init_adj_psi_V = adj_psi_V;
 238             _init_adj_psi_F = adj_psi_F;
 239             Regenerate();
 240         }
 241
 242         /// Initializes belief adjoints and with zero initial factor adjoints and regenerates
 243         void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
 244             init( adj_b_V, adj_b_F, getZeroAdjV(*_fg), getZeroAdjF(*_fg) );
 245         }
 246
 247         /// Initializes variable belief adjoints (and sets factor belief adjoints to zero) and with zero initial factor adjoints and regenerates
 248         void init( const std::vector<Prob> &adj_b_V ) {
 249             init(adj_b_V, getZeroAdjF(*_fg));
 250         }
 251
 252         /// Run until change is less than given tolerance
 253         void run();
 254
 255         /// Return number of iterations done so far
 256         size_t doneIters() { return _iters; }
 257
 258         /// Returns reference to variable factor adjoint
 259         Prob& adj_psi_V(size_t i) { return _adj_psi_V[i]; }
 260         /// Returns constant reference to variable factor adjoint
 261         const Prob& adj_psi_V(size_t i) const { return _adj_psi_V[i]; }
 262         /// Returns reference to factor adjoint
 263         Prob& adj_psi_F(size_t I) { return _adj_psi_F[I]; }
 264         /// Returns constant reference to factor adjoint
 265         const Prob& adj_psi_F(size_t I) const { return _adj_psi_F[I]; }
 266         /// Returns reference to variable belief adjoint
 267         Prob& adj_b_V(size_t i) { return _adj_b_V[i]; }
 268         /// Returns constant reference to variable belief adjoint
 269         const Prob& adj_b_V(size_t i) const { return _adj_b_V[i]; }
 270         /// Returns reference to factor belief adjoint
 271         Prob& adj_b_F(size_t I) { return _adj_b_F[I]; }
 272         /// Returns constant reference to factor belief adjoint
 273         const Prob& adj_b_F(size_t I) const { return _adj_b_F[I]; }
 274
 275      protected:
 276         /// Returns reference to variable->factor message adjoint
 277         Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
 278         /// Returns constant reference to variable->factor message adjoint
 279         const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
 280         /// Returns reference to factor->variable message adjoint
 281         Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
 282         /// Returns constant reference to factor->variable message adjoint
 283         const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
 284
 285      public:
 286         /// Parameters of this algorithm
 287         /* PROPERTIES(props,BBP) {
 288            /// Enumeration of possible update schedules
 289            DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
 290
 291            /// Verbosity
 292            size_t verbose;
 293
 294            /// Maximum number of iterations
 295            size_t maxiter;
 296
 297            /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
 298            Real tol;
 299
 300            /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
 301            Real damping;
 302
 303            /// Update schedule
 304            UpdateType updates;
 305
 306            // DISABLED BECAUSE IT IS BUGGY:
 307            // bool clean_updates;
 308         }
 309         */
 310 /* {{{ GENERATED CODE: DO NOT EDIT. Created by
 311     ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp
 312 */
 313         struct Properties {
 314             /// Enumeration of possible update schedules
 315             DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
 316             /// Verbosity
 317             size_t verbose;
 318             /// Maximum number of iterations
 319             size_t maxiter;
 320             /// Tolerance (not used for updates = SEQ_BP_REV, SEQ_BP_FWD)
 321             Real tol;
 322             /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
 323             Real damping;
 324             /// Update schedule
 325             UpdateType updates;
 326
 327             /// Set members from PropertySet
 328             void set(const PropertySet &opts);
 329             /// Get members into PropertySet
 330             PropertySet get() const;
 331             /// Convert to a string which can be parsed as a PropertySet
 332             std::string toString() const;
 333         } props;
 334 /* }}} END OF GENERATED CODE */
 335 };
 336
 337
 338 /// Enumeration of several cost functions that can be used with BBP.
 339 DAI_ENUM(bbp_cfn_t,CFN_GIBBS_B,CFN_GIBBS_B2,CFN_GIBBS_EXP,CFN_GIBBS_B_FACTOR,CFN_GIBBS_B2_FACTOR,CFN_GIBBS_EXP_FACTOR,CFN_VAR_ENT,CFN_FACTOR_ENT,CFN_BETHE_ENT);
 340
 341 /// Initialise BBP using InfAlg, cost function, and stateP
 342 /** Calls bbp.init with adjoints calculated from ia.beliefV and
 343  *  ia.beliefF. stateP is a Gibbs state and can be NULL, it will be
 344  *  initialised using a Gibbs run of 2*fg.Iterations() iterations.
 345  */
 346 void initBBPCostFnAdj( BBP &bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
 347
 348 /// Answers question: does the given cost function depend on having a Gibbs state?
 349 bool needGibbsState( bbp_cfn_t cfn );
 350
 351 /// Calculate actual value of cost function (cfn_type, stateP)
 352 /** This function returns the actual value of the cost function whose
 353  *  gradient with respect to singleton beliefs is given by
 354  *  gibbsToB1Adj on the same arguments
 355  */
 356 Real getCostFn( const InfAlg &fg, bbp_cfn_t cfn_type, const std::vector<size_t> *stateP );
 357
 358 /// Function to test the validity of adjoints computed by BBP given a state for each variable using numerical derivatives.
 359 /** Factors containing a variable are multiplied by psi_1 adjustments to verify accuracy of _adj_psi_V.
 360  *  \param bp BP object.
 361  *  \param state Global state of all variables.
 362  *  \param bbp_props BBP Properties.
 363  *  \param cfn Cost function to be used.
 364  *  \param h controls size of perturbation.
 365  */
 366 Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, bbp_cfn_t cfn, Real h );
 367
 368
 369 } // end of namespace dai
 370
 371
 372 #endif