Fixed bug (found by cax): when building MatLab MEX files, GMP libraries were not...
[libdai.git] / include / dai / bbp.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 BBP, which implements Back-Belief-Propagation
11 /// \todo Clean up code
12
13
14 #ifndef ___defined_libdai_bbp_h
15 #define ___defined_libdai_bbp_h
16
17
18 #include <vector>
19 #include <utility>
20
21 #include <dai/prob.h>
22 #include <dai/daialg.h>
23 #include <dai/factorgraph.h>
24 #include <dai/enum.h>
25 #include <dai/bp_dual.h>
26
27
28 namespace dai {
29
30
31 /// Enumeration of several cost functions that can be used with BBP
32 /** \note This class is meant as a base class for BBPCostFunction, which provides additional functionality.
33 */
34 DAI_ENUM(BBPCostFunctionBase,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);
35
36
37 /// Predefined cost functions that can be used with BBP
38 class BBPCostFunction : public BBPCostFunctionBase {
39 public:
40 /// Default constructor
41 BBPCostFunction() : BBPCostFunctionBase() {}
42
43 /// Construct from BBPCostFunctionBase \a x
44 BBPCostFunction( const BBPCostFunctionBase &x ) : BBPCostFunctionBase(x) {}
45
46 /// Returns whether this cost function depends on having a Gibbs state
47 bool needGibbsState() const;
48
49 /// Evaluates cost function in state \a stateP using the information in inference algorithm \a ia
50 Real evaluate( const InfAlg &ia, const std::vector<size_t> *stateP ) const;
51
52 /// Assignment operator
53 BBPCostFunction& operator=( const BBPCostFunctionBase &x ) {
54 BBPCostFunctionBase::operator=( x );
55 return *this;
56 }
57 };
58
59
60 /// Implements BBP (Back-Belief-Propagation) [\ref EaG09]
61 /** \author Frederik Eaton
62 */
63 class BBP {
64 private:
65 /// \name Input variables
66 //@{
67 /// Stores a BP_dual helper object
68 BP_dual _bp_dual;
69 /// Pointer to the factor graph
70 const FactorGraph *_fg;
71 /// Pointer to the approximate inference algorithm (currently, only BP objects are supported)
72 const InfAlg *_ia;
73 //@}
74
75 /// \name Output variables
76 //@{
77 /// Variable factor adjoints
78 std::vector<Prob> _adj_psi_V;
79 /// Factor adjoints
80 std::vector<Prob> _adj_psi_F;
81 /// Variable->factor message adjoints (indexed [i][_I])
82 std::vector<std::vector<Prob> > _adj_n;
83 /// Factor->variable message adjoints (indexed [i][_I])
84 std::vector<std::vector<Prob> > _adj_m;
85 /// Normalized variable belief adjoints
86 std::vector<Prob> _adj_b_V;
87 /// Normalized factor belief adjoints
88 std::vector<Prob> _adj_b_F;
89 //@}
90
91 /// \name Internal state variables
92 //@{
93 /// Initial variable factor adjoints
94 std::vector<Prob> _init_adj_psi_V;
95 /// Initial factor adjoints
96 std::vector<Prob> _init_adj_psi_F;
97
98 /// Unnormalized variable->factor message adjoint (indexed [i][_I])
99 std::vector<std::vector<Prob> > _adj_n_unnorm;
100 /// Unnormalized factor->variable message adjoint (indexed [i][_I])
101 std::vector<std::vector<Prob> > _adj_m_unnorm;
102 /// Updated normalized variable->factor message adjoint (indexed [i][_I])
103 std::vector<std::vector<Prob> > _new_adj_n;
104 /// Updated normalized factor->variable message adjoint (indexed [i][_I])
105 std::vector<std::vector<Prob> > _new_adj_m;
106 /// Unnormalized variable belief adjoints
107 std::vector<Prob> _adj_b_V_unnorm;
108 /// Unnormalized factor belief adjoints
109 std::vector<Prob> _adj_b_F_unnorm;
110
111 /// _Tmsg[i][_I] (see eqn. (41) in [\ref EaG09])
112 std::vector<std::vector<Prob > > _Tmsg;
113 /// _Umsg[I][_i] (see eqn. (42) in [\ref EaG09])
114 std::vector<std::vector<Prob > > _Umsg;
115 /// _Smsg[i][_I][_j] (see eqn. (43) in [\ref EaG09])
116 std::vector<std::vector<std::vector<Prob > > > _Smsg;
117 /// _Rmsg[I][_i][_J] (see eqn. (44) in [\ref EaG09])
118 std::vector<std::vector<std::vector<Prob > > > _Rmsg;
119
120 /// Number of iterations done
121 size_t _iters;
122 //@}
123
124 /// \name Index cache management (for performance)
125 //@{
126 /// Index type
127 typedef std::vector<size_t> _ind_t;
128 /// Cached indices (indexed [i][_I])
129 std::vector<std::vector<_ind_t> > _indices;
130 /// Prepares index cache _indices
131 /** \see bp.cpp
132 */
133 void RegenerateInds();
134 /// Returns an index from the cache
135 const _ind_t& _index(size_t i, size_t _I) const { return _indices[i][_I]; }
136 //@}
137
138 /// \name Initialization helper functions
139 //@{
140 /// Calculate T values; see eqn. (41) in [\ref EaG09]
141 void RegenerateT();
142 /// Calculate U values; see eqn. (42) in [\ref EaG09]
143 void RegenerateU();
144 /// Calculate S values; see eqn. (43) in [\ref EaG09]
145 void RegenerateS();
146 /// Calculate R values; see eqn. (44) in [\ref EaG09]
147 void RegenerateR();
148 /// Calculate _adj_b_V_unnorm and _adj_b_F_unnorm from _adj_b_V and _adj_b_F
149 void RegenerateInputs();
150 /// Initialise members for factor adjoints
151 /** \pre RegenerateInputs() should be called first
152 */
153 void RegeneratePsiAdjoints();
154 /// Initialise members for message adjoints for parallel algorithm
155 /** \pre RegenerateInputs() should be called first
156 */
157 void RegenerateParMessageAdjoints();
158 /// Initialise members for message adjoints for sequential algorithm
159 /** Same as RegenerateMessageAdjoints, but calls sendSeqMsgN rather
160 * than updating _adj_n (and friends) which are unused in the sequential algorithm.
161 * \pre RegenerateInputs() should be called first
162 */
163 void RegenerateSeqMessageAdjoints();
164 /// Called by \a init, recalculates intermediate values
165 void Regenerate();
166 //@}
167
168 /// \name Accessors/mutators
169 //@{
170 /// Returns reference to T value; see eqn. (41) in [\ref EaG09]
171 Prob & T(size_t i, size_t _I) { return _Tmsg[i][_I]; }
172 /// Returns constant reference to T value; see eqn. (41) in [\ref EaG09]
173 const Prob & T(size_t i, size_t _I) const { return _Tmsg[i][_I]; }
174 /// Returns reference to U value; see eqn. (42) in [\ref EaG09]
175 Prob & U(size_t I, size_t _i) { return _Umsg[I][_i]; }
176 /// Returns constant reference to U value; see eqn. (42) in [\ref EaG09]
177 const Prob & U(size_t I, size_t _i) const { return _Umsg[I][_i]; }
178 /// Returns reference to S value; see eqn. (43) in [\ref EaG09]
179 Prob & S(size_t i, size_t _I, size_t _j) { return _Smsg[i][_I][_j]; }
180 /// Returns constant reference to S value; see eqn. (43) in [\ref EaG09]
181 const Prob & S(size_t i, size_t _I, size_t _j) const { return _Smsg[i][_I][_j]; }
182 /// Returns reference to R value; see eqn. (44) in [\ref EaG09]
183 Prob & R(size_t I, size_t _i, size_t _J) { return _Rmsg[I][_i][_J]; }
184 /// Returns constant reference to R value; see eqn. (44) in [\ref EaG09]
185 const Prob & R(size_t I, size_t _i, size_t _J) const { return _Rmsg[I][_i][_J]; }
186
187 /// Returns reference to variable->factor message adjoint
188 Prob& adj_n(size_t i, size_t _I) { return _adj_n[i][_I]; }
189 /// Returns constant reference to variable->factor message adjoint
190 const Prob& adj_n(size_t i, size_t _I) const { return _adj_n[i][_I]; }
191 /// Returns reference to factor->variable message adjoint
192 Prob& adj_m(size_t i, size_t _I) { return _adj_m[i][_I]; }
193 /// Returns constant reference to factor->variable message adjoint
194 const Prob& adj_m(size_t i, size_t _I) const { return _adj_m[i][_I]; }
195 //@}
196
197 /// \name Parallel algorithm
198 //@{
199 /// Calculates new variable->factor message adjoint
200 /** Increases variable factor adjoint according to eqn. (27) in [\ref EaG09] and
201 * calculates the new variable->factor message adjoint according to eqn. (29) in [\ref EaG09].
202 */
203 void calcNewN( size_t i, size_t _I );
204 /// Calculates new factor->variable message adjoint
205 /** Increases factor adjoint according to eqn. (28) in [\ref EaG09] and
206 * calculates the new factor->variable message adjoint according to the r.h.s. of eqn. (30) in [\ref EaG09].
207 */
208 void calcNewM( size_t i, size_t _I );
209 /// Calculates unnormalized variable->factor message adjoint from the normalized one
210 void calcUnnormMsgN( size_t i, size_t _I );
211 /// Calculates unnormalized factor->variable message adjoint from the normalized one
212 void calcUnnormMsgM( size_t i, size_t _I );
213 /// Updates (un)normalized variable->factor message adjoints
214 void upMsgN( size_t i, size_t _I );
215 /// Updates (un)normalized factor->variable message adjoints
216 void upMsgM( size_t i, size_t _I );
217 /// Do one parallel update of all message adjoints
218 void doParUpdate();
219 //@}
220
221 /// \name Sequential algorithm
222 //@{
223 /// Helper function for sendSeqMsgM(): increases factor->variable message adjoint by \a p and calculates the corresponding unnormalized adjoint
224 void incrSeqMsgM( size_t i, size_t _I, const Prob& p );
225 // DISABLED BECAUSE IT IS BUGGY:
226 // void updateSeqMsgM( size_t i, size_t _I );
227 /// Sets normalized factor->variable message adjoint and calculates the corresponding unnormalized adjoint
228 void setSeqMsgM( size_t i, size_t _I, const Prob &p );
229 /// Implements routine Send-n in Figure 5 in [\ref EaG09]
230 void sendSeqMsgN( size_t i, size_t _I, const Prob &f );
231 /// Implements routine Send-m in Figure 5 in [\ref EaG09]
232 void sendSeqMsgM( size_t i, size_t _I );
233 //@}
234
235 /// Computes the adjoint of the unnormed probability vector from the normalizer and the adjoint of the normalized probability vector
236 /** \see eqn. (13) in [\ref EaG09]
237 */
238 Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w );
239
240 /// Calculates averaged L1 norm of unnormalized message adjoints
241 Real getUnMsgMag();
242 /// Calculates averaged L1 norms of current and new normalized message adjoints
243 void getMsgMags( Real &s, Real &new_s );
244 /// Returns indices and magnitude of the largest normalized factor->variable message adjoint
245 void getArgmaxMsgM( size_t &i, size_t &_I, Real &mag );
246 /// Returns magnitude of the largest (in L1-norm) normalized factor->variable message adjoint
247 Real getMaxMsgM();
248
249 /// Calculates sum of L1 norms of all normalized factor->variable message adjoints
250 Real getTotalMsgM();
251 /// Calculates sum of L1 norms of all updated normalized factor->variable message adjoints
252 Real getTotalNewMsgM();
253 /// Calculates sum of L1 norms of all normalized variable->factor message adjoints
254 Real getTotalMsgN();
255
256 /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the factors in the factor graph \a fg
257 std::vector<Prob> getZeroAdjF( const FactorGraph &fg );
258 /// Returns a vector of Probs (filled with zeroes) with state spaces corresponding to the variables in the factor graph \a fg
259 std::vector<Prob> getZeroAdjV( const FactorGraph &fg );
260
261 public:
262 /// \name Constructors/destructors
263 //@{
264 /// Construct BBP object from a InfAlg \a ia and a PropertySet \a opts
265 /** \param ia should be a BP object or something compatible
266 * \param opts Parameters @see Properties
267 */
268 BBP( const InfAlg *ia, const PropertySet &opts ) : _bp_dual(ia), _fg(&(ia->fg())), _ia(ia) {
269 props.set(opts);
270 }
271 //@}
272
273 /// \name Initialization
274 //@{
275 /// Initializes from given belief adjoints \a adj_b_V, \a adj_b_F and initial factor adjoints \a adj_psi_V, \a adj_psi_F
276 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 ) {
277 _adj_b_V = adj_b_V;
278 _adj_b_F = adj_b_F;
279 _init_adj_psi_V = adj_psi_V;
280 _init_adj_psi_F = adj_psi_F;
281 Regenerate();
282 }
283
284 /// Initializes from given belief adjoints \a adj_b_V and \a adj_b_F (setting initial factor adjoints to zero)
285 void init( const std::vector<Prob> &adj_b_V, const std::vector<Prob> &adj_b_F ) {
286 init( adj_b_V, adj_b_F, getZeroAdjV(*_fg), getZeroAdjF(*_fg) );
287 }
288
289 /// Initializes variable belief adjoints \a adj_b_V (and sets factor belief adjoints and initial factor adjoints to zero)
290 void init_V( const std::vector<Prob> &adj_b_V ) {
291 init( adj_b_V, getZeroAdjF(*_fg) );
292 }
293
294 /// Initializes factor belief adjoints \a adj_b_F (and sets variable belief adjoints and initial factor adjoints to zero)
295 void init_F( const std::vector<Prob> &adj_b_F ) {
296 init( getZeroAdjV(*_fg), adj_b_F );
297 }
298
299 /// Initializes with adjoints calculated from cost function \a cfn, and state \a stateP
300 /** Uses the internal pointer to an inference algorithm in combination with the cost function and state for initialization.
301 * \param cfn Cost function used for initialization;
302 * \param stateP is a Gibbs state and can be NULL; it will be initialised using a Gibbs run.
303 */
304 void initCostFnAdj( const BBPCostFunction &cfn, const std::vector<size_t> *stateP );
305 //@}
306
307 /// \name BBP Algorithm
308 //@{
309 /// Perform iterative updates until change is less than given tolerance
310 void run();
311 //@}
312
313 /// \name Query results
314 //@{
315 /// Returns reference to variable factor adjoint
316 Prob& adj_psi_V(size_t i) { return _adj_psi_V[i]; }
317 /// Returns constant reference to variable factor adjoint
318 const Prob& adj_psi_V(size_t i) const { return _adj_psi_V[i]; }
319 /// Returns reference to factor adjoint
320 Prob& adj_psi_F(size_t I) { return _adj_psi_F[I]; }
321 /// Returns constant reference to factor adjoint
322 const Prob& adj_psi_F(size_t I) const { return _adj_psi_F[I]; }
323 /// Returns reference to variable belief adjoint
324 Prob& adj_b_V(size_t i) { return _adj_b_V[i]; }
325 /// Returns constant reference to variable belief adjoint
326 const Prob& adj_b_V(size_t i) const { return _adj_b_V[i]; }
327 /// Returns reference to factor belief adjoint
328 Prob& adj_b_F(size_t I) { return _adj_b_F[I]; }
329 /// Returns constant reference to factor belief adjoint
330 const Prob& adj_b_F(size_t I) const { return _adj_b_F[I]; }
331 /// Return number of iterations done so far
332 size_t Iterations() { return _iters; }
333 //@}
334
335 public:
336 /// Parameters for BBP
337 /* PROPERTIES(props,BBP) {
338 /// Enumeration of possible update schedules
339 /// The following update schedules are defined:
340 /// - SEQ_FIX fixed sequential updates
341 /// - SEQ_MAX maximum residual updates (inspired by [\ref EMK06])
342 /// - SEQ_BP_REV schedule used by BP, but reversed
343 /// - SEQ_BP_FWD schedule used by BP
344 /// - PAR parallel updates
345 DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
346
347 /// Verbosity (amount of output sent to stderr)
348 size_t verbose = 0;
349
350 /// Maximum number of iterations
351 size_t maxiter;
352
353 /// Tolerance for convergence test
354 /// \note Not used for updates = SEQ_BP_REV, SEQ_BP_FWD
355 Real tol;
356
357 /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
358 Real damping;
359
360 /// Update schedule
361 UpdateType updates;
362
363 // DISABLED BECAUSE IT IS BUGGY:
364 // bool clean_updates;
365 }
366 */
367 /* {{{ GENERATED CODE: DO NOT EDIT. Created by
368 ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp
369 */
370 struct Properties {
371 /// Enumeration of possible update schedules
372 /** The following update schedules are defined:
373 * - SEQ_FIX fixed sequential updates
374 * - SEQ_MAX maximum residual updates (inspired by [\ref EMK06])
375 * - SEQ_BP_REV schedule used by BP, but reversed
376 * - SEQ_BP_FWD schedule used by BP
377 * - PAR parallel updates
378 */
379 DAI_ENUM(UpdateType,SEQ_FIX,SEQ_MAX,SEQ_BP_REV,SEQ_BP_FWD,PAR);
380 /// Verbosity (amount of output sent to stderr)
381 size_t verbose;
382 /// Maximum number of iterations
383 size_t maxiter;
384 /// Tolerance for convergence test
385 /** \note Not used for updates = SEQ_BP_REV, SEQ_BP_FWD
386 */
387 Real tol;
388 /// Damping constant (0 for none); damping = 1 - lambda where lambda is the damping constant used in [\ref EaG09]
389 Real damping;
390 /// Update schedule
391 UpdateType updates;
392
393 /// Set members from PropertySet
394 /** \throw UNKNOWN_PROPERTY if a Property key is not recognized
395 * \throw NOT_ALL_PROPERTIES_SPECIFIED if an expected Property is missing
396 */
397 void set(const PropertySet &opts);
398 /// Get members into PropertySet
399 PropertySet get() const;
400 /// Convert to a string which can be parsed as a PropertySet
401 std::string toString() const;
402 } props;
403 /* }}} END OF GENERATED CODE */
404 };
405
406
407 /// Function to verify the validity of adjoints computed by BBP using numerical differentiation.
408 /** Factors containing a variable are multiplied by small adjustments to verify accuracy of calculated variable factor adjoints.
409 * \param bp BP object;
410 * \param state Global state of all variables;
411 * \param bbp_props BBP parameters;
412 * \param cfn Cost function to be used;
413 * \param h Size of perturbation.
414 * \relates BBP
415 */
416 Real numericBBPTest( const InfAlg &bp, const std::vector<size_t> *state, const PropertySet &bbp_props, const BBPCostFunction &cfn, Real h );
417
418
419 } // end of namespace dai
420
421
422 #endif