Added DAG class and various minor improvements
[libdai.git] / src / graph.cpp
1 /* This file is part of libDAI - http://www.libdai.org/
2 *
3 * libDAI is licensed under the terms of the GNU General Public License version
4 * 2, or (at your option) any later version. libDAI is distributed without any
5 * warranty. See the file COPYING for more details.
6 *
7 * Copyright (C) 2006-2010 Joris Mooij [joris dot mooij at libdai dot org]
8 * Copyright (C) 2006-2007 Radboud University Nijmegen, The Netherlands
9 */
10
11
12 #include <dai/graph.h>
13
14
15 namespace dai {
16
17
18 using namespace std;
19
20
21 GraphAL& GraphAL::addEdge( size_t n1, size_t n2, bool check ) {
22 DAI_ASSERT( n1 < nrNodes() );
23 DAI_ASSERT( n2 < nrNodes() );
24 bool exists = false;
25 if( check ) {
26 // Check whether the edge already exists
27 foreach( const Neighbor &n, nb(n1) )
28 if( n == n2 ) {
29 exists = true;
30 break;
31 }
32 }
33 if( !exists && n1 != n2 ) { // Add edge
34 Neighbor nb_1( nb(n1).size(), n2, nb(n2).size() );
35 Neighbor nb_2( nb_1.dual, n1, nb_1.iter );
36 nb(n1).push_back( nb_1 );
37 nb(n2).push_back( nb_2 );
38 }
39 return *this;
40 }
41
42
43 void GraphAL::eraseNode( size_t n ) {
44 DAI_ASSERT( n < nrNodes() );
45 // Erase neighbor entry of node n
46 _nb.erase( _nb.begin() + n );
47 // Adjust neighbor entries of nodes
48 for( size_t n2 = 0; n2 < nrNodes(); n2++ ) {
49 for( size_t iter = 0; iter < nb(n2).size(); ) {
50 Neighbor &m = nb(n2, iter);
51 if( m.node == n ) {
52 // delete this entry, because it points to the deleted node
53 nb(n2).erase( nb(n2).begin() + iter );
54 } else {
55 // update this entry and the corresponding dual of the neighboring node
56 if( m.node > n )
57 m.node--;
58 nb( m.node, m.dual ).dual = iter;
59 m.iter = iter++;
60 }
61 }
62 }
63 }
64
65
66 void GraphAL::eraseEdge( size_t n1, size_t n2 ) {
67 DAI_ASSERT( n1 < nrNodes() );
68 DAI_ASSERT( n2 < nrNodes() );
69 size_t iter;
70 // Search for edge among neighbors of n1
71 for( iter = 0; iter < nb(n1).size(); iter++ )
72 if( nb(n1, iter).node == n2 ) {
73 // Remove it
74 nb(n1).erase( nb(n1).begin() + iter );
75 break;
76 }
77 // Change the iter and dual values of the subsequent neighbors
78 for( ; iter < nb(n1).size(); iter++ ) {
79 Neighbor &m = nb( n1, iter );
80 m.iter = iter;
81 nb( m.node, m.dual ).dual = iter;
82 }
83 // Search for edge among neighbors of n2
84 for( iter = 0; iter < nb(n2).size(); iter++ )
85 if( nb(n2, iter).node == n1 ) {
86 // Remove it
87 nb(n2).erase( nb(n2).begin() + iter );
88 break;
89 }
90 // Change the iter and node values of the subsequent neighbors
91 for( ; iter < nb(n2).size(); iter++ ) {
92 Neighbor &m = nb( n2, iter );
93 m.iter = iter;
94 nb( m.node, m.dual ).dual = iter;
95 }
96 }
97
98
99 SmallSet<size_t> GraphAL::nbSet( size_t n ) const {
100 SmallSet<size_t> result;
101 foreach( const Neighbor &m, nb(n) )
102 result |= m;
103 return result;
104 }
105
106
107 bool GraphAL::isConnected() const {
108 if( nrNodes() == 0 ) {
109 return true;
110 } else {
111 std::vector<bool> incomponent( nrNodes(), false );
112
113 incomponent[0] = true;
114 bool found_new_nodes;
115 do {
116 found_new_nodes = false;
117
118 // For all nodes, check if they are connected with the (growing) component
119 for( size_t n1 = 0; n1 < nrNodes(); n1++ )
120 if( !incomponent[n1] ) {
121 foreach( const Neighbor &n2, nb(n1) ) {
122 if( incomponent[n2] ) {
123 found_new_nodes = true;
124 incomponent[n1] = true;
125 break;
126 }
127 }
128 }
129 } while( found_new_nodes );
130
131 // Check if there are remaining nodes (not in the component)
132 bool all_connected = true;
133 for( size_t n1 = 0; (n1 < nrNodes()) && all_connected; n1++ )
134 if( !incomponent[n1] )
135 all_connected = false;
136
137 return all_connected;
138
139 // BGL implementation is slower...
140 /* using namespace boost;
141 typedef adjacency_list< vecS, vecS, undirectedS, property<vertex_distance_t, int> > boostGraphAL;
142 typedef pair<size_t, size_t> E;
143
144 // Copy graph structure into boostGraphAL object
145 vector<E> edges;
146 edges.reserve( nrEdges() );
147 for( size_t n1 = 0; n1 < nrNodes(); n1++ )
148 foreach( const Neighbor &n2, nb(n1) )
149 if( n1 < n2 )
150 edges.push_back( E( n1, n2 ) );
151 boostGraphAL g( edges.begin(), edges.end(), nrNodes() );
152
153 // Construct connected components using Boost GraphAL Library
154 std::vector<int> component( num_vertices( g ) );
155 int num_comp = connected_components( g, make_iterator_property_map(component.begin(), get(vertex_index, g)) );
156
157 return (num_comp == 1);
158 */
159 }
160 }
161
162
163 bool GraphAL::isTree() const {
164 typedef vector<Edge> levelType; // first is node, second is its parent
165 vector<levelType> levels;
166
167 if( nrNodes() == 0 )
168 return true;
169 else {
170 // start with root node 0
171 levels.push_back( levelType( 1, Edge( 0, 0 ) ) );
172 size_t treeSize = 1;
173 bool foundCycle = false;
174 do {
175 levels.push_back( levelType() );
176 const levelType &prevLevel = levels[levels.size() - 2];
177 // build new level: add all neighbors of nodes in the previous level
178 // (without backtracking), aborting if a cycle is detected
179 for( size_t e = 0; e < prevLevel.size(); e++ ) {
180 size_t n2 = prevLevel[e].first; // for all nodes n2 in the previous level
181 foreach( const Neighbor &n1, nb(n2) ) { // for all neighbors n1 of n2
182 if( n1 != prevLevel[e].second ) { // no backtracking allowed
183 for( size_t l = 0; l < levels.size() && !foundCycle; l++ )
184 for( size_t f = 0; f < levels[l].size() && !foundCycle; f++ )
185 if( levels[l][f].first == n1 )
186 // n1 has been visited before -> found a cycle
187 foundCycle = true;
188 if( !foundCycle )
189 // add n1 (and its parent n2) to current level
190 levels.back().push_back( Edge( n1, n2 ) );
191 }
192 if( foundCycle )
193 break;
194 }
195 if( foundCycle )
196 break;
197 }
198 treeSize += levels.back().size();
199 } while( (levels.back().size() != 0) && !foundCycle );
200 if( treeSize == nrNodes() && !foundCycle )
201 return true;
202 else
203 return false;
204 }
205 }
206
207
208 void GraphAL::printDot( std::ostream& os ) const {
209 os << "graph GraphAL {" << endl;
210 os << "node[shape=circle,width=0.4,fixedsize=true];" << endl;
211 for( size_t n = 0; n < nrNodes(); n++ )
212 os << "\tx" << n << ";" << endl;
213 for( size_t n1 = 0; n1 < nrNodes(); n1++ )
214 foreach( const Neighbor &n2, nb(n1) )
215 if( n1 < n2 )
216 os << "\tx" << n1 << " -- x" << n2 << ";" << endl;
217 os << "}" << endl;
218 }
219
220
221 void GraphAL::checkConsistency() const {
222 size_t N = nrNodes();
223 for( size_t n1 = 0; n1 < N; n1++ ) {
224 size_t iter = 0;
225 foreach( const Neighbor &n2, nb(n1) ) {
226 DAI_ASSERT( n2.iter == iter );
227 DAI_ASSERT( n2.node < N );
228 DAI_ASSERT( n2.dual < nb(n2).size() );
229 DAI_ASSERT( nb(n2, n2.dual) == n1 );
230 iter++;
231 }
232 }
233 }
234
235
236 GraphAL createGraphFull( size_t N ) {
237 GraphAL result( N );
238 for( size_t i = 0; i < N; i++ )
239 for( size_t j = i+1; j < N; j++ )
240 result.addEdge( i, j, false );
241 return result;
242 }
243
244
245 GraphAL createGraphGrid( size_t N1, size_t N2, bool periodic ) {
246 GraphAL result( N1*N2 );
247 if( N1 == 1 && N2 == 1 )
248 return result;
249 for( size_t i1 = 0; i1 < N1; i1++ )
250 for( size_t i2 = 0; i2 < N2; i2++ ) {
251 if( i1+1 < N1 || periodic )
252 result.addEdge( i1*N2 + i2, ((i1+1)%N1)*N2 + i2, N1 <= 2 );
253 if( i2+1 < N2 || periodic )
254 result.addEdge( i1*N2 + i2, i1*N2 + ((i2+1)%N2), N2 <= 2 );
255 }
256 return result;
257 }
258
259
260 GraphAL createGraphGrid3D( size_t N1, size_t N2, size_t N3, bool periodic ) {
261 GraphAL result( N1*N2*N3 );
262 for( size_t i1 = 0; i1 < N1; i1++ )
263 for( size_t i2 = 0; i2 < N2; i2++ )
264 for( size_t i3 = 0; i3 < N3; i3++ ) {
265 if( i1+1 < N1 || periodic )
266 result.addEdge( i1*N2*N3 + i2*N3 + i3, ((i1+1)%N1)*N2*N3 + i2*N3 + i3, N1 <= 2 );
267 if( i2+1 < N2 || periodic )
268 result.addEdge( i1*N2*N3 + i2*N3 + i3, i1*N2*N3 + ((i2+1)%N2)*N3 + i3, N2 <= 2 );
269 if( i3+1 < N3 || periodic )
270 result.addEdge( i1*N2*N3 + i2*N3 + i3, i1*N2*N3 + i2*N3 + ((i3+1)%N3), N3 <= 2 );
271 }
272 return result;
273 }
274
275
276 GraphAL createGraphLoop( size_t N ) {
277 GraphAL result( N );
278 for( size_t i = 0; i < N; i++ )
279 result.addEdge( i, (i+1)%N, N <= 2 );
280 return result;
281 }
282
283
284 GraphAL createGraphTree( size_t N ) {
285 GraphAL result( N );
286 for( size_t i = 1; i < N; i++ ) {
287 size_t j = rnd_int( 0, i-1 );
288 result.addEdge( i, j, false );
289 }
290 return result;
291 }
292
293
294 GraphAL createGraphRegular( size_t N, size_t d ) {
295 DAI_ASSERT( (N * d) % 2 == 0 );
296 DAI_ASSERT( d < N );
297
298 GraphAL G( N );
299 if( d > 0 ) {
300 bool ready = false;
301 size_t tries = 0;
302 while( !ready ) {
303 tries++;
304
305 // Start with N*d points {0,1,...,N*d-1} (N*d even) in N groups.
306 // Put U = {0,1,...,N*d-1}. (U denotes the set of unpaired points.)
307 vector<size_t> U;
308 U.reserve( N * d );
309 for( size_t i = 0; i < N * d; i++ )
310 U.push_back( i );
311
312 // Repeat the following until no suitable pair can be found: Choose
313 // two random points i and j in U, and if they are suitable, pair
314 // i with j and delete i and j from U.
315 G = GraphAL( N );
316 bool finished = false;
317 while( !finished ) {
318 random_shuffle( U.begin(), U.end() );
319 size_t i1, i2;
320 bool suit_pair_found = false;
321 for( i1 = 0; i1 < U.size()-1 && !suit_pair_found; i1++ )
322 for( i2 = i1+1; i2 < U.size() && !suit_pair_found; i2++ )
323 if( ((U[i1] / d) != (U[i2] / d)) && !G.hasEdge( U[i1] / d, U[i2] / d ) ) {
324 // they are suitable
325 suit_pair_found = true;
326 G.addEdge( U[i1] / d, U[i2] / d, false );
327 U.erase( U.begin() + i2 ); // first remove largest
328 U.erase( U.begin() + i1 ); // remove smallest
329 }
330 if( !suit_pair_found || U.empty() )
331 finished = true;
332 }
333
334 if( U.empty() ) {
335 // G is a graph with edge from vertex r to vertex s if and only if
336 // there is a pair containing points in the r'th and s'th groups.
337 // If G is d-regular, output, otherwise return to Step 1.
338 ready = true;
339 for( size_t n = 0; n < N; n++ )
340 if( G.nb(n).size() != d ) {
341 ready = false;
342 break;
343 }
344 } else
345 ready = false;
346 }
347 }
348
349 return G;
350 }
351
352
353 } // end of namespace dai