Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_cdlp: community detection using label propagation
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 Gabor Szarnyas and Balint Hegyi, Budapest University of
15 : // Technology and Economics (with accented characters: G\'{a}bor Sz\'{a}rnyas
16 : // and B\'{a}lint Hegyi, using LaTeX syntax).
17 : // https://inf.mit.bme.hu/en/members/szarnyasg .
18 :
19 : // Modified by Pascal Costanza, Intel, Belgium
20 :
21 : // NOTE: the calloc/free below must be thread-safe,
22 : // so it cannot use LAGraph_Malloc/LAGraph_Free (which can use
23 : // the Rapids Memory Manager methods when using CUDA, and those
24 : // methods are not yet thread-safe).
25 :
26 : //------------------------------------------------------------------------------
27 :
28 : // ## Background
29 : //
30 : // This function was originally written for the LDBC Graphalytics benchmark.
31 : //
32 : // The community detection using label propagation (CDLP) algorithm is
33 : // defined both for directed and undirected graphs.
34 : //
35 : // The definition implemented here is described in the following document:
36 : // https://ldbc.github.io/ldbc_graphalytics_docs/graphalytics_spec.pdf
37 : //
38 : // The algorithm is based on the one given in the following paper:
39 : //
40 : // Usha Raghavan, Reka Albert, and Soundar Kumara. "Near linear time algorithm
41 : // to detect community structures in large-scale networks". In: Physical
42 : // Review E 76.3 (2007), p. 036106, https://arxiv.org/abs/0709.2938
43 : //
44 : // The key idea of the algorithm is that each vertex is assigned the label
45 : // that is most frequent among its neighbors. To allow reproducible
46 : // experiments, the algorithm is modified to guarantee deterministic behavior:
47 : // it always picks the smallest label in case of a tie:
48 : //
49 : // min ( argmax_{l} (#neighbors with label l) )
50 : //
51 : // In other words, we need to compute the *minimum mode value* (minmode) for
52 : // the labels among the neighbors.
53 : //
54 : // For directed graphs, a label on a neighbor that is connected through both
55 : // an outgoing and on an incoming edge counts twice:
56 : //
57 : // min ( argmax_{l} (#incoming neighbors with l + #outgoing neighbors with l) )
58 :
59 : #define LG_FREE_ALL \
60 : { \
61 : GrB_free (&S) ; \
62 : GrB_free (&T) ; \
63 : LAGraph_Free ((void **) &Sp, NULL) ; \
64 : LAGraph_Free ((void **) &Si, NULL) ; \
65 : LAGraph_Free ((void **) &Tp, NULL) ; \
66 : LAGraph_Free ((void **) &Ti, NULL) ; \
67 : LAGraph_Free ((void **) &L, NULL) ; \
68 : LAGraph_Free ((void **) &L_next, NULL) ; \
69 : ptable_pool_free (counts_pool, max_threads) ; \
70 : counts_pool = NULL ; \
71 : GrB_free (&CDLP) ; \
72 : }
73 :
74 : #include <LAGraph.h>
75 : #include <LAGraphX.h>
76 : #if LAGRAPH_HAS_STDALIGN_H
77 : #include <stdalign.h>
78 : #define ALIGNAS_64 alignas(64)
79 : #else
80 : // the alignas keyword/macro is not available
81 : #define ALIGNAS_64
82 : #endif
83 : #include "LG_internal.h"
84 :
85 : // A Go-style slice / Lisp-style property list
86 : typedef struct {
87 : GrB_Index* entries;
88 : size_t len, cap;
89 : } plist;
90 :
91 5120 : void plist_free(plist *list) {
92 5120 : free(list->entries); // NOTE: cannot be LAGraph_Free
93 5120 : }
94 :
95 161728512 : void plist_clear(plist *list) {
96 161728512 : list->len = 0;
97 161728512 : }
98 :
99 727057 : void plist_append(plist* list, GrB_Index key, GrB_Index value) {
100 727057 : if (list->len == list->cap) {
101 1177 : size_t new_size = list->cap == 0 ? 16 : 2*list->cap;
102 1177 : list->entries = (GrB_Index*)realloc(list->entries, new_size * sizeof(GrB_Index));
103 1177 : list->cap = new_size;
104 : }
105 727057 : list->entries[list->len] = key;
106 727057 : list->entries[list->len+1] = value;
107 727057 : list->len += 2;
108 727057 : }
109 :
110 17652992 : GrB_Index plist_add(plist* list, GrB_Index entry) {
111 17655030 : for (size_t i = 0; i < list->len; i += 2) {
112 16927973 : if (list->entries[i] == entry) {
113 16925935 : return ++list->entries[i+1];
114 : }
115 : }
116 727057 : plist_append(list, entry, 1);
117 727057 : return 1;
118 : }
119 :
120 : typedef void (*plist_reducer) (GrB_Index* entry1, GrB_Index* count1, GrB_Index entry2, GrB_Index count2);
121 :
122 161728512 : void plist_reduce(plist* list, GrB_Index* entry, GrB_Index* count, plist_reducer reducer) {
123 162455569 : for (size_t i = 0; i < list->len; i += 2) {
124 727057 : reducer(entry, count, list->entries[i], list->entries[i+1]);
125 : }
126 161728512 : }
127 :
128 727057 : void counts_reducer(GrB_Index* e1, GrB_Index* c1, GrB_Index e2, GrB_Index c2) {
129 727057 : if (*c1 > c2) {
130 187542 : return;
131 : }
132 539515 : if (c2 > *c1) {
133 428889 : *e1 = e2;
134 428889 : *c1 = c2;
135 428889 : return;
136 : }
137 110626 : if (*e1 < e2) {
138 93071 : return;
139 : }
140 17555 : *e1 = e2;
141 17555 : *c1 = c2;
142 : }
143 :
144 :
145 : #define bucket_bits 9llu
146 : #define nof_buckets (1llu << bucket_bits)
147 : #define bucket_shift (64llu - bucket_bits)
148 :
149 : typedef struct {
150 : ALIGNAS_64 plist buckets[nof_buckets];
151 : } ptable;
152 :
153 10 : void ptable_free(ptable* table) {
154 5130 : for (size_t i = 0; i < nof_buckets; i++) {
155 5120 : plist_free(&table->buckets[i]);
156 : }
157 10 : }
158 :
159 20 : void ptable_pool_free(ptable* table, size_t n) {
160 20 : if (table == NULL) {
161 10 : return;
162 : }
163 20 : for (size_t i = 0; i < n; i++) {
164 10 : ptable_free(&table[i]);
165 : }
166 10 : free(table); // NOTE: cannot be LAGraph_Free
167 : }
168 :
169 315876 : void ptable_clear(ptable* table) {
170 162044388 : for (size_t i = 0; i < nof_buckets; i++) {
171 161728512 : plist_clear(&table->buckets[i]);
172 : }
173 315876 : }
174 :
175 17652992 : GrB_Index fib_reduce(GrB_Index x)
176 : {
177 : // 2^64 / golden ratio = 11400714819323198485
178 17652992 : GrB_Index fibhash = x * 11400714819323198485llu;
179 : // fast reduce
180 17652992 : return fibhash >> bucket_shift;
181 : }
182 :
183 17652992 : void ptable_add(ptable* table, GrB_Index entry) {
184 17652992 : plist_add(&table->buckets[fib_reduce(entry)], entry);
185 17652992 : }
186 :
187 315876 : void ptable_reduce(ptable* table, GrB_Index* entry, GrB_Index* count, plist_reducer reducer) {
188 315876 : *entry = GrB_INDEX_MAX + 1;
189 315876 : *count = 0;
190 162044388 : for (GrB_Index i = 0; i < nof_buckets; i++) {
191 161728512 : plist_reduce(&table->buckets[i], entry, count, reducer);
192 : }
193 315876 : }
194 :
195 : //****************************************************************************
196 :
197 11 : int LAGraph_cdlp
198 : (
199 : GrB_Vector *CDLP_handle, // output vector
200 : LAGraph_Graph G, // input graph
201 : int itermax, // max number of iterations
202 : char *msg
203 : )
204 : {
205 : GrB_Info info;
206 11 : LG_CLEAR_MSG ;
207 :
208 11 : GrB_Matrix S = NULL, T = NULL ;
209 : GrB_Vector CDLP ;
210 11 : GrB_Index *Sp = NULL, *Si = NULL, *Tp = NULL, *Ti = NULL, *L = NULL, *L_next = NULL ;
211 11 : ptable *counts_pool = NULL ;
212 :
213 : #ifdef _OPENMP
214 : size_t max_threads = omp_get_max_threads();
215 : #else
216 11 : size_t max_threads = 1 ;
217 : #endif
218 :
219 : //--------------------------------------------------------------------------
220 : // check inputs
221 : //--------------------------------------------------------------------------
222 :
223 11 : if (CDLP_handle == NULL)
224 : {
225 1 : return GrB_NULL_POINTER;
226 : }
227 :
228 10 : GrB_Matrix A = G->A ;
229 :
230 : //--------------------------------------------------------------------------
231 : // ensure input is binary and has no self-edges
232 : //--------------------------------------------------------------------------
233 :
234 : GrB_Index n;
235 10 : GRB_TRY (GrB_Matrix_nrows(&n, A)) ;
236 :
237 10 : GRB_TRY (GrB_Matrix_new (&S, GrB_UINT64, n, n)) ;
238 10 : GRB_TRY (GrB_apply (S, GrB_NULL, GrB_NULL, GrB_ONEB_UINT64, A, 0, GrB_NULL)) ;
239 :
240 10 : if (G->kind == LAGraph_ADJACENCY_DIRECTED) {
241 10 : GRB_TRY (GrB_Matrix_new (&T, GrB_UINT64, n, n)) ;
242 10 : GRB_TRY (GrB_transpose (T, GrB_NULL, GrB_NULL, S, GrB_NULL)) ;
243 10 : void * Tx = NULL ;
244 : GrB_Index Tps, Tis, Txs ;
245 : #if LAGRAPH_SUITESPARSE
246 : bool Tiso, Tjumbled ;
247 10 : GRB_TRY (GxB_Matrix_unpack_CSR (T, &Tp, &Ti, &Tx, &Tps, &Tis, &Txs, &Tiso, &Tjumbled, GrB_NULL)) ;
248 : #else
249 : GRB_TRY (GrB_Matrix_exportSize (&Tps, &Tis, &Txs, GrB_CSR_FORMAT, T)) ;
250 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Tp, Tps, sizeof(GrB_Index), NULL)) ;
251 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Ti, Tis, sizeof(GrB_Index), NULL)) ;
252 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Tx, Txs, sizeof(GrB_UINT64), NULL)) ;
253 : GRB_TRY (GrB_Matrix_export (Tp, Ti, (uint64_t *)Tx, &Tps, &Tis, &Txs, GrB_CSR_FORMAT, T)) ;
254 : #endif
255 10 : LAGRAPH_TRY (LAGraph_Free ((void *)&Tx, NULL)) ;
256 10 : GRB_TRY (GrB_free (&T)) ;
257 : }
258 :
259 : {
260 10 : void * Sx = NULL ;
261 : GrB_Index Sps, Sis, Sxs ;
262 : #if LAGRAPH_SUITESPARSE
263 : bool Siso, Sjumbled ;
264 10 : GRB_TRY (GxB_Matrix_unpack_CSR (S, &Sp, &Si, &Sx, &Sps, &Sis, &Sxs, &Siso, &Sjumbled, GrB_NULL)) ;
265 : #else
266 : GRB_TRY (GrB_Matrix_exportSize (&Sps, &Sis, &Sxs, GrB_CSR_FORMAT, S)) ;
267 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Sp, Sps, sizeof(GrB_Index), NULL)) ;
268 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Si, Sis, sizeof(GrB_Index), NULL)) ;
269 : LAGRAPH_TRY (LAGraph_Malloc ((void *)&Sx, Sxs, sizeof(GrB_UINT64), NULL)) ;
270 : GRB_TRY (GrB_Matrix_export (Sp, Si, (uint64_t *)Sx, &Sps, &Sis, &Sxs, GrB_CSR_FORMAT, S)) ;
271 : #endif
272 10 : LAGRAPH_TRY (LAGraph_Free((void *)&Sx, NULL)) ;
273 10 : GRB_TRY (GrB_free (&S)) ;
274 : }
275 :
276 10 : LG_TRY (LAGraph_Malloc ((void **) &L, n, sizeof (GrB_Index), msg)) ;
277 :
278 3304 : for (GrB_Index i = 0; i < n; i++) {
279 3294 : L[i] = i ;
280 : }
281 10 : LG_TRY (LAGraph_Malloc ((void **) &L_next, n, sizeof (GrB_Index), msg)) ;
282 :
283 10 : counts_pool = calloc(max_threads, sizeof(ptable));
284 :
285 353 : for (int iteration = 0; iteration < itermax; iteration++) {
286 :
287 : int64_t i;
288 : #pragma omp parallel for schedule(dynamic)
289 316226 : for (i = 0; i < n; i++) {
290 : #ifdef _OPENMP
291 : int thread_id = omp_get_thread_num() ;
292 : #else
293 315876 : int thread_id = 0 ;
294 : #endif
295 315876 : ptable *counts = &counts_pool [thread_id] ;
296 315876 : GrB_Index* neighbors = Si + Sp[i] ;
297 315876 : GrB_Index sz = Sp[i+1] - Sp[i] ;
298 9142372 : for (GrB_Index j = 0; j < sz; j++) {
299 8826496 : ptable_add (counts, L[neighbors[j]]) ;
300 : }
301 315876 : if (G->kind == LAGraph_ADJACENCY_DIRECTED) {
302 315876 : neighbors = Ti + Tp[i] ;
303 315876 : sz = Tp[i+1] - Tp[i] ;
304 9142372 : for (GrB_Index j = 0; j < sz; j++) {
305 8826496 : ptable_add (counts, L[neighbors[j]]) ;
306 : }
307 : }
308 : GrB_Index best_label, best_count ;
309 315876 : ptable_reduce (counts, &best_label, &best_count, counts_reducer) ;
310 315876 : L_next[i] = best_label ;
311 315876 : ptable_clear (counts) ;
312 : }
313 :
314 350 : GrB_Index* tmp = L ; L = L_next ; L_next = tmp ;
315 350 : bool changed = false ;
316 47967 : for (GrB_Index i = 0; i < n; i++) {
317 47960 : if (L[i] != L_next[i]) {
318 343 : changed = true ;
319 343 : break ;
320 : }
321 : }
322 350 : if (!changed) {
323 7 : break ;
324 : }
325 : }
326 :
327 10 : ptable_pool_free(counts_pool, max_threads); counts_pool = NULL;
328 :
329 : //--------------------------------------------------------------------------
330 : // extract final labels to the result vector
331 : //--------------------------------------------------------------------------
332 :
333 10 : GRB_TRY (GrB_Vector_new(&CDLP, GrB_UINT64, n))
334 3304 : for (GrB_Index i = 0; i < n; i++)
335 : {
336 3294 : GrB_Index l = L[i];
337 : // if (l == GrB_INDEX_MAX + 1) { l = i ; }
338 3294 : l = (l == GrB_INDEX_MAX + 1) ? i : l ;
339 3294 : GRB_TRY (GrB_Vector_setElement(CDLP, l, i))
340 : }
341 :
342 : //--------------------------------------------------------------------------
343 : // free workspace and return result
344 : //--------------------------------------------------------------------------
345 :
346 10 : (*CDLP_handle) = CDLP;
347 10 : CDLP = NULL; // set to NULL so LG_FREE_ALL doesn't free it
348 10 : LG_FREE_ALL;
349 :
350 10 : return (GrB_SUCCESS);
351 : }
|
