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