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