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