Partial adoption of contributions by Giuseppe:
[libdai.git] / src / regiongraph.cpp
1 /* Copyright (C) 2006-2008 Joris Mooij [j dot mooij at science dot ru dot nl]
2 Radboud University Nijmegen, The Netherlands
3
4 This file is part of libDAI.
5
6 libDAI is free software; you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation; either version 2 of the License, or
9 (at your option) any later version.
10
11 libDAI is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with libDAI; if not, write to the Free Software
18 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
19 */
20
21
22 #include <algorithm>
23 #include <cmath>
24 #include <dai/regiongraph.h>
25 #include <dai/factorgraph.h>
26 #include <dai/clustergraph.h>
27
28
29 namespace dai {
30
31
32 using namespace std;
33
34
35 RegionGraph::RegionGraph(const FactorGraph & fg, const vector<Region> & ors, const vector<Region> & irs, const vector<R_edge_t> & edges) : FactorGraph(fg), BipRegGraph(), _fac2OR() {
36 // Copy outer regions (giving them counting number 1.0)
37 ORs().reserve( ors.size() );
38 for( vector<Region>::const_iterator alpha = ors.begin(); alpha != ors.end(); alpha++ )
39 ORs().push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
40
41 // Incorporate factors into outer regions
42 /* if( !_hasNegatives ) {
43 // For each factor, find the outer regions that subsume that factor.
44 // Then, "divide" the factor over its subsuming outer regions.
45 for( size_t I = 0; I < nrFactors(); I++ ) {
46 vector<size_t> subsuming_alphas;
47
48 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
49 if( OR(alpha).vars() >> factor(I).vars() )
50 subsuming_alphas.push_back( alpha );
51 }
52
53 assert( subsuming_alphas.size() >= 1 );
54
55 for( vector<size_t>::const_iterator alpha = subsuming_alphas.begin(); alpha != subsuming_alphas.end(); alpha++ )
56 OR(*alpha) *= (factor(I) ^ (1.0 / subsuming_alphas.size()));
57 }
58 } else {
59 cout << "Careful composition of outer regions" << endl;
60 */
61 // For each factor, find an outer regions that subsumes that factor.
62 // Then, multiply the outer region with that factor.
63 for( size_t I = 0; I < nrFactors(); I++ ) {
64 size_t alpha;
65 for( alpha = 0; alpha < nr_ORs(); alpha++ )
66 if( OR(alpha).vars() >> factor(I).vars() ) {
67 _fac2OR[I] = alpha;
68 // OR(alpha) *= factor(I);
69 break;
70 }
71 assert( alpha != nr_ORs() );
72 }
73 RecomputeORs();
74 // }
75
76 // Copy inner regions
77 IRs().reserve( irs.size() );
78 for( vector<Region>::const_iterator beta = irs.begin(); beta != irs.end(); beta++ )
79 IRs().push_back( *beta );
80
81 // Copy edges
82 Redges().reserve( edges.size() );
83 for( vector<R_edge_t>::const_iterator e = edges.begin(); e != edges.end(); e++ )
84 Redges().push_back( *e );
85
86 // Regenerate BipartiteGraph internals
87 BipRegGraph::Regenerate();
88
89 // Check counting numbers
90 Check_Counting_Numbers();
91 }
92
93
94 // CVM style
95 RegionGraph::RegionGraph(const FactorGraph & fg, const vector<VarSet> & cl) : FactorGraph(fg), BipRegGraph() {
96 // Retain only maximal clusters
97 ClusterGraph cg( cl );
98 cg.eraseNonMaximal();
99
100 // Create outer regions, giving them counting number 1.0
101 ORs().reserve( cg.size() );
102 for( ClusterGraph::const_iterator alpha = cg.begin(); alpha != cg.end(); alpha++ )
103 ORs().push_back( FRegion(Factor(*alpha, 1.0), 1.0) );
104
105 // Incorporate factors into outer regions
106 /* if( !_hasNegatives ) {
107 // For each factor, find the outer regions that subsume that factor.
108 // Then, "divide" the factor over its subsuming outer regions.
109 for( size_t I = 0; I < nrFactors(); I++ ) {
110 vector<size_t> subsuming_alphas;
111
112 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
113 if( OR(alpha).vars() >> factor(I).vars() )
114 subsuming_alphas.push_back( alpha );
115 }
116
117 assert( subsuming_alphas.size() >= 1 );
118
119 for( vector<size_t>::const_iterator alpha = subsuming_alphas.begin(); alpha != subsuming_alphas.end(); alpha++ )
120 OR(*alpha) *= factor(I) ^ (1.0 / subsuming_alphas.size());
121 }
122 } else {
123 cout << "Careful composition of outer regions" << endl;
124 */
125 // For each factor, find an outer regions that subsumes that factor.
126 // Then, multiply the outer region with that factor.
127 for( size_t I = 0; I < nrFactors(); I++ ) {
128 size_t alpha;
129 for( alpha = 0; alpha < nr_ORs(); alpha++ )
130 if( OR(alpha).vars() >> factor(I).vars() ) {
131 _fac2OR[I] = alpha;
132 // OR(alpha) *= factor(I);
133 break;
134 }
135 assert( alpha != nr_ORs() );
136 }
137 RecomputeORs();
138 // }
139
140 // Create inner regions - first pass
141 set<VarSet> betas;
142 for( ClusterGraph::const_iterator alpha = cg.begin(); alpha != cg.end(); alpha++ )
143 for( ClusterGraph::const_iterator alpha2 = alpha; (++alpha2) != cg.end(); ) {
144 VarSet intersect = (*alpha) & (*alpha2);
145 if( intersect.size() > 0 )
146 betas.insert( intersect );
147 }
148
149 // Create inner regions - subsequent passes
150 set<VarSet> new_betas;
151 do {
152 new_betas.clear();
153 for( set<VarSet>::const_iterator gamma = betas.begin(); gamma != betas.end(); gamma++ )
154 for( set<VarSet>::const_iterator gamma2 = gamma; (++gamma2) != betas.end(); ) {
155 VarSet intersect = (*gamma) & (*gamma2);
156 if( (intersect.size() > 0) && (betas.count(intersect) == 0) )
157 new_betas.insert( intersect );
158 }
159 betas.insert(new_betas.begin(), new_betas.end());
160 } while( new_betas.size() );
161
162 // Create inner regions - store them in the bipartite graph
163 IRs().reserve( betas.size() );
164 for( set<VarSet>::const_iterator beta = betas.begin(); beta != betas.end(); beta++ )
165 IRs().push_back( Region(*beta,NAN) );
166
167 // Create edges
168 for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
169 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ ) {
170 if( OR(alpha).vars() >> IR(beta) )
171 Redges().push_back(R_edge_t(alpha,beta));
172 }
173 }
174
175 // Regenerate BipartiteGraph internals
176 BipRegGraph::Regenerate();
177
178 // Calculate counting numbers
179 Calc_Counting_Numbers();
180
181 // Check counting numbers
182 Check_Counting_Numbers();
183 }
184
185
186 void RegionGraph::Calc_Counting_Numbers() {
187 // Calculates counting numbers of inner regions based upon counting numbers of outer regions
188
189 vector<vector<size_t> > ancestors(nr_IRs());
190 for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
191 IR(beta).c() = NAN;
192 for( size_t beta2 = 0; beta2 < nr_IRs(); beta2++ )
193 if( (beta2 != beta) && IR(beta2) >> IR(beta) )
194 ancestors[beta].push_back(beta2);
195 }
196
197 bool new_counting;
198 do {
199 new_counting = false;
200 for( size_t beta = 0; beta < nr_IRs(); beta++ ) {
201 if( isnan( IR(beta).c() ) ) {
202 bool has_nan_ancestor = false;
203 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_nan_ancestor; beta2++ )
204 if( isnan( IR(*beta2).c() ) )
205 has_nan_ancestor = true;
206 if( !has_nan_ancestor ) {
207 double c = 1.0;
208 for( R_nb_cit alpha = nbIR(beta).begin(); alpha != nbIR(beta).end(); alpha++ )
209 c -= OR(*alpha).c();
210 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ )
211 c -= IR(*beta2).c();
212 IR(beta).c() = c;
213 new_counting = true;
214 }
215 }
216 }
217 } while( new_counting );
218 }
219
220
221 bool RegionGraph::Check_Counting_Numbers() {
222 // Checks whether the counting numbers satisfy the fundamental relation
223
224 bool all_valid = true;
225 for( vector<Var>::const_iterator n = vars().begin(); n != vars().end(); n++ ) {
226 double c_n = 0.0;
227 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
228 if( OR(alpha).vars() && *n )
229 c_n += OR(alpha).c();
230 for( size_t beta = 0; beta < nr_IRs(); beta++ )
231 if( IR(beta) && *n )
232 c_n += IR(beta).c();
233 if( fabs(c_n - 1.0) > 1e-15 ) {
234 all_valid = false;
235 cout << "WARNING: counting numbers do not satisfy relation for " << *n << "(c_n = " << c_n << ")." << endl;
236 }
237 }
238
239 return all_valid;
240 }
241
242
243 void RegionGraph::RecomputeORs() {
244 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
245 OR(alpha).fill( 1.0 );
246 for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
247 OR( I->second ) *= factor( I->first );
248 }
249
250
251 void RegionGraph::RecomputeORs( const VarSet &ns ) {
252 for( size_t alpha = 0; alpha < nr_ORs(); alpha++ )
253 if( OR(alpha).vars() && ns )
254 OR(alpha).fill( 1.0 );
255 for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
256 if( OR( I->second ).vars() && ns )
257 OR( I->second ) *= factor( I->first );
258 }
259
260
261 void RegionGraph::RecomputeOR( size_t I ) {
262 if( _fac2OR.count(I) ) {
263 size_t alpha = _fac2OR[I];
264 OR(alpha).fill( 1.0 );
265 for( fac2OR_cit I = _fac2OR.begin(); I != _fac2OR.end(); I++ )
266 if( I->second == alpha )
267 OR(alpha) *= factor( I->first );
268 }
269 }
270
271
272 ostream & operator << (ostream & os, const RegionGraph & rg) {
273 os << "Outer regions" << endl;
274 for( size_t alpha = 0; alpha < rg.nr_ORs(); alpha++ )
275 os << rg.OR(alpha).vars() << ": c = " << rg.OR(alpha).c() << endl;
276
277 os << "Inner regions" << endl;
278 for( size_t beta = 0; beta < rg.nr_IRs(); beta++ )
279 os << (VarSet)rg.IR(beta) << ": c = " << rg.IR(beta).c() << endl;
280
281 return(os);
282 }
283
284
285 } // end of namespace dai