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