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