Line data Source code
1 : //------------------------------------------------------------------------------
2 : // LAGraph_VertexCentrality_triangle: vertex triangle-centrality
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 Tim Davis, Texas A&M University.
15 :
16 : //------------------------------------------------------------------------------
17 :
18 : // TODO: ready for src, except remove all methods except method 3.
19 :
20 : // LAGraph_VertexCentrality_Triangle: computes the TriangleCentrality of
21 : // an undirected graph. No self edges are allowed on the input graph.
22 : // Methods 2 and 3 can tolerate any edge weights (they are ignored; only the
23 : // structure of G->A is used). Methods 1 and 1.5 require unit edge weights
24 : // (this could be modified); results are undefined if this condition doesn't
25 : // hold.
26 :
27 : // Paul Burkhardt. 2024. Triangle Centrality. ACM Trans. Knowl. Discov. Data
28 : // 18, 9, Article 217 (November 2024), 34 pages.
29 : // https://doi.org/10.1145/3685677
30 :
31 : // Method 3 is by far the fastest.
32 :
33 : // This method uses pure GrB* methods from the v2.0 C API only.
34 : // It does not rely on any SuiteSparse:GraphBLAS extensions.
35 :
36 : // TC0: in python (called TC1 in the first draft of the paper)
37 : //
38 : // def triangle_centrality1(A):
39 : // T = A.mxm(A, mask=A)
40 : // y = T.reduce_vector()
41 : // k = y.reduce_float()
42 : // return(1/k)*(3*(A @ y) - 2*(T @ y) + y)
43 : // note: T@y is wrong. should be plus_second semiring
44 :
45 : // def TC1(A):
46 : // # this was "Method 1.5" in a draft, note the T.one@y is now correct:
47 : // T = A.mxm(A, mask=A, desc=ST1)
48 : // y = T.reduce_vector()
49 : // k = y.reduce_float()
50 : // return (3 * (A @ y) - 2 * (T.one() @ y) + y) / k
51 :
52 : // def TC2(A):
53 : // # this was TC2 in the first submission
54 : // T = A.plus_pair(A, mask=A, desc=ST1)
55 : // y = Vector.dense(FP64, A.nrows)
56 : // T.reduce_vector(out=y, accum=FP64.plus)
57 : // k = y.reduce_float()
58 : // return (3 * A.plus_second(y) - 2 * T.plus_second(y) + y) / k
59 :
60 : // def TC3(A):
61 : // M = A.tril(-1)
62 : // T = A.plus_pair(A, mask=M, desc=ST1)
63 : // y = T.reduce() + T.reduce(desc=ST0)
64 : // k = y.reduce_float()
65 : // return (
66 : // 3 * A.plus_second(y) -
67 : // (2 * (T.plus_second(y) + T.plus_second(y, desc=ST0))) + y
68 : // ) / k
69 :
70 : //------------------------------------------------------------------------------
71 :
72 : #define LG_FREE_WORK \
73 : { \
74 : GrB_free (&T) ; \
75 : GrB_free (&u) ; \
76 : GrB_free (&w) ; \
77 : GrB_free (&y) ; \
78 : GrB_free (&L) ; \
79 : }
80 :
81 : #define LG_FREE_ALL \
82 : { \
83 : LG_FREE_WORK ; \
84 : GrB_free (centrality) ; \
85 : }
86 :
87 : #include "LG_internal.h"
88 :
89 : //------------------------------------------------------------------------------
90 : // LAGraph_VertexCentrality_Triangle: vertex triangle-centrality
91 : //------------------------------------------------------------------------------
92 :
93 49 : int LAGraph_VertexCentrality_Triangle // vertex triangle-centrality
94 : (
95 : // outputs:
96 : GrB_Vector *centrality, // centrality(i): triangle centrality of i
97 : uint64_t *ntriangles, // # of triangles in the graph
98 : // inputs:
99 : int method, // 0, 1, 2, or 3
100 : LAGraph_Graph G, // input graph
101 : char *msg
102 : )
103 : {
104 :
105 : //--------------------------------------------------------------------------
106 : // check inputs
107 : //--------------------------------------------------------------------------
108 :
109 49 : LG_CLEAR_MSG ;
110 49 : GrB_Matrix T = NULL, L = NULL, A = NULL ;
111 49 : GrB_Vector y = NULL, u = NULL, w = NULL ;
112 :
113 49 : LG_ASSERT (centrality != NULL && ntriangles != NULL, GrB_NULL_POINTER) ;
114 48 : (*centrality) = NULL ;
115 48 : LG_TRY (LAGraph_CheckGraph (G, msg)) ;
116 :
117 47 : if (G->kind == LAGraph_ADJACENCY_UNDIRECTED ||
118 10 : (G->kind == LAGraph_ADJACENCY_DIRECTED &&
119 10 : G->is_symmetric_structure == LAGraph_TRUE))
120 : {
121 : // the structure of A is known to be symmetric
122 46 : A = G->A ;
123 : }
124 : else
125 : {
126 : // A is not known to be symmetric
127 1 : LG_ASSERT_MSG (false, -1005, "G->A must be symmetric") ;
128 : }
129 :
130 : // no self edges can be present
131 46 : LG_ASSERT_MSG (G->nself_edges == 0, -1004, "G->nself_edges must be zero") ;
132 :
133 : //--------------------------------------------------------------------------
134 : // create the T matrix
135 : //--------------------------------------------------------------------------
136 :
137 : GrB_Index n ;
138 45 : GRB_TRY (GrB_Matrix_nrows (&n, A)) ;
139 45 : GRB_TRY (GrB_Matrix_new (&T, GrB_FP64, n, n)) ;
140 45 : double k = 0 ;
141 :
142 : //--------------------------------------------------------------------------
143 : // compute the Triangle Centrality
144 : //--------------------------------------------------------------------------
145 :
146 45 : if (method == 0 || method == 1)
147 27 : {
148 :
149 : //----------------------------------------------------------------------
150 : // TC0, TC1: simplest method, requires that A has all entries equal to 1
151 : //----------------------------------------------------------------------
152 :
153 : // todo: remove this method when moving this code from experimental/
154 : // to src/
155 :
156 27 : if (method == 0)
157 : {
158 : // T<A> = A*A : method 0 (was TC1 in the first paper submission)
159 18 : GRB_TRY (GrB_mxm (T, A, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, A, A,
160 : NULL)) ;
161 : }
162 : else
163 : {
164 : // this is faster than method 0
165 : // T<A> = A*A' : method TC1 (was method TC1.5)
166 9 : GRB_TRY (GrB_mxm (T, A, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, A, A,
167 : GrB_DESC_T1)) ;
168 : }
169 :
170 : // y = sum (T), where y(i) = sum (T (i,:)) and y(i)=0 of T(i,:) is empty
171 27 : GRB_TRY (GrB_Vector_new (&y, GrB_FP64, n)) ;
172 27 : GRB_TRY (GrB_reduce (y, NULL, NULL, GrB_PLUS_MONOID_FP64, T, NULL)) ;
173 :
174 : // k = sum (y)
175 27 : GRB_TRY (GrB_reduce (&k, NULL, GrB_PLUS_MONOID_FP64, y, NULL)) ;
176 :
177 : // T = spones (T)
178 27 : GRB_TRY (GrB_assign (T, T, NULL, (double) 1, GrB_ALL, n, GrB_ALL, n,
179 : GrB_DESC_S)) ;
180 :
181 : // centrality = (3*A*y - 2*T*y + y) / k
182 :
183 : // w = T*y
184 27 : GRB_TRY (GrB_Vector_new (&w, GrB_FP64, n)) ;
185 27 : GRB_TRY (GrB_mxv (w, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, T, y,
186 : NULL)) ;
187 :
188 : // w = (-2)*w
189 27 : double minus_two = -2 ;
190 27 : GRB_TRY (GrB_apply (w, NULL, NULL, GrB_TIMES_FP64, minus_two, w,
191 : NULL)) ;
192 :
193 : // u = A*y
194 27 : GRB_TRY (GrB_Vector_new (&u, GrB_FP64, n)) ;
195 27 : GRB_TRY (GrB_mxv (u, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, A, y,
196 : NULL)) ;
197 :
198 : }
199 18 : else if (method == 2)
200 : {
201 :
202 : //----------------------------------------------------------------------
203 : // TC2: using LAGraph_plus_one_fp64 semiring
204 : //----------------------------------------------------------------------
205 :
206 : // todo: remove this method when moving this code from experimental/
207 : // to src/
208 :
209 : // T{A} = A*A' (each triangle is seen 6 times)
210 9 : GRB_TRY (GrB_mxm (T, A, NULL, LAGraph_plus_one_fp64, A, A,
211 : GrB_DESC_ST1)) ;
212 :
213 : // y = sum (T), where y(i) = sum (T (i,:)) and y(i)=0 of T(i,:) is empty
214 9 : GRB_TRY (GrB_Vector_new (&y, GrB_FP64, n)) ;
215 9 : GRB_TRY (GrB_assign (y, NULL, NULL, ((double) 0), GrB_ALL, n, NULL)) ;
216 9 : GRB_TRY (GrB_reduce (y, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, T,
217 : NULL)) ;
218 :
219 : // k = sum (y)
220 9 : GRB_TRY (GrB_reduce (&k, NULL, GrB_PLUS_MONOID_FP64, y, NULL)) ;
221 :
222 : // centrality = (3*A*y - 2*T*y + y) / k
223 :
224 : // w = T*y
225 9 : GRB_TRY (GrB_Vector_new (&w, GrB_FP64, n)) ;
226 9 : GRB_TRY (GrB_mxv (w, NULL, NULL, LAGraph_plus_second_fp64, T, y,
227 : NULL)) ;
228 :
229 : // w = (-2)*w
230 9 : double minus_two = -2 ;
231 9 : GRB_TRY (GrB_apply (w, NULL, NULL, GrB_TIMES_FP64, minus_two, w,
232 : NULL)) ;
233 :
234 : // u = A*y
235 9 : GRB_TRY (GrB_Vector_new (&u, GrB_FP64, n)) ;
236 9 : GRB_TRY (GrB_mxv (u, NULL, NULL, LAGraph_plus_second_fp64, A, y,
237 : NULL)) ;
238 :
239 : }
240 9 : else if (method == 3)
241 : {
242 :
243 : //----------------------------------------------------------------------
244 : // TC3: using tril. This is the fastest method.
245 : //----------------------------------------------------------------------
246 :
247 : // todo: When this method is moved to src/, keep this method only.
248 :
249 : // L = tril (A,-1)
250 9 : GRB_TRY (GrB_Matrix_new (&L, GrB_FP64, n, n)) ;
251 9 : GRB_TRY (GrB_select (L, NULL, NULL, GrB_TRIL, A, (int64_t) (-1),
252 : NULL)) ;
253 :
254 : // T{L}= A*A' (each triangle is seen 3 times; T is lower triangular)
255 9 : GRB_TRY (GrB_mxm (T, L, NULL, LAGraph_plus_one_fp64, A, A,
256 : GrB_DESC_ST1)) ;
257 9 : GRB_TRY (GrB_free (&L)) ;
258 :
259 : // y = sum (T'), where y(j) = sum (T (:,j)) and y(j)=0 if T(:,j) empty
260 9 : GRB_TRY (GrB_Vector_new (&y, GrB_FP64, n)) ;
261 9 : GRB_TRY (GrB_assign (y, NULL, NULL, ((double) 0), GrB_ALL, n, NULL)) ;
262 9 : GRB_TRY (GrB_reduce (y, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, T,
263 : GrB_DESC_T0)) ;
264 : // y += sum (T)
265 9 : GRB_TRY (GrB_reduce (y, NULL, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, T,
266 : NULL)) ;
267 :
268 : // k = sum (y). y is the same as the other methods, above, just
269 : // computed using the lower triangular matrix T. So k/6 is the total
270 : // number of triangles in the graph.
271 9 : GRB_TRY (GrB_reduce (&k, NULL, GrB_PLUS_MONOID_FP64, y, NULL)) ;
272 :
273 : // centrality = (3*A*y - 2* (T*y + T'*y) + y) / k
274 :
275 : // w = T*y
276 9 : GRB_TRY (GrB_Vector_new (&w, GrB_FP64, n)) ;
277 9 : GRB_TRY (GrB_mxv (w, NULL, NULL, LAGraph_plus_second_fp64, T, y,
278 : NULL)) ;
279 : // w += T'*y
280 9 : GRB_TRY (GrB_mxv (w, NULL, GrB_PLUS_FP64, LAGraph_plus_second_fp64,
281 : T, y, GrB_DESC_T0)) ;
282 :
283 : // w = (-2)*w
284 9 : double minus_two = -2 ;
285 9 : GRB_TRY (GrB_apply (w, NULL, NULL, GrB_TIMES_FP64, minus_two, w,
286 : NULL)) ;
287 :
288 : // u = A*y
289 9 : GRB_TRY (GrB_Vector_new (&u, GrB_FP64, n)) ;
290 9 : GRB_TRY (GrB_mxv (u, NULL, NULL, LAGraph_plus_second_fp64, A, y,
291 : NULL)) ;
292 :
293 : }
294 :
295 : //--------------------------------------------------------------------------
296 : // centrality = (3*u + w + y) / k for all 4 methods
297 : //--------------------------------------------------------------------------
298 :
299 : // centrality = 3*u
300 45 : GRB_TRY (GrB_Vector_new (centrality, GrB_FP64, n)) ;
301 45 : const double three = 3 ;
302 45 : GRB_TRY (GrB_apply (*centrality, NULL, NULL, GrB_TIMES_FP64, three, u,
303 : NULL)) ;
304 :
305 : // centrality += (w + y)
306 45 : GRB_TRY (GrB_eWiseAdd (*centrality, NULL, GrB_PLUS_FP64, GrB_PLUS_FP64,
307 : w, y, NULL)) ;
308 :
309 : // centrality = centrality / k
310 45 : GRB_TRY (GrB_apply (*centrality, NULL, NULL, GrB_TIMES_FP64,
311 : ((k == 0) ? 1.0 : (1.0/k)), *centrality, NULL)) ;
312 :
313 45 : (*ntriangles) = (uint64_t) (k/6) ; // # triangles is k/6 for all methods
314 :
315 : //--------------------------------------------------------------------------
316 : // free workspace and return result
317 : //--------------------------------------------------------------------------
318 :
319 45 : LG_FREE_WORK ;
320 45 : return (GrB_SUCCESS) ;
321 : }
|
