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