Improved ClusterGraph implementation and MaxSpanningTreePrims implementation.
[libdai.git] / include / dai / bipgraph.h
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 #ifndef __defined_libdai_bipgraph_h
23 #define __defined_libdai_bipgraph_h
24
25
26 #include <ostream>
27 #include <vector>
28 #include <cassert>
29 #include <dai/util.h>
30
31
32 namespace dai {
33
34
35 /// A BipartiteGraph represents a bipartite graph, with two types of nodes (both are numbered
36 /// as 0,1,2,...), with edges only between nodes of different type. The edges are stored as
37 /// lists of adjacent nodes for each node.
38 class BipartiteGraph {
39 public:
40 /// A Neighbor describes a neighboring node of some other node.
41 /** Iterating over all neighbors of some node i can be done in the following way:
42 * \code
43 * foreach( const BipartiteGraph::Neighbor &I, nb1(i) ) {
44 * size_t _I = I.iter;
45 * size_t _i = I.dual;
46 * // I == I.node;
47 * // The _I'th neighbor of i is I, and the _i'th neighbor of I is i:
48 * // nb1(i)[_I] == I, nb2(I)[_i] == i
49 * }
50 * \endcode
51 */
52 struct Neighbor {
53 /// iter corresponds to the index of this Neighbor entry in the list of neighbors
54 unsigned iter;
55 /// node contains the number of the neighboring node
56 unsigned node;
57 /// dual contains the "dual" iter
58 unsigned dual;
59 /// cast to unsigned returns node
60 operator unsigned () const { return node; };
61 };
62
63 /// Neighbors is a vector of Neigbor entries
64 typedef std::vector<Neighbor> Neighbors;
65
66 /// Edge is used as index of an edge
67 typedef std::pair<size_t,size_t> Edge;
68
69 private:
70 /// _nb1 contains for each node of the first kind a list of its neighbors
71 std::vector<Neighbors> _nb1;
72 /// _nb2 contains for each node of the second kind a list of its neighbors
73 std::vector<Neighbors> _nb2;
74
75 /// Used internally by isTree()
76 struct levelType {
77 std::vector<size_t> ind1; // indices of vertices of type 1
78 std::vector<size_t> ind2; // indices of vertices of type 2
79 };
80
81 public:
82 /// Default constructor
83 BipartiteGraph() : _nb1(), _nb2() {}
84
85 /// Copy constructor
86 BipartiteGraph( const BipartiteGraph & x ) : _nb1(x._nb1), _nb2(x._nb2) {}
87
88 /// Assignment operator
89 BipartiteGraph & operator=( const BipartiteGraph & x ) {
90 if( this != &x ) {
91 _nb1 = x._nb1;
92 _nb2 = x._nb2;
93 }
94 return *this;
95 }
96
97 /// Create bipartite graph from a range of edges, encoded as pairs of node numbers
98 /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
99 /// to the node of the first type, and the second integer corresponds to the node of the
100 /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
101 /// of the second type. The values between the iterators begin and end should be of type Edge.
102 template<typename EdgeInputIterator>
103 void create( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end );
104
105 /// Construct bipartite graph from a range of edges, encoded as pairs of node numbers
106 /// (more precisely, a std::pair<unsigned, unsigned> where the first integer corresponds
107 /// to the node of the first type, and the second integer corresponds to the node of the
108 /// second type). nr1 is the number of nodes of the first type, nr2 the number of nodes
109 /// of the second type. The values between the iterators begin and end should be of type Edge.
110 template<typename EdgeInputIterator>
111 BipartiteGraph( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) : _nb1( nr1 ), _nb2( nr2 ) {
112 create( nr1, nr2, begin, end );
113 }
114
115 /// Returns constant reference to the _i2'th neighbor of node i1 of first type
116 const Neighbor & nb1( size_t i1, size_t _i2 ) const { return _nb1[i1][_i2]; }
117 /// Returns reference to the _i2'th neighbor of node i1 of first type
118 Neighbor & nb1( size_t i1, size_t _i2 ) { return _nb1[i1][_i2]; }
119
120 /// Returns constant reference to the _i1'th neighbor of node i2 of second type
121 const Neighbor & nb2( size_t i2, size_t _i1 ) const { return _nb2[i2][_i1]; }
122 /// Returns reference to the _i1'th neighbor of node i2 of second type
123 Neighbor & nb2( size_t i2, size_t _i1 ) { return _nb2[i2][_i1]; }
124
125 /// Returns constant reference to all neighbors of node of first type
126 const Neighbors & nb1( size_t i1 ) const { return _nb1[i1]; }
127 /// Returns reference to all neighbors of node of first type
128 Neighbors & nb1( size_t i1 ) { return _nb1[i1]; }
129
130 /// Returns constant reference to all neighbors of node of second type
131 const Neighbors & nb2( size_t i2 ) const { return _nb2[i2]; }
132 /// Returns reference to all neighbors of node of second type
133 Neighbors & nb2( size_t i2 ) { return _nb2[i2]; }
134
135 /// Returns number of nodes of first type
136 size_t nr1() const { return _nb1.size(); }
137 /// Returns number of nodes of second type
138 size_t nr2() const { return _nb2.size(); }
139
140 /// Calculates the number of edges
141 size_t nrEdges() const {
142 size_t sum = 0;
143 for( size_t i1 = 0; i1 < nr1(); i1++ )
144 sum += nb1(i1).size();
145 return sum;
146 }
147
148 /// Add node of first type without neighbors.
149 void add1() {
150 _nb1.push_back( Neighbors() );
151 }
152
153 /// Add node of first type without neighbors.
154 void add2() {
155 _nb2.push_back( Neighbors() );
156 }
157
158 /// Add node of first type with neighbors specified by a range.
159 /// For improved efficiency, the size of the range may be specified by sizeHint.
160 /// *NodeInputIterator should be a size_t.
161 template <typename NodeInputIterator>
162 void add1( NodeInputIterator begin, NodeInputIterator end, size_t sizeHint = 0 ) {
163 Neighbors nbs1new;
164 nbs1new.reserve( sizeHint );
165 size_t iter = 0;
166 for( NodeInputIterator it = begin; it != end; ++it ) {
167 assert( *it < nr2() );
168
169 Neighbor nb1new;
170 nb1new.node = *it;
171 nb1new.iter = iter;
172 nb1new.dual = nb2(*it).size();
173
174 Neighbor nb2new;
175 nb2new.node = nr1();
176 nb2new.iter = nb2(*it).size();
177 nb2new.dual = iter++;
178
179 nbs1new.push_back( nb1new );
180 nb2( *it ).push_back( nb2new );
181 }
182 _nb1.push_back( nbs1new );
183 }
184
185 /// Add node of second type with neighbors specified by a range.
186 /// For improved efficiency, the size of the range may be specified by sizeHint.
187 /// *NodeInputIterator should be a size_t.
188 template <typename NodeInputIterator>
189 void add2( NodeInputIterator begin, NodeInputIterator end, size_t sizeHint = 0 ) {
190 Neighbors nbs2new;
191 nbs2new.reserve( sizeHint );
192 size_t iter = 0;
193 for( NodeInputIterator it = begin; it != end; ++it ) {
194 assert( *it < nr1() );
195
196 Neighbor nb2new;
197 nb2new.node = *it;
198 nb2new.iter = iter;
199 nb2new.dual = nb1(*it).size();
200
201 Neighbor nb1new;
202 nb1new.node = nr2();
203 nb1new.iter = nb1(*it).size();
204 nb1new.dual = iter++;
205
206 nbs2new.push_back( nb2new );
207 nb1( *it ).push_back( nb1new );
208 }
209 _nb2.push_back( nbs2new );
210 }
211
212 /// Remove node of first type
213 void erase1( size_t n1 ) {
214 assert( n1 < nr1() );
215 // Erase neighbor entry of node n1
216 _nb1.erase( _nb1.begin() + n1 );
217 // Adjust neighbor entries of nodes of second type
218 for( size_t n2 = 0; n2 < nr2(); n2++ )
219 for( size_t iter = 0; iter < nb2(n2).size(); ) {
220 if( nb2(n2, iter).node == n1 ) {
221 // delete this entry, because it points to the deleted node
222 nb2(n2).erase( nb2(n2).begin() + iter );
223 // adjust all subsequent entries:
224 // update their iter and the corresponding dual of the neighboring node of the first kind
225 for( size_t newiter = iter; newiter < nb2(n2).size(); newiter++ ) {
226 nb2( n2, newiter ).iter = newiter;
227 nb1( nb2(n2, newiter).node, nb2(n2, newiter).dual ).dual = newiter;
228 }
229 } else if( nb2(n2, iter).node > n1 ) {
230 nb2(n2, iter).node--;
231 iter++;
232 } else
233 iter++;
234 }
235 }
236
237 /// Remove node of second type
238 void erase2( size_t n2 ) {
239 assert( n2 < nr2() );
240 // Erase neighbor entry of node n2
241 _nb2.erase( _nb2.begin() + n2 );
242 // Adjust neighbor entries of nodes of first type
243 for( size_t n1 = 0; n1 < nr1(); n1++ )
244 for( size_t iter = 0; iter < nb1(n1).size(); ) {
245 if( nb1(n1, iter).node == n2 ) {
246 // delete this entry, because it points to the deleted node
247 nb1(n1).erase( nb1(n1).begin() + iter );
248 // adjust all subsequent entries:
249 // update their iter and the corresponding dual of the neighboring node of the first kind
250 for( size_t newiter = iter; newiter < nb1(n1).size(); newiter++ ) {
251 nb1( n1, newiter ).iter = newiter;
252 nb2( nb1(n1, newiter).node, nb1(n1, newiter).dual ).dual = newiter;
253 }
254 } else if( nb1(n1, iter).node > n2 ) {
255 nb1(n1, iter).node--;
256 iter++;
257 } else
258 iter++;
259 }
260 }
261
262 /// Calculate neighbors of neighbors of node n1 of first type.
263 /// If include == true, include n1 itself, otherwise exclude n1.
264 std::vector<size_t> delta1( size_t n1, bool include = false ) const {
265 std::vector<size_t> result;
266 foreach( const Neighbor &n2, nb1(n1) )
267 foreach( const Neighbor &m1, nb2(n2) )
268 if( include || (m1 != n1) )
269 result.push_back( m1 );
270 // remove duplicates
271 std::vector<size_t>::iterator it = unique( result.begin(), result.end() );
272 result.erase( it, result.end() );
273 return result;
274 }
275
276 /// Calculate neighbors of neighbors of node n2 of second type
277 /// If include == true, include n2 itself, otherwise exclude n2.
278 std::vector<size_t> delta2( size_t n2, bool include = false ) const {
279 std::vector<size_t> result;
280 foreach( const Neighbor &n1, nb2(n2) )
281 foreach( const Neighbor &m2, nb1(n1) )
282 if( include || (m2 != n2) )
283 result.push_back( m2 );
284 // remove duplicates
285 std::vector<size_t>::iterator it = unique( result.begin(), result.end() );
286 result.erase( it, result.end() );
287 return result;
288 }
289
290 /// Returns true if the graph is connected
291 bool isConnected() const {
292 if( nr1() == 0 ) {
293 return true;
294 } else {
295 std::vector<bool> incomponent1( nr1(), false );
296 std::vector<bool> incomponent2( nr2(), false );
297
298 incomponent1[0] = true;
299 bool found_new_nodes;
300 do {
301 found_new_nodes = false;
302
303 // For all nodes of second type, check if they are connected with the (growing) component
304 for( size_t n2 = 0; n2 < nr2(); n2++ )
305 if( !incomponent2[n2] ) {
306 foreach( const Neighbor &n1, nb2(n2) ) {
307 if( incomponent1[n1] ) {
308 found_new_nodes = true;
309 incomponent2[n2] = true;
310 break;
311 }
312 }
313 }
314
315 // For all nodes of first type, check if they are connected with the (growing) component
316 for( size_t n1 = 0; n1 < nr1(); n1++ )
317 if( !incomponent1[n1] ) {
318 foreach( const Neighbor &n2, nb1(n1) ) {
319 if( incomponent2[n2] ) {
320 found_new_nodes = true;
321 incomponent1[n1] = true;
322 break;
323 }
324 }
325 }
326 } while( found_new_nodes );
327
328 // Check if there are remaining nodes (not in the component)
329 bool all_connected = true;
330 for( size_t n1 = 0; (n1 < nr1()) && all_connected; n1++ )
331 if( !incomponent1[n1] )
332 all_connected = false;
333 for( size_t n2 = 0; (n2 < nr2()) && all_connected; n2++ )
334 if( !incomponent2[n2] )
335 all_connected = false;
336
337 return all_connected;
338 }
339 }
340
341 /// Returns true if the graph is a tree, i.e., if it is singly connected and connected.
342 /// This is equivalent to whether for each pair of vertices in the graph, there exists
343 /// a unique path in the graph that starts at the first and ends at the second vertex.
344 bool isTree() const {
345 using namespace std;
346 vector<levelType> levels;
347
348 bool foundCycle = false;
349 size_t nr_1 = 0;
350 size_t nr_2 = 0;
351
352 if( nr1() == 0 || nr2() == 0 )
353 return true;
354 else {
355 levelType newLevel;
356 do {
357 newLevel.ind1.clear();
358 newLevel.ind2.clear();
359 if( levels.size() == 0 ) {
360 size_t n1 = 0;
361 // add n1 to ind1
362 newLevel.ind1 = vector<size_t>( 1, n1 );
363 // add all neighbors of n1 to ind2
364 newLevel.ind2.reserve( nb1(n1).size() );
365 foreach( const Neighbor &n2, nb1(n1) )
366 newLevel.ind2.push_back( n2 );
367 } else {
368 const levelType &prevLevel = levels.back();
369 // build newLevel.ind1
370 foreach( size_t n2, prevLevel.ind2 ) { // for all n2 in the previous level
371 foreach( const Neighbor &n1, nb2(n2) ) { // for all neighbors n1 of n2
372 if( find( prevLevel.ind1.begin(), prevLevel.ind1.end(), n1 ) == prevLevel.ind1.end() ) { // n1 not in previous level
373 if( find( newLevel.ind1.begin(), newLevel.ind1.end(), n1 ) != newLevel.ind1.end() )
374 foundCycle = true; // n1 already in new level: we found a cycle
375 else
376 newLevel.ind1.push_back( n1 ); // add n1 to new level
377 }
378 if( foundCycle )
379 break;
380 }
381 if( foundCycle )
382 break;
383 }
384 // build newLevel.ind2
385 foreach( size_t n1, newLevel.ind1 ) { // for all n1 in this level
386 foreach( const Neighbor &n2, nb1(n1) ) { // for all neighbors n2 of n1
387 if( find( prevLevel.ind2.begin(), prevLevel.ind2.end(), n2 ) == prevLevel.ind2.end() ) { // n2 not in previous level
388 if( find( newLevel.ind2.begin(), newLevel.ind2.end(), n2 ) != newLevel.ind2.end() )
389 foundCycle = true; // n2 already in new level: we found a cycle
390 else
391 newLevel.ind2.push_back( n2 ); // add n2 to new level
392 }
393 if( foundCycle )
394 break;
395 }
396 if( foundCycle )
397 break;
398 }
399 }
400 levels.push_back( newLevel );
401 nr_1 += newLevel.ind1.size();
402 nr_2 += newLevel.ind2.size();
403 } while( ((newLevel.ind1.size() != 0) || (newLevel.ind2.size() != 0)) && !foundCycle );
404 if( nr_1 == nr1() && nr_2 == nr2() && !foundCycle )
405 return true;
406 else
407 return false;
408 }
409 }
410
411 /// Stream to output stream os in graphviz .dot syntax
412 void display( std::ostream& os ) const {
413 using namespace std;
414 os << "graph G {" << endl;
415 os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
416 for( size_t n1 = 0; n1 < nr1(); n1++ )
417 os << "\tx" << n1 << ";" << endl;
418 os << "node[shape=box,width=0.3,height=0.3,fixedsize=true];" << endl;
419 for( size_t n2 = 0; n2 < nr2(); n2++ )
420 os << "\ty" << n2 << ";" << endl;
421 for( size_t n1 = 0; n1 < nr1(); n1++ )
422 foreach( const Neighbor &n2, nb1(n1) )
423 os << "\tx" << n1 << " -- y" << n2 << ";" << endl;
424 os << "}" << endl;
425 }
426 };
427
428
429 template<typename EdgeInputIterator>
430 void BipartiteGraph::create( size_t nr1, size_t nr2, EdgeInputIterator begin, EdgeInputIterator end ) {
431 _nb1.clear();
432 _nb1.resize( nr1 );
433 _nb2.clear();
434 _nb2.resize( nr2 );
435
436 for( EdgeInputIterator e = begin; e != end; e++ ) {
437 // Each edge yields a neighbor pair
438 Neighbor nb_1;
439 nb_1.iter = _nb1[e->first].size();
440 nb_1.node = e->second;
441 nb_1.dual = _nb2[e->second].size();
442
443 Neighbor nb_2;
444 nb_2.iter = nb_1.dual;
445 nb_2.node = e->first;
446 nb_2.dual = nb_1.iter;
447
448 _nb1[e->first].push_back( nb_1 );
449 _nb2[e->second].push_back( nb_2 );
450 }
451 }
452
453
454 } // end of namespace dai
455
456
457 #endif