1 | // Copyright 2004 The Trustees of Indiana University. |
2 | |
3 | // Use, modification and distribution is subject to the Boost Software |
4 | // License, Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at |
5 | // http://www.boost.org/LICENSE_1_0.txt) |
6 | |
7 | // Authors: Douglas Gregor |
8 | // Andrew Lumsdaine |
9 | #include <boost/graph/betweenness_centrality.hpp> |
10 | |
11 | #include <boost/graph/adjacency_list.hpp> |
12 | #include <vector> |
13 | #include <stack> |
14 | #include <queue> |
15 | #include <boost/property_map/property_map.hpp> |
16 | #include <boost/core/lightweight_test.hpp> |
17 | #include <boost/random/uniform_01.hpp> |
18 | #include <boost/random/linear_congruential.hpp> |
19 | #include <boost/lexical_cast.hpp> |
20 | |
21 | using namespace boost; |
22 | |
23 | const double error_tolerance = 0.001; |
24 | |
25 | typedef property< edge_weight_t, double, property< edge_index_t, std::size_t > > |
26 | EdgeProperties; |
27 | |
28 | struct weighted_edge |
29 | { |
30 | int source, target; |
31 | double weight; |
32 | }; |
33 | |
34 | template < typename Graph > |
35 | void run_weighted_test(Graph*, int V, weighted_edge edge_init[], int E, |
36 | double correct_centrality[]) |
37 | { |
38 | Graph g(V); |
39 | typedef typename graph_traits< Graph >::vertex_descriptor Vertex; |
40 | typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator; |
41 | typedef typename graph_traits< Graph >::edge_descriptor Edge; |
42 | |
43 | std::vector< Vertex > vertices(V); |
44 | { |
45 | vertex_iterator v, v_end; |
46 | int index = 0; |
47 | for (boost::tie(v, v_end) = boost::vertices(g); v != v_end; |
48 | ++v, ++index) |
49 | { |
50 | put(vertex_index, g, *v, index); |
51 | vertices[index] = *v; |
52 | } |
53 | } |
54 | |
55 | std::vector< Edge > edges(E); |
56 | for (int e = 0; e < E; ++e) |
57 | { |
58 | edges[e] = add_edge( |
59 | vertices[edge_init[e].source], vertices[edge_init[e].target], g) |
60 | .first; |
61 | put(edge_weight, g, edges[e], 1.0); |
62 | } |
63 | |
64 | std::vector< double > centrality(V); |
65 | brandes_betweenness_centrality(g, |
66 | centrality_map(make_iterator_property_map( |
67 | centrality.begin(), get(vertex_index, g), double())) |
68 | .vertex_index_map(get(vertex_index, g)) |
69 | .weight_map(get(edge_weight, g))); |
70 | |
71 | for (int v = 0; v < V; ++v) |
72 | { |
73 | BOOST_TEST(centrality[v] == correct_centrality[v]); |
74 | } |
75 | } |
76 | |
77 | struct unweighted_edge |
78 | { |
79 | int source, target; |
80 | }; |
81 | |
82 | template < typename Graph > |
83 | void run_unweighted_test(Graph*, int V, unweighted_edge edge_init[], int E, |
84 | double correct_centrality[], double* correct_edge_centrality = 0) |
85 | { |
86 | Graph g(V); |
87 | typedef typename graph_traits< Graph >::vertex_descriptor Vertex; |
88 | typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator; |
89 | typedef typename graph_traits< Graph >::edge_descriptor Edge; |
90 | |
91 | std::vector< Vertex > vertices(V); |
92 | { |
93 | vertex_iterator v, v_end; |
94 | int index = 0; |
95 | for (boost::tie(v, v_end) = boost::vertices(g); v != v_end; |
96 | ++v, ++index) |
97 | { |
98 | put(vertex_index, g, *v, index); |
99 | vertices[index] = *v; |
100 | } |
101 | } |
102 | |
103 | std::vector< Edge > edges(E); |
104 | for (int e = 0; e < E; ++e) |
105 | { |
106 | edges[e] = add_edge( |
107 | vertices[edge_init[e].source], vertices[edge_init[e].target], g) |
108 | .first; |
109 | put(edge_weight, g, edges[e], 1.0); |
110 | put(edge_index, g, edges[e], e); |
111 | } |
112 | |
113 | std::vector< double > centrality(V); |
114 | std::vector< double > edge_centrality1(E); |
115 | |
116 | brandes_betweenness_centrality(g, |
117 | centrality_map(make_iterator_property_map( |
118 | centrality.begin(), get(vertex_index, g), double())) |
119 | .edge_centrality_map(make_iterator_property_map( |
120 | edge_centrality1.begin(), get(edge_index, g), double())) |
121 | .vertex_index_map(get(vertex_index, g))); |
122 | |
123 | std::vector< double > centrality2(V); |
124 | std::vector< double > edge_centrality2(E); |
125 | brandes_betweenness_centrality(g, |
126 | vertex_index_map(get(vertex_index, g)) |
127 | .weight_map(get(edge_weight, g)) |
128 | .centrality_map(make_iterator_property_map( |
129 | centrality2.begin(), get(vertex_index, g), double())) |
130 | .edge_centrality_map(make_iterator_property_map( |
131 | edge_centrality2.begin(), get(edge_index, g), double()))); |
132 | |
133 | std::vector< double > edge_centrality3(E); |
134 | brandes_betweenness_centrality(g, |
135 | edge_centrality_map(make_iterator_property_map( |
136 | edge_centrality3.begin(), get(edge_index, g), double()))); |
137 | |
138 | for (int v = 0; v < V; ++v) |
139 | { |
140 | BOOST_TEST(centrality[v] == centrality2[v]); |
141 | |
142 | double relative_error = correct_centrality[v] == 0.0 |
143 | ? centrality[v] |
144 | : (centrality[v] - correct_centrality[v]) / correct_centrality[v]; |
145 | if (relative_error < 0) |
146 | relative_error = -relative_error; |
147 | BOOST_TEST(relative_error < error_tolerance); |
148 | } |
149 | |
150 | for (int e = 0; e < E; ++e) |
151 | { |
152 | BOOST_TEST(edge_centrality1[e] == edge_centrality2[e]); |
153 | BOOST_TEST(edge_centrality1[e] == edge_centrality3[e]); |
154 | |
155 | if (correct_edge_centrality) |
156 | { |
157 | double relative_error = correct_edge_centrality[e] == 0.0 |
158 | ? edge_centrality1[e] |
159 | : (edge_centrality1[e] - correct_edge_centrality[e]) |
160 | / correct_edge_centrality[e]; |
161 | if (relative_error < 0) |
162 | relative_error = -relative_error; |
163 | BOOST_TEST(relative_error < error_tolerance); |
164 | |
165 | if (relative_error >= error_tolerance) |
166 | { |
167 | std::cerr << "Edge " << e << " has edge centrality " |
168 | << edge_centrality1[e] << ", should be " |
169 | << correct_edge_centrality[e] << std::endl; |
170 | } |
171 | } |
172 | } |
173 | } |
174 | |
175 | template < typename Graph > void run_wheel_test(Graph*, int V) |
176 | { |
177 | typedef typename graph_traits< Graph >::vertex_descriptor Vertex; |
178 | typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator; |
179 | typedef typename graph_traits< Graph >::edge_descriptor Edge; |
180 | |
181 | Graph g(V); |
182 | Vertex center = *boost::vertices(g).first; |
183 | |
184 | std::vector< Vertex > vertices(V); |
185 | { |
186 | vertex_iterator v, v_end; |
187 | int index = 0; |
188 | for (boost::tie(v, v_end) = boost::vertices(g); v != v_end; |
189 | ++v, ++index) |
190 | { |
191 | put(vertex_index, g, *v, index); |
192 | vertices[index] = *v; |
193 | if (*v != center) |
194 | { |
195 | Edge e = add_edge(*v, center, g).first; |
196 | put(edge_weight, g, e, 1.0); |
197 | } |
198 | } |
199 | } |
200 | |
201 | std::vector< double > centrality(V); |
202 | brandes_betweenness_centrality(g, |
203 | make_iterator_property_map( |
204 | centrality.begin(), get(vertex_index, g), double())); |
205 | |
206 | std::vector< double > centrality2(V); |
207 | brandes_betweenness_centrality(g, |
208 | centrality_map(make_iterator_property_map( |
209 | centrality2.begin(), get(vertex_index, g), double())) |
210 | .vertex_index_map(get(vertex_index, g)) |
211 | .weight_map(get(edge_weight, g))); |
212 | |
213 | relative_betweenness_centrality(g, |
214 | make_iterator_property_map( |
215 | centrality.begin(), get(vertex_index, g), double())); |
216 | |
217 | relative_betweenness_centrality(g, |
218 | make_iterator_property_map( |
219 | centrality2.begin(), get(vertex_index, g), double())); |
220 | |
221 | for (int v = 0; v < V; ++v) |
222 | { |
223 | BOOST_TEST(centrality[v] == centrality2[v]); |
224 | BOOST_TEST( |
225 | (v == 0 && centrality[v] == 1) || (v != 0 && centrality[v] == 0)); |
226 | } |
227 | |
228 | double dominance = central_point_dominance(g, |
229 | make_iterator_property_map( |
230 | centrality2.begin(), get(vertex_index, g), double())); |
231 | BOOST_TEST(dominance == 1.0); |
232 | } |
233 | |
234 | template < typename MutableGraph > |
235 | void randomly_add_edges(MutableGraph& g, double edge_probability) |
236 | { |
237 | typedef typename graph_traits< MutableGraph >::directed_category |
238 | directed_category; |
239 | |
240 | minstd_rand gen; |
241 | uniform_01< minstd_rand, double > rand_gen(gen); |
242 | |
243 | typedef typename graph_traits< MutableGraph >::vertex_descriptor vertex; |
244 | typename graph_traits< MutableGraph >::vertex_iterator vi, vi_end; |
245 | for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) |
246 | { |
247 | vertex v = *vi; |
248 | typename graph_traits< MutableGraph >::vertex_iterator wi |
249 | = is_same< directed_category, undirected_tag >::value |
250 | ? vi |
251 | : vertices(g).first; |
252 | while (wi != vi_end) |
253 | { |
254 | vertex w = *wi++; |
255 | if (v != w) |
256 | { |
257 | if (rand_gen() < edge_probability) |
258 | add_edge(v, w, g); |
259 | } |
260 | } |
261 | } |
262 | } |
263 | |
264 | template < typename Graph, typename VertexIndexMap, typename CentralityMap > |
265 | void simple_unweighted_betweenness_centrality( |
266 | const Graph& g, VertexIndexMap index, CentralityMap centrality) |
267 | { |
268 | typedef typename boost::graph_traits< Graph >::vertex_descriptor vertex; |
269 | typedef |
270 | typename boost::graph_traits< Graph >::vertex_iterator vertex_iterator; |
271 | typedef typename boost::graph_traits< Graph >::adjacency_iterator |
272 | adjacency_iterator; |
273 | typedef typename boost::graph_traits< Graph >::vertices_size_type |
274 | vertices_size_type; |
275 | typedef typename boost::property_traits< CentralityMap >::value_type |
276 | centrality_type; |
277 | |
278 | vertex_iterator vi, vi_end; |
279 | for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) |
280 | put(centrality, *vi, 0); |
281 | |
282 | vertex_iterator si, si_end; |
283 | for (boost::tie(si, si_end) = vertices(g); si != si_end; ++si) |
284 | { |
285 | vertex s = *si; |
286 | |
287 | // S <-- empty stack |
288 | std::stack< vertex > S; |
289 | |
290 | // P[w] <-- empty list, w \in V |
291 | typedef std::vector< vertex > Predecessors; |
292 | std::vector< Predecessors > predecessors(num_vertices(g)); |
293 | |
294 | // sigma[t] <-- 0, t \in V |
295 | std::vector< vertices_size_type > sigma(num_vertices(g), 0); |
296 | |
297 | // sigma[s] <-- 1 |
298 | sigma[get(index, s)] = 1; |
299 | |
300 | // d[t] <-- -1, t \in V |
301 | std::vector< int > d(num_vertices(g), -1); |
302 | |
303 | // d[s] <-- 0 |
304 | d[get(index, s)] = 0; |
305 | |
306 | // Q <-- empty queue |
307 | std::queue< vertex > Q; |
308 | |
309 | // enqueue s --> Q |
310 | Q.push(s); |
311 | |
312 | while (!Q.empty()) |
313 | { |
314 | // dequeue v <-- Q |
315 | vertex v = Q.front(); |
316 | Q.pop(); |
317 | |
318 | // push v --> S |
319 | S.push(v); |
320 | |
321 | adjacency_iterator wi, wi_end; |
322 | for (boost::tie(wi, wi_end) = adjacent_vertices(v, g); wi != wi_end; |
323 | ++wi) |
324 | { |
325 | vertex w = *wi; |
326 | |
327 | // w found for the first time? |
328 | if (d[get(index, w)] < 0) |
329 | { |
330 | // enqueue w --> Q |
331 | Q.push(w); |
332 | |
333 | // d[w] <-- d[v] + 1 |
334 | d[get(index, w)] = d[get(index, v)] + 1; |
335 | } |
336 | |
337 | // shortest path to w via v? |
338 | if (d[get(index, w)] == d[get(index, v)] + 1) |
339 | { |
340 | // sigma[w] = sigma[w] + sigma[v] |
341 | sigma[get(index, w)] += sigma[get(index, v)]; |
342 | |
343 | // append v --> P[w] |
344 | predecessors[get(index, w)].push_back(v); |
345 | } |
346 | } |
347 | } |
348 | |
349 | // delta[v] <-- 0, v \in V |
350 | std::vector< centrality_type > delta(num_vertices(g), 0); |
351 | |
352 | // S returns vertices in order of non-increasing distance from s |
353 | while (!S.empty()) |
354 | { |
355 | // pop w <-- S |
356 | vertex w = S.top(); |
357 | S.pop(); |
358 | |
359 | const Predecessors& w_preds = predecessors[get(index, w)]; |
360 | for (typename Predecessors::const_iterator vi = w_preds.begin(); |
361 | vi != w_preds.end(); ++vi) |
362 | { |
363 | vertex v = *vi; |
364 | // delta[v] <-- delta[v] + (sigma[v]/sigma[w])*(1 + delta[w]) |
365 | delta[get(index, v)] += ((centrality_type)sigma[get(index, v)] |
366 | / sigma[get(index, w)]) |
367 | * (1 + delta[get(index, w)]); |
368 | } |
369 | |
370 | if (w != s) |
371 | { |
372 | // C_B[w] <-- C_B[w] + delta[w] |
373 | centrality[w] += delta[get(index, w)]; |
374 | } |
375 | } |
376 | } |
377 | |
378 | typedef typename graph_traits< Graph >::directed_category directed_category; |
379 | const bool is_undirected |
380 | = is_same< directed_category, undirected_tag >::value; |
381 | if (is_undirected) |
382 | { |
383 | vertex_iterator v, v_end; |
384 | for (boost::tie(v, v_end) = vertices(g); v != v_end; ++v) |
385 | { |
386 | put(centrality, *v, get(centrality, *v) / centrality_type(2)); |
387 | } |
388 | } |
389 | } |
390 | |
391 | template < typename Graph > void random_unweighted_test(Graph*, int n) |
392 | { |
393 | Graph g(n); |
394 | |
395 | { |
396 | typename graph_traits< Graph >::vertex_iterator v, v_end; |
397 | int index = 0; |
398 | for (boost::tie(v, v_end) = boost::vertices(g); v != v_end; |
399 | ++v, ++index) |
400 | { |
401 | put(vertex_index, g, *v, index); |
402 | } |
403 | } |
404 | |
405 | randomly_add_edges(g, 0.20); |
406 | |
407 | std::cout << "Random graph with " << n << " vertices and " << num_edges(g) |
408 | << " edges.\n" ; |
409 | |
410 | std::cout << " Direct translation of Brandes' algorithm..." ; |
411 | std::vector< double > centrality(n); |
412 | simple_unweighted_betweenness_centrality(g, get(vertex_index, g), |
413 | make_iterator_property_map( |
414 | centrality.begin(), get(vertex_index, g), double())); |
415 | std::cout << "DONE.\n" ; |
416 | |
417 | std::cout << " Real version, unweighted..." ; |
418 | std::vector< double > centrality2(n); |
419 | brandes_betweenness_centrality(g, |
420 | make_iterator_property_map( |
421 | centrality2.begin(), get(vertex_index, g), double())); |
422 | std::cout << "DONE.\n" ; |
423 | |
424 | if (!std::equal(first1: centrality.begin(), last1: centrality.end(), first2: centrality2.begin())) |
425 | { |
426 | for (std::size_t v = 0; v < centrality.size(); ++v) |
427 | { |
428 | double relative_error = centrality[v] == 0.0 |
429 | ? centrality2[v] |
430 | : (centrality2[v] - centrality[v]) / centrality[v]; |
431 | if (relative_error < 0) |
432 | relative_error = -relative_error; |
433 | BOOST_TEST(relative_error < error_tolerance); |
434 | } |
435 | } |
436 | |
437 | std::cout << " Real version, weighted..." ; |
438 | std::vector< double > centrality3(n); |
439 | |
440 | for (typename graph_traits< Graph >::edge_iterator ei = edges(g).first; |
441 | ei != edges(g).second; ++ei) |
442 | put(edge_weight, g, *ei, 1); |
443 | |
444 | brandes_betweenness_centrality(g, |
445 | weight_map(get(edge_weight, g)) |
446 | .centrality_map(make_iterator_property_map( |
447 | centrality3.begin(), get(vertex_index, g), double()))); |
448 | std::cout << "DONE.\n" ; |
449 | |
450 | if (!std::equal(first1: centrality.begin(), last1: centrality.end(), first2: centrality3.begin())) |
451 | { |
452 | for (std::size_t v = 0; v < centrality.size(); ++v) |
453 | { |
454 | double relative_error = centrality[v] == 0.0 |
455 | ? centrality3[v] |
456 | : (centrality3[v] - centrality[v]) / centrality[v]; |
457 | if (relative_error < 0) |
458 | relative_error = -relative_error; |
459 | BOOST_TEST(relative_error < error_tolerance); |
460 | } |
461 | } |
462 | } |
463 | |
464 | int main(int argc, char* argv[]) |
465 | { |
466 | int random_test_num_vertices = 300; |
467 | if (argc >= 2) |
468 | random_test_num_vertices = boost::lexical_cast< int >(arg: argv[1]); |
469 | typedef adjacency_list< listS, listS, undirectedS, |
470 | property< vertex_index_t, int >, EdgeProperties > |
471 | Graph; |
472 | typedef adjacency_list< listS, listS, directedS, |
473 | property< vertex_index_t, int >, EdgeProperties > |
474 | Digraph; |
475 | |
476 | struct unweighted_edge ud_edge_init1[5] |
477 | = { { .source: 0, .target: 1 }, { .source: 0, .target: 3 }, { .source: 1, .target: 2 }, { .source: 3, .target: 2 }, { .source: 2, .target: 4 } }; |
478 | double ud_centrality1[5] = { 0.5, 1.0, 3.5, 1.0, 0.0 }; |
479 | run_unweighted_test((Graph*)0, V: 5, edge_init: ud_edge_init1, E: 5, correct_centrality: ud_centrality1); |
480 | |
481 | // Example borrowed from the JUNG test suite |
482 | struct unweighted_edge ud_edge_init2[10] = { |
483 | { .source: 0, .target: 1 }, |
484 | { .source: 0, .target: 6 }, |
485 | { .source: 1, .target: 2 }, |
486 | { .source: 1, .target: 3 }, |
487 | { .source: 2, .target: 4 }, |
488 | { .source: 3, .target: 4 }, |
489 | { .source: 4, .target: 5 }, |
490 | { .source: 5, .target: 8 }, |
491 | { .source: 7, .target: 8 }, |
492 | { .source: 6, .target: 7 }, |
493 | }; |
494 | double ud_centrality2[9] |
495 | = { 0.2142 * 28, 0.2797 * 28, 0.0892 * 28, 0.0892 * 28, 0.2797 * 28, |
496 | 0.2142 * 28, 0.1666 * 28, 0.1428 * 28, 0.1666 * 28 }; |
497 | double ud_edge_centrality2[10] = { 10.66666, 9.33333, 6.5, 6.5, 6.5, 6.5, |
498 | 10.66666, 9.33333, 8.0, 8.0 }; |
499 | |
500 | run_unweighted_test( |
501 | (Graph*)0, V: 9, edge_init: ud_edge_init2, E: 10, correct_centrality: ud_centrality2, correct_edge_centrality: ud_edge_centrality2); |
502 | |
503 | weighted_edge dw_edge_init1[6] = { { .source: 0, .target: 1, .weight: 1.0 }, { .source: 0, .target: 3, .weight: 1.0 }, |
504 | { .source: 1, .target: 2, .weight: 0.5 }, { .source: 3, .target: 1, .weight: 1.0 }, { .source: 3, .target: 4, .weight: 1.0 }, { .source: 4, .target: 2, .weight: 0.5 } }; |
505 | double dw_centrality1[5] = { 0.0, 1.5, 0.0, 1.0, 0.5 }; |
506 | run_weighted_test((Digraph*)0, V: 5, edge_init: dw_edge_init1, E: 6, correct_centrality: dw_centrality1); |
507 | |
508 | run_wheel_test((Graph*)0, V: 15); |
509 | |
510 | random_unweighted_test((Graph*)0, n: random_test_num_vertices); |
511 | |
512 | return boost::report_errors(); |
513 | } |
514 | |