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