Cleanup of BBP code
[libdai.git] / src / bbp.cpp
1 /* Copyright (C) 2009 Frederik Eaton [frederik at ofb dot net]
2
3 This file is part of libDAI.
4
5 libDAI is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License, or
8 (at your option) any later version.
9
10 libDAI is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with libDAI; if not, write to the Free Software
17 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18 */
19
20
21 #include <dai/bp.h>
22 #include <dai/bbp.h>
23 #include <dai/gibbs.h>
24 #include <dai/util.h>
25 #include <dai/bipgraph.h>
26
27
28 namespace dai {
29
30
31 using namespace std;
32
33
34 typedef BipartiteGraph::Neighbor Neighbor;
35
36
37 Prob unnormAdjoint( const Prob &w, Real Z_w, const Prob &adj_w ) {
38 assert( w.size() == adj_w.size() );
39 Prob adj_w_unnorm( w.size(), 0.0 );
40 Real s = 0.0;
41 for( size_t i = 0; i < w.size(); i++ )
42 s += w[i] * adj_w[i];
43 for( size_t i = 0; i < w.size(); i++ )
44 adj_w_unnorm[i] = (adj_w[i] - s) / Z_w;
45 return adj_w_unnorm;
46 // THIS WOULD BE ABOUT 50% SLOWER: return (adj_w - (w * adj_w).sum()) / Z_w;
47 }
48
49
50 size_t getFactorEntryForState( const FactorGraph &fg, size_t I, const vector<size_t> &state ) {
51 size_t f_entry = 0;
52 for( int _j = fg.nbF(I).size() - 1; _j >= 0; _j-- ) {
53 // note that iterating over nbF(I) yields the same ordering
54 // of variables as iterating over factor(I).vars()
55 size_t j = fg.nbF(I)[_j];
56 f_entry *= fg.var(j).states();
57 f_entry += state[j];
58 }
59 return f_entry;
60 }
61
62
63 void initBBPCostFnAdj( BBP& bbp, const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
64 const FactorGraph &fg = ia.fg();
65
66 switch( (size_t)cfn_type ) {
67 case bbp_cfn_t::cfn_bethe_ent: {
68 vector<Prob> b1_adj;
69 vector<Prob> b2_adj;
70 vector<Prob> psi1_adj;
71 vector<Prob> psi2_adj;
72 b1_adj.reserve(fg.nrVars());
73 psi1_adj.reserve(fg.nrVars());
74 b2_adj.reserve(fg.nrFactors());
75 psi2_adj.reserve(fg.nrFactors());
76 for( size_t i = 0; i < fg.nrVars(); i++ ) {
77 size_t dim = fg.var(i).states();
78 int c = fg.nbV(i).size();
79 Prob p(dim,0.0);
80 for( size_t xi = 0; xi < dim; xi++ )
81 p[xi] = (1-c)*(1+log(ia.beliefV(i)[xi]));
82 b1_adj.push_back(p);
83
84 for( size_t xi = 0; xi < dim; xi++ )
85 p[xi] = -ia.beliefV(i)[xi];
86 psi1_adj.push_back(p);
87 }
88 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
89 size_t dim = fg.factor(I).states();
90 Prob p(dim,0.0);
91 for( size_t xI = 0; xI < dim; xI++ )
92 p[xI] = 1 + log(ia.beliefF(I)[xI]/fg.factor(I).p()[xI]);
93 b2_adj.push_back(p);
94
95 for( size_t xI = 0; xI < dim; xI++ )
96 p[xI] = -ia.beliefF(I)[xI]/fg.factor(I).p()[xI];
97 psi2_adj.push_back(p);
98 }
99 bbp.init(b1_adj, b2_adj, psi1_adj, psi2_adj);
100 break;
101 }
102 case bbp_cfn_t::cfn_factor_ent: {
103 vector<Prob> b2_adj;
104 b2_adj.reserve(fg.nrFactors());
105 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
106 size_t dim = fg.factor(I).states();
107 Prob p(dim,0.0);
108 for( size_t xI = 0; xI < dim; xI++ ) {
109 double bIxI = ia.beliefF(I)[xI];
110 if( bIxI < 1.0e-15 )
111 p[xI] = -1.0e10;
112 else
113 p[xI] = 1+log(bIxI);
114 }
115 b2_adj.push_back(p);
116 }
117 bbp.init(get_zero_adj_V(fg), b2_adj);
118 break;
119 }
120 case bbp_cfn_t::cfn_var_ent: {
121 vector<Prob> b1_adj;
122 b1_adj.reserve(fg.nrVars());
123 for(size_t i=0; i<fg.nrVars(); i++) {
124 size_t dim = fg.var(i).states();
125 Prob p(dim,0.0);
126 for( size_t xi = 0; xi < fg.var(i).states(); xi++ ) {
127 double bixi = ia.beliefV(i)[xi];
128 if( bixi < 1.0e-15 )
129 p[xi] = -1.0e10;
130 else
131 p[xi] = 1+log(bixi);
132 }
133 b1_adj.push_back(p);
134 }
135 bbp.init(b1_adj);
136 break;
137 }
138 case bbp_cfn_t::cfn_gibbs_b:
139 case bbp_cfn_t::cfn_gibbs_b2:
140 case bbp_cfn_t::cfn_gibbs_exp: {
141 // cost functions that use Gibbs sample, summing over variable
142 // marginals
143 vector<size_t> state;
144 if(stateP==NULL)
145 state = getGibbsState(ia,2*ia.Iterations());
146 else
147 state = *stateP;
148 assert(state.size()==fg.nrVars());
149
150 vector<Prob> b1_adj;
151 b1_adj.reserve(fg.nrVars());
152 for( size_t i = 0; i < state.size(); i++ ) {
153 size_t n = fg.var(i).states();
154 Prob delta(n,0.0);
155 assert(/*0<=state[i] &&*/ state[i] < n);
156 double b = ia.beliefV(i)[state[i]];
157 switch( (size_t)cfn_type ) {
158 case bbp_cfn_t::cfn_gibbs_b:
159 delta[state[i]] = 1.0; break;
160 case bbp_cfn_t::cfn_gibbs_b2:
161 delta[state[i]] = b; break;
162 case bbp_cfn_t::cfn_gibbs_exp:
163 delta[state[i]] = exp(b); break;
164 default: DAI_THROW(INTERNAL_ERROR);
165 }
166 b1_adj.push_back(delta);
167 }
168 bbp.init(b1_adj);
169 break;
170 }
171 case bbp_cfn_t::cfn_gibbs_b_factor:
172 case bbp_cfn_t::cfn_gibbs_b2_factor:
173 case bbp_cfn_t::cfn_gibbs_exp_factor: {
174 // cost functions that use Gibbs sample, summing over factor
175 // marginals
176 vector<size_t> state;
177 if(stateP==NULL)
178 state = getGibbsState(ia,2*ia.Iterations());
179 else
180 state = *stateP;
181 assert(state.size()==fg.nrVars());
182
183 vector<Prob> b2_adj;
184 b2_adj.reserve(fg.nrVars());
185 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
186 size_t n = fg.factor(I).states();
187 Prob delta(n,0.0);
188
189 size_t x_I = getFactorEntryForState(fg, I, state);
190 assert(/*0<=x_I &&*/ x_I < n);
191
192 double b = ia.beliefF(I)[x_I];
193 switch( (size_t)cfn_type ) {
194 case bbp_cfn_t::cfn_gibbs_b_factor:
195 delta[x_I] = 1.0; break;
196 case bbp_cfn_t::cfn_gibbs_b2_factor:
197 delta[x_I] = b; break;
198 case bbp_cfn_t::cfn_gibbs_exp_factor:
199 delta[x_I] = exp(b); break;
200 default: DAI_THROW(INTERNAL_ERROR);
201 }
202 b2_adj.push_back(delta);
203 }
204 bbp.init(get_zero_adj_V(fg), b2_adj);
205 break;
206 }
207 default: DAI_THROW(UNKNOWN_ENUM_VALUE);
208 }
209 }
210
211
212 Real getCostFn( const InfAlg &ia, bbp_cfn_t cfn_type, const vector<size_t> *stateP ) {
213 double cf = 0.0;
214 const FactorGraph &fg = ia.fg();
215
216 switch( (size_t)cfn_type ) {
217 case bbp_cfn_t::cfn_bethe_ent: // ignores state
218 cf = -ia.logZ();
219 break;
220 case bbp_cfn_t::cfn_var_ent: // ignores state
221 for(size_t i=0; i<fg.nrVars(); i++)
222 cf += -ia.beliefV(i).entropy();
223 break;
224 case bbp_cfn_t::cfn_factor_ent: // ignores state
225 for(size_t I=0; I<fg.nrFactors(); I++)
226 cf += -ia.beliefF(I).entropy();
227 break;
228 case bbp_cfn_t::cfn_gibbs_b:
229 case bbp_cfn_t::cfn_gibbs_b2:
230 case bbp_cfn_t::cfn_gibbs_exp: {
231 assert(stateP != NULL);
232 vector<size_t> state = *stateP;
233 assert(state.size()==fg.nrVars());
234 for( size_t i = 0; i < fg.nrVars(); i++ ) {
235 double b = ia.beliefV(i)[state[i]];
236 switch( (size_t)cfn_type ) {
237 case bbp_cfn_t::cfn_gibbs_b: cf += b; break;
238 case bbp_cfn_t::cfn_gibbs_b2: cf += b*b/2; break;
239 case bbp_cfn_t::cfn_gibbs_exp: cf += exp(b); break;
240 default: DAI_THROW(INTERNAL_ERROR);
241 }
242 }
243 break;
244 }
245 case bbp_cfn_t::cfn_gibbs_b_factor:
246 case bbp_cfn_t::cfn_gibbs_b2_factor:
247 case bbp_cfn_t::cfn_gibbs_exp_factor: {
248 assert(stateP != NULL);
249 vector<size_t> state = *stateP;
250 assert(state.size()==fg.nrVars());
251 for( size_t I = 0; I < fg.nrFactors(); I++ ) {
252 size_t x_I = getFactorEntryForState(fg, I, state);
253 double b = ia.beliefF(I)[x_I];
254 switch((size_t)cfn_type) {
255 case bbp_cfn_t::cfn_gibbs_b_factor: cf += b; break;
256 case bbp_cfn_t::cfn_gibbs_b2_factor: cf += b*b/2; break;
257 case bbp_cfn_t::cfn_gibbs_exp_factor: cf += exp(b); break;
258 default: DAI_THROW(INTERNAL_ERROR);
259 }
260 }
261 break;
262 }
263 default:
264 cerr << "Unknown cost function " << cfn_type << endl;
265 DAI_THROW(UNKNOWN_ENUM_VALUE);
266 }
267 return cf;
268 }
269
270
271 vector<Prob> get_zero_adj_F( const FactorGraph &fg ) {
272 vector<Prob> adj_2;
273 adj_2.reserve( fg.nrFactors() );
274 for( size_t I = 0; I < fg.nrFactors(); I++ )
275 adj_2.push_back( Prob( fg.factor(I).states(), Real(0.0) ) );
276 return adj_2;
277 }
278
279
280 vector<Prob> get_zero_adj_V( const FactorGraph &fg ) {
281 vector<Prob> adj_1;
282 adj_1.reserve( fg.nrVars() );
283 for( size_t i = 0; i < fg.nrVars(); i++ )
284 adj_1.push_back( Prob( fg.var(i).states(), Real(0.0) ) );
285 return adj_1;
286 }
287
288
289 #define LOOP_ij(body) { \
290 size_t i_states = _fg->var(i).states(); \
291 size_t j_states = _fg->var(j).states(); \
292 if(_fg->var(i) > _fg->var(j)) { \
293 size_t xij=0; \
294 for(size_t xi=0; xi<i_states; xi++) { \
295 for(size_t xj=0; xj<j_states; xj++) { \
296 body; \
297 xij++; \
298 } \
299 } \
300 } else { \
301 size_t xij=0; \
302 for(size_t xj=0; xj<j_states; xj++) { \
303 for(size_t xi=0; xi<i_states; xi++) { \
304 body; \
305 xij++; \
306 } \
307 } \
308 } \
309 }
310
311
312 void BBP::RegenerateInds() {
313 // initialise _indices
314 // typedef std::vector<size_t> _ind_t;
315 // std::vector<std::vector<_ind_t> > _indices;
316 _indices.resize(_fg->nrVars());
317 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
318 _indices[i].reserve(_fg->nbV(i).size());
319 foreach(const Neighbor &I, _fg->nbV(i)) {
320 _ind_t index;
321 index.reserve(_fg->factor(I).states());
322 for(IndexFor k( _fg->var(i), _fg->factor(I).vars() ); k >= 0; ++k) {
323 index.push_back(k);
324 }
325 _indices[i].push_back(index);
326 }
327 }
328 }
329
330
331 void BBP::RegenerateT() {
332 // _T[i][_I]
333 _T.clear();
334 _T.resize( _fg->nrVars() );
335 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
336 _T[i].resize( _fg->nbV(i).size() );
337 foreach( const Neighbor &I, _fg->nbV(i) ) {
338 Prob prod( _fg->var(i).states(), 1.0 );
339 foreach( const Neighbor &J, _fg->nbV(i) )
340 if( J.node != I.node )
341 prod *= _bp_dual.msgM( i, J.iter );
342 _T[i][I.iter] = prod;
343 }
344 }
345 }
346
347
348 void BBP::RegenerateU() {
349 // _U[I][_i]
350 _U.clear();
351 _U.resize( _fg->nrFactors() );
352 for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
353 _U[I].resize( _fg->nbF(I).size() );
354 foreach( const Neighbor &i, _fg->nbF(I) ) {
355 Prob prod( _fg->factor(I).states(), 1.0 );
356 foreach( const Neighbor &j, _fg->nbF(I) )
357 if( i.node != j.node ) {
358 Prob n_jI( _bp_dual.msgN( j, j.dual ) );
359 const _ind_t &ind = _index( j, j.dual );
360 // multiply prod by n_jI
361 for( size_t x_I = 0; x_I < prod.size(); x_I++ )
362 prod[x_I] *= n_jI[ind[x_I]];
363 }
364 _U[I][i.iter] = prod;
365 }
366 }
367 }
368
369
370 void BBP::RegenerateS() {
371 // _S[i][_I][_j]
372 _S.clear();
373 _S.resize( _fg->nrVars() );
374 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
375 _S[i].resize( _fg->nbV(i).size() );
376 foreach( const Neighbor &I, _fg->nbV(i) ) {
377 _S[i][I.iter].resize( _fg->nbF(I).size() );
378 foreach( const Neighbor &j, _fg->nbF(I) )
379 if( i != j ) {
380 Factor prod( _fg->factor(I) );
381 foreach( const Neighbor &k, _fg->nbF(I) ) {
382 if( k != i && k.node != j.node ) {
383 const _ind_t &ind = _index( k, k.dual );
384 Prob p( _bp_dual.msgN( k, k.dual ) );
385 for( size_t x_I = 0; x_I < prod.states(); x_I++ )
386 prod.p()[x_I] *= p[ind[x_I]];
387 }
388 }
389 // "Marginalize" onto i|j (unnormalized)
390 // XXX improve efficiency?
391 Prob marg;
392 marg = prod.marginal( VarSet(_fg->var(i), _fg->var(j)), false ).p();
393 _S[i][I.iter][j.iter] = marg;
394 }
395 }
396 }
397 }
398
399
400 void BBP::RegenerateR() {
401 // _R[I][_i][_J]
402 _R.clear();
403 _R.resize( _fg->nrFactors() );
404 for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
405 _R[I].resize( _fg->nbF(I).size() );
406 foreach( const Neighbor &i, _fg->nbF(I) ) {
407 _R[I][i.iter].resize( _fg->nbV(i).size() );
408 foreach( const Neighbor &J, _fg->nbV(i) ) {
409 if( I != J ) {
410 Prob prod( _fg->var(i).states(), 1.0 );
411 foreach( const Neighbor &K, _fg->nbV(i) )
412 if( K.node != I && K.node != J.node )
413 prod *= _bp_dual.msgM(i,K.iter);
414 _R[I][i.iter][J.iter] = prod;
415 }
416 }
417 }
418 }
419 }
420
421
422 void BBP::RegenerateInputs() {
423 _adj_b_V_unnorm.clear();
424 _adj_b_V_unnorm.reserve(_fg->nrVars());
425 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
426 _adj_b_V_unnorm.push_back(
427 unnormAdjoint(_bp_dual.beliefV(i).p(),
428 _bp_dual.beliefVZ(i), _adj_b_V[i]));
429 }
430 _adj_b_F_unnorm.clear();
431 _adj_b_F_unnorm.reserve(_fg->nrFactors());
432 for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
433 _adj_b_F_unnorm.push_back(
434 unnormAdjoint(_bp_dual.beliefF(I).p(),
435 _bp_dual.beliefFZ(I), _adj_b_F[I]));
436 }
437 }
438
439
440 void BBP::RegeneratePsiAdjoints() {
441 _adj_psi_V.clear();
442 _adj_psi_V.reserve(_fg->nrVars());
443 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
444 Prob p(_adj_b_V_unnorm[i]);
445 assert(p.size()==_fg->var(i).states());
446 foreach(const Neighbor& I, _fg->nbV(i)) {
447 p *= _bp_dual.msgM(i,I.iter);
448 }
449 p += _init_adj_psi_V[i];
450 _adj_psi_V.push_back(p);
451 }
452 _adj_psi_F.clear();
453 _adj_psi_F.reserve(_fg->nrFactors());
454 for( size_t I = 0; I < _fg->nrFactors(); I++ ) {
455 Prob p(_adj_b_F_unnorm[I]);
456 assert(p.size()==_fg->factor(I).states());
457 foreach(const Neighbor& i, _fg->nbF(I)) {
458 Prob n_iI(_bp_dual.msgN(i,i.dual));
459 const _ind_t& ind = _index(i,i.dual);
460 // multiply prod with n_jI
461 for(size_t x_I = 0; x_I < p.size(); x_I++)
462 p[x_I] *= n_iI[ind[x_I]];
463 }
464 p += _init_adj_psi_F[I];
465 _adj_psi_F.push_back(p);
466 }
467 }
468
469
470 void BBP::RegenerateParMessageAdjoints() {
471 size_t nv = _fg->nrVars();
472 _adj_n.resize(nv);
473 _adj_m.resize(nv);
474 _adj_n_unnorm.resize(nv);
475 _adj_m_unnorm.resize(nv);
476 _new_adj_n.clear();
477 _new_adj_n.resize(nv);
478 _new_adj_m.clear();
479 _new_adj_m.resize(nv);
480 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
481 size_t n_i = _fg->nbV(i).size();
482 _adj_n[i].resize(n_i);
483 _adj_m[i].resize(n_i);
484 _adj_n_unnorm[i].resize(n_i);
485 _adj_m_unnorm[i].resize(n_i);
486 _new_adj_n[i].resize(n_i);
487 _new_adj_m[i].resize(n_i);
488 foreach(const Neighbor& I, _fg->nbV(i)) {
489 // calculate adj_n
490 {
491 Prob prod(_fg->factor(I).p());
492 prod *= _adj_b_F_unnorm[I];
493 foreach(const Neighbor& j, _fg->nbF(I)) {
494 if(i != j) {
495 Prob n_jI(_bp_dual.msgN(j,j.dual));
496 const _ind_t& ind = _index(j,j.dual);
497 // multiply prod with n_jI
498 for( size_t x_I = 0; x_I < prod.size(); x_I++ )
499 prod[x_I] *= n_jI[ind[x_I]];
500 }
501 }
502 Prob marg(_fg->var(i).states(), 0.0);
503 const _ind_t &ind = _index(i,I.iter);
504 for( size_t r = 0; r < prod.size(); r++ )
505 marg[ind[r]] += prod[r];
506 _new_adj_n[i][I.iter] = marg;
507 upMsgN(i,I.iter);
508 }
509
510 // calculate adj_m
511 {
512 Prob prod(_adj_b_V_unnorm[i]);
513 assert(prod.size()==_fg->var(i).states());
514 foreach(const Neighbor& J, _fg->nbV(i)) {
515 if(size_t(J) != size_t(I)) {
516 prod *= _bp_dual.msgM(i,J.iter);
517 }
518 }
519 _new_adj_m[i][I.iter] = prod;
520 upMsgM(i,I.iter);
521 }
522 }
523 }
524 }
525
526
527 void BBP::incrSeqMsgM( size_t i, size_t _I, const Prob &p ) {
528 if( props.clean_updates )
529 _new_adj_m[i][_I] += p;
530 else {
531 _adj_m[i][_I] += p;
532 calcUnnormMsgM(i, _I);
533 }
534 }
535
536
537 #if 0
538 Real pv_thresh=1000;
539 #endif
540
541
542 void BBP::updateSeqMsgM( size_t i, size_t _I ) {
543 if( props.clean_updates ) {
544 #if 0
545 if(_new_adj_m[i][_I].sumAbs() > pv_thresh ||
546 _adj_m[i][_I].sumAbs() > pv_thresh) {
547
548 DAI_DMSG("in updateSeqMsgM");
549 DAI_PV(i);
550 DAI_PV(_I);
551 DAI_PV(_adj_m[i][_I]);
552 DAI_PV(_new_adj_m[i][_I]);
553 }
554 #endif
555 _adj_m[i][_I] += _new_adj_m[i][_I];
556 calcUnnormMsgM( i, _I );
557 _new_adj_m[i][_I].fill( 0.0 );
558 }
559 }
560
561
562 void BBP::setSeqMsgM( size_t i, size_t _I, const Prob &p ) {
563 _adj_m[i][_I] = p;
564 calcUnnormMsgM(i, _I);
565 }
566
567
568 void BBP::sendSeqMsgN( size_t i, size_t _I, const Prob &f ) {
569 Prob f_unnorm = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), f );
570 const Neighbor &I = _fg->nbV(i)[_I];
571 assert( I.iter == _I );
572 _adj_psi_V[i] += f_unnorm * T(i,_I);
573 #if 0
574 if(f_unnorm.sumAbs() > pv_thresh) {
575 DAI_DMSG("in sendSeqMsgN");
576 DAI_PV(i);
577 DAI_PV(I);
578 DAI_PV(_I);
579 DAI_PV(_bp_dual.msgN(i,_I));
580 DAI_PV(_bp_dual.zN(i,_I));
581 DAI_PV(_bp_dual.msgM(i,_I));
582 DAI_PV(_bp_dual.zM(i,_I));
583 DAI_PV(_fg->factor(I).p());
584 }
585 #endif
586 foreach( const Neighbor &J, _fg->nbV(i) ) {
587 if( J.node != I.node ) {
588 #if 0
589 if(f_unnorm.sumAbs() > pv_thresh) {
590 DAI_DMSG("in sendSeqMsgN loop");
591 DAI_PV(J);
592 DAI_PV(f_unnorm);
593 DAI_PV(_R[J][J.dual][_I]);
594 DAI_PV(f_unnorm * _R[J][J.dual][_I]);
595 }
596 #endif
597 incrSeqMsgM( i, J.iter, f_unnorm * R(J, J.dual, _I) );
598 }
599 }
600 }
601
602
603 void BBP::sendSeqMsgM( size_t j, size_t _I ) {
604 const Neighbor &I = _fg->nbV(j)[_I];
605
606 // DAI_PV(j);
607 // DAI_PV(I);
608 // DAI_PV(_adj_m_unnorm_jI);
609 // DAI_PV(_adj_m[j][_I]);
610 // DAI_PV(_bp_dual.zM(j,_I));
611
612 size_t _j = I.dual;
613 const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][_I];
614 Prob um( U(I, _j) );
615 const _ind_t &ind = _index(j, _I);
616 for( size_t x_I = 0; x_I < um.size(); x_I++ )
617 um[x_I] *= _adj_m_unnorm_jI[ind[x_I]];
618 um *= 1 - props.damping;
619 _adj_psi_F[I] += um;
620
621 /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
622 _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, _j) ) * Factor( _fg->var(j), _adj_m_unnorm[j][_I] )).p() * (1.0 - props.damping);
623 */
624
625 // DAI_DMSG("in sendSeqMsgM");
626 // DAI_PV(j);
627 // DAI_PV(I);
628 // DAI_PV(_I);
629 // DAI_PV(_fg->nbF(I).size());
630 foreach( const Neighbor &i, _fg->nbF(I) ) {
631 if( i.node != j ) {
632 const Prob &S = _S[i][i.dual][_j];
633 Prob msg( _fg->var(i).states(), 0.0 );
634 LOOP_ij(
635 msg[xi] += S[xij] * _adj_m_unnorm_jI[xj];
636 );
637 msg *= 1.0 - props.damping;
638 /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
639 Var vi = _fg->var(i);
640 Var vj = _fg->var(j);
641 msg = (Factor(VarSet(vi,vj), S) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p() * (1.0 - props.damping);
642 */
643 #if 0
644 if(msg.sumAbs() > pv_thresh) {
645 DAI_DMSG("in sendSeqMsgM loop");
646
647 DAI_PV(j);
648 DAI_PV(I);
649 DAI_PV(_I);
650 DAI_PV(_fg->nbF(I).size());
651 DAI_PV(_fg->factor(I).p());
652 DAI_PV(_S[i][i.dual][_j]);
653
654 DAI_PV(i);
655 DAI_PV(i.dual);
656 DAI_PV(msg);
657 DAI_PV(_fg->nbV(i).size());
658 }
659 #endif
660 assert( _fg->nbV(i)[i.dual].node == I );
661 sendSeqMsgN( i, i.dual, msg );
662 }
663 }
664 setSeqMsgM( j, _I, _adj_m[j][_I] * props.damping );
665 }
666
667
668 void BBP::RegenerateSeqMessageAdjoints() {
669 size_t nv = _fg->nrVars();
670 _adj_m.resize(nv);
671 _adj_m_unnorm.resize(nv);
672 _new_adj_m.resize(nv);
673 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
674 size_t n_i = _fg->nbV(i).size();
675 _adj_m[i].resize(n_i);
676 _adj_m_unnorm[i].resize(n_i);
677 _new_adj_m[i].resize(n_i);
678 foreach(const Neighbor& I, _fg->nbV(i)) {
679 // calculate adj_m
680 Prob prod(_adj_b_V_unnorm[i]);
681 assert(prod.size()==_fg->var(i).states());
682 foreach(const Neighbor& J, _fg->nbV(i)) {
683 if(size_t(J) != size_t(I)) {
684 prod *= _bp_dual.msgM(i,J.iter);
685 }
686 }
687 _adj_m[i][I.iter] = prod;
688 calcUnnormMsgM(i, I.iter);
689 _new_adj_m[i][I.iter] = Prob(_fg->var(i).states(),0.0);
690 }
691 }
692 for( size_t i = 0; i < _fg->nrVars(); i++ ) {
693 foreach(const Neighbor& I, _fg->nbV(i)) {
694 // calculate adj_n
695 Prob prod(_fg->factor(I).p());
696 prod *= _adj_b_F_unnorm[I];
697 foreach(const Neighbor& j, _fg->nbF(I)) {
698 if(i != j) {
699 Prob n_jI(_bp_dual.msgN(j,j.dual));
700 const _ind_t& ind = _index(j,j.dual);
701 // multiply prod with n_jI
702 for( size_t x_I = 0; x_I < prod.size(); x_I++ )
703 prod[x_I] *= n_jI[ind[x_I]];
704 }
705 }
706 Prob marg(_fg->var(i).states(), 0.0);
707 const _ind_t &ind = _index(i,I.iter);
708 for( size_t r = 0; r < prod.size(); r++ )
709 marg[ind[r]] += prod[r];
710 sendSeqMsgN(i,I.iter,marg);
711 }
712 }
713 }
714
715
716 void BBP::Regenerate() {
717 RegenerateInds();
718 RegenerateT();
719 RegenerateU();
720 RegenerateS();
721 RegenerateR();
722 RegenerateInputs();
723 RegeneratePsiAdjoints();
724 if( props.updates == Properties::UpdateType::PAR )
725 RegenerateParMessageAdjoints();
726 else
727 RegenerateSeqMessageAdjoints();
728 _iters = 0;
729 }
730
731
732 void BBP::calcNewN( size_t i, size_t _I ) {
733 _adj_psi_V[i] += T(i,_I) * _adj_n_unnorm[i][_I];
734 Prob &new_adj_n_iI = _new_adj_n[i][_I];
735 new_adj_n_iI = Prob( _fg->var(i).states(), 0.0 );
736 size_t I = _fg->nbV(i)[_I];
737 foreach( const Neighbor &j, _fg->nbF(I) )
738 if( j != i ) {
739 const Prob &p = _S[i][_I][j.iter];
740 const Prob &_adj_m_unnorm_jI = _adj_m_unnorm[j][j.dual];
741 LOOP_ij(
742 new_adj_n_iI[xi] += p[xij] * _adj_m_unnorm_jI[xj];
743 );
744 /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
745 Var vi = _fg->var(i);
746 Var vj = _fg->var(j);
747 new_adj_n_iI = (Factor(VarSet(vi, vj), p) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p();
748 */
749 }
750 }
751
752
753 void BBP::calcNewM( size_t i, size_t _I ) {
754 const Neighbor &I = _fg->nbV(i)[_I];
755 Prob p( U(I, I.dual) );
756 const Prob &adj = _adj_m_unnorm[i][_I];
757 const _ind_t &ind = _index(i,_I);
758 for( size_t x_I = 0; x_I < p.size(); x_I++ )
759 p[x_I] *= adj[ind[x_I]];
760 _adj_psi_F[I] += p;
761 /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
762 _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
763 */
764
765 _new_adj_m[i][_I] = Prob( _fg->var(i).states(), 0.0 );
766 foreach( const Neighbor &J, _fg->nbV(i) )
767 if( J != I )
768 _new_adj_m[i][_I] += _R[I][I.dual][J.iter] * _adj_n_unnorm[i][J.iter];
769 }
770
771
772 void BBP::calcUnnormMsgN( size_t i, size_t _I ) {
773 _adj_n_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgN(i,_I), _bp_dual.zN(i,_I), _adj_n[i][_I] );
774 }
775
776
777 void BBP::calcUnnormMsgM( size_t i, size_t _I ) {
778 _adj_m_unnorm[i][_I] = unnormAdjoint( _bp_dual.msgM(i,_I), _bp_dual.zM(i,_I), _adj_m[i][_I] );
779 }
780
781
782 void BBP::upMsgN( size_t i, size_t _I ) {
783 _adj_n[i][_I] = _new_adj_n[i][_I];
784 calcUnnormMsgN( i, _I );
785 }
786
787
788 void BBP::upMsgM( size_t i, size_t _I ) {
789 _adj_m[i][_I] = _new_adj_m[i][_I];
790 calcUnnormMsgM( i, _I );
791 }
792
793
794 void BBP::doParUpdate() {
795 for( size_t i = 0; i < _fg->nrVars(); i++ )
796 foreach( const Neighbor &I, _fg->nbV(i) ) {
797 calcNewM( i, I.iter );
798 calcNewN( i, I.iter );
799 }
800 for( size_t i = 0; i < _fg->nrVars(); i++ )
801 foreach( const Neighbor &I, _fg->nbV(i) ) {
802 upMsgM( i, I.iter );
803 upMsgN( i, I.iter );
804 }
805 }
806
807
808 Real BBP::getUnMsgMag() {
809 Real s = 0.0;
810 size_t e = 0;
811 for( size_t i = 0; i < _fg->nrVars(); i++ )
812 foreach( const Neighbor &I, _fg->nbV(i) ) {
813 s += _adj_m_unnorm[i][I.iter].sumAbs();
814 s += _adj_n_unnorm[i][I.iter].sumAbs();
815 e++;
816 }
817 return s / e;
818 }
819
820
821 void BBP::getMsgMags( Real &s, Real &new_s ) {
822 s = 0.0;
823 new_s = 0.0;
824 size_t e = 0;
825 for( size_t i = 0; i < _fg->nrVars(); i++ )
826 foreach( const Neighbor &I, _fg->nbV(i) ) {
827 s += _adj_m[i][I.iter].sumAbs();
828 s += _adj_n[i][I.iter].sumAbs();
829 new_s += _new_adj_m[i][I.iter].sumAbs();
830 new_s += _new_adj_n[i][I.iter].sumAbs();
831 e++;
832 }
833 s /= e;
834 new_s /= e;
835 }
836
837 // tuple<size_t,size_t,Real> BBP::getArgMaxPsi1Adj() {
838 // size_t argmax_var=0;
839 // size_t argmax_var_state=0;
840 // Real max_var=0;
841 // for( size_t i = 0; i < _fg->nrVars(); i++ ) {
842 // pair<size_t,Real> argmax_state = adj_psi_V(i).argmax();
843 // if(i==0 || argmax_state.second>max_var) {
844 // argmax_var = i;
845 // max_var = argmax_state.second;
846 // argmax_var_state = argmax_state.first;
847 // }
848 // }
849 // assert(/*0 <= argmax_var_state &&*/
850 // argmax_var_state < _fg->var(argmax_var).states());
851 // return tuple<size_t,size_t,Real>(argmax_var,argmax_var_state,max_var);
852 // }
853
854
855 void BBP::getArgmaxMsgM( size_t &out_i, size_t &out__I, Real &mag ) {
856 bool found = false;
857 for( size_t i = 0; i < _fg->nrVars(); i++ )
858 foreach( const Neighbor &I, _fg->nbV(i) ) {
859 Real thisMag = _adj_m[i][I.iter].sumAbs();
860 if( !found || mag < thisMag ) {
861 found = true;
862 mag = thisMag;
863 out_i = i;
864 out__I = I.iter;
865 }
866 }
867 assert( found );
868 }
869
870
871 Real BBP::getMaxMsgM() {
872 size_t dummy;
873 Real mag;
874 getArgmaxMsgM( dummy, dummy, mag );
875 return mag;
876 }
877
878
879 Real BBP::getTotalMsgM() {
880 Real mag = 0.0;
881 for( size_t i = 0; i < _fg->nrVars(); i++ )
882 foreach( const Neighbor &I, _fg->nbV(i) )
883 mag += _adj_m[i][I.iter].sumAbs();
884 return mag;
885 }
886
887
888 Real BBP::getTotalNewMsgM() {
889 Real mag = 0.0;
890 for( size_t i = 0; i < _fg->nrVars(); i++ )
891 foreach( const Neighbor &I, _fg->nbV(i) )
892 mag += _new_adj_m[i][I.iter].sumAbs();
893 return mag;
894 }
895
896
897 Real BBP::getTotalMsgN() {
898 Real mag = 0.0;
899 for( size_t i = 0; i < _fg->nrVars(); i++ )
900 foreach( const Neighbor &I, _fg->nbV(i) )
901 mag += _adj_n[i][I.iter].sumAbs();
902 return mag;
903 }
904
905
906 void BBP::run() {
907 typedef BBP::Properties::UpdateType UT;
908 Real tol = props.tol;
909 UT updates = props.updates;
910 Real tic=toc();
911 switch((size_t)updates) {
912 case UT::SEQ_MAX: {
913 size_t i, _I;
914 Real mag;
915 do {
916 _iters++;
917 getArgmaxMsgM(i,_I,mag);
918 sendSeqMsgM(i,_I);
919 } while(mag > tol && _iters < props.maxiter);
920
921 if(_iters >= props.maxiter) {
922 cerr << "Warning: BBP didn't converge in " << _iters
923 << " iterations (greatest message magnitude = " << mag << ")"
924 << endl;
925 }
926 break;
927 }
928 case UT::SEQ_FIX: {
929 Real mag;
930 do {
931 _iters++;
932 mag = getTotalMsgM();
933 if(mag < tol) break;
934
935 for( size_t i = 0; i < _fg->nrVars(); i++ )
936 foreach(const Neighbor &I, _fg->nbV(i))
937 sendSeqMsgM(i, I.iter);
938 for( size_t i = 0; i < _fg->nrVars(); i++ )
939 foreach(const Neighbor &I, _fg->nbV(i))
940 updateSeqMsgM(i, I.iter);
941 } while(mag > tol && _iters < props.maxiter);
942
943 if(_iters >= props.maxiter) {
944 cerr << "Warning: BBP didn't converge in " << _iters
945 << " iterations (greatest message magnitude = " << mag << ")"
946 << endl;
947 break;
948 }
949 break;
950 }
951 case UT::SEQ_BP_REV:
952 case UT::SEQ_BP_FWD: {
953 const BP *bp = static_cast<const BP*>(_ia);
954 vector<pair<size_t, size_t> > sentMessages = bp->getSentMessages();
955 size_t totalMessages = sentMessages.size();
956 if(totalMessages==0) {
957 cerr << "Asked for updates = " << updates << " but no BP messages; did you forget to set recordSentMessages?" << endl;
958 DAI_THROW(INTERNAL_ERROR);
959 }
960 if(updates==UT::SEQ_BP_FWD) {
961 reverse(sentMessages.begin(), sentMessages.end());
962 }
963 // DAI_PV(sentMessages.size());
964 // DAI_PV(_iters);
965 // DAI_PV(props.maxiter);
966 while(sentMessages.size()>0 && _iters < props.maxiter) {
967 // DAI_PV(sentMessages.size());
968 // DAI_PV(_iters);
969 _iters++;
970 pair<size_t, size_t> e = sentMessages.back();
971 sentMessages.pop_back();
972 size_t i = e.first, _I = e.second;
973 sendSeqMsgM(i,_I);
974 }
975 if(_iters >= props.maxiter) {
976 cerr << "Warning: BBP updates limited to " << props.maxiter
977 << " iterations, but using UpdateType " << updates
978 << " with " << totalMessages << " messages"
979 << endl;
980 }
981 break;
982 }
983 case UT::PAR: {
984 do {
985 _iters++;
986 doParUpdate();
987 } while((_iters < 2 || getUnMsgMag() > tol) && _iters < props.maxiter);
988 if(_iters==props.maxiter) {
989 Real s, new_s;
990 getMsgMags(s,new_s);
991 cerr << "Warning: BBP didn't converge in " << _iters
992 << " iterations (unnorm message magnitude = " << getUnMsgMag()
993 << ", norm message mags = " << s << " -> " << new_s
994 << ")" << endl;
995 }
996 break;
997 }
998 }
999 if(props.verbose >= 3) {
1000 cerr << "BBP::run() took " << toc()-tic << " seconds "
1001 << doneIters() << " iterations" << endl;
1002 }
1003 }
1004
1005
1006 bool needGibbsState(bbp_cfn_t cfn) {
1007 switch((size_t)cfn) {
1008 case bbp_cfn_t::cfn_gibbs_b:
1009 case bbp_cfn_t::cfn_gibbs_b2:
1010 case bbp_cfn_t::cfn_gibbs_exp:
1011 case bbp_cfn_t::cfn_gibbs_b_factor:
1012 case bbp_cfn_t::cfn_gibbs_b2_factor:
1013 case bbp_cfn_t::cfn_gibbs_exp_factor:
1014 return true;
1015 default:
1016 return false;
1017 }
1018 }
1019
1020
1021 vector<size_t> getGibbsState(const InfAlg& ia, size_t iters) {
1022 PropertySet gibbsProps;
1023 gibbsProps.Set("iters", iters);
1024 gibbsProps.Set("verbose", size_t(0));
1025 Gibbs gibbs(ia.fg(), gibbsProps);
1026 gibbs.run();
1027 return gibbs.state();
1028 }
1029
1030
1031 double numericBBPTest(const InfAlg& bp, const vector<size_t> *state, const PropertySet& bbp_props, bbp_cfn_t cfn, double h) {
1032 // calculate the value of the unperturbed cost function
1033 Real cf0 = getCostFn(bp, cfn, state);
1034
1035 // run BBP to estimate adjoints
1036 BBP bbp(&bp,bbp_props);
1037 initBBPCostFnAdj(bbp, bp, cfn, state);
1038 bbp.run();
1039
1040 Real d=0;
1041 const FactorGraph& fg = bp.fg();
1042
1043 if(1) {
1044 // verify bbp.adj_psi_V
1045
1046 // for each variable i
1047 for(size_t i=0; i<fg.nrVars(); i++) {
1048 vector<double> adj_est;
1049 // for each value xi
1050 for(size_t xi=0; xi<fg.var(i).states(); xi++) {
1051 // Clone 'bp' (which may be any InfAlg)
1052 InfAlg *bp_prb = bp.clone();
1053
1054 // perturb it
1055 size_t n = bp_prb->fg().var(i).states();
1056 Prob psi_1_prb(n,1.0);
1057 psi_1_prb[xi] += h;
1058 // psi_1_prb.normalize();
1059 size_t I = bp_prb->fg().nbV(i)[0]; // use first factor in list of neighbors of i
1060 bp_prb->fg().factor(I) *= Factor(bp_prb->fg().var(i), psi_1_prb);
1061
1062 // call 'init' on the perturbed variables
1063 bp_prb->init(bp_prb->fg().var(i));
1064
1065 // run copy to convergence
1066 bp_prb->run();
1067
1068 // calculate new value of cost function
1069 Real cf_prb = getCostFn(*bp_prb, cfn, state);
1070
1071 // use to estimate adjoint for i
1072 adj_est.push_back((cf_prb-cf0)/h);
1073
1074 // free cloned InfAlg
1075 delete bp_prb;
1076 }
1077 Prob p_adj_est(adj_est.begin(), adj_est.end());
1078 // compare this numerical estimate to the BBP estimate; sum the distances
1079 cerr << "i: " << i
1080 << ", p_adj_est: " << p_adj_est
1081 << ", bbp.adj_psi_V(i): " << bbp.adj_psi_V(i) << endl;
1082 d += dist(p_adj_est, bbp.adj_psi_V(i), Prob::DISTL1);
1083 }
1084 }
1085 /* if(1) {
1086 // verify bbp.adj_n and bbp.adj_m
1087
1088 // We actually want to check the responsiveness of objective
1089 // function to changes in the final messages. But at the end of a
1090 // BBP run, the message adjoints are for the initial messages (and
1091 // they should be close to zero, see paper). So this resets the
1092 // BBP adjoints to the refer to the desired final messages
1093 bbp.RegenerateMessageAdjoints();
1094
1095 // for each variable i
1096 for(size_t i=0; i<bp_dual.nrVars(); i++) {
1097 // for each factor I ~ i
1098 foreach(size_t I, bp_dual.nbV(i)) {
1099 vector<double> adj_n_est;
1100 // for each value xi
1101 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1102 BP_dual bp_dual_prb(bp_dual);
1103 // make h-sized change to newMsgN
1104 bp_dual_prb.newMsgN(i,I)[xi] += h;
1105 // recalculate beliefs
1106 bp_dual_prb.CalcBeliefs();
1107 // get cost function value
1108 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1109 // add it to list of adjoints
1110 adj_n_est.push_back((cf_prb-cf0)/h);
1111 }
1112
1113 vector<double> adj_m_est;
1114 // for each value xi
1115 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1116 BP_dual bp_dual_prb(bp_dual);
1117 // make h-sized change to newMsgM
1118 bp_dual_prb.newMsgM(I,i)[xi] += h;
1119 // recalculate beliefs
1120 bp_dual_prb.CalcBeliefs();
1121 // get cost function value
1122 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1123 // add it to list of adjoints
1124 adj_m_est.push_back((cf_prb-cf0)/h);
1125 }
1126
1127 Prob p_adj_n_est(adj_n_est.begin(), adj_n_est.end());
1128 // compare this numerical estimate to the BBP estimate; sum the distances
1129 cerr << "i: " << i << ", I: " << I
1130 << ", adj_n_est: " << p_adj_n_est
1131 << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
1132 d += dist(p_adj_n_est, bbp.adj_n(i,I), Prob::DISTL1);
1133
1134 Prob p_adj_m_est(adj_m_est.begin(), adj_m_est.end());
1135 // compare this numerical estimate to the BBP estimate; sum the distances
1136 cerr << "i: " << i << ", I: " << I
1137 << ", adj_m_est: " << p_adj_m_est
1138 << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
1139 d += dist(p_adj_m_est, bbp.adj_m(I,i), Prob::DISTL1);
1140 }
1141 }
1142 }
1143 */
1144 /* if(0) {
1145 // verify bbp.adj_b_V
1146 for(size_t i=0; i<bp_dual.nrVars(); i++) {
1147 vector<double> adj_b_V_est;
1148 // for each value xi
1149 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1150 BP_dual bp_dual_prb(bp_dual);
1151
1152 // make h-sized change to b_1(i)[x_i]
1153 bp_dual_prb._beliefs.b1[i][xi] += h;
1154
1155 // get cost function value
1156 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1157
1158 // add it to list of adjoints
1159 adj_b_V_est.push_back((cf_prb-cf0)/h);
1160 }
1161 Prob p_adj_b_V_est(adj_b_V_est.begin(), adj_b_V_est.end());
1162 // compare this numerical estimate to the BBP estimate; sum the distances
1163 cerr << "i: " << i
1164 << ", adj_b_V_est: " << p_adj_b_V_est
1165 << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
1166 d += dist(p_adj_b_V_est, bbp.adj_b_V(i), Prob::DISTL1);
1167 }
1168 }
1169 */
1170
1171 // return total of distances
1172 return d;
1173 }
1174
1175
1176 } // end of namespace dai
1177
1178
1179 /* {{{ GENERATED CODE: DO NOT EDIT. Created by
1180 ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp
1181 */
1182 namespace dai {
1183
1184 void BBP::Properties::set(const PropertySet &opts)
1185 {
1186 const std::set<PropertyKey> &keys = opts.allKeys();
1187 std::set<PropertyKey>::const_iterator i;
1188 bool die=false;
1189 for(i=keys.begin(); i!=keys.end(); i++) {
1190 if(*i == "verbose") continue;
1191 if(*i == "maxiter") continue;
1192 if(*i == "tol") continue;
1193 if(*i == "damping") continue;
1194 if(*i == "updates") continue;
1195 if(*i == "clean_updates") continue;
1196 cerr << "BBP: Unknown property " << *i << endl;
1197 die=true;
1198 }
1199 if(die) {
1200 DAI_THROW(UNKNOWN_PROPERTY_TYPE);
1201 }
1202 if(!opts.hasKey("verbose")) {
1203 cerr << "BBP: Missing property \"verbose\" for method \"BBP\"" << endl;
1204 die=true;
1205 }
1206 if(!opts.hasKey("maxiter")) {
1207 cerr << "BBP: Missing property \"maxiter\" for method \"BBP\"" << endl;
1208 die=true;
1209 }
1210 if(!opts.hasKey("tol")) {
1211 cerr << "BBP: Missing property \"tol\" for method \"BBP\"" << endl;
1212 die=true;
1213 }
1214 if(!opts.hasKey("damping")) {
1215 cerr << "BBP: Missing property \"damping\" for method \"BBP\"" << endl;
1216 die=true;
1217 }
1218 if(!opts.hasKey("updates")) {
1219 cerr << "BBP: Missing property \"updates\" for method \"BBP\"" << endl;
1220 die=true;
1221 }
1222 if(!opts.hasKey("clean_updates")) {
1223 cerr << "BBP: Missing property \"clean_updates\" for method \"BBP\"" << endl;
1224 die=true;
1225 }
1226 if(die) {
1227 DAI_THROW(NOT_ALL_PROPERTIES_SPECIFIED);
1228 }
1229 verbose = opts.getStringAs<size_t>("verbose");
1230 maxiter = opts.getStringAs<size_t>("maxiter");
1231 tol = opts.getStringAs<double>("tol");
1232 damping = opts.getStringAs<double>("damping");
1233 updates = opts.getStringAs<UpdateType>("updates");
1234 clean_updates = opts.getStringAs<bool>("clean_updates");
1235 }
1236 PropertySet BBP::Properties::get() const {
1237 PropertySet opts;
1238 opts.Set("verbose", verbose);
1239 opts.Set("maxiter", maxiter);
1240 opts.Set("tol", tol);
1241 opts.Set("damping", damping);
1242 opts.Set("updates", updates);
1243 opts.Set("clean_updates", clean_updates);
1244 return opts;
1245 }
1246 string BBP::Properties::toString() const {
1247 stringstream s(stringstream::out);
1248 s << "[";
1249 s << "verbose=" << verbose << ",";
1250 s << "maxiter=" << maxiter << ",";
1251 s << "tol=" << tol << ",";
1252 s << "damping=" << damping << ",";
1253 s << "updates=" << updates << ",";
1254 s << "clean_updates=" << clean_updates;
1255 s << "]";
1256 return s.str();
1257 }
1258 } // end of namespace dai
1259 /* }}} END OF GENERATED CODE */