Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LG_check_lcc: compute local clustering coefficients of a graph (simple)
3 : //------------------------------------------------------------------------------
4 :
5 : // LAGraph, (c) 2023 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 Pascal Costanza, Intel, Belgium
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // An implementation of local clustering coefficient based on the descripion at
19 : // https://en.wikipedia.org/wiki/Clustering_coefficient
20 : // This method is for testing only, to check the result of other, faster methods.
21 : // Do not benchmark this method; it is simple by design.
22 :
23 : #define LG_FREE_ALL \
24 : { \
25 : GrB_free (&S) ; \
26 : GrB_free (&W) ; \
27 : GrB_free (&LCC) ; \
28 : LAGraph_Free ((void **) &Sp, msg) ; \
29 : LAGraph_Free ((void **) &Si, msg) ; \
30 : LAGraph_Free ((void **) &Tp, msg) ; \
31 : LAGraph_Free ((void **) &Ti, msg) ; \
32 : LAGraph_Free ((void **) &vb, msg) ; \
33 : LAGraph_Free ((void **) &vx, msg) ; \
34 : }
35 :
36 : #include <stdlib.h>
37 : #include "LG_internal.h"
38 : #include "LG_test.h"
39 : #include "LG_test.h"
40 : #include "LG_Xtest.h"
41 :
42 : // assumes that the indices array is sorted
43 2496168 : GrB_Index find(const GrB_Index* indices, GrB_Index n, GrB_Index index) {
44 2496168 : GrB_Index i = 0, j = n, h;
45 14259134 : while (i < j) {
46 11762966 : h = (i+j) >> 1;
47 11762966 : if (indices[h] < index) {
48 1617048 : i = h + 1;
49 : } else {
50 10145918 : j = h;
51 : }
52 : }
53 2496168 : return i;
54 : }
55 :
56 : // computes how many elements are in the intersection of both sets
57 : // assumes that both arrays are sorted and don't have any duplicates
58 178144 : GrB_Index intersection_size(
59 : GrB_Index* x, GrB_Index nx,
60 : GrB_Index* y, GrB_Index ny
61 : ) {
62 178144 : GrB_Index n = 0;
63 2497316 : while ((nx > 0) && (ny > 0)) {
64 2319172 : if (y[0] > x[0]) {
65 413706 : GrB_Index* tmp = x; x = y; y = tmp;
66 413706 : GrB_Index ntmp = nx; nx = ny; ny = ntmp;
67 : }
68 2319172 : GrB_Index index = find(y, ny, x[0]);
69 2319172 : if ((index < ny) && (y[index] == x[0])) {
70 2067108 : n++;
71 2067108 : y += index+1;
72 2067108 : ny -= index+1;
73 : } else {
74 252064 : y += index;
75 252064 : ny -= index;
76 : }
77 2319172 : x += 1;
78 2319172 : nx -= 1;
79 : }
80 178144 : return n;
81 : }
82 :
83 : // This is a slow version.
84 20 : int LG_check_lcc(
85 : // outputs:
86 : GrB_Vector *coefficients, // the local clustering coefficients
87 : // inputs
88 : LAGraph_Graph G, // input graph
89 : char *msg
90 : )
91 : {
92 : #if LAGRAPH_SUITESPARSE
93 :
94 : //--------------------------------------------------------------------------
95 : // check inputs
96 : //--------------------------------------------------------------------------
97 :
98 20 : LG_CLEAR_MSG ;
99 :
100 40 : bool undirected = (G->kind == LAGraph_ADJACENCY_UNDIRECTED) ||
101 20 : ((G->kind == LAGraph_ADJACENCY_DIRECTED) &&
102 20 : (G->is_symmetric_structure == LAGraph_TRUE)) ;
103 20 : bool directed = !undirected ;
104 :
105 20 : GrB_Matrix S = NULL, T = NULL, W = NULL ;
106 20 : GrB_Index *Sp = NULL, *Si = NULL, *Tp = NULL, *Ti = NULL ;
107 20 : void *Sx = NULL , *Tx = NULL ;
108 20 : int8_t *vb = NULL ;
109 20 : double *vx = NULL ;
110 20 : GrB_Vector LCC = NULL ;
111 : GrB_Index n, ncols ;
112 20 : LG_ASSERT (coefficients != NULL, GrB_NULL_POINTER) ;
113 20 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
114 20 : GRB_TRY (GrB_Matrix_nrows (&n, G->A)) ; //set n to number of rows
115 20 : GRB_TRY (GrB_Matrix_ncols (&ncols, G->A)) ;
116 20 : LG_ASSERT (n == ncols, GrB_INVALID_OBJECT) ;
117 :
118 20 : GRB_TRY (GrB_Vector_new (&LCC, GrB_FP64, n)) ;
119 :
120 20 : GrB_Matrix A = G->A ;
121 :
122 20 : GRB_TRY(GrB_Matrix_new(&S, GrB_BOOL, n, n)) ;
123 20 : GRB_TRY(GrB_assign(S, A, GrB_NULL, true, GrB_ALL, n, GrB_ALL, n, GrB_DESC_S)) ;
124 20 : if (G->nself_edges != 0) {
125 6 : GRB_TRY(GrB_select(S, GrB_NULL, GrB_NULL, GrB_OFFDIAG, S, 0, GrB_NULL)) ;
126 : }
127 :
128 20 : if (directed)
129 : {
130 2 : GRB_TRY(GrB_Matrix_new(&W, GrB_BOOL, n, n)) ;
131 2 : GRB_TRY(GrB_eWiseAdd(W, GrB_NULL, GrB_NULL, GrB_ONEB_BOOL, S, S, GrB_DESC_T1)) ;
132 2 : T = W ;
133 : } else {
134 18 : T = S ;
135 : }
136 :
137 : GrB_Index Sp_size, Si_size, Sx_size ;
138 : bool Siso ;
139 20 : GRB_TRY(GxB_Matrix_unpack_CSR(S, &Sp, &Si, &Sx, &Sp_size, &Si_size, &Sx_size, &Siso, NULL, GrB_NULL));
140 20 : LAGraph_Free (&Sx, msg) ;
141 :
142 : GrB_Index Tp_size, Ti_size, Tx_size ;
143 : bool Tiso ;
144 20 : if (directed) {
145 2 : GRB_TRY(GxB_Matrix_unpack_CSR(T, &Tp, &Ti, &Tx, &Tp_size, &Ti_size, &Tx_size, &Tiso, NULL, GrB_NULL));
146 2 : LAGraph_Free (&Tx, msg) ;
147 : } else {
148 18 : Tp = Sp; Tp_size = Sp_size;
149 18 : Ti = Si; Ti_size = Si_size;
150 : }
151 :
152 20 : LAGraph_Calloc ((void **) &vb, n, sizeof (int8_t), msg) ;
153 20 : LAGraph_Malloc ((void **) &vx, n, sizeof (double), msg) ;
154 :
155 : int64_t i;
156 20 : GrB_Index nvals = 0 ;
157 :
158 : #pragma omp parallel for schedule(dynamic) reduction(+ : nvals)
159 6608 : for (i = 0; i < n; i++) {
160 6588 : GrB_Index *neighbors = Ti + Tp[i];
161 6588 : GrB_Index k = Tp[i + 1] - Tp[i];
162 6588 : if (k < 2) continue;
163 :
164 6576 : GrB_Index j, esum = 0;
165 184720 : for (j = 0; j < k; j++) {
166 178144 : GrB_Index e = neighbors[j];
167 178144 : GrB_Index *links = Si + Sp[e];
168 178144 : GrB_Index l = Sp[e + 1] - Sp[e];
169 178144 : if (l == 0) continue;
170 178144 : if (undirected) {
171 176996 : l = find(links, l, e);
172 : }
173 178144 : esum += intersection_size(neighbors, k, links, l);
174 : }
175 :
176 6576 : if (esum == 0) continue;
177 :
178 6558 : if (undirected) {
179 6428 : esum *= 2;
180 : }
181 6558 : vb[i] = true;
182 6558 : vx[i] = ((double) esum) / ((double) (k * (k - 1)));
183 6558 : nvals++;
184 : }
185 :
186 20 : if (!directed)
187 : {
188 : // do not free these; they are aliases to Sp and Si:
189 18 : Tp = NULL ;
190 18 : Ti = NULL ;
191 : }
192 :
193 20 : GRB_TRY(GxB_Vector_pack_Bitmap(LCC, &vb, (void**)&vx, n*sizeof(int8_t), n*sizeof(double), false, nvals, GrB_NULL)) ;
194 :
195 20 : *coefficients = LCC ; LCC = NULL ;
196 20 : LG_FREE_ALL;
197 20 : return (GrB_SUCCESS);
198 : #else
199 : return (GrB_NOT_IMPLEMENTED) ;
200 : #endif
201 : }
|
