Finished release 0.2.4
[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-2010 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 void RegionGraph::construct( const FactorGraph &fg, const std::vector<VarSet> &ors, const std::vector<Region> &irs, const std::vector<std::pair<size_t,size_t> > &edges ) {
27 // Copy factor graph structure
28 FactorGraph::operator=( fg );
29
30 // Copy inner regions
31 IRs = irs;
32
33 // Construct outer regions (giving them counting number 1.0)
34 ORs.clear();
35 ORs.reserve( ors.size() );
36 foreach( const VarSet &alpha, ors )
37 ORs.push_back( FRegion(Factor(alpha, 1.0), 1.0) );
38
39 // For each factor, find an outer region that subsumes that factor.
40 // Then, multiply the outer region with that factor.
41 fac2OR.clear();
42 fac2OR.reserve( nrFactors() );
43 for( size_t I = 0; I < nrFactors(); I++ ) {
44 size_t alpha;
45 for( alpha = 0; alpha < nrORs(); alpha++ )
46 if( OR(alpha).vars() >> factor(I).vars() ) {
47 fac2OR.push_back( alpha );
48 break;
49 }
50 DAI_ASSERT( alpha != nrORs() );
51 }
52 RecomputeORs();
53
54 // Create bipartite graph
55 G.construct( nrORs(), nrIRs(), edges.begin(), edges.end() );
56 }
57
58
59 void RegionGraph::constructCVM( const FactorGraph &fg, const std::vector<VarSet> &cl ) {
60 // Retain only maximal clusters
61 ClusterGraph cg( cl );
62 cg.eraseNonMaximal();
63
64 // Create inner regions - first pass
65 set<VarSet> betas;
66 for( size_t alpha = 0; alpha < cg.clusters.size(); alpha++ )
67 for( size_t alpha2 = alpha; (++alpha2) != cg.clusters.size(); ) {
68 VarSet intersection = cg.clusters[alpha] & cg.clusters[alpha2];
69 if( intersection.size() > 0 )
70 betas.insert( intersection );
71 }
72
73 // Create inner regions - subsequent passes
74 set<VarSet> new_betas;
75 do {
76 new_betas.clear();
77 for( set<VarSet>::const_iterator gamma = betas.begin(); gamma != betas.end(); gamma++ )
78 for( set<VarSet>::const_iterator gamma2 = gamma; (++gamma2) != betas.end(); ) {
79 VarSet intersection = (*gamma) & (*gamma2);
80 if( (intersection.size() > 0) && (betas.count(intersection) == 0) )
81 new_betas.insert( intersection );
82 }
83 betas.insert(new_betas.begin(), new_betas.end());
84 } while( new_betas.size() );
85
86 // Create inner regions - final phase
87 vector<Region> irs;
88 irs.reserve( betas.size() );
89 for( set<VarSet>::const_iterator beta = betas.begin(); beta != betas.end(); beta++ )
90 irs.push_back( Region(*beta,0.0) );
91
92 // Create edges
93 vector<pair<size_t,size_t> > edges;
94 for( size_t beta = 0; beta < irs.size(); beta++ )
95 for( size_t alpha = 0; alpha < cg.clusters.size(); alpha++ )
96 if( cg.clusters[alpha] >> irs[beta] )
97 edges.push_back( pair<size_t,size_t>(alpha,beta) );
98
99 // Construct region graph
100 construct( fg, cg.clusters, irs, edges );
101
102 // Calculate counting numbers
103 calcCountingNumbers();
104 }
105
106
107 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 ) {
108 vector<VarSet> ors_alt;
109 ors_alt.reserve( ors.size() );
110 for( size_t alpha = 0; alpha < ors.size(); alpha++ ) {
111 ors_alt.push_back( ors[alpha] );
112 DAI_ASSERT( ors[alpha].c() == 1.0 );
113 }
114 construct( fg, ors_alt, irs, edges );
115
116 // Check counting numbers
117 #ifdef DAI_DEBUG
118 checkCountingNumbers();
119 #endif
120 }
121
122
123 void RegionGraph::calcCountingNumbers() {
124 // Calculates counting numbers of inner regions based upon counting numbers of outer regions
125
126 vector<vector<size_t> > ancestors(nrIRs());
127 boost::dynamic_bitset<> assigned(nrIRs());
128 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
129 IR(beta).c() = 0.0;
130 for( size_t beta2 = 0; beta2 < nrIRs(); beta2++ )
131 if( (beta2 != beta) && IR(beta2) >> IR(beta) )
132 ancestors[beta].push_back(beta2);
133 }
134
135 bool new_counting;
136 do {
137 new_counting = false;
138 for( size_t beta = 0; beta < nrIRs(); beta++ ) {
139 if( !assigned[beta] ) {
140 bool has_unassigned_ancestor = false;
141 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); (beta2 != ancestors[beta].end()) && !has_unassigned_ancestor; beta2++ )
142 if( !assigned[*beta2] )
143 has_unassigned_ancestor = true;
144 if( !has_unassigned_ancestor ) {
145 Real c = 1.0;
146 foreach( const Neighbor &alpha, nbIR(beta) )
147 c -= OR(alpha).c();
148 for( vector<size_t>::const_iterator beta2 = ancestors[beta].begin(); beta2 != ancestors[beta].end(); beta2++ )
149 c -= IR(*beta2).c();
150 IR(beta).c() = c;
151 assigned.set(beta, true);
152 new_counting = true;
153 }
154 }
155 }
156 } while( new_counting );
157 }
158
159
160 bool RegionGraph::checkCountingNumbers() const {
161 // Checks whether the counting numbers satisfy the fundamental relation
162
163 bool all_valid = true;
164 for( vector<Var>::const_iterator n = vars().begin(); n != vars().end(); n++ ) {
165 Real c_n = 0.0;
166 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
167 if( OR(alpha).vars().contains( *n ) )
168 c_n += OR(alpha).c();
169 for( size_t beta = 0; beta < nrIRs(); beta++ )
170 if( IR(beta).contains( *n ) )
171 c_n += IR(beta).c();
172 if( fabs(c_n - 1.0) > 1e-15 ) {
173 all_valid = false;
174 cerr << "WARNING: counting numbers do not satisfy relation for " << *n << "(c_n = " << c_n << ")." << endl;
175 }
176 }
177
178 return all_valid;
179 }
180
181
182 void RegionGraph::RecomputeORs() {
183 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
184 OR(alpha).fill( 1.0 );
185 for( size_t I = 0; I < nrFactors(); I++ )
186 if( fac2OR[I] != -1U )
187 OR( fac2OR[I] ) *= factor( I );
188 }
189
190
191 void RegionGraph::RecomputeORs( const VarSet &ns ) {
192 for( size_t alpha = 0; alpha < nrORs(); alpha++ )
193 if( OR(alpha).vars().intersects( ns ) )
194 OR(alpha).fill( 1.0 );
195 for( size_t I = 0; I < nrFactors(); I++ )
196 if( fac2OR[I] != -1U )
197 if( OR( fac2OR[I] ).vars().intersects( ns ) )
198 OR( fac2OR[I] ) *= factor( I );
199 }
200
201
202 void RegionGraph::RecomputeOR( size_t I ) {
203 DAI_ASSERT( I < nrFactors() );
204 if( fac2OR[I] != -1U ) {
205 size_t alpha = fac2OR[I];
206 OR(alpha).fill( 1.0 );
207 for( size_t J = 0; J < nrFactors(); J++ )
208 if( fac2OR[J] == alpha )
209 OR(alpha) *= factor( J );
210 }
211 }
212
213
214 /// Send RegionGraph to output stream
215 ostream & operator << (ostream & os, const RegionGraph & rg) {
216 os << "Outer regions" << endl;
217 for( size_t alpha = 0; alpha < rg.nrORs(); alpha++ )
218 os << alpha << ": " << rg.OR(alpha).vars() << ": c = " << rg.OR(alpha).c() << endl;
219
220 os << "Inner regions" << endl;
221 for( size_t beta = 0; beta < rg.nrIRs(); beta++ )
222 os << beta << ": " << (VarSet)rg.IR(beta) << ": c = " << rg.IR(beta).c() << endl;
223
224 os << "Edges" << endl;
225 for( size_t alpha = 0; alpha < rg.nrORs(); alpha++ )
226 foreach( const RegionGraph::Neighbor &beta, rg.nbOR(alpha) )
227 os << alpha << "->" << beta << endl;
228
229 return(os);
230 }
231
232
233 } // end of namespace dai