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