Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_RichClubCoefficient: rich club coefficient of a graph
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved.
6 : // SPDX-License-Identifier: BSD-2-Clause
7 : //
8 : // For additional details (including references to third party source code and
9 : // other files) see the LICENSE 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 Gabriel Gomez, Texas A&M University
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // Get the rich club coefficient of a graph.
19 :
20 : // Given a Symetric Graph with no self edges, LAGraph_RichClubCoefficient will
21 : // calculate the rich club coefficients of the graph.
22 :
23 : // The values will be output as a sparse GrB_Vector, the rich club coefficient
24 : // of k will be found at the closest entry at or above k.
25 :
26 : // The G->out_degree cached property must be defined for this method.
27 :
28 : // References:
29 :
30 : // Julian J. McAuley, Luciano da Fontoura Costa, and Tibério S. Caetano, “The
31 : // rich-club phenomenon across complex network hierarchies”, Applied Physics
32 : // Letters Vol 91 Issue 8, August 2007. https://arxiv.org/abs/physics/0701290
33 :
34 : #define LG_FREE_WORK \
35 : { \
36 : /* free any workspace used here */ \
37 : GrB_free(&D) ; \
38 : GrB_free(&P) ; \
39 : GrB_free(&A_deg) ; \
40 : GrB_free(°rees) ; \
41 : GrB_free(°_x) ; \
42 : GrB_free(&node_edges) ; \
43 : GrB_free(&node_edges_x) ; \
44 : GrB_free(&ones_v) ; \
45 : GrB_free(&edges_per_deg) ; \
46 : GrB_free(&verts_per_deg) ; \
47 : GrB_free(&iseq_2lt) ; \
48 : GrB_free(&plus_2le) ; \
49 : GrB_free(&rcCalculation) ; \
50 : GrB_free(&ramp_v) ; \
51 : LAGraph_Free(&a_space, NULL) ; \
52 : }
53 :
54 :
55 : #define LG_FREE_ALL \
56 : { \
57 : /* free any workspace used here */ \
58 : LG_FREE_WORK ; \
59 : /* free all the output variable(s) */ \
60 : GrB_free (rccs) ; \
61 : }
62 :
63 : #include "LG_internal.h"
64 : #include "LAGraphX.h"
65 :
66 : typedef void (*LAGraph_binary_function) (void *, const void *, const void *) ;
67 :
68 : #define ISEQ_2ISLT \
69 : "void iseq_2islt(int64_t *z, const int64_t *x, const int64_t *y) \n"\
70 : "{ \n"\
71 : "(*z) = (int64_t)((*x < *y) + (*x <= *y)) ; \n"\
72 : "}"
73 5816980 : void iseq_2islt(int64_t *z, const int64_t *x, const int64_t *y)
74 : {
75 5816980 : (*z) = (int64_t)((*x < *y) + (*x <= *y)) ;
76 5816980 : }
77 :
78 : #define RICH_CLUB_FORMULA \
79 : "void rich_club_formula(double *z, const int64_t *x, const int64_t *y) \n"\
80 : "{ \n"\
81 : " (*z) = ((double)(*x)) / (((double)(*y)) * (((double)(*y)) - 1.0)) ; \n"\
82 : "}"
83 126 : void rich_club_formula(double *z, const int64_t *x, const int64_t *y)
84 : {
85 126 : (*z) = ((double)(*x)) / (((double)(*y)) * (((double)(*y)) - 1.0));
86 126 : }
87 402 : int LAGraph_RichClubCoefficient
88 : (
89 : // output:
90 : //rccs(i): rich club coefficent of i
91 : GrB_Vector *rccs,
92 :
93 : // input:
94 : LAGraph_Graph G, //input graph
95 : char *msg
96 : )
97 : {
98 : //--------------------------------------------------------------------------
99 : // Declorations
100 : //--------------------------------------------------------------------------
101 402 : LG_CLEAR_MSG ;
102 :
103 : // n x n Adjacency Matrix
104 : // With values cooresponding to the degree of its column
105 402 : GrB_Matrix A_deg = NULL;
106 :
107 : // n x n Diagonal Matrix
108 : // entries corresponding to degrees.
109 402 : GrB_Matrix D = NULL;
110 :
111 : // n degrees vector
112 402 : GrB_Vector degrees = NULL, deg_x = NULL;
113 :
114 : // n x 1
115 : // contains the number of edges for which the ith node is
116 : // the smallest degree node * 2 + # edges w/ same degree as the other node
117 : // to account for double counting of edges w/ same degree as the other node.
118 402 : GrB_Vector node_edges = NULL, node_edges_x = NULL;
119 :
120 : // max_degree x 1
121 : // the ith entry contains the number of edges whose lowest degree is i.
122 402 : GrB_Vector edges_per_deg = NULL;
123 :
124 : // max_degree x 1
125 : // the ith entry contains the number of verticies whose degree is i.
126 402 : GrB_Vector verts_per_deg = NULL;
127 :
128 : // edge_vec_nvals x 1
129 : // Vector of ones
130 402 : GrB_Vector ones_v = NULL;
131 :
132 : // Ramp vector
133 402 : GrB_Vector ramp_v = NULL;
134 :
135 : // 2 * (x < y) + (x == y)
136 402 : GrB_BinaryOp iseq_2lt = NULL;
137 :
138 : // [+].[iseq_2lt]
139 402 : GrB_Semiring plus_2le = NULL;
140 :
141 : // 2E_K / (N_k (N_k -1))
142 402 : GrB_BinaryOp rcCalculation = NULL;
143 :
144 402 : GrB_Matrix A = NULL; // G->A, the adjacency matrix
145 :
146 : // Matrix used for row reduction
147 402 : GrB_Matrix P = NULL;
148 :
149 : GrB_Index n ;
150 :
151 : GrB_Index edge_vec_nvals;
152 : int64_t max_deg;
153 402 : bool iso = false;
154 :
155 402 : void *a_space = NULL;
156 :
157 402 : int64_t *node_edges_arr = NULL, *deg_arr = NULL,
158 402 : *epd_arr = NULL, *ones = NULL,
159 402 : *vpd_arr = NULL;
160 402 : GrB_Type epd_type = NULL, vpd_type = NULL;
161 402 : uint64_t epd_n = 0, vpd_n = 0, epd_size = 0, vpd_size = 0;
162 402 : int epd_h = 0, vpd_h = 0;
163 402 : GrB_Index *epd_index = NULL, *vpd_index = NULL;
164 :
165 : //--------------------------------------------------------------------------
166 : // Check inputs
167 : //--------------------------------------------------------------------------
168 402 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
169 402 : LG_ASSERT (rccs != NULL, GrB_NULL_POINTER);
170 :
171 402 : LG_ASSERT_MSG(
172 : G->kind == LAGraph_ADJACENCY_UNDIRECTED, GrB_INVALID_VALUE,
173 : "G->A must be symmetric") ;
174 402 : LG_ASSERT_MSG (G->out_degree != NULL, GrB_EMPTY_OBJECT,
175 : "G->out_degree must be defined") ;
176 402 : LG_ASSERT_MSG (G->nself_edges == 0, GrB_INVALID_VALUE,
177 : "G->nself_edges must be zero") ;
178 :
179 : //--------------------------------------------------------------------------
180 : // Initializations
181 : //--------------------------------------------------------------------------
182 402 : A = G->A ;
183 402 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
184 402 : GRB_TRY (GrB_Matrix_new(&A_deg, GrB_INT64,n,n)) ;
185 :
186 390 : GRB_TRY (GrB_Vector_new(°rees, GrB_INT64, n)) ;
187 382 : GRB_TRY (GrB_Vector_new(&node_edges, GrB_INT64, n)) ;
188 : #if LAGRAPH_SUITESPARSE
189 374 : GRB_TRY (GxB_BinaryOp_new(
190 : &iseq_2lt, (LAGraph_binary_function) (&iseq_2islt),
191 : GrB_INT64, GrB_INT64, GrB_INT64, "iseq_2islt", ISEQ_2ISLT)) ;
192 366 : GRB_TRY (GxB_BinaryOp_new(
193 : &rcCalculation, (LAGraph_binary_function) (&rich_club_formula),
194 : GrB_FP64, GrB_INT64, GrB_INT64,
195 : "rich_club_formula", RICH_CLUB_FORMULA)) ;
196 : #else
197 : GRB_TRY (GrB_BinaryOp_new(
198 : &iseq_2lt, (LAGraph_binary_function) (&iseq_2islt),
199 : GrB_INT64, GrB_INT64, GrB_INT64)) ;
200 : GRB_TRY (GrB_BinaryOp_new(
201 : &rcCalculation, (LAGraph_binary_function) (&rich_club_formula),
202 : GrB_FP64, GrB_INT64, GrB_INT64 )) ;
203 : #endif
204 :
205 358 : GRB_TRY (GrB_Semiring_new(&plus_2le, GrB_PLUS_MONOID_INT64, iseq_2lt)) ;
206 :
207 350 : GRB_TRY (GrB_Vector_reduce_INT64(
208 : &max_deg, NULL, GrB_MAX_MONOID_INT64, G->out_degree, NULL)) ;
209 350 : GRB_TRY (GrB_Vector_new(&edges_per_deg, GrB_INT64, max_deg)) ;
210 342 : GRB_TRY (GrB_Vector_new(&verts_per_deg, GrB_INT64, max_deg)) ;
211 334 : GRB_TRY (GrB_Vector_new(rccs, GrB_FP64, max_deg)) ;
212 :
213 : //--------------------------------------------------------------------------
214 : // Calculations
215 : //--------------------------------------------------------------------------
216 :
217 : // degrees = G->out_degree - 1
218 : // Fill out degree vector, to target col_scale mxm on graphs
219 : // with singletons, scalar value irrelevant.
220 326 : GRB_TRY (GrB_Vector_assign_INT64(
221 : degrees, NULL, NULL, (int64_t) -1, GrB_ALL, 0, NULL)) ;
222 322 : GRB_TRY (GrB_Vector_assign(
223 : degrees, NULL, GrB_PLUS_INT64, G->out_degree, GrB_ALL, 0, NULL)) ;
224 318 : GRB_TRY (GrB_Matrix_diag(&D, degrees, 0)) ;
225 :
226 : // Each edge in the graph gets the value of the degree of its row node
227 : #if LAGRAPH_SUITESPARSE
228 294 : GRB_TRY (GrB_mxm(
229 : A_deg, NULL, NULL, GxB_ANY_FIRST_INT64, D, A, NULL)) ;
230 : #else
231 : GRB_TRY (GrB_mxm(
232 : A_deg, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_INT64, D, A, NULL)) ;
233 : #endif
234 : // Sum the number of edges each node is "responsible" for.
235 282 : GRB_TRY (GrB_mxv(
236 : node_edges, NULL, GrB_PLUS_INT64, plus_2le, A_deg, degrees, NULL)) ;
237 :
238 : // The rest of this is indexing the number of edges and number of nodes at
239 : // each degree and then doing a cummulative sum to know the amount of edges
240 : // and nodes at degree geq k.
241 276 : GRB_TRY (GrB_Vector_nvals (&edge_vec_nvals, node_edges)) ;
242 : #if LG_SUITESPARSE_GRAPHBLAS_V10
243 276 : if(n == edge_vec_nvals)
244 : {
245 199 : deg_x = degrees;
246 199 : degrees = NULL;
247 199 : node_edges_x = node_edges;
248 199 : node_edges = NULL;
249 : }
250 : else
251 : {
252 77 : GRB_TRY (GrB_Vector_assign(
253 : degrees, G->out_degree, NULL, degrees, GrB_ALL, 0, GrB_DESC_RS
254 : )) ;
255 75 : GRB_TRY (GrB_Vector_new(°_x, GrB_BOOL, 0)) ;
256 73 : GRB_TRY (GrB_Vector_new(&node_edges_x, GrB_BOOL, 0)) ;
257 71 : GRB_TRY (GxB_Vector_extractTuples_Vector(
258 : NULL, deg_x, degrees, NULL
259 : )) ;
260 69 : GRB_TRY (GxB_Vector_extractTuples_Vector(
261 : NULL, node_edges_x, node_edges, NULL
262 : )) ;
263 : }
264 266 : GRB_TRY (GrB_Vector_nvals(&edge_vec_nvals, node_edges_x))
265 266 : GRB_TRY (GrB_Vector_new(&ones_v, GrB_INT64, edge_vec_nvals)) ;
266 :
267 :
268 259 : GRB_TRY (GrB_Vector_assign_INT64(
269 : edges_per_deg, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ;
270 255 : GRB_TRY (GrB_Vector_assign_INT64(
271 : verts_per_deg, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ;
272 251 : GRB_TRY (GrB_Vector_assign_INT64(
273 : ones_v, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ;
274 :
275 : #ifndef COVERAGE
276 : GRB_TRY (GrB_Vector_new(&ramp_v, GrB_INT64, edge_vec_nvals + 1)) ;
277 : GRB_TRY (GrB_Vector_assign_INT64(
278 : ramp_v, NULL, NULL, (int64_t) 0, GrB_ALL, 0, NULL)) ;
279 : GRB_TRY (GrB_apply (
280 : ramp_v, NULL, NULL, GrB_ROWINDEX_INT64, ramp_v, 0, NULL)) ;
281 : #endif
282 :
283 247 : LG_TRY (LAGraph_FastAssign_Semiring (
284 : edges_per_deg, NULL, GrB_PLUS_INT64, deg_x, node_edges_x, ramp_v,
285 : GxB_PLUS_SECOND_INT64, NULL, msg
286 : )) ;
287 135 : LG_TRY (LAGraph_FastAssign_Semiring (
288 : verts_per_deg, NULL, GrB_PLUS_INT64, deg_x, ones_v, ramp_v,
289 : GxB_PLUS_PAIR_INT64, NULL, msg
290 : )) ;
291 :
292 23 : GRB_TRY (GxB_Vector_unload(
293 : edges_per_deg, (void **) &epd_arr, &epd_type,
294 : &epd_n, &epd_size, &epd_h, NULL)) ;
295 23 : GRB_TRY (GxB_Vector_unload(
296 : verts_per_deg, (void **) &vpd_arr, &vpd_type,
297 : &vpd_n, &vpd_size, &vpd_h, NULL)) ;
298 :
299 23 : LG_ASSERT (max_deg == vpd_n && max_deg == epd_n, GrB_INVALID_VALUE) ;
300 : //run a cummulative sum (backwards) on vpd_arr
301 2501 : for(int64_t i = max_deg - 1; i > 0; --i)
302 : {
303 2478 : vpd_arr[i-1] += vpd_arr[i] ;
304 2478 : epd_arr[i-1] += epd_arr[i] ;
305 : }
306 23 : GRB_TRY(GxB_Vector_load(
307 : edges_per_deg, (void **) &epd_arr, epd_type,
308 : epd_n, epd_size, epd_h, NULL)) ;
309 23 : GRB_TRY(GxB_Vector_load(
310 : verts_per_deg, (void **) &vpd_arr, vpd_type,
311 : vpd_n, vpd_size, vpd_h, NULL)) ;
312 : #else
313 : LG_TRY (LAGraph_Malloc(
314 : &a_space, edge_vec_nvals * 3 + max_deg * 4, sizeof(int64_t), NULL
315 : )) ;
316 : int64_t *T = a_space;
317 : deg_arr = T; T += edge_vec_nvals;
318 : node_edges_arr = T; T += edge_vec_nvals;
319 : ones = T; T += edge_vec_nvals;
320 : epd_arr = T; T += max_deg;
321 : vpd_arr = T; T += max_deg;
322 : epd_index = T; T += max_deg;
323 : vpd_index = T; T += max_deg;
324 :
325 : #pragma omp parallel for schedule(static)
326 : for(uint64_t i = 0; i < edge_vec_nvals; ++i)
327 : {
328 : ones[i] = 1ll;
329 : }
330 : GRB_TRY (GrB_Vector_apply_BinaryOp2nd_INT64(
331 : degrees, NULL, NULL, GrB_MINUS_INT64, G->out_degree, 1, NULL)) ;
332 : //TODO: remove NULL for Vanilla GB
333 : GRB_TRY (GrB_Vector_extractTuples_INT64(
334 : NULL, deg_arr, &edge_vec_nvals, degrees
335 : )) ;
336 : GRB_TRY (GrB_Vector_extractTuples_INT64(
337 : NULL, node_edges_arr, &edge_vec_nvals, node_edges
338 : )) ;
339 :
340 : // Build with degrees as indecies and handle duplicates via adition
341 : GRB_TRY (GrB_Vector_build_INT64 (
342 : edges_per_deg, deg_arr, node_edges_arr, edge_vec_nvals,
343 : GrB_PLUS_INT64)) ;
344 : GRB_TRY (GrB_Vector_build_INT64 (
345 : verts_per_deg, deg_arr, ones, edge_vec_nvals, GrB_PLUS_INT64)) ;
346 : GRB_TRY (GrB_Vector_assign_INT64(
347 : edges_per_deg, edges_per_deg, NULL, (int64_t) 0,
348 : GrB_ALL, 0, GrB_DESC_SC)) ;
349 : GRB_TRY (GrB_Vector_assign_INT64(
350 : verts_per_deg, verts_per_deg, NULL, (int64_t) 0,
351 : GrB_ALL, 0, GrB_DESC_SC)) ;
352 :
353 : // Extract into arrays
354 : GRB_TRY (GrB_Vector_extractTuples_INT64(
355 : epd_index, epd_arr, &max_deg, edges_per_deg
356 : )) ;
357 : GRB_TRY (GrB_Vector_extractTuples_INT64(
358 : vpd_index, vpd_arr, &max_deg, verts_per_deg
359 : )) ;
360 : //run a cummulative sum (backwards) on vpd_arr
361 : for(int64_t i = max_deg - 1; i > 0; --i)
362 : {
363 : vpd_arr[i-1] += vpd_arr[i] ;
364 : epd_arr[i-1] += epd_arr[i] ;
365 : }
366 : GRB_TRY (GrB_Vector_clear(edges_per_deg)) ;
367 : GRB_TRY (GrB_Vector_clear(verts_per_deg)) ;
368 : GRB_TRY (GrB_Vector_build_INT64(
369 : edges_per_deg, epd_index, epd_arr, max_deg, NULL
370 : )) ;
371 : GRB_TRY (GrB_Vector_build_INT64(
372 : verts_per_deg, vpd_index, vpd_arr, max_deg, NULL
373 : )) ;
374 : T = deg_arr = node_edges_arr = ones = NULL ;
375 : epd_index = vpd_index = epd_arr = vpd_arr = NULL ;
376 : #endif
377 :
378 : /**
379 : * Cumulative sum (TODO: should be a GBLAS method!)
380 : *
381 : * GrB_cumsum(GrB_Matrix C, const GrB_Matrix mask, const GrB_BinaryOp accum,
382 : * const GrB_BinaryOp plus, GrB_Matrix A, const GrB_Descriptor desc)
383 : *
384 : * By default sums rows. Returns a nearly full matrix:
385 : * [., ., 1, 1, 1, 1, ., ., 1] --> [., ., 1, 2, 3, 4, 4, 4, 5]
386 : * Mask can be A, then returns a matrix with the same pattern.
387 : * [., ., 1, 1, 1, 1, ., ., 1] --> [., ., 1, 2, 3, 4, ., ., 5]
388 : *
389 : * Should we be able to sum in the opposite direction?
390 : * Yes since not all monoids have inverse operations.
391 : *
392 : * If plus biop is not a monoid, this method should still work?
393 : */
394 :
395 : //Computes the RCC of a matrix
396 23 : GRB_TRY(GrB_eWiseMult(
397 : *rccs, NULL, NULL, rcCalculation,
398 : edges_per_deg, verts_per_deg, NULL
399 : )) ;
400 :
401 15 : LG_FREE_WORK ;
402 15 : return (GrB_SUCCESS) ;
403 : }
|
