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>
16 #define DAI_FBP_FAST 1
25 // This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
26 Real
FBP::logZ() const {
28 for( size_t I
= 0; I
< nrFactors(); I
++ ) {
29 sum
+= (beliefF(I
) * factor(I
).log(true)).sum(); // FBP
30 sum
+= Weight(I
) * beliefF(I
).entropy(); // FBP
32 for( size_t i
= 0; i
< nrVars(); ++i
) {
34 bforeach( const Neighbor
&I
, nbV(i
) )
37 sum
+= (1.0 - c_i
) * beliefV(i
).entropy(); // FBP
43 // This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
44 Prob
FBP::calcIncomingMessageProduct( size_t I
, bool without_i
, size_t i
) const {
45 Real c_I
= Weight(I
); // FBP: c_I
47 Factor
Fprod( factor(I
) );
48 Prob
&prod
= Fprod
.p();
50 if( props
.logdomain
) {
54 prod
^= (1.0 / c_I
); // FBP
56 // Calculate product of incoming messages and factor I
57 bforeach( const Neighbor
&j
, nbF(I
) )
58 if( !(without_i
&& (j
== i
)) ) {
59 // prod_j will be the product of messages coming into j
60 // FBP: corresponds to messages n_jI
61 Prob
prod_j( var(j
).states(), props
.logdomain
? 0.0 : 1.0 );
62 bforeach( const Neighbor
&J
, nbV(j
) )
63 if( J
!= I
) { // for all J in nb(j) \ I
65 prod_j
+= message( j
, J
.iter
);
67 prod_j
*= message( j
, J
.iter
);
68 } else if( c_I
!= 1.0 ) {
69 // FBP: multiply by m_Ij^(1-1/c_I)
71 prod_j
+= newMessage( j
, J
.iter
) * (1.0 - 1.0 / c_I
);
73 prod_j
*= newMessage( j
, J
.iter
) ^ (1.0 - 1.0 / c_I
);
76 // multiply prod with prod_j
78 // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
80 Fprod
+= Factor( var(j
), prod_j
);
82 Fprod
*= Factor( var(j
), prod_j
);
86 // ind is the precalculated IndexFor(j,I) i.e. to x_I == k corresponds x_j == ind[k]
87 const ind_t
&ind
= index(j
, _I
);
89 for( size_t r
= 0; r
< prod
.size(); ++r
)
91 prod
.set( r
, prod
[r
] + prod_j
[ind
[r
]] );
93 prod
.set( r
, prod
[r
] * prod_j
[ind
[r
]] );
100 // This code has been copied from bp.cpp, except where comments indicate FBP-specific behaviour
101 void FBP::calcNewMessage( size_t i
, size_t _I
) {
102 // calculate updated message I->i
103 size_t I
= nbV(i
,_I
);
105 Real c_I
= Weight(I
); // FBP: c_I
107 Factor
Fprod( factor(I
) );
108 Prob
&prod
= Fprod
.p();
109 prod
= calcIncomingMessageProduct( I
, true, i
);
111 if( props
.logdomain
) {
116 // Marginalize onto i
118 if( !DAI_FBP_FAST
) {
119 // UNOPTIMIZED (SIMPLE TO READ, BUT SLOW) VERSION
120 if( props
.inference
== Properties::InfType::SUMPROD
)
121 marg
= Fprod
.marginal( var(i
) ).p();
123 marg
= Fprod
.maxMarginal( var(i
) ).p();
126 marg
= Prob( var(i
).states(), 0.0 );
127 // ind is the precalculated IndexFor(i,I) i.e. to x_I == k corresponds x_i == ind[k]
128 const ind_t ind
= index(i
,_I
);
129 if( props
.inference
== Properties::InfType::SUMPROD
)
130 for( size_t r
= 0; r
< prod
.size(); ++r
)
131 marg
.set( ind
[r
], marg
[ind
[r
]] + prod
[r
] );
133 for( size_t r
= 0; r
< prod
.size(); ++r
)
134 if( prod
[r
] > marg
[ind
[r
]] )
135 marg
.set( ind
[r
], prod
[r
] );
143 if( props
.logdomain
)
144 newMessage(i
,_I
) = marg
.log();
146 newMessage(i
,_I
) = marg
;
148 // Update the residual if necessary
149 if( props
.updates
== Properties::UpdateType::SEQMAX
)
150 updateResidual( i
, _I
, dist( newMessage( i
, _I
), message( i
, _I
), DISTLINF
) );
154 void FBP::construct() {
156 _weight
.resize( nrFactors(), 1.0 );
160 } // end of namespace dai