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>
17 #include <dai/bp_dual.h>
19 #include <dai/bipgraph.h>
28 void BP_dual::init() {
36 void BP_dual::regenerateMessages() {
37 size_t nv
= fg().nrVars();
42 for( size_t i
= 0; i
< nv
; i
++ ) {
43 size_t nvf
= fg().nbV(i
).size();
44 _msgs
.Zn
[i
].resize(nvf
, 1.0);
45 _msgs
.Zm
[i
].resize(nvf
, 1.0);
46 size_t states
= fg().var(i
).states();
47 _msgs
.n
[i
].resize(nvf
, Prob(states
));
48 _msgs
.m
[i
].resize(nvf
, Prob(states
));
53 void BP_dual::regenerateBeliefs() {
55 _beliefs
.b1
.reserve(fg().nrVars());
56 _beliefs
.Zb1
.resize(fg().nrVars(), 1.0);
58 _beliefs
.b2
.reserve(fg().nrFactors());
59 _beliefs
.Zb2
.resize(fg().nrFactors(), 1.0);
61 for( size_t i
= 0; i
< fg().nrVars(); i
++ )
62 _beliefs
.b1
.push_back( Prob( fg().var(i
).states() ) );
63 for( size_t I
= 0; I
< fg().nrFactors(); I
++ )
64 _beliefs
.b2
.push_back( Prob( fg().factor(I
).nrStates() ) );
68 void BP_dual::calcMessages() {
69 // calculate 'n' messages from "factor marginal / factor"
70 for( size_t I
= 0; I
< fg().nrFactors(); I
++ ) {
71 Factor f
= _ia
->beliefF(I
) / fg().factor(I
);
72 bforeach( const Neighbor
&i
, fg().nbF(I
) )
73 msgN(i
, i
.dual
) = f
.marginal( fg().var(i
) ).p();
75 // calculate 'm' messages and normalizers from 'n' messages
76 for( size_t i
= 0; i
< fg().nrVars(); i
++ )
77 bforeach( const Neighbor
&I
, fg().nbV(i
) )
78 calcNewM( i
, I
.iter
);
79 // recalculate 'n' messages and normalizers from 'm' messages
80 for( size_t i
= 0; i
< fg().nrVars(); i
++ )
81 bforeach( const Neighbor
&I
, fg().nbV(i
) )
86 void BP_dual::calcNewM( size_t i
, size_t _I
) {
87 // calculate updated message I->i
88 const Neighbor
&I
= fg().nbV(i
)[_I
];
89 Prob
prod( fg().factor(I
).p() );
90 bforeach( const Neighbor
&j
, fg().nbF(I
) )
91 if( j
!= i
) { // for all j in I \ i
92 Prob
&n
= msgN(j
,j
.dual
);
93 IndexFor
ind( fg().var(j
), fg().factor(I
).vars() );
94 for( size_t x
= 0; ind
.valid(); x
++, ++ind
)
95 prod
.set( x
, prod
[x
] * n
[ind
] );
98 Prob
marg( fg().var(i
).states(), 0.0 );
99 // ind is the precalculated Index(i,I) i.e. to x_I == k corresponds x_i == ind[k]
100 IndexFor
ind( fg().var(i
), fg().factor(I
).vars() );
101 for( size_t x
= 0; ind
.valid(); x
++, ++ind
)
102 marg
.set( ind
, marg
[ind
] + prod
[x
] );
104 _msgs
.Zm
[i
][_I
] = marg
.normalize();
105 _msgs
.m
[i
][_I
] = marg
;
109 void BP_dual::calcNewN( size_t i
, size_t _I
) {
110 // calculate updated message i->I
111 const Neighbor
&I
= fg().nbV(i
)[_I
];
112 Prob
prod( fg().var(i
).states(), 1.0 );
113 bforeach( const Neighbor
&J
, fg().nbV(i
) )
114 if( J
.node
!= I
.node
) // for all J in i \ I
115 prod
*= msgM(i
,J
.iter
);
116 _msgs
.Zn
[i
][_I
] = prod
.normalize();
117 _msgs
.n
[i
][_I
] = prod
;
121 void BP_dual::calcBeliefs() {
122 for( size_t i
= 0; i
< fg().nrVars(); i
++ )
123 calcBeliefV(i
); // calculate b_i
124 for( size_t I
= 0; I
< fg().nrFactors(); I
++ )
125 calcBeliefF(I
); // calculate b_I
129 void BP_dual::calcBeliefV( size_t i
) {
130 Prob
prod( fg().var(i
).states(), 1.0 );
131 bforeach( const Neighbor
&I
, fg().nbV(i
) )
132 prod
*= msgM(i
,I
.iter
);
133 _beliefs
.Zb1
[i
] = prod
.normalize();
134 _beliefs
.b1
[i
] = prod
;
138 void BP_dual::calcBeliefF( size_t I
) {
139 Prob
prod( fg().factor(I
).p() );
140 bforeach( const Neighbor
&j
, fg().nbF(I
) ) {
141 IndexFor
ind( fg().var(j
), fg().factor(I
).vars() );
142 Prob
n( msgN(j
,j
.dual
) );
143 for( size_t x
= 0; ind
.valid(); x
++, ++ind
)
144 prod
.set( x
, prod
[x
] * n
[ind
] );
146 _beliefs
.Zb2
[I
] = prod
.normalize();
147 _beliefs
.b2
[I
] = prod
;
151 } // end of namespace dai