Multiple changes: changes in build system, one workaround and one bug fix
[libdai.git] / src / bp_dual.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2011, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <dai/dai_config.h>
10 #ifdef DAI_WITH_CBP
11
12
13 #include <iostream>
14 #include <sstream>
15 #include <algorithm>
16
17 #include <dai/bp_dual.h>
18 #include <dai/util.h>
19 #include <dai/bipgraph.h>
20
21
22 namespace dai {
23
24
25 using namespace std;
26
27
28 void BP_dual::init() {
29 regenerateMessages();
30 regenerateBeliefs();
31 calcMessages();
32 calcBeliefs();
33 }
34
35
36 void BP_dual::regenerateMessages() {
37 size_t nv = fg().nrVars();
38 _msgs.Zn.resize(nv);
39 _msgs.Zm.resize(nv);
40 _msgs.m.resize(nv);
41 _msgs.n.resize(nv);
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));
49 }
50 }
51
52
53 void BP_dual::regenerateBeliefs() {
54 _beliefs.b1.clear();
55 _beliefs.b1.reserve(fg().nrVars());
56 _beliefs.Zb1.resize(fg().nrVars(), 1.0);
57 _beliefs.b2.clear();
58 _beliefs.b2.reserve(fg().nrFactors());
59 _beliefs.Zb2.resize(fg().nrFactors(), 1.0);
60
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() ) );
65 }
66
67
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();
74 }
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) )
82 calcNewN(i, I.iter);
83 }
84
85
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] );
96 }
97 // Marginalize onto i
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] );
103
104 _msgs.Zm[i][_I] = marg.normalize();
105 _msgs.m[i][_I] = marg;
106 }
107
108
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;
118 }
119
120
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
126 }
127
128
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;
135 }
136
137
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] );
145 }
146 _beliefs.Zb2[I] = prod.normalize();
147 _beliefs.b2[I] = prod;
148 }
149
150
151 } // end of namespace dai
152
153
154 #endif