Various fixes to make libDAI build successfully under modern OSX versions
[libdai.git] / src / cobwebgraph.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * Copyright (c) 2006-2012, The libDAI authors. All rights reserved.
4 *
5 * Use of this source code is governed by a BSD-style license that can be found in the LICENSE file.
6 */
7
8
9 #include <dai/cobwebgraph.h>
10 #include <algorithm>
11 #include <map>
12 #include <set>
13
14
15 namespace dai {
16
17
18 using namespace std;
19
20
21 void CobwebGraph::setRgn( vector<SmallSet<size_t> > partition, NeighborType neighbors, bool debugging ) {
22 if( checkPartition( partition ) )
23 isPartition = true;
24 else
25 isPartition = false;
26
27 eraseNonMaximal( partition );
28 if( debugging )
29 cerr << "isPartition:" << isPartition << endl;
30 _INRs = partition;
31 if( debugging )
32 cerr << "setting external factors and regional factors" << endl;
33 setExtnFact();
34 if( debugging )
35 cerr << "setting msgs" << endl;
36 setMSGs( neighbors );
37 // if(!isPartition)
38 if( debugging )
39 cerr << "setting the counting numbers" << endl;
40 if( !isPartition )
41 setCountingNumbers( debugging );
42 if( debugging )
43 cerr << "setting var2dr" << endl;
44 setVar2CW();
45 }
46
47
48 void CobwebGraph::eraseNonMaximal( vector<SmallSet<size_t> >& partition ) {
49 vector<SmallSet<size_t> >::iterator del = partition.begin();
50 while( del != partition.end() ){
51 bool isNonMax = false;
52 for( vector<SmallSet<size_t> >::iterator it = partition.begin(); it != partition.end(); it++ )
53 if( (*it) >> (*del) && it != del ) {
54 isNonMax = true;
55 break;
56 }
57 if( isNonMax )
58 del = partition.erase( del );
59 else
60 del++;
61 }
62 }
63
64
65 void CobwebGraph::setCountingNumbers( bool debugging ) {
66 // should handle the case when one region in the msg-region graph is subset of another or even to top regions are the same
67 _cn.reserve( nrCWs() );
68 for( size_t R = 0; R < _INRs.size(); R++ ) {
69 vector<VarSet> topR;
70 // finding the intersection of messages
71 SmallSet<VarSet> betas;
72 for( size_t m1 = 0; m1 < M(R).size(); m1++ ) {
73 topR.push_back( M(R,m1).msg.vars() );
74 for( size_t m2 = 0; m2 < M(R).size(); m2++ ) {
75 if( m1 == m2 )
76 continue;
77 VarSet b = M(R,m1).msg.vars() & M(R,m2).msg.vars();
78 if( b.size() == 0 )
79 continue;
80 betas.insert( b );
81 }
82 }
83 if( debugging )
84 cerr << "topR: " << topR << endl;
85
86 // finding the intersections of intersections
87 bool somechange = true;
88 while( somechange ) {
89 SmallSet<VarSet> newbetas;
90 somechange = false;
91 bforeach( const VarSet &b1, betas ) {
92 bforeach( const VarSet &b2, betas ) {
93 if( b1 == b2 )
94 continue;
95 VarSet b3 = b1 & b2;
96 if( betas.contains(b3) || b3.size() == 0 )
97 continue;
98 newbetas |= b3;
99 somechange = true;
100 }
101 }
102 betas |= newbetas;
103 }
104 if( debugging )
105 cerr << "betas: " << betas << endl;
106
107 // set the counting number
108 _cn.push_back( map<VarSet, pair<int, vector<size_t> > >() );
109 // adding sub-regions of every message
110 for( size_t i = 0; i < topR.size(); i++ ) {
111 M(R,i).subregions.clear();
112 bforeach( const VarSet& b, betas )
113 if( b << topR[i] )
114 M(R,i).subregions.push_back( b );
115 }
116 SmallSet<VarSet> subVisited;
117 SmallSet<VarSet> topRSet( topR.begin(), topR.end(), topR.size() );
118 // looks to see if all parents of a sub-region got their counting number
119 while( !betas.empty() ) {
120 bforeach( const VarSet &beta, betas ) {
121 bool allparentsset = true;
122 bforeach( const VarSet &beta2, betas )
123 if( beta2 >> beta && beta2 != beta ) {
124 allparentsset = false;
125 break;
126 }
127 if( allparentsset ) {
128 // the first in the pair is cn and the second the index of top regions containing it
129 _cn[R][beta] = make_pair( 1, vector<size_t>() );
130 for( size_t TR = 0; TR < topR.size(); TR++ ) {
131 if( topR[TR] >> beta ) {
132 _cn[R][beta].first--;
133 _cn[R][beta].second.push_back( TR );
134 }
135 }
136 bforeach( const VarSet& possibleparent, subVisited )
137 if( possibleparent >> beta )
138 _cn[R][beta].first -= _cn[R][possibleparent].first;
139
140 if( debugging )
141 cerr << "cn[" << R << "][" << beta << "] <- " << _cn[R][beta] << endl;
142 subVisited.insert( beta );
143 betas /= beta;
144 break; // since betas has changed we need to enter the loop again
145 }
146 }
147 }
148 // if( debugging )
149 // cerr << "cn[" << R << "] " << _cn[R] << endl;
150 }
151 }
152
153
154 bool CobwebGraph::checkPartition( const vector<SmallSet<size_t> >& regions ) const {
155 for( size_t i = 0; i < regions.size(); i++ )
156 for( size_t j = 0; j < regions.size(); j++ ) {
157 if( j == i )
158 continue;
159 if( regions[i].intersects( regions[j] ) )
160 return false;
161 }
162 return true;
163 }
164
165
166 void CobwebGraph::setVar2CW() {
167 _var2CW.clear();
168 _var2CW.resize( nrVars() );
169 for( size_t R = 0; R < _INRs.size(); R++ ) {
170 bforeach( size_t i, _INRs[R] )
171 _var2CW[i].push_back( R );
172 }
173 }
174
175
176 // assumes availablity of externals and internals
177 void CobwebGraph::setMSGs( NeighborType neighbors ) {
178 // initialize
179 _outM.clear();
180 _outM.reserve( nrCWs() );
181 for( size_t R = 0; R < nrCWs(); R++ )
182 _outM.push_back( vector<VarSet>() );
183 _M.clear();
184 _M.reserve( _INRs.size() );
185 for( size_t R = 0; R < _INRs.size(); R++ ) {
186 _M.push_back( vector<Connection>() );
187 for( size_t R2 = 0; R2 < _INRs.size(); R2++ ) {
188 if( R2 == R )
189 continue;
190 SmallSet<size_t> tmpSet = _EXRs[R] & (_INRs[R2]); // if R2 has boundary of R1 inside it
191 if( tmpSet.size() == 0 )
192 continue;
193
194 Connection tmp;
195 tmp.fc = (_Rfs[R] & _Rfs[R2]).elements();
196 tmp.my = R;
197 tmp.iter = _M[R].size();
198 tmp.his = R2;
199 tmp.dual = _outM[R2].size();
200 VarSet tmpVarSet;
201 bforeach( const size_t &i, tmpSet )
202 tmpVarSet.insert( var(i) );
203 tmp.msg = Factor( tmpVarSet, (Real)1.0 ); //*
204 _outM[R2].push_back( tmpVarSet ); //*
205 // tmp.varinds = tmpSet.elements();
206 _M[R].push_back( tmp );
207 }
208 if( neighbors == NeighborType::CLOSEST || neighbors == NeighborType::TOP ) {
209 bool foundone = true;
210 while( foundone ) {
211 for( size_t Rfirst = 0; Rfirst < _M[R].size(); Rfirst++ ) {
212 foundone = false;
213 for( size_t Rsec = 0; Rsec < _M[R].size(); Rsec++ ) {
214 if( Rfirst == Rsec )
215 continue;
216 if( _M[R][Rfirst].msg.vars() >> _M[R][Rsec].msg.vars() ) { // one message is subseteq of the other
217 // only if one exclusively contains the other remove it
218 if( ((neighbors == NeighborType::TOP || neighbors == NeighborType::CLOSEST) && _M[R][Rfirst].msg.vars().size() > _M[R][Rsec].msg.vars().size()) ) {
219 _M[R].erase( _M[R].begin() + Rsec );
220 foundone = true;
221 break;
222 } else if( neighbors == NeighborType::CLOSEST && // if both have the same size then in case of CLOSEST delete the second one if it has smaller inner region overlap
223 (INRs(_M[R][Rfirst].his) & INRs(R)).size() > (INRs(_M[R][Rsec].his) & INRs(R)).size() ) {
224 _M[R].erase( _M[R].begin() + Rsec );
225 foundone = true;
226 break;
227 }
228 }
229 }
230 if( foundone )
231 break;
232 }
233 }
234 }
235 }
236 }
237
238
239 void CobwebGraph::setExtnFact() {
240 // finds and sets all the externals using _INRs
241 SmallSet<size_t> allpots;
242 _EXRs.clear();
243 _Rfs.clear();
244 _EXRs.reserve( _INRs.size() );
245 _Rfs.reserve( _INRs.size() );
246 _Rifs.reserve( _INRs.size() );
247 _Rxfs.reserve( _INRs.size() );
248 for( size_t R = 0; R < _INRs.size(); R++ ) {
249 _EXRs.push_back( SmallSet<size_t>() );
250 _Rfs.push_back( SmallSet<size_t>() );
251 _Rifs.push_back( SmallSet<size_t>() );
252 _Rxfs.push_back( SmallSet<size_t>() );
253 bforeach( const size_t &i, _INRs[R] ) {
254 _EXRs[R] |= deltai(i) / _INRs[R];
255 bforeach( const Neighbor &I, nbV(i) ) {
256 _Rfs[R] |= I.node;
257 if( factor(I).vars() << inds2vars(INRs(R).elements()) )
258 _Rifs[R] |= I.node;
259 else
260 _Rxfs[R] |= I.node;
261 }
262 }
263 }
264 }
265
266
267 } // end of namespace dai