Improved documentation of include/dai/bp_dual.h
[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 */
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