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