Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LG_check_edgeBetweennessCentrality: reference implementation for edge
3 : // betweenness centrality
4 : //------------------------------------------------------------------------------
5 :
6 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
7 : // SPDX-License-Identifier: BSD-2-Clause
8 : //
9 : // For additional details (including references to third party source code and
10 : // other files) see the LICENSE file or contact permission@sei.cmu.edu. See
11 : // Contributors.txt for a full list of contributors. Created, in part, with
12 : // funding and support from the U.S. Government (see Acknowledgments.txt file).
13 : // DM22-0790
14 :
15 : // Contributed by Casey Pei, Texas A&M University
16 :
17 : //------------------------------------------------------------------------------
18 :
19 : #define LG_FREE_WORK \
20 : { \
21 : LAGraph_Free ((void **) &queue, NULL) ; \
22 : LAGraph_Free ((void **) &depth, NULL) ; \
23 : LAGraph_Free ((void **) &bc_vertex_flow, NULL) ;\
24 : LAGraph_Free ((void **) &S, NULL) ; \
25 : LAGraph_Free ((void **) &paths, NULL) ; \
26 : LAGraph_Free ((void **) &Pj, NULL) ; \
27 : LAGraph_Free ((void **) &Ptail, NULL) ; \
28 : if(AT != G->AT) \
29 : { \
30 : GrB_free (&AT) ; \
31 : } \
32 : }
33 :
34 : #define LG_FREE_ALL \
35 : { \
36 : LG_FREE_WORK ; \
37 : LAGraph_Free ((void **) &Ap, NULL) ; \
38 : LAGraph_Free ((void **) &Aj, NULL) ; \
39 : LAGraph_Free ((void **) &Ax, NULL) ; \
40 : LAGraph_Free ((void **) &ATp, NULL) ; \
41 : LAGraph_Free ((void **) &ATj, NULL) ; \
42 : LAGraph_Free ((void **) &ATx, NULL) ; \
43 : LAGraph_Free ((void **) &result, NULL) ; \
44 : }
45 :
46 : #include "LG_internal.h"
47 : #include <LAGraphX.h>
48 : #include "LG_Xtest.h"
49 :
50 :
51 : //------------------------------------------------------------------------------
52 : // test the results from a Edge Betweenness Centrality
53 : //------------------------------------------------------------------------------
54 :
55 17 : int LG_check_edgeBetweennessCentrality
56 : (
57 : // output
58 : GrB_Matrix *C, // centrality matrix
59 : // input
60 : LAGraph_Graph G,
61 : GrB_Vector sources, // source vertices to compute shortest paths (if NULL or empty, use all vertices)
62 : char *msg
63 : )
64 : {
65 : #if LAGRAPH_SUITESPARSE
66 :
67 : //--------------------------------------------------------------------------
68 : // initialize workspace variables
69 : //--------------------------------------------------------------------------
70 :
71 17 : double tt = LAGraph_WallClockTime ( ) ;
72 : GrB_Info info ;
73 :
74 : // Array storing shortest path distances from source to each vertex
75 17 : int64_t *depth = NULL ;
76 :
77 : // Array storing dependency scores during accumulation phase
78 17 : double *bc_vertex_flow = NULL ;
79 :
80 : // Stack used for backtracking phase in dependency accumulation
81 17 : int64_t *S = NULL ;
82 :
83 : // Queue used for BFS traversal
84 17 : int64_t *queue = NULL ;
85 :
86 : // Predecessor list components:
87 : // Pj: array of predecessor vertices
88 : // Ptail: end indices for each vertex's predecessor list
89 : // Phead: start indices for each vertex's predecessor list
90 17 : GrB_Index *Pj = NULL ;
91 17 : GrB_Index *Ptail = NULL ;
92 17 : GrB_Index *Phead = NULL ;
93 :
94 : // Array storing number of shortest paths to each vertex
95 17 : double *paths = NULL ;
96 :
97 : // Temporary array for centrality results
98 17 : double *result = NULL;
99 :
100 17 : GrB_Vector internal_sources = NULL;
101 17 : bool created_sources = false;
102 :
103 17 : GrB_Matrix AT = NULL ;
104 :
105 : //--------------------------------------------------------------------------
106 : // check inputs
107 : //--------------------------------------------------------------------------
108 :
109 17 : GrB_Index *Ap = NULL, *Aj = NULL, *neighbors = NULL, *ATp = NULL, *ATj = NULL ;
110 17 : void *Ax = NULL, *ATx = NULL ;
111 : GrB_Index Ap_size, Aj_size, Ax_size, n, nvals, ATp_size, ATj_size, ATx_size ;
112 17 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
113 17 : GRB_TRY (GrB_Matrix_nrows (&n, G->A)) ;
114 17 : GRB_TRY (GrB_Matrix_nvals (&nvals, G->A)) ;
115 : #if defined ( COVERAGE )
116 17 : bool print_timings = true ;
117 : #else
118 : bool print_timings = (n >= 2000) ;
119 : #endif
120 :
121 17 : LG_TRY (LAGraph_DeleteSelfEdges (G, msg)) ;
122 :
123 17 : GrB_Matrix A = G->A ;
124 :
125 17 : LG_TRY (LAGraph_Cached_AT (G, msg)) ;
126 :
127 17 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
128 15 : G->is_symmetric_structure == LAGraph_TRUE)
129 : {
130 : // A and A' have the same structure
131 : // AT = A;
132 : // GrB_Matrix_new (&AT, GrB_FP64, n, n) ;
133 2 : GrB_Matrix_dup (&AT, A) ;
134 : }
135 : else
136 : {
137 : // A and A' differ
138 15 : AT = G->AT ;
139 15 : LG_ASSERT_MSG (AT != NULL, LAGRAPH_NOT_CACHED, "G->AT is required") ;
140 : }
141 :
142 :
143 : //--------------------------------------------------------------------------
144 :
145 17 : LG_CLEAR_MSG ;
146 :
147 : //--------------------------------------------------------------------------
148 : // allocate workspace
149 : //--------------------------------------------------------------------------
150 :
151 17 : LG_TRY(LAGraph_Malloc((void **)&depth, n, sizeof(int64_t), msg));
152 :
153 17 : LG_TRY(LAGraph_Calloc((void **)&bc_vertex_flow, n, sizeof(double), msg));
154 :
155 17 : LG_TRY(LAGraph_Malloc((void **)&S, n, sizeof(int64_t), msg));
156 :
157 17 : LG_TRY(LAGraph_Malloc((void **)&queue, n, sizeof(int64_t), msg));
158 :
159 : //--------------------------------------------------------------------------
160 : // bfs on the A
161 : //--------------------------------------------------------------------------
162 :
163 17 : if (print_timings)
164 : {
165 17 : tt = LAGraph_WallClockTime ( ) - tt ;
166 17 : printf ("LG_check_bfs init time: %g sec\n", tt) ;
167 17 : tt = LAGraph_WallClockTime ( ) ;
168 : }
169 :
170 : // Initialize centrality matrix result to 0
171 : // 1. result [(v, w)] ← 0, ∀(v, w) ∈ E
172 : // A temporary result centrality matrix initialized to 0 for all vertice,
173 : // -- further changes would need to be made to make it a dictionary of edges.
174 17 : GrB_Index result_size = n * n ;
175 17 : LG_TRY(LAGraph_Calloc((void **)&result, result_size, sizeof(double), msg));
176 :
177 : // result (v,w) is held in result (INDEX(v,w)):
178 : #define INDEX(i,j) ((i)*n+(j))
179 :
180 : //--------------------------------------------------------------------------
181 : // unpack the A matrix in CSR form for SuiteSparse:GraphBLAS
182 : //--------------------------------------------------------------------------
183 :
184 : #if LAGRAPH_SUITESPARSE
185 : bool iso, AT_iso ;
186 17 : GRB_TRY (GxB_Matrix_unpack_CSR (A,
187 : &Ap, &Aj, &Ax, &Ap_size, &Aj_size, &Ax_size, &iso, NULL, NULL)) ;
188 :
189 17 : GRB_TRY (GxB_Matrix_unpack_CSR (AT,
190 : &ATp, &ATj, &ATx, &ATp_size, &ATj_size, &ATx_size, &AT_iso, NULL, NULL)) ;
191 : #endif
192 :
193 17 : Phead = ATp ;
194 :
195 : //--------------------------------------------------------------------------
196 :
197 17 : LG_TRY(LAGraph_Malloc((void **)&Pj, nvals, sizeof(GrB_Index), msg));
198 17 : LG_TRY(LAGraph_Malloc((void **)&Ptail, n, sizeof(GrB_Index), msg)); // might need to be + 1
199 :
200 17 : LAGraph_Calloc ((void **) &paths, n, sizeof (double), msg) ;
201 :
202 17 : if (sources == NULL)
203 : {
204 : // Create a vector with all nodes as sources
205 8 : GRB_TRY (GrB_Vector_new (&internal_sources, GrB_INT64, n)) ;
206 :
207 : // internal_sources (0:n-1) = 0
208 8 : GRB_TRY (GrB_assign (internal_sources, NULL, NULL, 0, GrB_ALL, n, NULL)) ;
209 :
210 : // internal_sources (0:n-1) = 0:n-1
211 8 : GRB_TRY (GrB_apply (internal_sources, NULL, NULL, GrB_ROWINDEX_INT64,
212 : internal_sources, 0, NULL)) ;
213 :
214 : // Use this vector instead
215 8 : sources = internal_sources;
216 8 : created_sources = true;
217 : }
218 :
219 : // Extract number of source nodes
220 : GrB_Index nvals_sources;
221 17 : GRB_TRY (GrB_Vector_nvals (&nvals_sources, sources));
222 :
223 17 : if (nvals_sources == 0)
224 : {
225 : // If sources vector is empty, return an empty centrality matrix
226 : // (Create an empty C matrix or set it to zeros as needed)
227 1 : GRB_TRY (GrB_Matrix_new(C, GrB_FP64, n, n));
228 :
229 : // Clean up resources
230 1 : if (created_sources) GRB_TRY (GrB_free(&internal_sources));
231 1 : LG_FREE_ALL ;
232 1 : return (GrB_SUCCESS);
233 : }
234 :
235 : // =========================================================================
236 : // === Main computation loop ==============================================
237 : // =========================================================================
238 :
239 : GrB_Index s;
240 :
241 : // Process each source vertex
242 3126 : for (GrB_Index i = 0; i < nvals_sources; i++) {
243 3110 : GRB_TRY (GrB_Vector_extractElement(&s, sources, i)) ;
244 :
245 : // check for invalid indices
246 3110 : LG_ASSERT (s < n, GrB_INVALID_VALUE) ;
247 :
248 3110 : size_t sp = 0; // stack pointer
249 3110 : memcpy(Ptail, ATp, n * sizeof(GrB_Index));
250 :
251 : // Initialize path counts
252 2594714 : for (int64_t i = 0; i < n; i++) {
253 2591604 : paths[i] = 0;
254 : }
255 3110 : paths[s] = 1;
256 :
257 : // Initialize distances
258 2594714 : for (size_t t = 0; t < n; t++) {
259 2591604 : depth[t] = -1;
260 : }
261 3110 : depth[s] = 0;
262 :
263 : //----------------------------------------------------------------------
264 : // BFS phase to compute shortest paths
265 : //----------------------------------------------------------------------
266 :
267 3110 : int64_t qh = 0, qt = 0; // queue head and tail
268 3110 : queue[qt++] = s; // enqueue source
269 :
270 2594673 : while (qh < qt) {
271 2591563 : int64_t v = queue[qh++];
272 2591563 : S[sp++] = v;
273 :
274 : // Process neighbors of current vertex
275 44269633 : for (int64_t p = Ap[v]; p < Ap[v+1]; p++) {
276 41678070 : int64_t w = Aj[p];
277 :
278 : // Handle unvisited vertices
279 41678070 : if (depth[w] < 0) {
280 2588453 : queue[qt++] = w;
281 2588453 : depth[w] = depth[v] + 1;
282 : }
283 :
284 : // Update path counts for vertices at next level
285 41678070 : if (depth[w] == depth[v] + 1) {
286 19588517 : paths[w] += paths[v];
287 :
288 : #if 1
289 19588517 : LG_ASSERT (Ptail [w] < Phead [w+1], GrB_INVALID_VALUE) ;
290 19588517 : LG_ASSERT (Ptail [w] >= Phead [w], GrB_INVALID_VALUE) ;
291 : #else
292 : if (Ptail [w] >= Phead [w+1] || Ptail [w] < Phead [w])
293 : {
294 : printf ("Ack! w=%ld Ptail [w]=%ld, Phead [w]=%ld Phead[w+1]=%ld\n",
295 : w, Ptail [w], Phead [w], Phead [w+1]) ;
296 : fflush (stdout) ; abort ( ) ;
297 : }
298 : #endif
299 :
300 19588517 : Pj[Ptail[w]++] = v;
301 :
302 : }
303 : }
304 : }
305 :
306 : //----------------------------------------------------------------------
307 : // Dependency accumulation phase
308 : //----------------------------------------------------------------------
309 :
310 : // Initialize dependency scores
311 2594714 : for (size_t v = 0; v < n; v++) {
312 2591604 : bc_vertex_flow[v] = 0;
313 : }
314 :
315 : // Process vertices in reverse order of discovery
316 2594673 : while (sp > 0) {
317 2591563 : int64_t w = S[--sp];
318 :
319 : // Update dependencies through predecessors
320 22180080 : for (int64_t p = Phead[w]; p < Ptail[w]; p++) {
321 19588517 : int64_t v = Pj[p];
322 :
323 : // Compute and accumulate dependency
324 19588517 : double centrality = paths[v] * ((bc_vertex_flow[w] + 1) / paths[w]);
325 19588517 : bc_vertex_flow[v] += centrality;
326 :
327 19588517 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED) {
328 1976 : result[INDEX(v,w)] += centrality / 2;
329 1976 : result[INDEX(w,v)] += centrality / 2;
330 : }
331 : else {
332 19586541 : result[INDEX(v,w)] += centrality;
333 : }
334 : }
335 : }
336 : }
337 :
338 16 : if (print_timings)
339 : {
340 16 : tt = LAGraph_WallClockTime ( ) - tt ;
341 16 : printf ("LG_check_edgeBetweenessCentrality time: %g sec\n", tt) ;
342 16 : tt = LAGraph_WallClockTime ( ) ;
343 : }
344 :
345 : //--------------------------------------------------------------------------
346 : // repack the A matrix in CSR form for SuiteSparse:GraphBLAS
347 : //--------------------------------------------------------------------------
348 :
349 : #if LAGRAPH_SUITESPARSE
350 16 : GRB_TRY (GxB_Matrix_pack_CSR (A,
351 : &Ap, &Aj, &Ax, Ap_size, Aj_size, Ax_size, iso, false, NULL)) ;
352 16 : GRB_TRY (GxB_Matrix_pack_CSR (AT,
353 : &ATp, &ATj, &ATx, ATp_size, ATj_size, ATx_size, AT_iso, false, NULL)) ;
354 : #endif
355 :
356 : #if 0
357 : GrB_Info GxB_Matrix_pack_FullR // pack a full matrix, held by row
358 : (
359 : GrB_Matrix A, // matrix to create (type, nrows, ncols unchanged)
360 : void **Ax, // values, Ax_size >= nrows*ncols * (type size)
361 : // or Ax_size >= (type size), if iso is true
362 : GrB_Index Ax_size, // size of Ax in bytes
363 : bool iso, // if true, A is iso
364 : const GrB_Descriptor desc
365 : ) ;
366 : #endif
367 :
368 : GrB_Matrix C_temp;
369 16 : LG_TRY (GrB_Matrix_new(&C_temp, GrB_FP64, n, n)) ;
370 16 : LG_TRY (GxB_Matrix_pack_FullR(C_temp, (void **) &result, result_size * sizeof(double), false, NULL) ) ;
371 :
372 16 : LG_TRY (GrB_assign(C_temp, A, NULL, C_temp, GrB_ALL, n, GrB_ALL, n, GrB_DESC_RS)) ;
373 :
374 16 : *C = C_temp;
375 :
376 16 : if (created_sources) {
377 8 : GRB_TRY (GrB_free(&internal_sources));
378 : }
379 :
380 : //--------------------------------------------------------------------------
381 : // free workspace and return result
382 : //--------------------------------------------------------------------------
383 16 : LG_FREE_WORK ;
384 :
385 16 : if (print_timings)
386 : {
387 16 : tt = LAGraph_WallClockTime ( ) - tt ;
388 16 : printf ("LG_check_edgeBetweennessCentrality check time: %g sec\n", tt) ;
389 : }
390 16 : return (GrB_SUCCESS) ;
391 : #else
392 : return (GrB_NOT_IMPLEMENTED) ;
393 : #endif
394 : }
|
