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
5 This file is part of libDAI.
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.
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.
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
25 #include <boost/dynamic_bitset.hpp>
26 #include <dai/regiongraph.h>
27 #include <dai/factorgraph.h>
28 #include <dai/clustergraph.h>
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) );
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
++ ) {
48 for( alpha
= 0; alpha
< nrORs(); alpha
++ )
49 if( OR(alpha
).vars() >> factor(I
).vars() ) {
50 fac2OR
.push_back( alpha
);
53 assert( alpha
!= nrORs() );
57 // create bipartite graph
58 G
.construct( nrORs(), nrIRs(), edges
.begin(), edges
.end() );
60 // Check counting numbers
62 Check_Counting_Numbers();
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
);
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) );
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
++ ) {
83 for( alpha
= 0; alpha
< nrORs(); alpha
++ )
84 if( OR(alpha
).vars() >> factor(I
).vars() ) {
85 fac2OR
.push_back( alpha
);
88 assert( alpha
!= nrORs() );
92 // Create inner regions - first pass
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
);
101 // Create inner regions - subsequent passes
102 set
<VarSet
> new_betas
;
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
);
111 betas
.insert(new_betas
.begin(), new_betas
.end());
112 } while( new_betas
.size() );
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) );
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
) );
128 // create bipartite graph
129 G
.construct( nrORs(), nrIRs(), edges
.begin(), edges
.end() );
131 // Calculate counting numbers
132 Calc_Counting_Numbers();
134 // Check counting numbers
136 Check_Counting_Numbers();
141 void RegionGraph::Calc_Counting_Numbers() {
142 // Calculates counting numbers of inner regions based upon counting numbers of outer regions
144 vector
<vector
<size_t> > ancestors(nrIRs());
145 boost::dynamic_bitset
<> assigned(nrIRs());
146 for( size_t beta
= 0; beta
< nrIRs(); beta
++ ) {
148 for( size_t beta2
= 0; beta2
< nrIRs(); beta2
++ )
149 if( (beta2
!= beta
) && IR(beta2
) >> IR(beta
) )
150 ancestors
[beta
].push_back(beta2
);
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
) {
164 foreach( const Neighbor
&alpha
, nbIR(beta
) )
166 for( vector
<size_t>::const_iterator beta2
= ancestors
[beta
].begin(); beta2
!= ancestors
[beta
].end(); beta2
++ )
169 assigned
.set(beta
, true);
174 } while( new_counting
);
178 bool RegionGraph::Check_Counting_Numbers() {
179 // Checks whether the counting numbers satisfy the fundamental relation
181 bool all_valid
= true;
182 for( vector
<Var
>::const_iterator n
= vars().begin(); n
!= vars().end(); n
++ ) {
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
) )
190 if( fabs(c_n
- 1.0) > 1e-15 ) {
192 cerr
<< "WARNING: counting numbers do not satisfy relation for " << *n
<< "(c_n = " << c_n
<< ")." << endl
;
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
);
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
);
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
);
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
;
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
;
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
;
251 } // end of namespace dai