Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGr_EdgeBetweennessCentrality: edge betweenness-centrality
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-Licene-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICEnE file or contact permission@sei.cmu.edu. See
10 : // Contributors.txt for a full list of contributors. Created, in part, with
11 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
12 : // DM22-0790
13 :
14 : // Contributed by Casey Pei and Tim Davis, Texas A&M University;
15 : // Adapted and revised from GraphBLAS C API Spec, Appendix B.4.
16 :
17 : //------------------------------------------------------------------------------
18 :
19 : // LAGr_EdgeBetweennessCentrality: Exact algorithm for computing
20 : // betweeness centrality.
21 :
22 : // This is an Advanced algorithm (no self edges allowed)
23 :
24 : //------------------------------------------------------------------------------
25 :
26 : #define useAssign
27 : // #define debug
28 :
29 : #define LG_FREE_WORK \
30 : { \
31 : GrB_free (&frontier) ; \
32 : GrB_free (&J_vec) ; \
33 : GrB_free (&I_vec) ; \
34 : GrB_free (&J_matrix) ; \
35 : GrB_free (&I_matrix) ; \
36 : GrB_free (&Fd1A) ; \
37 : GrB_free (&paths) ; \
38 : GrB_free (&bc_vertex_flow) ; \
39 : GrB_free (&temp_update) ; \
40 : GrB_free (&Add_One_Divide) ; \
41 : GrB_free (&Update) ; \
42 : GrB_free (&HalfUpdate) ; \
43 : GrB_free (&HalfUpdateT) ; \
44 : GrB_free (&SymmetricUpdate) ; \
45 : GrB_free (&internal_sources) ; \
46 : if (Search != NULL) \
47 : { \
48 : for (int64_t i = 0 ; i <= n ; i++) \
49 : { \
50 : GrB_free (&(Search [i])) ; \
51 : } \
52 : LAGraph_Free ((void **) &Search, NULL) ; \
53 : } \
54 : }
55 :
56 : #define LG_FREE_ALL \
57 : { \
58 : LG_FREE_WORK ; \
59 : GrB_free (centrality) ; \
60 : }
61 :
62 : #include "LG_internal.h"
63 : #include <LAGraphX.h>
64 :
65 : //------------------------------------------------------------------------------
66 : // (1+x)/y function for double: z = (1 + x) / y
67 : //------------------------------------------------------------------------------
68 :
69 24048 : void add_one_divide_function (double *z, const double *x, const double *y)
70 : {
71 24048 : double a = (*(x)) ;
72 24048 : double b = (*(y)) ;
73 24048 : (*(z)) = (1 + a) / b ;
74 24048 : }
75 :
76 : #define ADD_ONE_DIVIDE_FUNCTION_DEFN \
77 : "void add_one_divide_function (double *z, const double *x, const double *y)\n" \
78 : "{ \n" \
79 : " double a = (*(x)) ; \n" \
80 : " double b = (*(y)) ; \n" \
81 : " (*(z)) = (1 + a) / b ; \n" \
82 : "}"
83 :
84 : //------------------------------------------------------------------------------
85 : // LAGr_EdgeBetweennessCentrality: edge betweenness-centrality
86 : //------------------------------------------------------------------------------
87 :
88 23 : int LAGr_EdgeBetweennessCentrality
89 : (
90 : // output:
91 : GrB_Matrix *centrality, // centrality(i): betweeness centrality of i
92 : // input:
93 : LAGraph_Graph G, // input graph
94 : GrB_Vector sources, // source vertices to compute shortest paths (if NULL or empty, use all vertices)
95 : char *msg
96 : )
97 : {
98 :
99 : #if LAGRAPH_SUITESPARSE
100 :
101 : //--------------------------------------------------------------------------
102 : // check inputs
103 : //--------------------------------------------------------------------------
104 :
105 23 : LG_CLEAR_MSG ;
106 :
107 : // Array of BFS search matrices.
108 : // Search[i] is a sparse matrix that stores the depth at which each vertex is
109 : // first seen thus far in each BFS at the current depth i. Each column
110 : // corresponds to a BFS traversal starting from a source node.
111 23 : GrB_Vector *Search = NULL ;
112 :
113 : // Frontier vector, a sparse matrix.
114 : // Stores # of shortest paths to vertices at current BFS depth
115 23 : GrB_Vector frontier = NULL ;
116 :
117 : // Paths matrix holds the number of shortest paths for each node and
118 : // starting node discovered so far. A dense vector that is updated with
119 : // sparse updates, and also used as a mask.
120 23 : GrB_Vector paths = NULL ;
121 :
122 : // The betweenness centrality for each vertex. A dense vector that
123 : // accumulates flow values during backtracking.
124 23 : GrB_Vector bc_vertex_flow = NULL ;
125 :
126 : // Update matrix for betweenness centrality for each edge. A sparse matrix
127 : // that holds intermediate centrality updates.
128 23 : GrB_Matrix Update = NULL ;
129 :
130 : // Binary operator for computing (1+x)/y in centrality calculations
131 23 : GrB_BinaryOp Add_One_Divide = NULL ;
132 :
133 : // Temporary vectors and matrices for intermediate calculations
134 : // Diagonal values for J_matrix
135 23 : GrB_Vector J_vec = NULL ;
136 :
137 : // Diagonal values for I_matrix
138 23 : GrB_Vector I_vec = NULL ;
139 :
140 : // Matrix for previous level contributions
141 23 : GrB_Matrix I_matrix = NULL ;
142 :
143 : // Matrix for current level contributions
144 23 : GrB_Matrix J_matrix = NULL ;
145 :
146 : // Intermediate product matrix
147 23 : GrB_Matrix Fd1A = NULL ;
148 :
149 : // Temporary vector for centrality updates
150 23 : GrB_Vector temp_update = NULL ;
151 :
152 : // Temporary matrices for doing updates on
153 : // approximate and undirected graphs
154 23 : GrB_Matrix HalfUpdate = NULL ;
155 23 : GrB_Matrix HalfUpdateT = NULL ;
156 23 : GrB_Matrix SymmetricUpdate = NULL ;
157 :
158 : // Source nodes vector (will be created if NULL is passed)
159 23 : GrB_Vector internal_sources = NULL;
160 23 : bool created_sources = false;
161 :
162 23 : GrB_Index n = 0 ; // # nodes in the graph
163 :
164 23 : double t1_total = 0;
165 23 : double t2_total = 0;
166 23 : double t3_total = 0;
167 :
168 23 : LG_ASSERT (centrality != NULL, GrB_NULL_POINTER) ;
169 23 : (*centrality) = NULL ;
170 23 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
171 :
172 23 : GrB_Matrix A = G->A ;
173 : #if 0
174 : GrB_Matrix AT ;
175 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
176 : G->is_symmetric_structure == LAGraph_TRUE)
177 : {
178 : // A and A' have the same structure
179 : AT = A ;
180 : }
181 : else
182 : {
183 : // A and A' differ
184 : AT = G->AT ;
185 : LG_ASSERT_MSG (AT != NULL, LAGRAPH_NOT_CACHED, "G->AT is required") ;
186 : }
187 : #endif
188 :
189 23 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
190 : GrB_Index nsources ;
191 23 : if (sources == NULL)
192 : {
193 8 : nsources = n ;
194 : }
195 : else
196 : {
197 15 : GRB_TRY (GrB_Vector_nvals (&nsources, sources)) ;
198 : }
199 23 : LG_ASSERT (nsources > 0, GrB_INVALID_VALUE) ;
200 :
201 : // =========================================================================
202 : // === initialization =====================================================
203 : // =========================================================================
204 :
205 22 : GRB_TRY (GxB_BinaryOp_new (&Add_One_Divide,
206 : (GxB_binary_function) add_one_divide_function,
207 : GrB_FP64, GrB_FP64, GrB_FP64,
208 : "add_one_divide_function", ADD_ONE_DIVIDE_FUNCTION_DEFN)) ;
209 :
210 : // Initialize the frontier, paths, Update, and bc_vertex_flow
211 22 : GRB_TRY (GrB_Vector_new (&paths, GrB_FP64, n)) ;
212 22 : GRB_TRY (GrB_Vector_new (&frontier, GrB_FP64, n)) ;
213 22 : GRB_TRY (GrB_Matrix_new (&Update, GrB_FP64, n, n)) ;
214 22 : GRB_TRY (GrB_Vector_new (&bc_vertex_flow, GrB_FP64, n)) ;
215 :
216 :
217 : // Initialize centrality matrix with zeros using A as structural mask
218 22 : LG_TRY (GrB_Matrix_new(centrality, GrB_FP64, n, n)) ;
219 22 : GRB_TRY (GrB_assign (*centrality, A, NULL, 0.0, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S)) ;
220 :
221 : // Allocate memory for the array of S vectors
222 22 : LG_TRY (LAGraph_Calloc ((void **) &Search, n + 1, sizeof (GrB_Vector), msg)) ;
223 :
224 : // =========================================================================
225 : // === Process source nodes ================================================
226 : // =========================================================================
227 :
228 : // If sources is NULL, create a dense vector with all vertices
229 22 : if (sources == NULL)
230 : {
231 : // Create a vector with all nodes as sources
232 8 : GRB_TRY (GrB_Vector_new (&internal_sources, GrB_INT64, n)) ;
233 :
234 : // internal_sources (0:n-1) = 0
235 8 : GRB_TRY (GrB_assign (internal_sources, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
236 :
237 : // internal_sources (0:n-1) = 0:n-1
238 8 : GRB_TRY (GrB_apply (internal_sources, NULL, NULL, GrB_ROWINDEX_INT64,
239 : internal_sources, 0, NULL)) ;
240 :
241 : /*
242 : int64_t ns = n;
243 : for (GrB_Index i = 0; i < ns; i++)
244 : {
245 : GRB_TRY (GrB_Vector_setElement_INT64 (internal_sources, i, i)) ;
246 : }
247 : */
248 :
249 : // Use this vector instead
250 8 : sources = internal_sources;
251 8 : created_sources = true;
252 : }
253 :
254 : // =========================================================================
255 : // === Breadth-first search stage ==========================================
256 : // =========================================================================
257 :
258 22 : GrB_Index frontier_size, last_frontier_size = 0 ;
259 22 : GRB_TRY (GrB_Vector_nvals (&frontier_size, frontier)) ;
260 :
261 : int64_t depth;
262 : GrB_Index root;
263 :
264 22 : GRB_TRY (GrB_Vector_new(&J_vec, GrB_FP64, n)) ;
265 22 : GRB_TRY (GrB_Vector_new (&I_vec, GrB_FP64, n)) ;
266 22 : GRB_TRY (GrB_Matrix_new (&Fd1A, GrB_FP64, n, n)) ;
267 22 : GRB_TRY (GrB_Vector_new(&temp_update, GrB_FP64, n)) ; // Create a temporary vector
268 :
269 22 : GRB_TRY (GrB_Matrix_new(&HalfUpdate, GrB_FP64, n, n)) ;
270 22 : GRB_TRY (GrB_Matrix_new(&HalfUpdateT, GrB_FP64, n, n)) ;
271 22 : GRB_TRY (GrB_Matrix_new(&SymmetricUpdate, GrB_FP64, n, n)) ;
272 :
273 : // Iterate through source nodes
274 3180 : for (GrB_Index i = 0; i < nsources; i++)
275 : {
276 3158 : GRB_TRY (GrB_Vector_extractElement(&root, sources, i)) ;
277 :
278 : // Verify the root index is valid
279 3158 : LG_ASSERT (root < n, GrB_INVALID_VALUE) ;
280 :
281 3158 : depth = 0 ;
282 :
283 : // root frontier: Search [0](root) = true
284 3158 : GrB_free (&(Search [0])) ;
285 3158 : GRB_TRY (GrB_Vector_new(&(Search [0]), GrB_BOOL, n)) ;
286 3158 : GRB_TRY (GrB_Vector_setElement_BOOL(Search [0], (bool) true, root)) ;
287 :
288 : // clear paths, and then set paths (root) = 1
289 3158 : GRB_TRY (GrB_Vector_clear (paths)) ;
290 3158 : GRB_TRY (GrB_Vector_setElement (paths, (double) 1.0, root)) ;
291 :
292 3158 : GRB_TRY (GrB_Matrix_clear (Update)) ;
293 :
294 : // Extract row root from A into frontier vector: frontier = A(root,:)
295 3158 : GRB_TRY (GrB_Col_extract (frontier, NULL, NULL, A, GrB_ALL, n, root,
296 : GrB_DESC_T0)) ;
297 :
298 3158 : GRB_TRY (GrB_Vector_nvals (&frontier_size, frontier)) ;
299 3158 : GRB_TRY (GrB_assign (frontier, frontier, NULL, 1.0, GrB_ALL, n, GrB_DESC_S)) ;
300 :
301 129801 : while (frontier_size != 0)
302 : {
303 126643 : depth++ ;
304 :
305 : //----------------------------------------------------------------------
306 : // paths += frontier
307 : // Accumulate path counts for vertices at current depth
308 : //----------------------------------------------------------------------
309 :
310 126643 : GRB_TRY (GrB_assign (paths, NULL, GrB_PLUS_FP64, frontier, GrB_ALL, n,
311 : NULL)) ;
312 :
313 : //----------------------------------------------------------------------
314 : // Search[depth] = structure(frontier)
315 : // Record the frontier structure at current depth
316 : //----------------------------------------------------------------------
317 :
318 126643 : GrB_free (&(Search [depth])) ;
319 126643 : LG_TRY (LAGraph_Vector_Structure (&(Search [depth]), frontier, msg)) ;
320 :
321 : //----------------------------------------------------------------------
322 : // frontier<!paths> = frontier * A
323 : //----------------------------------------------------------------------
324 :
325 126643 : GRB_TRY (LG_SET_FORMAT_HINT (frontier, LG_SPARSE)) ;
326 126643 : GRB_TRY (GrB_vxm (frontier, paths, NULL, /* LAGraph_plus_first_fp64 */
327 : GxB_PLUS_FIRST_FP64, frontier,
328 : A, GrB_DESC_RSC )) ;
329 :
330 : //----------------------------------------------------------------------
331 : // Get size of current frontier: frontier_size = nvals(frontier)
332 : //----------------------------------------------------------------------
333 :
334 126643 : last_frontier_size = frontier_size ;
335 126643 : GRB_TRY (GrB_Vector_nvals (&frontier_size, frontier)) ;
336 : }
337 :
338 :
339 : // =========================================================================
340 : // === Betweenness centrality computation phase ============================
341 : // =========================================================================
342 :
343 : // bc_vertex_flow = ones (n, n) ; a full matrix (and stays full)
344 3158 : GRB_TRY (GrB_assign(bc_vertex_flow, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ;
345 :
346 3158 : GRB_TRY (GrB_Matrix_clear (HalfUpdate)) ;
347 3158 : GRB_TRY (GrB_Matrix_clear (HalfUpdateT)) ;
348 3158 : GRB_TRY (GrB_Matrix_clear (SymmetricUpdate)) ;
349 3158 : GRB_TRY (GrB_Matrix_clear (Fd1A)) ;
350 3158 : GRB_TRY (GrB_Vector_clear (J_vec)) ;
351 3158 : GRB_TRY (GrB_Vector_clear (I_vec)) ;
352 3158 : GRB_TRY (GrB_Vector_clear (temp_update)) ;
353 :
354 :
355 :
356 :
357 : // Backtrack through the BFS and compute centrality updates for each vertex
358 : // GrB_Index fd1_size;
359 :
360 129801 : while (depth >= 1)
361 : {
362 126643 : GrB_Vector f_d = Search [depth] ;
363 126643 : GrB_Vector f_d1 = Search [depth - 1] ;
364 :
365 : //----------------------------------------------------------------------
366 : // j<S(depth, :)> = (1 + v) / p
367 : // J = diag(j)
368 : // Compute weighted contributions from current level
369 : //----------------------------------------------------------------------
370 :
371 126643 : GRB_TRY (GrB_eWiseMult(J_vec, f_d, NULL, Add_One_Divide, bc_vertex_flow, paths, GrB_DESC_RS)) ;
372 126643 : GRB_TRY (GrB_Matrix_diag(&J_matrix, J_vec, 0)) ;
373 :
374 : //----------------------------------------------------------------------
375 : // i<S(depth-1, :)> = p
376 : // I = diag(i)
377 : // Compute weighted contributions from previous level
378 : //----------------------------------------------------------------------
379 :
380 126643 : GRB_TRY (GrB_Vector_extract (I_vec, f_d1, NULL, paths, GrB_ALL, n, GrB_DESC_RS)) ;
381 126643 : GRB_TRY (GrB_Matrix_diag(&I_matrix, I_vec, 0)) ;
382 :
383 : //----------------------------------------------------------------------
384 : // Update = I × A × J
385 : // Compute edge updates based on current level weights
386 : //----------------------------------------------------------------------
387 :
388 126643 : double t1 = LAGraph_WallClockTime();
389 126643 : GRB_TRY(GrB_mxm(Fd1A, NULL, NULL, LAGraph_plus_first_fp64,
390 : I_matrix, A, NULL));
391 126643 : t1 = LAGraph_WallClockTime() - t1;
392 126643 : t1_total += t1;
393 :
394 126643 : double t2 = LAGraph_WallClockTime();
395 126643 : GRB_TRY(GrB_mxm(Update, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64,
396 : Fd1A, J_matrix, NULL));
397 126643 : t2 = LAGraph_WallClockTime() - t2;
398 126643 : t2_total += t2;
399 126643 : GRB_TRY (GrB_free (&I_matrix)) ;
400 126643 : GRB_TRY (GrB_free (&J_matrix)) ;
401 : //----------------------------------------------------------------------
402 : // centrality<A> += Update
403 : // Accumulate centrality values for edges
404 : //----------------------------------------------------------------------
405 :
406 : #ifdef useAssign
407 : // centrality{A} += Update, using assign
408 126643 : double t3 = LAGraph_WallClockTime();
409 :
410 126643 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED) {
411 : // First divide the Update matrix by 2 for symmetric distribution
412 153 : GrB_apply(HalfUpdate, NULL, NULL, GrB_DIV_FP64, Update, 2.0, NULL);
413 :
414 : // Create a transposed version of the update
415 153 : GrB_transpose(HalfUpdateT, NULL, NULL, HalfUpdate, NULL);
416 :
417 : // Add the original and transposed matrices to create a symmetric update
418 153 : GrB_eWiseAdd(SymmetricUpdate, NULL, NULL, GrB_PLUS_FP64, HalfUpdate, HalfUpdateT, NULL);
419 :
420 : // Apply the symmetric update to the centrality
421 153 : GRB_TRY(GrB_assign(*centrality, A, GrB_PLUS_FP64, SymmetricUpdate, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S));
422 :
423 : }
424 : else {
425 126490 : GRB_TRY (GrB_assign(*centrality, A, GrB_PLUS_FP64, Update, GrB_ALL, n, GrB_ALL, n,
426 : GrB_DESC_S));
427 : }
428 :
429 126643 : t3 = LAGraph_WallClockTime() - t3;
430 126643 : t3_total += t3;
431 : #else
432 : // TODO: Approx update using ewise add not implemented
433 : // centrality = centrality + Update using eWiseAdd
434 : double t3 = LAGraph_WallClockTime();
435 : GRB_TRY (GrB_eWiseAdd (*centrality, NULL, NULL, GrB_PLUS_FP64, *centrality, Update, NULL));
436 : t3 = LAGraph_WallClockTime() - t3;
437 : t3_total += t3;
438 : #endif
439 :
440 :
441 : //----------------------------------------------------------------------
442 : // v = Update +.
443 : // Reduce update matrix to vector for next iteration
444 : //----------------------------------------------------------------------
445 :
446 126643 : GRB_TRY (GrB_reduce(temp_update, NULL, NULL, GrB_PLUS_MONOID_FP64, Update, NULL)) ;
447 126643 : GRB_TRY (GrB_eWiseAdd(bc_vertex_flow, NULL, NULL, GrB_PLUS_FP64, bc_vertex_flow, temp_update, NULL)) ;
448 :
449 : // 24 d = d − 1
450 126643 : depth-- ;
451 : }
452 :
453 : }
454 :
455 : #ifdef debug
456 : printf(" I*A time: %g\n", t1_total);
457 :
458 : printf(" (I*A)*J time: %g\n", t2_total);
459 :
460 : #ifdef useAssign
461 : printf(" Centrality update using assign time: %g\n", t3_total);
462 : #else
463 : printf(" Centrality update using eWiseAdd time: %g\n", t3_total);
464 : #endif
465 :
466 : GxB_print (*centrality, GxB_FULL) ;
467 :
468 : #endif
469 :
470 :
471 : // =========================================================================
472 : // === finalize the centrality =============================================
473 : // =========================================================================
474 :
475 9164 : LG_FREE_WORK ;
476 22 : return (GrB_SUCCESS) ;
477 : #else
478 : return (GrB_NOT_IMPLEMENTED) ;
479 : #endif
480 : }
|
