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