1 /* This file is part of libDAI - http://www.libdai.org/
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
9 #include <dai/dai_config.h>
15 #include <dai/gibbs.h>
17 #include <dai/bipgraph.h>
26 /// Returns the entry of the I'th factor corresponding to a global state
27 size_t getFactorEntryForState( const FactorGraph
&fg
, size_t I
, const vector
<size_t> &state
) {
29 for( int _j
= fg
.nbF(I
).size() - 1; _j
>= 0; _j
-- ) {
30 // note that iterating over nbF(I) yields the same ordering
31 // of variables as iterating over factor(I).vars()
32 size_t j
= fg
.nbF(I
)[_j
];
33 f_entry
*= fg
.var(j
).states();
40 bool BBPCostFunction::needGibbsState() const {
41 switch( (size_t)(*this) ) {
45 case CFN_GIBBS_B_FACTOR
:
46 case CFN_GIBBS_B2_FACTOR
:
47 case CFN_GIBBS_EXP_FACTOR
:
55 Real
BBPCostFunction::evaluate( const InfAlg
&ia
, const vector
<size_t> *stateP
) const {
57 const FactorGraph
&fg
= ia
.fg();
59 switch( (size_t)(*this) ) {
60 case CFN_BETHE_ENT
: // ignores state
63 case CFN_VAR_ENT
: // ignores state
64 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
65 cf
+= -ia
.beliefV(i
).entropy();
67 case CFN_FACTOR_ENT
: // ignores state
68 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
69 cf
+= -ia
.beliefF(I
).entropy();
74 DAI_ASSERT( stateP
!= NULL
);
75 vector
<size_t> state
= *stateP
;
76 DAI_ASSERT( state
.size() == fg
.nrVars() );
77 for( size_t i
= 0; i
< fg
.nrVars(); i
++ ) {
78 Real b
= ia
.beliefV(i
)[state
[i
]];
79 switch( (size_t)(*this) ) {
90 DAI_THROW(UNKNOWN_ENUM_VALUE
);
94 } case CFN_GIBBS_B_FACTOR
:
95 case CFN_GIBBS_B2_FACTOR
:
96 case CFN_GIBBS_EXP_FACTOR
: {
97 DAI_ASSERT( stateP
!= NULL
);
98 vector
<size_t> state
= *stateP
;
99 DAI_ASSERT( state
.size() == fg
.nrVars() );
100 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
101 size_t x_I
= getFactorEntryForState( fg
, I
, state
);
102 Real b
= ia
.beliefF(I
)[x_I
];
103 switch( (size_t)(*this) ) {
104 case CFN_GIBBS_B_FACTOR
:
107 case CFN_GIBBS_B2_FACTOR
:
110 case CFN_GIBBS_EXP_FACTOR
:
114 DAI_THROW(UNKNOWN_ENUM_VALUE
);
119 DAI_THROWE(UNKNOWN_ENUM_VALUE
, "Unknown cost function " + std::string(*this));
125 #define LOOP_ij(body) { \
126 size_t i_states = _fg->var(i).states(); \
127 size_t j_states = _fg->var(j).states(); \
128 if(_fg->var(i) > _fg->var(j)) { \
130 for(size_t xi=0; xi<i_states; xi++) { \
131 for(size_t xj=0; xj<j_states; xj++) { \
138 for(size_t xj=0; xj<j_states; xj++) { \
139 for(size_t xi=0; xi<i_states; xi++) { \
148 void BBP::RegenerateInds() {
149 // initialise _indices
150 // typedef std::vector<size_t> _ind_t;
151 // std::vector<std::vector<_ind_t> > _indices;
152 _indices
.resize( _fg
->nrVars() );
153 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
155 _indices
[i
].reserve( _fg
->nbV(i
).size() );
156 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
158 index
.reserve( _fg
->factor(I
).nrStates() );
159 for( IndexFor
k(_fg
->var(i
), _fg
->factor(I
).vars()); k
.valid(); ++k
)
160 index
.push_back( k
);
161 _indices
[i
].push_back( index
);
167 void BBP::RegenerateT() {
169 _Tmsg
.resize( _fg
->nrVars() );
170 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
171 _Tmsg
[i
].resize( _fg
->nbV(i
).size() );
172 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
173 Prob
prod( _fg
->var(i
).states(), 1.0 );
174 bforeach( const Neighbor
&J
, _fg
->nbV(i
) )
175 if( J
.node
!= I
.node
)
176 prod
*= _bp_dual
.msgM( i
, J
.iter
);
177 _Tmsg
[i
][I
.iter
] = prod
;
183 void BBP::RegenerateU() {
185 _Umsg
.resize( _fg
->nrFactors() );
186 for( size_t I
= 0; I
< _fg
->nrFactors(); I
++ ) {
187 _Umsg
[I
].resize( _fg
->nbF(I
).size() );
188 bforeach( const Neighbor
&i
, _fg
->nbF(I
) ) {
189 Prob
prod( _fg
->factor(I
).nrStates(), 1.0 );
190 bforeach( const Neighbor
&j
, _fg
->nbF(I
) )
191 if( i
.node
!= j
.node
) {
192 Prob
n_jI( _bp_dual
.msgN( j
, j
.dual
) );
193 const _ind_t
&ind
= _index( j
, j
.dual
);
194 // multiply prod by n_jI
195 for( size_t x_I
= 0; x_I
< prod
.size(); x_I
++ )
196 prod
.set( x_I
, prod
[x_I
] * n_jI
[ind
[x_I
]] );
198 _Umsg
[I
][i
.iter
] = prod
;
204 void BBP::RegenerateS() {
206 _Smsg
.resize( _fg
->nrVars() );
207 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
208 _Smsg
[i
].resize( _fg
->nbV(i
).size() );
209 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
210 _Smsg
[i
][I
.iter
].resize( _fg
->nbF(I
).size() );
211 bforeach( const Neighbor
&j
, _fg
->nbF(I
) )
213 Factor
prod( _fg
->factor(I
) );
214 bforeach( const Neighbor
&k
, _fg
->nbF(I
) ) {
215 if( k
!= i
&& k
.node
!= j
.node
) {
216 const _ind_t
&ind
= _index( k
, k
.dual
);
217 Prob
p( _bp_dual
.msgN( k
, k
.dual
) );
218 for( size_t x_I
= 0; x_I
< prod
.nrStates(); x_I
++ )
219 prod
.set( x_I
, prod
[x_I
] * p
[ind
[x_I
]] );
222 // "Marginalize" onto i|j (unnormalized)
224 marg
= prod
.marginal( VarSet(_fg
->var(i
), _fg
->var(j
)), false ).p();
225 _Smsg
[i
][I
.iter
][j
.iter
] = marg
;
232 void BBP::RegenerateR() {
234 _Rmsg
.resize( _fg
->nrFactors() );
235 for( size_t I
= 0; I
< _fg
->nrFactors(); I
++ ) {
236 _Rmsg
[I
].resize( _fg
->nbF(I
).size() );
237 bforeach( const Neighbor
&i
, _fg
->nbF(I
) ) {
238 _Rmsg
[I
][i
.iter
].resize( _fg
->nbV(i
).size() );
239 bforeach( const Neighbor
&J
, _fg
->nbV(i
) ) {
241 Prob
prod( _fg
->var(i
).states(), 1.0 );
242 bforeach( const Neighbor
&K
, _fg
->nbV(i
) )
243 if( K
.node
!= I
&& K
.node
!= J
.node
)
244 prod
*= _bp_dual
.msgM( i
, K
.iter
);
245 _Rmsg
[I
][i
.iter
][J
.iter
] = prod
;
253 void BBP::RegenerateInputs() {
254 _adj_b_V_unnorm
.clear();
255 _adj_b_V_unnorm
.reserve( _fg
->nrVars() );
256 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
257 _adj_b_V_unnorm
.push_back( unnormAdjoint( _bp_dual
.beliefV(i
).p(), _bp_dual
.beliefVZ(i
), _adj_b_V
[i
] ) );
258 _adj_b_F_unnorm
.clear();
259 _adj_b_F_unnorm
.reserve( _fg
->nrFactors() );
260 for( size_t I
= 0; I
< _fg
->nrFactors(); I
++ )
261 _adj_b_F_unnorm
.push_back( unnormAdjoint( _bp_dual
.beliefF(I
).p(), _bp_dual
.beliefFZ(I
), _adj_b_F
[I
] ) );
265 void BBP::RegeneratePsiAdjoints() {
267 _adj_psi_V
.reserve( _fg
->nrVars() );
268 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
269 Prob
p( _adj_b_V_unnorm
[i
] );
270 DAI_ASSERT( p
.size() == _fg
->var(i
).states() );
271 bforeach( const Neighbor
&I
, _fg
->nbV(i
) )
272 p
*= _bp_dual
.msgM( i
, I
.iter
);
273 p
+= _init_adj_psi_V
[i
];
274 _adj_psi_V
.push_back( p
);
277 _adj_psi_F
.reserve( _fg
->nrFactors() );
278 for( size_t I
= 0; I
< _fg
->nrFactors(); I
++ ) {
279 Prob
p( _adj_b_F_unnorm
[I
] );
280 DAI_ASSERT( p
.size() == _fg
->factor(I
).nrStates() );
281 bforeach( const Neighbor
&i
, _fg
->nbF(I
) ) {
282 Prob
n_iI( _bp_dual
.msgN( i
, i
.dual
) );
283 const _ind_t
& ind
= _index( i
, i
.dual
);
284 // multiply prod with n_jI
285 for( size_t x_I
= 0; x_I
< p
.size(); x_I
++ )
286 p
.set( x_I
, p
[x_I
] * n_iI
[ind
[x_I
]] );
288 p
+= _init_adj_psi_F
[I
];
289 _adj_psi_F
.push_back( p
);
294 void BBP::RegenerateParMessageAdjoints() {
295 size_t nv
= _fg
->nrVars();
298 _adj_n_unnorm
.resize( nv
);
299 _adj_m_unnorm
.resize( nv
);
300 _new_adj_n
.resize( nv
);
301 _new_adj_m
.resize( nv
);
302 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
303 size_t n_i
= _fg
->nbV(i
).size();
304 _adj_n
[i
].resize( n_i
);
305 _adj_m
[i
].resize( n_i
);
306 _adj_n_unnorm
[i
].resize( n_i
);
307 _adj_m_unnorm
[i
].resize( n_i
);
308 _new_adj_n
[i
].resize( n_i
);
309 _new_adj_m
[i
].resize( n_i
);
310 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
312 Prob
prod( _fg
->factor(I
).p() );
313 prod
*= _adj_b_F_unnorm
[I
];
314 bforeach( const Neighbor
&j
, _fg
->nbF(I
) )
316 Prob
n_jI( _bp_dual
.msgN( j
, j
.dual
) );
317 const _ind_t
&ind
= _index( j
, j
.dual
);
318 // multiply prod with n_jI
319 for( size_t x_I
= 0; x_I
< prod
.size(); x_I
++ )
320 prod
.set( x_I
, prod
[x_I
] * n_jI
[ind
[x_I
]] );
322 Prob
marg( _fg
->var(i
).states(), 0.0 );
323 const _ind_t
&ind
= _index( i
, I
.iter
);
324 for( size_t r
= 0; r
< prod
.size(); r
++ )
325 marg
.set( ind
[r
], marg
[ind
[r
]] + prod
[r
] );
326 _new_adj_n
[i
][I
.iter
] = marg
;
331 Prob
prod( _adj_b_V_unnorm
[i
] );
332 DAI_ASSERT( prod
.size() == _fg
->var(i
).states() );
333 bforeach( const Neighbor
&J
, _fg
->nbV(i
) )
334 if( J
.node
!= I
.node
)
335 prod
*= _bp_dual
.msgM(i
,J
.iter
);
336 _new_adj_m
[i
][I
.iter
] = prod
;
344 void BBP::RegenerateSeqMessageAdjoints() {
345 size_t nv
= _fg
->nrVars();
347 _adj_m_unnorm
.resize( nv
);
348 _new_adj_m
.resize( nv
);
349 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
350 size_t n_i
= _fg
->nbV(i
).size();
351 _adj_m
[i
].resize( n_i
);
352 _adj_m_unnorm
[i
].resize( n_i
);
353 _new_adj_m
[i
].resize( n_i
);
354 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
356 Prob
prod( _adj_b_V_unnorm
[i
] );
357 DAI_ASSERT( prod
.size() == _fg
->var(i
).states() );
358 bforeach( const Neighbor
&J
, _fg
->nbV(i
) )
359 if( J
.node
!= I
.node
)
360 prod
*= _bp_dual
.msgM( i
, J
.iter
);
361 _adj_m
[i
][I
.iter
] = prod
;
362 calcUnnormMsgM( i
, I
.iter
);
363 _new_adj_m
[i
][I
.iter
] = Prob( _fg
->var(i
).states(), 0.0 );
366 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ ) {
367 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
369 Prob
prod( _fg
->factor(I
).p() );
370 prod
*= _adj_b_F_unnorm
[I
];
371 bforeach( const Neighbor
&j
, _fg
->nbF(I
) )
373 Prob
n_jI( _bp_dual
.msgN( j
, j
.dual
) );
374 const _ind_t
& ind
= _index( j
, j
.dual
);
375 // multiply prod with n_jI
376 for( size_t x_I
= 0; x_I
< prod
.size(); x_I
++ )
377 prod
.set( x_I
, prod
[x_I
] * n_jI
[ind
[x_I
]] );
379 Prob
marg( _fg
->var(i
).states(), 0.0 );
380 const _ind_t
&ind
= _index( i
, I
.iter
);
381 for( size_t r
= 0; r
< prod
.size(); r
++ )
382 marg
.set( ind
[r
], marg
[ind
[r
]] + prod
[r
] );
383 sendSeqMsgN( i
, I
.iter
,marg
);
389 void BBP::Regenerate() {
396 RegeneratePsiAdjoints();
397 if( props
.updates
== Properties::UpdateType::PAR
)
398 RegenerateParMessageAdjoints();
400 RegenerateSeqMessageAdjoints();
405 void BBP::calcNewN( size_t i
, size_t _I
) {
406 _adj_psi_V
[i
] += T(i
,_I
) * _adj_n_unnorm
[i
][_I
];
407 Prob
&new_adj_n_iI
= _new_adj_n
[i
][_I
];
408 new_adj_n_iI
= Prob( _fg
->var(i
).states(), 0.0 );
409 size_t I
= _fg
->nbV(i
)[_I
];
410 bforeach( const Neighbor
&j
, _fg
->nbF(I
) )
412 const Prob
&p
= _Smsg
[i
][_I
][j
.iter
];
413 const Prob
&_adj_m_unnorm_jI
= _adj_m_unnorm
[j
][j
.dual
];
415 new_adj_n_iI
.set( xi
, new_adj_n_iI
[xi
] + p
[xij
] * _adj_m_unnorm_jI
[xj
] );
417 /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
418 Var vi = _fg->var(i);
419 Var vj = _fg->var(j);
420 new_adj_n_iI = (Factor(VarSet(vi, vj), p) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p();
426 void BBP::calcNewM( size_t i
, size_t _I
) {
427 const Neighbor
&I
= _fg
->nbV(i
)[_I
];
428 Prob
p( U(I
, I
.dual
) );
429 const Prob
&adj
= _adj_m_unnorm
[i
][_I
];
430 const _ind_t
&ind
= _index(i
,_I
);
431 for( size_t x_I
= 0; x_I
< p
.size(); x_I
++ )
432 p
.set( x_I
, p
[x_I
] * adj
[ind
[x_I
]] );
434 /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
435 _adj_psi_F[I] += (Factor( _fg->factor(I).vars(), U(I, I.dual) ) * Factor( _fg->var(i), _adj_m_unnorm[i][_I] )).p();
438 _new_adj_m
[i
][_I
] = Prob( _fg
->var(i
).states(), 0.0 );
439 bforeach( const Neighbor
&J
, _fg
->nbV(i
) )
441 _new_adj_m
[i
][_I
] += _Rmsg
[I
][I
.dual
][J
.iter
] * _adj_n_unnorm
[i
][J
.iter
];
445 void BBP::calcUnnormMsgN( size_t i
, size_t _I
) {
446 _adj_n_unnorm
[i
][_I
] = unnormAdjoint( _bp_dual
.msgN(i
,_I
), _bp_dual
.zN(i
,_I
), _adj_n
[i
][_I
] );
450 void BBP::calcUnnormMsgM( size_t i
, size_t _I
) {
451 _adj_m_unnorm
[i
][_I
] = unnormAdjoint( _bp_dual
.msgM(i
,_I
), _bp_dual
.zM(i
,_I
), _adj_m
[i
][_I
] );
455 void BBP::upMsgN( size_t i
, size_t _I
) {
456 _adj_n
[i
][_I
] = _new_adj_n
[i
][_I
];
457 calcUnnormMsgN( i
, _I
);
461 void BBP::upMsgM( size_t i
, size_t _I
) {
462 _adj_m
[i
][_I
] = _new_adj_m
[i
][_I
];
463 calcUnnormMsgM( i
, _I
);
467 void BBP::doParUpdate() {
468 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
469 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
470 calcNewM( i
, I
.iter
);
471 calcNewN( i
, I
.iter
);
473 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
474 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
481 void BBP::incrSeqMsgM( size_t i
, size_t _I
, const Prob
&p
) {
482 /* if( props.clean_updates )
483 _new_adj_m[i][_I] += p;
486 calcUnnormMsgM(i
, _I
);
497 void BBP::updateSeqMsgM( size_t i, size_t _I ) {
498 if( props.clean_updates ) {
500 if(_new_adj_m[i][_I].sumAbs() > pv_thresh ||
501 _adj_m[i][_I].sumAbs() > pv_thresh) {
503 DAI_DMSG("in updateSeqMsgM");
506 DAI_PV(_adj_m[i][_I]);
507 DAI_PV(_new_adj_m[i][_I]);
510 _adj_m[i][_I] += _new_adj_m[i][_I];
511 calcUnnormMsgM( i, _I );
512 _new_adj_m[i][_I].fill( 0.0 );
517 void BBP::setSeqMsgM( size_t i
, size_t _I
, const Prob
&p
) {
519 calcUnnormMsgM( i
, _I
);
523 void BBP::sendSeqMsgN( size_t i
, size_t _I
, const Prob
&f
) {
524 Prob f_unnorm
= unnormAdjoint( _bp_dual
.msgN(i
,_I
), _bp_dual
.zN(i
,_I
), f
);
525 const Neighbor
&I
= _fg
->nbV(i
)[_I
];
526 DAI_ASSERT( I
.iter
== _I
);
527 _adj_psi_V
[i
] += f_unnorm
* T( i
, _I
);
529 if(f_unnorm
.sumAbs() > pv_thresh
) {
530 DAI_DMSG("in sendSeqMsgN");
534 DAI_PV(_bp_dual
.msgN(i
,_I
));
535 DAI_PV(_bp_dual
.zN(i
,_I
));
536 DAI_PV(_bp_dual
.msgM(i
,_I
));
537 DAI_PV(_bp_dual
.zM(i
,_I
));
538 DAI_PV(_fg
->factor(I
).p());
541 bforeach( const Neighbor
&J
, _fg
->nbV(i
) ) {
542 if( J
.node
!= I
.node
) {
544 if(f_unnorm
.sumAbs() > pv_thresh
) {
545 DAI_DMSG("in sendSeqMsgN loop");
548 DAI_PV(_Rmsg
[J
][J
.dual
][_I
]);
549 DAI_PV(f_unnorm
* _Rmsg
[J
][J
.dual
][_I
]);
552 incrSeqMsgM( i
, J
.iter
, f_unnorm
* R(J
, J
.dual
, _I
) );
558 void BBP::sendSeqMsgM( size_t j
, size_t _I
) {
559 const Neighbor
&I
= _fg
->nbV(j
)[_I
];
563 // DAI_PV(_adj_m_unnorm_jI);
564 // DAI_PV(_adj_m[j][_I]);
565 // DAI_PV(_bp_dual.zM(j,_I));
568 const Prob
&_adj_m_unnorm_jI
= _adj_m_unnorm
[j
][_I
];
570 const _ind_t
&ind
= _index(j
, _I
);
571 for( size_t x_I
= 0; x_I
< um
.size(); x_I
++ )
572 um
.set( x_I
, um
[x_I
] * _adj_m_unnorm_jI
[ind
[x_I
]] );
573 um
*= 1 - props
.damping
;
576 /* THE FOLLOWING WOULD BE SLIGHTLY SLOWER:
577 _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);
580 // DAI_DMSG("in sendSeqMsgM");
584 // DAI_PV(_fg->nbF(I).size());
585 bforeach( const Neighbor
&i
, _fg
->nbF(I
) ) {
587 const Prob
&S
= _Smsg
[i
][i
.dual
][_j
];
588 Prob
msg( _fg
->var(i
).states(), 0.0 );
590 msg
.set( xi
, msg
[xi
] + S
[xij
] * _adj_m_unnorm_jI
[xj
] );
592 msg
*= 1.0 - props
.damping
;
593 /* THE FOLLOWING WOULD BE ABOUT TWICE AS SLOW:
594 Var vi = _fg->var(i);
595 Var vj = _fg->var(j);
596 msg = (Factor(VarSet(vi,vj), S) * Factor(vj,_adj_m_unnorm_jI)).marginal(vi,false).p() * (1.0 - props.damping);
599 if(msg
.sumAbs() > pv_thresh
) {
600 DAI_DMSG("in sendSeqMsgM loop");
605 DAI_PV(_fg
->nbF(I
).size());
606 DAI_PV(_fg
->factor(I
).p());
607 DAI_PV(_Smsg
[i
][i
.dual
][_j
]);
612 DAI_PV(_fg
->nbV(i
).size());
615 DAI_ASSERT( _fg
->nbV(i
)[i
.dual
].node
== I
);
616 sendSeqMsgN( i
, i
.dual
, msg
);
619 setSeqMsgM( j
, _I
, _adj_m
[j
][_I
] * props
.damping
);
623 Prob
BBP::unnormAdjoint( const Prob
&w
, Real Z_w
, const Prob
&adj_w
) {
624 DAI_ASSERT( w
.size() == adj_w
.size() );
625 Prob
adj_w_unnorm( w
.size(), 0.0 );
627 for( size_t i
= 0; i
< w
.size(); i
++ )
628 s
+= w
[i
] * adj_w
[i
];
629 for( size_t i
= 0; i
< w
.size(); i
++ )
630 adj_w_unnorm
.set( i
, (adj_w
[i
] - s
) / Z_w
);
632 // THIS WOULD BE ABOUT 50% SLOWER: return (adj_w - (w * adj_w).sum()) / Z_w;
636 Real
BBP::getUnMsgMag() {
639 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
640 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
641 s
+= _adj_m_unnorm
[i
][I
.iter
].sumAbs();
642 s
+= _adj_n_unnorm
[i
][I
.iter
].sumAbs();
649 void BBP::getMsgMags( Real
&s
, Real
&new_s
) {
653 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
654 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
655 s
+= _adj_m
[i
][I
.iter
].sumAbs();
656 s
+= _adj_n
[i
][I
.iter
].sumAbs();
657 new_s
+= _new_adj_m
[i
][I
.iter
].sumAbs();
658 new_s
+= _new_adj_n
[i
][I
.iter
].sumAbs();
665 // tuple<size_t,size_t,Real> BBP::getArgMaxPsi1Adj() {
666 // size_t argmax_var=0;
667 // size_t argmax_var_state=0;
669 // for( size_t i = 0; i < _fg->nrVars(); i++ ) {
670 // pair<size_t,Real> argmax_state = adj_psi_V(i).argmax();
671 // if(i==0 || argmax_state.second>max_var) {
673 // max_var = argmax_state.second;
674 // argmax_var_state = argmax_state.first;
677 // DAI_ASSERT(/*0 <= argmax_var_state &&*/
678 // argmax_var_state < _fg->var(argmax_var).states());
679 // return tuple<size_t,size_t,Real>(argmax_var,argmax_var_state,max_var);
683 void BBP::getArgmaxMsgM( size_t &out_i
, size_t &out__I
, Real
&mag
) {
685 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
686 bforeach( const Neighbor
&I
, _fg
->nbV(i
) ) {
687 Real thisMag
= _adj_m
[i
][I
.iter
].sumAbs();
688 if( !found
|| mag
< thisMag
) {
699 Real
BBP::getMaxMsgM() {
702 getArgmaxMsgM( dummy
, dummy
, mag
);
707 Real
BBP::getTotalMsgM() {
709 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
710 bforeach( const Neighbor
&I
, _fg
->nbV(i
) )
711 mag
+= _adj_m
[i
][I
.iter
].sumAbs();
716 Real
BBP::getTotalNewMsgM() {
718 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
719 bforeach( const Neighbor
&I
, _fg
->nbV(i
) )
720 mag
+= _new_adj_m
[i
][I
.iter
].sumAbs();
725 Real
BBP::getTotalMsgN() {
727 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
728 bforeach( const Neighbor
&I
, _fg
->nbV(i
) )
729 mag
+= _adj_n
[i
][I
.iter
].sumAbs();
734 std::vector
<Prob
> BBP::getZeroAdjF( const FactorGraph
&fg
) {
736 adj_2
.reserve( fg
.nrFactors() );
737 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ )
738 adj_2
.push_back( Prob( fg
.factor(I
).nrStates(), 0.0 ) );
743 std::vector
<Prob
> BBP::getZeroAdjV( const FactorGraph
&fg
) {
745 adj_1
.reserve( fg
.nrVars() );
746 for( size_t i
= 0; i
< fg
.nrVars(); i
++ )
747 adj_1
.push_back( Prob( fg
.var(i
).states(), 0.0 ) );
752 void BBP::initCostFnAdj( const BBPCostFunction
&cfn
, const vector
<size_t> *stateP
) {
753 const FactorGraph
&fg
= _ia
->fg();
755 switch( (size_t)cfn
) {
756 case BBPCostFunction::CFN_BETHE_ENT
: {
759 vector
<Prob
> psi1_adj
;
760 vector
<Prob
> psi2_adj
;
761 b1_adj
.reserve( fg
.nrVars() );
762 psi1_adj
.reserve( fg
.nrVars() );
763 b2_adj
.reserve( fg
.nrFactors() );
764 psi2_adj
.reserve( fg
.nrFactors() );
765 for( size_t i
= 0; i
< fg
.nrVars(); i
++ ) {
766 size_t dim
= fg
.var(i
).states();
767 int c
= fg
.nbV(i
).size();
769 for( size_t xi
= 0; xi
< dim
; xi
++ )
770 p
.set( xi
, (1 - c
) * (1 + log( _ia
->beliefV(i
)[xi
] )) );
771 b1_adj
.push_back( p
);
773 for( size_t xi
= 0; xi
< dim
; xi
++ )
774 p
.set( xi
, -_ia
->beliefV(i
)[xi
] );
775 psi1_adj
.push_back( p
);
777 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
778 size_t dim
= fg
.factor(I
).nrStates();
780 for( size_t xI
= 0; xI
< dim
; xI
++ )
781 p
.set( xI
, 1 + log( _ia
->beliefF(I
)[xI
] / fg
.factor(I
).p()[xI
] ) );
782 b2_adj
.push_back( p
);
784 for( size_t xI
= 0; xI
< dim
; xI
++ )
785 p
.set( xI
, -_ia
->beliefF(I
)[xI
] / fg
.factor(I
).p()[xI
] );
786 psi2_adj
.push_back( p
);
788 init( b1_adj
, b2_adj
, psi1_adj
, psi2_adj
);
790 } case BBPCostFunction::CFN_FACTOR_ENT
: {
792 b2_adj
.reserve( fg
.nrFactors() );
793 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
794 size_t dim
= fg
.factor(I
).nrStates();
796 for( size_t xI
= 0; xI
< dim
; xI
++ ) {
797 Real bIxI
= _ia
->beliefF(I
)[xI
];
799 p
.set( xI
, -1.0e10
);
801 p
.set( xI
, 1 + log( bIxI
) );
807 } case BBPCostFunction::CFN_VAR_ENT
: {
809 b1_adj
.reserve( fg
.nrVars() );
810 for( size_t i
= 0; i
< fg
.nrVars(); i
++ ) {
811 size_t dim
= fg
.var(i
).states();
813 for( size_t xi
= 0; xi
< fg
.var(i
).states(); xi
++ ) {
814 Real bixi
= _ia
->beliefV(i
)[xi
];
816 p
.set( xi
, -1.0e10
);
818 p
.set( xi
, 1 + log( bixi
) );
820 b1_adj
.push_back( p
);
824 } case BBPCostFunction::CFN_GIBBS_B
:
825 case BBPCostFunction::CFN_GIBBS_B2
:
826 case BBPCostFunction::CFN_GIBBS_EXP
: {
827 // cost functions that use Gibbs sample, summing over variable marginals
828 vector
<size_t> state
;
830 state
= getGibbsState( _ia
->fg(), 2*_ia
->Iterations() );
833 DAI_ASSERT( state
.size() == fg
.nrVars() );
836 b1_adj
.reserve(fg
.nrVars());
837 for( size_t i
= 0; i
< state
.size(); i
++ ) {
838 size_t n
= fg
.var(i
).states();
839 Prob
delta( n
, 0.0 );
840 DAI_ASSERT(/*0<=state[i] &&*/ state
[i
] < n
);
841 Real b
= _ia
->beliefV(i
)[state
[i
]];
842 switch( (size_t)cfn
) {
843 case BBPCostFunction::CFN_GIBBS_B
:
844 delta
.set( state
[i
], 1.0 );
846 case BBPCostFunction::CFN_GIBBS_B2
:
847 delta
.set( state
[i
], b
);
849 case BBPCostFunction::CFN_GIBBS_EXP
:
850 delta
.set( state
[i
], exp(b
) );
853 DAI_THROW(UNKNOWN_ENUM_VALUE
);
855 b1_adj
.push_back( delta
);
859 } case BBPCostFunction::CFN_GIBBS_B_FACTOR
:
860 case BBPCostFunction::CFN_GIBBS_B2_FACTOR
:
861 case BBPCostFunction::CFN_GIBBS_EXP_FACTOR
: {
862 // cost functions that use Gibbs sample, summing over factor marginals
863 vector
<size_t> state
;
865 state
= getGibbsState( _ia
->fg(), 2*_ia
->Iterations() );
868 DAI_ASSERT( state
.size() == fg
.nrVars() );
871 b2_adj
.reserve( fg
.nrVars() );
872 for( size_t I
= 0; I
< fg
.nrFactors(); I
++ ) {
873 size_t n
= fg
.factor(I
).nrStates();
874 Prob
delta( n
, 0.0 );
876 size_t x_I
= getFactorEntryForState( fg
, I
, state
);
877 DAI_ASSERT(/*0<=x_I &&*/ x_I
< n
);
879 Real b
= _ia
->beliefF(I
)[x_I
];
880 switch( (size_t)cfn
) {
881 case BBPCostFunction::CFN_GIBBS_B_FACTOR
:
882 delta
.set( x_I
, 1.0 );
884 case BBPCostFunction::CFN_GIBBS_B2_FACTOR
:
887 case BBPCostFunction::CFN_GIBBS_EXP_FACTOR
:
888 delta
.set( x_I
, exp( b
) );
891 DAI_THROW(UNKNOWN_ENUM_VALUE
);
893 b2_adj
.push_back( delta
);
898 DAI_THROW(UNKNOWN_ENUM_VALUE
);
904 typedef BBP::Properties::UpdateType UT
;
905 Real tol
= props
.tol
;
906 UT
&updates
= props
.updates
;
909 switch( (size_t)updates
) {
915 getArgmaxMsgM( i
, _I
, mag
);
916 sendSeqMsgM( i
, _I
);
917 } while( mag
> tol
&& _iters
< props
.maxiter
);
919 if( _iters
>= props
.maxiter
)
920 if( props
.verbose
>= 1 )
921 cerr
<< "Warning: BBP didn't converge in " << _iters
<< " iterations (greatest message magnitude = " << mag
<< ")" << endl
;
923 } case UT::SEQ_FIX
: {
927 mag
= getTotalMsgM();
931 for( size_t i
= 0; i
< _fg
->nrVars(); i
++ )
932 bforeach( const Neighbor
&I
, _fg
->nbV(i
) )
933 sendSeqMsgM( i
, I
.iter
);
934 /* for( size_t i = 0; i < _fg->nrVars(); i++ )
935 bforeach( const Neighbor &I, _fg->nbV(i) )
936 updateSeqMsgM( i, I.iter );*/
937 } while( mag
> tol
&& _iters
< props
.maxiter
);
939 if( _iters
>= props
.maxiter
)
940 if( props
.verbose
>= 1 )
941 cerr
<< "Warning: BBP didn't converge in " << _iters
<< " iterations (greatest message magnitude = " << mag
<< ")" << endl
;
943 } case UT::SEQ_BP_REV
:
944 case UT::SEQ_BP_FWD
: {
945 const BP
*bp
= static_cast<const BP
*>(_ia
);
946 vector
<pair
<size_t, size_t> > sentMessages
= bp
->getSentMessages();
947 size_t totalMessages
= sentMessages
.size();
948 if( totalMessages
== 0 )
949 DAI_THROWE(INTERNAL_ERROR
, "Asked for updates=" + std::string(updates
) + " but no BP messages; did you forget to set recordSentMessages?");
950 if( updates
==UT::SEQ_BP_FWD
)
951 reverse( sentMessages
.begin(), sentMessages
.end() );
952 // DAI_PV(sentMessages.size());
954 // DAI_PV(props.maxiter);
955 while( sentMessages
.size() > 0 && _iters
< props
.maxiter
) {
956 // DAI_PV(sentMessages.size());
959 pair
<size_t, size_t> e
= sentMessages
.back();
960 sentMessages
.pop_back();
961 size_t i
= e
.first
, _I
= e
.second
;
962 sendSeqMsgM( i
, _I
);
964 if( _iters
>= props
.maxiter
)
965 if( props
.verbose
>= 1 )
966 cerr
<< "Warning: BBP updates limited to " << props
.maxiter
<< " iterations, but using UpdateType " << updates
<< " with " << totalMessages
<< " messages" << endl
;
972 } while( (_iters
< 2 || getUnMsgMag() > tol
) && _iters
< props
.maxiter
);
973 if( _iters
== props
.maxiter
) {
975 getMsgMags( s
, new_s
);
976 if( props
.verbose
>= 1 )
977 cerr
<< "Warning: BBP didn't converge in " << _iters
<< " iterations (unnorm message magnitude = " << getUnMsgMag() << ", norm message mags = " << s
<< " -> " << new_s
<< ")" << endl
;
982 if( props
.verbose
>= 3 )
983 cerr
<< "BBP::run() took " << toc()-tic
<< " seconds " << Iterations() << " iterations" << endl
;
987 Real
numericBBPTest( const InfAlg
&bp
, const std::vector
<size_t> *state
, const PropertySet
&bbp_props
, const BBPCostFunction
&cfn
, Real h
) {
988 BBP
bbp( &bp
, bbp_props
);
989 // calculate the value of the unperturbed cost function
990 Real cf0
= cfn
.evaluate( bp
, state
);
991 // run BBP to estimate adjoints
992 bbp
.initCostFnAdj( cfn
, state
);
996 const FactorGraph
& fg
= bp
.fg();
999 // verify bbp.adj_psi_V
1001 // for each variable i
1002 for( size_t i
= 0; i
< fg
.nrVars(); i
++ ) {
1003 vector
<Real
> adj_est
;
1004 // for each value xi
1005 for( size_t xi
= 0; xi
< fg
.var(i
).states(); xi
++ ) {
1006 // Clone 'bp' (which may be any InfAlg)
1007 InfAlg
*bp_prb
= bp
.clone();
1010 size_t n
= bp_prb
->fg().var(i
).states();
1011 Prob
psi_1_prb( n
, 1.0 );
1012 psi_1_prb
.set( xi
, psi_1_prb
[xi
] + h
);
1013 // psi_1_prb.normalize();
1014 size_t I
= bp_prb
->fg().nbV(i
)[0]; // use first factor in list of neighbors of i
1015 Factor tmp
= bp_prb
->fg().factor(I
) * Factor( bp_prb
->fg().var(i
), psi_1_prb
);
1016 bp_prb
->fg().setFactor( I
, tmp
);
1018 // call 'init' on the perturbed variables
1019 bp_prb
->init( bp_prb
->fg().var(i
) );
1021 // run copy to convergence
1024 // calculate new value of cost function
1025 Real cf_prb
= cfn
.evaluate( *bp_prb
, state
);
1027 // use to estimate adjoint for i
1028 adj_est
.push_back( (cf_prb
- cf0
) / h
);
1030 // free cloned InfAlg
1033 Prob
p_adj_est( adj_est
);
1034 // compare this numerical estimate to the BBP estimate; sum the distances
1036 << ", p_adj_est: " << p_adj_est
1037 << ", bbp.adj_psi_V(i): " << bbp
.adj_psi_V(i
) << endl
;
1038 d
+= dist( p_adj_est
, bbp
.adj_psi_V(i
), DISTL1
);
1042 // verify bbp.adj_n and bbp.adj_m
1044 // We actually want to check the responsiveness of objective
1045 // function to changes in the final messages. But at the end of a
1046 // BBP run, the message adjoints are for the initial messages (and
1047 // they should be close to zero, see paper). So this resets the
1048 // BBP adjoints to the refer to the desired final messages
1049 bbp.RegenerateMessageAdjoints();
1051 // for each variable i
1052 for(size_t i=0; i<bp_dual.nrVars(); i++) {
1053 // for each factor I ~ i
1054 bforeach(size_t I, bp_dual.nbV(i)) {
1055 vector<Real> adj_n_est;
1056 // for each value xi
1057 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1058 BP_dual bp_dual_prb(bp_dual);
1059 // make h-sized change to newMsgN
1060 bp_dual_prb.newMsgN(i,I)[xi] += h;
1061 // recalculate beliefs
1062 bp_dual_prb.CalcBeliefs();
1063 // get cost function value
1064 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1065 // add it to list of adjoints
1066 adj_n_est.push_back((cf_prb-cf0)/h);
1069 vector<Real> adj_m_est;
1070 // for each value xi
1071 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1072 BP_dual bp_dual_prb(bp_dual);
1073 // make h-sized change to newMsgM
1074 bp_dual_prb.newMsgM(I,i)[xi] += h;
1075 // recalculate beliefs
1076 bp_dual_prb.CalcBeliefs();
1077 // get cost function value
1078 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1079 // add it to list of adjoints
1080 adj_m_est.push_back((cf_prb-cf0)/h);
1083 Prob p_adj_n_est( adj_n_est );
1084 // compare this numerical estimate to the BBP estimate; sum the distances
1085 cerr << "i: " << i << ", I: " << I
1086 << ", adj_n_est: " << p_adj_n_est
1087 << ", bbp.adj_n(i,I): " << bbp.adj_n(i,I) << endl;
1088 d += dist(p_adj_n_est, bbp.adj_n(i,I), DISTL1);
1090 Prob p_adj_m_est( adj_m_est );
1091 // compare this numerical estimate to the BBP estimate; sum the distances
1092 cerr << "i: " << i << ", I: " << I
1093 << ", adj_m_est: " << p_adj_m_est
1094 << ", bbp.adj_m(I,i): " << bbp.adj_m(I,i) << endl;
1095 d += dist(p_adj_m_est, bbp.adj_m(I,i), DISTL1);
1101 // verify bbp.adj_b_V
1102 for(size_t i=0; i<bp_dual.nrVars(); i++) {
1103 vector<Real> adj_b_V_est;
1104 // for each value xi
1105 for(size_t xi=0; xi<bp_dual.var(i).states(); xi++) {
1106 BP_dual bp_dual_prb(bp_dual);
1108 // make h-sized change to b_1(i)[x_i]
1109 bp_dual_prb._beliefs.b1[i][xi] += h;
1111 // get cost function value
1112 Real cf_prb = getCostFn(bp_dual_prb, cfn, &state);
1114 // add it to list of adjoints
1115 adj_b_V_est.push_back((cf_prb-cf0)/h);
1117 Prob p_adj_b_V_est( adj_b_V_est );
1118 // compare this numerical estimate to the BBP estimate; sum the distances
1120 << ", adj_b_V_est: " << p_adj_b_V_est
1121 << ", bbp.adj_b_V(i): " << bbp.adj_b_V(i) << endl;
1122 d += dist(p_adj_b_V_est, bbp.adj_b_V(i), DISTL1);
1127 // return total of distances
1132 } // end of namespace dai
1135 /* {{{ GENERATED CODE: DO NOT EDIT. Created by
1136 ./scripts/regenerate-properties include/dai/bbp.h src/bbp.cpp
1140 void BBP::Properties::set(const PropertySet
&opts
)
1142 const std::set
<PropertyKey
> &keys
= opts
.keys();
1143 std::string errormsg
;
1144 for( std::set
<PropertyKey
>::const_iterator i
= keys
.begin(); i
!= keys
.end(); i
++ ) {
1145 if( *i
== "verbose" ) continue;
1146 if( *i
== "maxiter" ) continue;
1147 if( *i
== "tol" ) continue;
1148 if( *i
== "damping" ) continue;
1149 if( *i
== "updates" ) continue;
1150 errormsg
= errormsg
+ "BBP: Unknown property " + *i
+ "\n";
1152 if( !errormsg
.empty() )
1153 DAI_THROWE(UNKNOWN_PROPERTY
, errormsg
);
1154 if( !opts
.hasKey("maxiter") )
1155 errormsg
= errormsg
+ "BBP: Missing property \"maxiter\" for method \"BBP\"\n";
1156 if( !opts
.hasKey("tol") )
1157 errormsg
= errormsg
+ "BBP: Missing property \"tol\" for method \"BBP\"\n";
1158 if( !opts
.hasKey("damping") )
1159 errormsg
= errormsg
+ "BBP: Missing property \"damping\" for method \"BBP\"\n";
1160 if( !opts
.hasKey("updates") )
1161 errormsg
= errormsg
+ "BBP: Missing property \"updates\" for method \"BBP\"\n";
1162 if( !errormsg
.empty() )
1163 DAI_THROWE(NOT_ALL_PROPERTIES_SPECIFIED
,errormsg
);
1164 if( opts
.hasKey("verbose") ) {
1165 verbose
= opts
.getStringAs
<size_t>("verbose");
1169 maxiter
= opts
.getStringAs
<size_t>("maxiter");
1170 tol
= opts
.getStringAs
<Real
>("tol");
1171 damping
= opts
.getStringAs
<Real
>("damping");
1172 updates
= opts
.getStringAs
<UpdateType
>("updates");
1174 PropertySet
BBP::Properties::get() const {
1176 opts
.set("verbose", verbose
);
1177 opts
.set("maxiter", maxiter
);
1178 opts
.set("tol", tol
);
1179 opts
.set("damping", damping
);
1180 opts
.set("updates", updates
);
1183 string
BBP::Properties::toString() const {
1184 stringstream
s(stringstream::out
);
1186 s
<< "verbose=" << verbose
<< ",";
1187 s
<< "maxiter=" << maxiter
<< ",";
1188 s
<< "tol=" << tol
<< ",";
1189 s
<< "damping=" << damping
<< ",";
1190 s
<< "updates=" << updates
;
1194 } // end of namespace dai
1195 /* }}} END OF GENERATED CODE */