diff --git a/experimental/algorithm/LAGr_ClosenessCentrality.c b/experimental/algorithm/LAGr_ClosenessCentrality.c new file mode 100644 index 0000000000..409a977ca3 --- /dev/null +++ b/experimental/algorithm/LAGr_ClosenessCentrality.c @@ -0,0 +1,302 @@ +//------------------------------------------------------------------------------ +// LAGr_ClosenessCentrality: Closeness centrality algorithm +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Karan Bhalla and Timothy A. Davis, Texas A&M University + +//------------------------------------------------------------------------------ + +#define LG_FREE_WORK \ + { \ + GrB_free(&distance_vector); \ + GrB_free(&bfs_level); \ + GrB_free(&inf_scalar); \ + GrB_free(&dist_sums); \ + GrB_free(&reachable_counts); \ + GrB_free(&D_APSP); \ + GrB_free(&A_work); \ + GrB_free(&x); \ + LAGraph_Delete(&G_AT, NULL); \ + LAGraph_Free((void **)&source_indices, NULL); \ + } + +#define LG_FREE_ALL \ + { \ + LG_FREE_WORK; \ + GrB_free(¢rality_vector); \ + } + +#include +#include "LG_internal.h" + +//------------------------------------------------------------------------------ +// LAGr_ClosenessCentrality +//------------------------------------------------------------------------------ + +int LAGr_ClosenessCentrality( + // output: + GrB_Vector *centrality, + // input: + LAGraph_Graph G, + GrB_Vector sources, // nodes to score; NULL or empty => all nodes + bool use_weights, // if true, use edge weights in shortest paths + cc_algo_t algorithm, // shortest-path algorithm to use + GrB_Scalar Delta, // delta for SSSP; if NULL, derived from G->emin + char *msg) +{ + + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + LG_CLEAR_MSG; + GrB_Info info; + GrB_Vector centrality_vector = NULL; + GrB_Vector distance_vector = NULL; // per-node BF / SSSP result + GrB_Vector bfs_level = NULL; // BFS: level (distance) of each node + GrB_Scalar inf_scalar = NULL; // INFINITY threshold for SSSP pruning + GrB_Vector dist_sums = NULL; // FW: row sums of D + GrB_Vector reachable_counts = NULL; // FW: row non-zero counts of D + GrB_Matrix D_APSP = NULL; // FW: APSP result + GrB_Matrix A_work = NULL; // AT copy for G_AT + LAGraph_Graph G_AT = NULL; // temporary graph wrapping AT (directed only) + GrB_Vector x = NULL; // FW: dense all-zeros vector for row counting + GrB_Index *source_indices = NULL; + GrB_Index n = 0; + GrB_Index source_count = 0; + + LG_ASSERT(centrality != NULL, GrB_NULL_POINTER); + (*centrality) = NULL; + LG_TRY(LAGraph_CheckGraph(G, msg)); + + GRB_TRY(GrB_Matrix_nrows(&n, G->A)); + + bool use_all_nodes = (sources == NULL); + if (!use_all_nodes) + { + GRB_TRY(GrB_Vector_nvals(&source_count, sources)); + use_all_nodes = (source_count == 0); + } + + cc_algo_t algo = algorithm; + + // set up G_AT: graph over the incoming adjacency. Running any shortest-path + // algorithm from v on G_AT traverses incoming edges, giving distances to + // v in G. + + bool is_directed = (G->kind != LAGraph_ADJACENCY_UNDIRECTED && + G->is_symmetric_structure != LAGraph_TRUE); + + if (is_directed) + { + LG_ASSERT_MSG(G->AT != NULL, LAGRAPH_NOT_CACHED, "G->AT is required"); + GRB_TRY(GrB_Matrix_dup(&A_work, G->AT)); + LG_TRY(LAGraph_New(&G_AT, &A_work, + LAGraph_ADJACENCY_DIRECTED, msg)); + + if (G->emin != NULL) + { + GRB_TRY(GrB_Scalar_dup(&G_AT->emin, G->emin)); + G_AT->emin_state = G->emin_state; + } + } + + // G_trav->A is the incoming adjacency used by all algorithm paths. + LAGraph_Graph G_trav = (G_AT != NULL) ? G_AT : G; + + //-------------------------------------------------------------------------- + // build source node list + //-------------------------------------------------------------------------- + + if (!use_all_nodes) + { + LG_TRY(LAGraph_Malloc((void **)&source_indices, + source_count, sizeof(GrB_Index), msg)); + GRB_TRY(GrB_Vector_extractTuples_UINT64(source_indices, NULL, + &source_count, sources)); + for (GrB_Index k = 0; k < source_count; k++) + { + LG_ASSERT(source_indices[k] < n, GrB_INVALID_INDEX); + } + } + else + { + source_count = n; + } + + GRB_TRY(GrB_Vector_new(¢rality_vector, GrB_FP64, n)); + + //========================================================================== + // Floyd-Warshall path + //========================================================================== + + if (algo == CC_FLOYD_WARSHALL) + { + GrB_Type D_type; + GRB_TRY(LAGraph_FW(G_trav->A, &D_APSP, &D_type)); + + // Remove self-loops so they do not bias sums or counts. + GRB_TRY(GrB_select(D_APSP, NULL, NULL, + GrB_OFFDIAG, D_APSP, (int64_t)0, NULL)); + + // dist_sums[v] = sum of row v + GRB_TRY(GrB_Vector_new(&dist_sums, GrB_FP64, n)); + GRB_TRY(GrB_reduce(dist_sums, NULL, NULL, + GrB_PLUS_MONOID_FP64, D_APSP, NULL)); + + // reachable_counts[v] = number of reachable nodes from v (excluding self). + // Uses plus-one semiring to count structural non-zeros per row. + GRB_TRY(GrB_Vector_new(&x, GrB_FP64, n)); + GRB_TRY(GrB_assign(x, NULL, NULL, (double)0, GrB_ALL, n, NULL)); + GRB_TRY(GrB_Vector_new(&reachable_counts, GrB_FP64, n)); + GRB_TRY(GrB_mxv(reachable_counts, NULL, NULL, LAGraph_plus_one_fp64, + D_APSP, x, NULL)); + GRB_TRY(GrB_free(&D_APSP)); + + // centrality[v] = R(v) / dist_sums[v]. + GRB_TRY(GrB_eWiseMult(centrality_vector, NULL, NULL, GrB_DIV_FP64, + reachable_counts, dist_sums, NULL)); + + (*centrality) = centrality_vector; + LG_FREE_WORK; + return GrB_SUCCESS; + } + + // allocate output vector + if (use_all_nodes) + { + // Dense output: isolated nodes keep score 0. + GRB_TRY(GrB_assign(centrality_vector, NULL, NULL, + 0.0, GrB_ALL, n, NULL)); + } + + if (algo == CC_SSSP) + { + GRB_TRY(GrB_Scalar_new(&inf_scalar, GrB_FP64)); + GRB_TRY(GrB_Scalar_setElement_FP64(inf_scalar, (double)INFINITY)); + } + + //========================================================================== + // per-node shortest-path loop + //========================================================================== + + for (GrB_Index k = 0; k < source_count; k++) + { + GrB_Index node_to_score = + use_all_nodes ? k : source_indices[k]; + + double distance_sum = 0; + GrB_Index reachable_count = 0; + + //---------------------------------------------------------------------- + // BFS (unweighted shortest paths) + //---------------------------------------------------------------------- + + if (algo == CC_BFS) + { + // IMPORTANT: LAGr_BreadthFirstSearch sets *level = NULL without + // freeing the prior handle. We must free bfs_level before every + // call to avoid a memory leak. + GRB_TRY(GrB_free(&bfs_level)); + LG_TRY(LAGr_BreadthFirstSearch(&bfs_level, NULL, + G_trav, node_to_score, msg)); + + GRB_TRY(GrB_Vector_removeElement(bfs_level, node_to_score)); + GRB_TRY(GrB_Vector_nvals(&reachable_count, bfs_level)); + if (reachable_count > 0) + { + GRB_TRY(GrB_reduce(&distance_sum, NULL, + GrB_PLUS_MONOID_FP64, bfs_level, NULL)); + } + } + + //---------------------------------------------------------------------- + // Delta-stepping SSSP (Dijkstra-like, non-negative weights) + //---------------------------------------------------------------------- + + else if (algo == CC_SSSP) + { + GRB_TRY(GrB_free(&distance_vector)); + LG_TRY(LAGr_SingleSourceShortestPath(&distance_vector, + G_trav, + node_to_score, + Delta, msg)); + + // Prune unreachable nodes (value == INFINITY). + GRB_TRY(GrB_select(distance_vector, NULL, NULL, + GrB_VALUELT_FP64, + distance_vector, inf_scalar, NULL)); + + // Exclude self (distance = 0). + GRB_TRY(GrB_Vector_removeElement(distance_vector, + node_to_score)); + GRB_TRY(GrB_Vector_nvals(&reachable_count, distance_vector)); + + if (reachable_count > 0) + { + GRB_TRY(GrB_reduce(&distance_sum, NULL, + GrB_PLUS_MONOID_FP64, + distance_vector, NULL)); + } + } + + //---------------------------------------------------------------------- + // Bellman-Ford + //---------------------------------------------------------------------- + + else // CC_BELLMAN_FORD + { + GRB_TRY(GrB_free(&distance_vector)); + + info = LAGraph_BF_basic(&distance_vector, + G_trav->A, node_to_score); + if (info == GrB_NO_VALUE) + { + // TODO: revisit + // Negative-weight cycle reachable from this node; + // Set centrality to NaN to indicate an invalid score, and continue. + GRB_TRY(GrB_Vector_setElement(centrality_vector, + (double)NAN, node_to_score)); + continue; + } + GRB_TRY(info); + + // Exclude self (distance = 0). + GRB_TRY(GrB_Vector_removeElement(distance_vector, + node_to_score)); + GRB_TRY(GrB_Vector_nvals(&reachable_count, distance_vector)); + + if (reachable_count > 0) + { + GRB_TRY(GrB_reduce(&distance_sum, NULL, + GrB_PLUS_MONOID_FP64, + distance_vector, NULL)); + } + } + + // closeness(v) = R(v) / sum_{u->v} d(u, v) + if (reachable_count > 0 && distance_sum > 0) + { + double closeness = ((double)reachable_count) / distance_sum; + GRB_TRY(GrB_Vector_setElement(centrality_vector, + closeness, node_to_score)); + } + } + + //-------------------------------------------------------------------------- + // free workspace and return result + //-------------------------------------------------------------------------- + + (*centrality) = centrality_vector; + LG_FREE_WORK; + return GrB_SUCCESS; +} \ No newline at end of file diff --git a/experimental/algorithm/LAGr_KatzCentrality.c b/experimental/algorithm/LAGr_KatzCentrality.c new file mode 100644 index 0000000000..b4a1f64466 --- /dev/null +++ b/experimental/algorithm/LAGr_KatzCentrality.c @@ -0,0 +1,167 @@ +//------------------------------------------------------------------------------ +// LAGr_KatzCentrality: Katz centrality algorithm +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Karan Bhalla and Timothy A. Davis, Texas A&M University; + +//------------------------------------------------------------------------------ + +#define LG_FREE_WORK \ +{ \ + GrB_free(&x_prev); \ + GrB_free(&b); \ + GrB_free(&t); \ +} + +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK; \ + GrB_free(&x); \ +} + +#include "LG_internal.h" +#include + +int LAGr_KatzCentrality +( + // output: + GrB_Vector *centrality, + int64_t *iters, + // input: + LAGraph_Graph G, + double alpha, + double beta, + int64_t max_iter, + double tol, + bool normalize, + bool use_weights, // uses weights if true, else treats all weights as 1 + char *msg +) +{ + //-------------------------------------------------------------------------- + // check inputs + //-------------------------------------------------------------------------- + + LG_CLEAR_MSG; + GrB_Index n = 0; + GrB_Vector x = NULL, x_prev = NULL, b = NULL, t = NULL; + + LG_ASSERT (centrality != NULL && iters != NULL, GrB_NULL_POINTER) ; + (*centrality) = NULL; + LG_TRY (LAGraph_CheckGraph (G, msg)) ; + + GrB_Matrix AT; + if (G->kind == LAGraph_ADJACENCY_UNDIRECTED || + G->is_symmetric_structure == LAGraph_TRUE) + { + AT = G->A ; + } + else + { + AT = G->AT ; + LG_ASSERT_MSG (AT != NULL, LAGRAPH_NOT_CACHED, "G->AT is required") ; + } + + // compute correct semiring based on whether edge weights are used + // TODO: add FP32 support + GrB_Semiring semiring = use_weights ? GrB_PLUS_TIMES_SEMIRING_FP64 + : GxB_PLUS_SECOND_FP64 ; + + if (use_weights) + { + // ensure min edge weight is available and nonnegative + LG_TRY (LAGraph_Cached_EMin (G, msg)) ; + LG_ASSERT_MSG ( + G->emin != NULL && + (G->emin_state == LAGraph_VALUE || + G->emin_state == LAGraph_BOUND), + LAGRAPH_NOT_CACHED, "G->emin is required") ; + + double emin = 0 ; + GRB_TRY (GrB_Scalar_extractElement_FP64 (&emin, G->emin)) ; + LG_ASSERT_MSG (emin >= 0, GrB_INVALID_VALUE, + "use_weights=true requires nonnegative edge weights") ; + } + + //-------------------------------------------------------------------------- + // initializations + //-------------------------------------------------------------------------- + + GRB_TRY (GrB_Matrix_nrows (&n, AT)) ; + + GRB_TRY (GrB_Vector_new (&x, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&x_prev, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&b, GrB_FP64, n)) ; + GRB_TRY (GrB_Vector_new (&t, GrB_FP64, n)) ; + + GRB_TRY (GrB_assign (x, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_assign (x_prev, NULL, NULL, 0.0, GrB_ALL, n, NULL)) ; + GRB_TRY (GrB_assign (b, NULL, NULL, beta, GrB_ALL, n, NULL)) ; + + // first iteration is always done + double rdiff = 1 ; + + // TODO: determine best way to stop + for ((*iters) = 0 ; ; (*iters)++) + { + // check for convergence failure + LG_ASSERT_MSGF ((*iters) < max_iter, LAGRAPH_CONVERGENCE_FAILURE, + "katz centrality failed to converge in %" PRId64 " iterations", max_iter) ; + + // swap x and x_prev + GrB_Vector temp = x_prev ; + x_prev = x ; + x = temp ; + + // x = A' * x_prev + GRB_TRY (GrB_mxv (x, NULL, NULL, semiring, AT, x_prev, NULL)) ; + + // x = alpha * x + beta + GRB_TRY (GrB_apply (x, NULL, NULL, GrB_TIMES_FP64, alpha, x, NULL)) ; + GRB_TRY (GrB_eWiseAdd (x, NULL, NULL, GrB_PLUS_FP64, x, b, NULL)) ; + + // t = x - x_prev + GRB_TRY (GrB_eWiseAdd (t, NULL, NULL, GrB_MINUS_FP64, x, x_prev, + NULL)) ; + // t = abs (t) + GRB_TRY (GrB_apply (t, NULL, NULL, GrB_ABS_FP64, t, NULL)) ; + // rdiff = sum (t) + GRB_TRY (GrB_reduce (&rdiff, NULL, GrB_PLUS_MONOID_FP64, t, NULL)) ; + + if (rdiff < tol) break ; // TODO: revisit + } + + // normalize using the L2 norm if flag is set + if (normalize) + { + double sumsq = 0 ; + + // sumsq = x' * x = sum (x .* x) + GRB_TRY (GrB_eWiseMult (t, NULL, NULL, GrB_TIMES_FP64, x, x, NULL)) ; + GRB_TRY (GrB_reduce (&sumsq, NULL, GrB_PLUS_MONOID_FP64, t, NULL)) ; + + if (sumsq > 0) + { + GRB_TRY (GrB_apply (x, NULL, NULL, GrB_DIV_FP64, x, sqrt (sumsq), + NULL)) ; + } + } + + //-------------------------------------------------------------------------- + // free workspace and return result + //-------------------------------------------------------------------------- + + (*centrality) = x; + LG_FREE_WORK; + + return GrB_SUCCESS; +} diff --git a/experimental/algorithm/LAGraph_FW.c b/experimental/algorithm/LAGraph_FW.c new file mode 100644 index 0000000000..2e20151b28 --- /dev/null +++ b/experimental/algorithm/LAGraph_FW.c @@ -0,0 +1,92 @@ +//------------------------------------------------------------------------------ +// LAGraph_FW: Floyd-Warshall method: all pairs shortest paths +//------------------------------------------------------------------------------ + +// The input is a square unsymmetric matrix G, for a directed graph. G can be +// of any type. The output is always of type GrB_FP64, regardless of the input +// type, since shortest path distances are naturally floating point and callers +// such as LAGr_ClosenessCentrality require FP64 output. +// TODO: revisit in future; consider giving the user control over the data type of D to use. + +// G(i,j) is the edge weight for edge (i,j). D(i,j) on output is the length of +// the shortest path from node i to j, if the entry is present. If D(i,j) is +// not present then there is no path from i to j. The shortest path itself +// is not returned. + +// Negative weights are OK, unless there is a negative weight cycle. In +// that case, the output is undefined. + +#define LG_FREE_WORK \ +{ \ + GrB_free (&A) ; \ + GrB_free (&B) ; \ +} + +#define LG_FREE_ALL \ +{ \ + LG_FREE_WORK ; \ + GrB_free (&D_matrix) ; \ +} + +#include "LG_internal.h" +#include "LAGraphX.h" + +GrB_Info LAGraph_FW +( + const GrB_Matrix G, // input graph, with edge weights + GrB_Matrix *D, // output graph, created on output + GrB_Type *D_type // type of D (always GrB_FP64) +) +{ + GrB_Info info ; + char *msg = NULL ; + GrB_Matrix D_matrix = NULL, A = NULL, B = NULL ; + + // make sure outputs and input are valid + LG_ASSERT (G != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (D != NULL, GrB_NULL_POINTER) ; + LG_ASSERT (D_type != NULL, GrB_NULL_POINTER) ; + (*D) = NULL ; + (*D_type) = NULL ; + + // output is always FP64 + GrB_BinaryOp op = GrB_MIN_FP64 ; + GrB_UnaryOp idop = GrB_IDENTITY_FP64 ; + GrB_Semiring semiring = GxB_MIN_PLUS_FP64 ; + + GrB_Index n, ncols ; + GRB_TRY (GrB_Matrix_nrows (&n, G)) ; + GRB_TRY (GrB_Matrix_ncols (&ncols, G)) ; + LG_ASSERT (n == ncols, GrB_INVALID_VALUE) ; + + (*D_type) = GrB_FP64 ; + GRB_TRY (GrB_Matrix_new (&D_matrix, GrB_FP64, n, n)) ; + GRB_TRY (LG_SET_FORMAT_HINT (D_matrix, LG_BITMAP)) ; + GRB_TRY (GrB_Matrix_new (&A, GrB_FP64, n, 1)) ; + GRB_TRY (GrB_Matrix_new (&B, GrB_FP64, 1, n)) ; + + // copy G into D, casting to FP64 + GRB_TRY (GrB_apply (D_matrix, GrB_NULL, GrB_NULL, idop, G, GrB_NULL)) ; + + // Set D(i,i) = 0 for all i. This overrides any self-edge weights in G + for (GrB_Index i = 0 ; i < n ; i++) + { + GRB_TRY (GrB_Matrix_setElement_FP64 (D_matrix, 0.0, i, i)) ; + } + + for (GrB_Index k = 0 ; k < n ; k++) + { + // A = D(:,k), the kth column + GRB_TRY (GrB_extract (A, GrB_NULL, GrB_NULL, D_matrix, + GrB_ALL, n, &k, 1, GrB_NULL)) ; + // B = D(k,:), the kth row + GRB_TRY (GrB_extract (B, GrB_NULL, GrB_NULL, D_matrix, + &k, 1, GrB_ALL, n, GrB_NULL)) ; + // D = min (D, A*B) using the min-plus semiring + GRB_TRY (GrB_mxm (D_matrix, GrB_NULL, op, semiring, A, B, GrB_NULL)) ; + } + + LG_FREE_WORK ; + (*D) = D_matrix ; + return (GrB_SUCCESS) ; +} \ No newline at end of file diff --git a/experimental/benchmark/closenessCentrality_demo.c b/experimental/benchmark/closenessCentrality_demo.c new file mode 100644 index 0000000000..403f757339 --- /dev/null +++ b/experimental/benchmark/closenessCentrality_demo.c @@ -0,0 +1,202 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/benchmark/closenessCentrality_demo.c: demo for Closeness +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2022 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Karan Bhalla and Timothy A. Davis, Texas A&M University + +//------------------------------------------------------------------------------ + +#include "../../src/benchmark/LAGraph_demo.h" +#include "LAGraphX.h" +#include +#include + +#define LG_FREE_ALL \ + { \ + GrB_free (¢rality) ; \ + GrB_free (&sources) ; \ + GrB_free (&Delta) ; \ + LAGraph_Delete (&G, msg) ; \ + } + +int main (int argc, char **argv) +{ + //-------------------------------------------------------------------------- + // startup LAGraph and GraphBLAS + //-------------------------------------------------------------------------- + + char msg [LAGRAPH_MSG_LEN] ; + LAGraph_Graph G = NULL ; + GrB_Vector centrality = NULL, sources = NULL ; + GrB_Scalar Delta = NULL ; + + bool burble = false ; + demo_init (burble) ; + + //-------------------------------------------------------------------------- + // parse command-line arguments + //-------------------------------------------------------------------------- + + // Usage: + // closenessCentrality_demo + // [num_sources] [use_weights] [delta] + // + // : required; path to .mtx file + // : required int; shortest-path algorithm: + // 0 = CC_BFS (unweighted BFS) + // 1 = CC_SSSP (delta-stepping SSSP) + // 2 = CC_BELLMAN_FORD (Bellman-Ford) + // 3 = CC_FLOYD_WARSHALL (all-pairs, no sources) + // [num_sources] : optional; number of random source nodes to score + // (0 or omitted => score all nodes) + // [use_weights] : optional 0/1 (default 0); use edge weights + // [delta] : optional float; delta for CC_SSSP + // (ignored for other algorithms) + + if (argc < 3 || argc > 6) + { + printf ("Usage: %s " + " [num_sources] [use_weights] [delta]\n" + " algorithm: 0=CC_BFS 1=CC_SSSP 2=CC_BELLMAN_FORD" + " 3=CC_FLOYD_WARSHALL\n", + argv [0]) ; + return (GrB_INVALID_VALUE) ; + } + + cc_algo_t algorithm = (cc_algo_t) atoi (argv [2]) ; + bool use_weights = false ; + double delta_val = -1 ; + int num_sources = 0 ; + + if (argc > 3) num_sources = atoi (argv [3]) ; + if (argc > 4) use_weights = (atoi (argv [4]) != 0) ; + if (argc > 5) delta_val = strtod (argv [5], NULL) ; + + //-------------------------------------------------------------------------- + // read in the graph + //-------------------------------------------------------------------------- + + // char *matrix_name = (argc > 1) ? argv [1] : "stdin" ; + double t = LAGraph_WallClockTime () ; + LAGRAPH_TRY( + readproblem(&G, NULL, false, false, false, NULL, false, argc, argv)) ; + + GrB_Index n ; + GRB_TRY (GrB_Matrix_nrows (&n, G->A)) ; + + LAGRAPH_TRY (LAGraph_Cached_AT (G, msg)) ; + + if (use_weights) + { + LAGRAPH_TRY (LAGraph_Cached_EMin (G, msg)) ; + } + + GrB_Index nvals ; + GRB_TRY (GrB_Matrix_nvals (&nvals, G->A)) ; + + t = LAGraph_WallClockTime () - t ; + printf ("Time to read the graph: %g sec\n", t) ; + printf ("Nodes: %" PRIu64 " Edges: %" PRIu64 "\n", + (uint64_t) n, (uint64_t) nvals) ; + + printf ("\n==========================\nThe input graph matrix G:\n") ; + // LAGRAPH_TRY (LAGraph_Graph_Print (G, LAGraph_SHORT, stdout, msg)) ; + + //-------------------------------------------------------------------------- + // build optional source-node vector + //-------------------------------------------------------------------------- + + if (num_sources > 0) + { + if ((uint64_t) num_sources >= n) + { + printf ("Error: num_sources (%" PRId32 ") must be less than" + " n (%" PRIu64 ").\n", num_sources, n) ; + LG_FREE_ALL ; + return (GrB_INVALID_VALUE) ; + } + + GRB_TRY (GrB_Vector_new (&sources, GrB_UINT64, num_sources)) ; + + double t_seed = LAGraph_WallClockTime () ; + srand ((int) t_seed) ; + + bool *used = NULL ; + LAGRAPH_TRY (LAGraph_Calloc ((void **) &used, n, sizeof (bool), msg)) ; + + for (int i = 0 ; i < num_sources ; i++) + { + GrB_Index idx ; + do { idx = rand () % n ; } while (used [idx]) ; + used [idx] = true ; + GRB_TRY (GrB_Vector_setElement (sources, idx, i)) ; + } + LAGRAPH_TRY (LAGraph_Free ((void **) &used, msg)) ; + + printf ("Using %d random source nodes.\n", num_sources) ; + } + + //-------------------------------------------------------------------------- + // compute closeness centrality + //-------------------------------------------------------------------------- + + printf ("\nCloseness params: use_weights=%d algorithm=%d" + " num_sources=%d\n\n", + (int) use_weights, (int) algorithm, num_sources) ; + + if (algorithm == CC_SSSP && delta_val > 0) + { + GRB_TRY (GrB_Scalar_new (&Delta, GrB_FP64)) ; + GRB_TRY (GrB_Scalar_setElement_FP64 (Delta, delta_val)) ; + } + + t = LAGraph_WallClockTime () ; + LAGRAPH_TRY (LAGr_ClosenessCentrality (¢rality, G, sources, + use_weights, algorithm, Delta, + msg)) ; + t = LAGraph_WallClockTime () - t ; + printf ("Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + //-------------------------------------------------------------------------- + // print results + //-------------------------------------------------------------------------- + + GrB_Index scoredvals = 0 ; + GRB_TRY (GrB_Vector_nvals (&scoredvals, centrality)) ; + printf ("Number of scored nodes: %" PRIu64 "\n", (uint64_t) scoredvals) ; + + // print first few scores + // GrB_Index print_limit = (n < 20) ? n : 20 ; + // printf ("Closeness centrality (first %" PRIu64 " nodes):\n", + // (uint64_t) print_limit) ; + // for (GrB_Index i = 0 ; i < print_limit ; i++) + // { + // double score = 0 ; + // GrB_Info info = GrB_Vector_extractElement_FP64 (&score, centrality, i) ; + // if (info == GrB_SUCCESS) + // { + // printf (" node %" PRIu64 ": %g\n", (uint64_t) i, score) ; + // } + // else + // { + // printf (" node %" PRIu64 ": (unreachable)\n", (uint64_t) i) ; + // } + // } + + //-------------------------------------------------------------------------- + // free everything and finish + //-------------------------------------------------------------------------- + + LG_FREE_ALL ; + LAGRAPH_TRY (LAGraph_Finalize (msg)) ; + return (GrB_SUCCESS) ; +} diff --git a/experimental/benchmark/katzCentrality_demo.c b/experimental/benchmark/katzCentrality_demo.c new file mode 100644 index 0000000000..d021bdca84 --- /dev/null +++ b/experimental/benchmark/katzCentrality_demo.c @@ -0,0 +1,157 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/benchmark/katzCentrality_demo.c: benchmark for Katz +//------------------------------------------------------------------------------ + +#include "../../src/benchmark/LAGraph_demo.h" +#include "LAGraphX.h" +#include + +#define NTHREAD_LIST 1 +#define THREAD_LIST 0 + +#define LG_FREE_ALL \ +{ \ + LAGraph_Delete (&G, NULL) ; \ + GrB_free (&c) ; \ +} + +int main (int argc, char **argv) +{ + //-------------------------------------------------------------------------- + // initialize LAGraph and GraphBLAS + //-------------------------------------------------------------------------- + + char msg [LAGRAPH_MSG_LEN] ; + LAGraph_Graph G = NULL ; + GrB_Vector c = NULL ; + + // start GraphBLAS and LAGraph + bool burble = false ; + demo_init (burble) ; + + int ntrials = 3 ; + printf ("# of trials: %d\n", ntrials) ; + + int nt = NTHREAD_LIST ; + int Nthreads [20] = { 0, THREAD_LIST } ; + + int nthreads_max, nthreads_outer, nthreads_inner ; + LAGRAPH_TRY (LAGraph_GetNumThreads (&nthreads_outer, &nthreads_inner, msg)) ; + nthreads_max = nthreads_outer * nthreads_inner ; + + if (Nthreads [1] == 0) + { + Nthreads [1] = nthreads_max ; + for (int t = 2 ; t <= nt ; t++) + { + Nthreads [t] = Nthreads [t-1] / 2 ; + if (Nthreads [t] == 0) + nt = t-1 ; + } + } + + printf ("threads to test:") ; + for (int t = 1 ; t <= nt ; t++) + { + int nthreads = Nthreads [t] ; + if (nthreads > nthreads_max) + continue ; + printf (" %d", nthreads) ; + } + printf ("\n") ; + + //-------------------------------------------------------------------------- + // read in the graph and parse optional Katz parameters + //-------------------------------------------------------------------------- + + // Usage: + // katzCentrality_demo + // [beta] [max_iters] [tol] [normalize] [use_weights] + // where normalize and use_weights are 0/1. + if (argc != 1 && (argc < 3 || argc > 8)) + { + printf ("Usage: %s " + " [beta] [max_iters] [tol] [normalize] [use_weights]\n", + argv [0]) ; + return (GrB_INVALID_VALUE) ; + } + + char *matrix_name = (argc > 1) ? argv [1] : "stdin" ; + + double alpha = 0.01 ; + double beta = 1.00 ; + int64_t iters = 0 ; + int64_t max_iter = 1000 ; + double tol = 1e-6 ; + bool normalize = false ; + bool use_weights = false ; + + if (argc > 2) alpha = strtod (argv [2], NULL) ; + if (argc > 3) beta = strtod (argv [3], NULL) ; + if (argc > 4) max_iter = (int64_t) strtoll (argv [4], NULL, 10) ; + if (argc > 5) tol = strtod (argv [5], NULL) ; + if (argc > 6) normalize = (atoi (argv [6]) != 0) ; + if (argc > 7) use_weights = (atoi (argv [7]) != 0) ; + + LAGRAPH_TRY( + readproblem(&G, NULL, false, false, false, NULL, false, argc, argv)); + + // Katz uses incoming edges; AT is required for directed graphs. + int cache_result = LAGraph_Cached_AT (G, msg) ; + LG_ASSERT_MSG (cache_result == GrB_SUCCESS || + cache_result == LAGRAPH_CACHE_NOT_NEEDED, + cache_result, "LAGraph_Cached_AT failed") ; + + //-------------------------------------------------------------------------- + // benchmark Katz centrality + //-------------------------------------------------------------------------- + + printf ("\nKatz params: alpha=%g beta=%g max_iter=%lld tol=%g" + " normalize=%d use_weights=%d\n\n", + alpha, beta, (long long) max_iter, tol, + (int) normalize, (int) use_weights) ; + + // warmup for more accurate timing + double tt = LAGraph_WallClockTime () ; + LAGRAPH_TRY (LAGr_KatzCentrality (&c, &iters, G, alpha, beta, + max_iter, tol, normalize, use_weights, msg)) ; + tt = LAGraph_WallClockTime () - tt ; + GRB_TRY (GrB_free (&c)) ; + printf ("warmup time %g sec\n", tt) ; + + for (int t = 1 ; t <= nt ; t++) + { + int nthreads = Nthreads [t] ; + if (nthreads > nthreads_max) + continue ; + LAGRAPH_TRY (LAGraph_SetNumThreads (1, nthreads, msg)) ; + + double ttot = 0, ttrial [100] ; + for (int trial = 0 ; trial < ntrials ; trial++) + { + double t1 = LAGraph_WallClockTime () ; + LAGRAPH_TRY (LAGr_KatzCentrality (&c, &iters, G, alpha, beta, + max_iter, tol, normalize, use_weights, msg)) ; + GRB_TRY (GrB_free (&c)) ; + ttrial [trial] = LAGraph_WallClockTime () - t1 ; + ttot += ttrial [trial] ; + + printf ("threads %2d trial %2d: %12.6f sec\n", + nthreads, trial, ttrial [trial]) ; + fprintf (stderr, "threads %2d trial %2d: %12.6f sec\n", + nthreads, trial, ttrial [trial]) ; + } + + ttot = ttot / ntrials ; + printf ( "Avg: KatzCentrality nthreads: %3d time: %12.6f sec" + " (%" PRId64 " iterations) matrix: %s\n", + nthreads, ttot, iters, matrix_name) ; + fprintf (stderr, "Avg: KatzCentrality nthreads: %3d time: %12.6f sec" + " (%" PRId64 " iterations) matrix: %s\n", + nthreads, ttot, iters, matrix_name) ; + } + + LG_FREE_ALL ; + LAGRAPH_TRY (LAGraph_Finalize (msg)) ; + return (GrB_SUCCESS) ; +} diff --git a/experimental/test/test_ClosenessCentrality.c b/experimental/test/test_ClosenessCentrality.c new file mode 100644 index 0000000000..0d8df400b2 --- /dev/null +++ b/experimental/test/test_ClosenessCentrality.c @@ -0,0 +1,398 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/test/test_ClosenessCentrality.c: tests for Closeness +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Karan Bhalla and Timothy A. Davis, Texas A&M University + +//------------------------------------------------------------------------------ + +#include +#include +#include + +#include "LAGraphX.h" +#include "LAGraph_test.h" +#include "LG_internal.h" + +#define LEN 512 +char msg [LAGRAPH_MSG_LEN] ; +char filename [LEN+1] ; + +//------------------------------------------------------------------------------ +// reference results from NetworkX +// note: WF_improved is set to FALSE for all tests +//------------------------------------------------------------------------------ + +double diamonds_closeness[8] = { + 0.0, 1.0, 1.0, 1.0, 0.8, 0.5, 0.5, 0.4117647058823529 +} ; + +double karate_closeness[34] = { + 0.5689655172413793, 0.4852941176470588, 0.559322033898305, 0.4647887323943662, + 0.3793103448275862, 0.38372093023255816, 0.38372093023255816, 0.44, + 0.515625, 0.4342105263157895, 0.3793103448275862, 0.36666666666666664, + 0.3707865168539326, 0.515625, 0.3707865168539326, 0.3707865168539326, + 0.28448275862068967, 0.375, 0.3707865168539326, 0.5, + 0.3707865168539326, 0.375, 0.3707865168539326, 0.39285714285714285, + 0.375, 0.375, 0.3626373626373626, 0.4583333333333333, + 0.4520547945205479, 0.38372093023255816, 0.4583333333333333, 0.5409836065573771, + 0.515625, 0.55 +} ; + +double ldbc_directed_example_closeness[10] = { + 1.6891891892, 0.0000000000, 1.5151515152, 1.3937282230, + 1.8115942029, 0.0000000000, 0.0000000000, 3.2258064516, + 0.0000000000, 1.1928429423 +} ; + +//------------------------------------------------------------------------------ +// difference: compare closeness vector result with reference values +//------------------------------------------------------------------------------ + +double difference (GrB_Vector c, double *reference_c, GrB_Index n) ; + +double difference (GrB_Vector c, double *reference_c, GrB_Index n) +{ + GrB_Vector diff = NULL, reference_c_vector = NULL ; + OK (GrB_Vector_new (&reference_c_vector, GrB_FP64, n)) ; + + for (GrB_Index i = 0 ; i < n ; i++) + { + OK (GrB_Vector_setElement_FP64 (reference_c_vector, reference_c [i], i)) ; + } + + OK (GrB_Vector_new (&diff, GrB_FP64, n)) ; + OK (GrB_eWiseAdd (diff, NULL, NULL, GrB_MINUS_FP64, reference_c_vector, c, + NULL)) ; + OK (GrB_apply (diff, NULL, NULL, GrB_ABS_FP64, diff, NULL)) ; + + double err = 0 ; + OK (GrB_reduce (&err, NULL, GrB_MAX_MONOID_FP64, diff, NULL)) ; + + OK (GrB_free (&diff)) ; + OK (GrB_free (&reference_c_vector)) ; + + return err ; +} + +//------------------------------------------------------------------------------ +// directed and unweighted graph, BFS, all sources +//------------------------------------------------------------------------------ +void test_closeness_diamonds_bfs (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "diamonds.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nDiamonds graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + false, CC_BFS, NULL, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, diamonds_closeness, 8) ; + printf (" diamonds: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// undirected and unweighted graph, BFS, all sources +//------------------------------------------------------------------------------ +void test_closeness_karate_bfs (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nKarate graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + false, CC_BFS, NULL, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, karate_closeness, 34) ; + printf (" karate: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// undirected and unweighted graph, Floyd-Warshall, all sources +//------------------------------------------------------------------------------ +void test_closeness_karate_fw (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nKarate graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + false, CC_FLOYD_WARSHALL, NULL, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, karate_closeness, 34) ; + printf (" karate: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// directed and weighted graph, Floyd-Warshall, all sources +//------------------------------------------------------------------------------ +void test_closeness_ldbc_directed_example_fw (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "ldbc-directed-example.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nldbc directed example graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + true, CC_FLOYD_WARSHALL, NULL, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, ldbc_directed_example_closeness, 10) ; + printf (" ldbc directed example: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// directed and weighted graph, Bellman-Ford, all sources +//------------------------------------------------------------------------------ +void test_closeness_ldbc_directed_example_bf (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "ldbc-directed-example.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nldbc directed example graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + true, CC_BELLMAN_FORD, NULL, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, ldbc_directed_example_closeness, 10) ; + printf (" ldbc directed example: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// directed and weighted graph, SSSP, all sources +//------------------------------------------------------------------------------ +void test_closeness_ldbc_directed_example_sssp (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + + snprintf (filename, LEN, LG_DATA_DIR "%s", "ldbc-directed-example.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + uint64_t n, nedges ; + OK (GrB_Matrix_nrows (&n, G->A)) ; + OK (GrB_Matrix_nvals (&nedges, G->A)) ; + printf ("\n\nldbc directed example graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", + n, nedges) ; + + GrB_Scalar Delta = NULL ; + OK (GrB_Scalar_new (&Delta, GrB_FP64)) ; + OK (GrB_Scalar_setElement_FP64 (Delta, 1.6)) ; + + double t = LAGraph_WallClockTime () ; + OK (LAGr_ClosenessCentrality (¢rality, G, NULL, + true, CC_SSSP, Delta, msg)) ; + t = LAGraph_WallClockTime () - t ; + printf (" Time for LAGr_ClosenessCentrality: %g sec\n", t) ; + + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + + double err = difference (centrality, ldbc_directed_example_closeness, 10) ; + printf (" ldbc directed example: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (&Delta)) ; + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; +#endif +} + +//------------------------------------------------------------------------------ +// list of tests +//------------------------------------------------------------------------------ + +TEST_LIST = { + {"test_closeness_diamonds_bfs", test_closeness_diamonds_bfs}, + {"test_closeness_karate_bfs", test_closeness_karate_bfs}, + {"test_closeness_karate_fw", test_closeness_karate_fw}, + {"test_closeness_ldbc_directed_example_fw", test_closeness_ldbc_directed_example_fw}, + {"test_closeness_ldbc_directed_example_bf", test_closeness_ldbc_directed_example_bf}, + {"test_closeness_ldbc_directed_example_sssp", test_closeness_ldbc_directed_example_sssp}, + {NULL, NULL} +} ; diff --git a/experimental/test/test_KatzCentrality.c b/experimental/test/test_KatzCentrality.c new file mode 100644 index 0000000000..c76d71868a --- /dev/null +++ b/experimental/test/test_KatzCentrality.c @@ -0,0 +1,324 @@ +//------------------------------------------------------------------------------ +// LAGraph/src/test/test_KatzCentrality.c: testing for Katz centrality +// ----------------------------------------------------------------------------- + +// LAGraph, (c) 2019-2025 by The LAGraph Contributors, All Rights Reserved. +// SPDX-License-Identifier: BSD-2-Clause +// +// For additional details (including references to third party source code and +// other files) see the LICENSE file or contact permission@sei.cmu.edu. See +// Contributors.txt for a full list of contributors. Created, in part, with +// funding and support from the U.S. Government (see Acknowledgments.txt file). +// DM22-0790 + +// Contributed by Karan Bhalla and Timothy A. Davis, Texas A&M University + +//------------------------------------------------------------------------------ + +#include +#include +#include + +#include "LAGraphX.h" +#include "LAGraph_test.h" +#include "LG_internal.h" + +#define LEN 512 +char msg [LAGRAPH_MSG_LEN] ; +char filename [LEN+1] ; + +//------------------------------------------------------------------------------ +// reference results from NetworkX +//------------------------------------------------------------------------------ + +double karate_katz [34] = { + 4.9829901679, 3.6518079520, 4.1214049703, 3.0226461989, 1.8904511529, + 2.0310725862, 2.0310725862, 2.5778843937, 3.1126620496, 1.9260736983, + 1.8904511529, 1.4982988513, 1.8005633729, 3.0918177427, 1.9405256847, + 1.9405256847, 1.4062144433, 1.8634795227, 1.9405256847, 2.3774128717, + 1.9405256847, 1.8634795227, 1.9405256847, 2.5865313998, 1.7091149570, + 1.7301052384, 1.7513648061, 2.3556382376, 2.2266143974, 2.3743152002, + 2.6169724552, 3.0054078802, 4.2659247944, 5.1393352269 +} ; + +double diamonds_katz [8] = { + 1.0000000000, 1.1000000000, 1.3200000000, 1.1000000000, 1.3520000000, + 1.1352000000, 1.1352000000, 1.2270400000 +} ; + +double ldbc_directed_example_katz [10] = { + 1.1016778324, 1.0000000000, 1.1563813542, 1.2611415319, 1.1347459707, + 1.0000000000, 1.0000000000, 1.0356314652, 1.0000000000, 1.0721318253 +} ; + +double rand_katz [20] = { + 0.1835353701, 0.2018889071, 0.2018889071, 0.2018889071, 0.2018889071, + 0.2239131515, 0.2239131515, 0.2239131515, 0.2283180003, 0.2283180003, + 0.2283180003, 0.2291989701, 0.2291989701, 0.2291989701, 0.2293751641, + 0.2293751641, 0.2293751641, 0.2294104029, 0.2294104029, 0.2752924834 +} ; + +//------------------------------------------------------------------------------ +// difference: compare Katz vector result with reference values +//------------------------------------------------------------------------------ + +double difference (GrB_Vector c, double *reference_c, GrB_Index n) ; + +double difference (GrB_Vector c, double *reference_c, GrB_Index n) +{ + GrB_Vector diff = NULL, reference_c_vector = NULL ; + OK (GrB_Vector_new (&reference_c_vector, GrB_FP64, n)) ; + + for (GrB_Index i = 0 ; i < n ; i++) + { + OK (GrB_Vector_setElement_FP64 (reference_c_vector, reference_c [i], i)) ; + } + + OK (GrB_Vector_new (&diff, GrB_FP64, n)) ; + OK (GrB_eWiseAdd (diff, NULL, NULL, GrB_MINUS_FP64, reference_c_vector, c, + NULL)) ; + OK (GrB_apply (diff, NULL, NULL, GrB_ABS_FP64, diff, NULL)) ; + + double err = 0 ; + OK (GrB_reduce (&err, NULL, GrB_MAX_MONOID_FP64, diff, NULL)) ; + + OK (GrB_free (&diff)) ; + OK (GrB_free (&reference_c_vector)) ; + + return err ; +} + + +void test_katz_diamonds (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + int64_t niters = 0 ; + + // Create diamonds graph (directed) + snprintf (filename, LEN, LG_DATA_DIR "%s", "diamonds.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + // Check that AT is cached + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + // Print graph stats + uint64_t n, nedges ; + OK (GrB_Matrix_nrows(&n, G->A)) ; + OK (GrB_Matrix_nvals(&nedges, G->A)) ; + printf ("\n\nDiamonds graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", n, nedges) ; + + double alpha = 0.1 ; + + // Compute katz centrality + double t = LAGraph_WallClockTime() ; + OK (LAGr_KatzCentrality (¢rality, &niters, G, alpha, 1.0, 1000, 1e-6, false, false, msg)) ; + t = LAGraph_WallClockTime() - t ; + printf (" Time for LAGr_KatzCentrality: %g sec\n", t) ; + printf (" Iterations for LAGr_KatzCentrality: %" PRId64 "\n", niters) ; + + // Compare with reference values. + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + TEST_CHECK (cnvals == n) ; + + double err = difference (centrality, diamonds_katz, 8) ; + printf (" diamonds: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; + +#endif +} + +void test_katz_karate (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + int64_t niters = 0 ; + + // Create karate graph (undirected), treated as unweighted + snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_UNDIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + // Check that AT is cached + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + // Print graph stats + uint64_t n, nedges ; + OK (GrB_Matrix_nrows(&n, G->A)) ; + OK (GrB_Matrix_nvals(&nedges, G->A)) ; + printf ("\n\nKarate graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", n, nedges) ; + + // alpha must be less than 1/lambda_max, where lambda_max is 6.725697727631724 for karate graph + double alpha = 0.1 ; + + // Compute katz centrality + double t = LAGraph_WallClockTime() ; + OK (LAGr_KatzCentrality (¢rality, &niters, G, alpha, 1.0, 1000, 1e-6, false, false, msg)) ; + t = LAGraph_WallClockTime() - t ; + printf (" Time for LAGr_KatzCentrality: %g sec\n", t) ; + printf (" Iterations for LAGr_KatzCentrality: %" PRId64 "\n", niters) ; + + // Compare with reference values. + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + TEST_CHECK (cnvals == n) ; + + double err = difference (centrality, karate_katz, 34) ; + printf (" karate: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; + +#endif +} + +void test_katz_ldbc_directed_example (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + int64_t niters = 0 ; + + // Create ldbc-directed-example graph (weighted and directed) + snprintf (filename, LEN, LG_DATA_DIR "%s", "ldbc-directed-example.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + // Check that AT is cached + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + // Print graph stats + uint64_t n, nedges ; + OK (GrB_Matrix_nrows(&n, G->A)) ; + OK (GrB_Matrix_nvals(&nedges, G->A)) ; + printf ("\n\nldbc directed example graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", n, nedges) ; + + double alpha = 0.1 ; + + // Compute katz centrality with edge weights considered + double t = LAGraph_WallClockTime() ; + OK (LAGr_KatzCentrality (¢rality, &niters, G, alpha, 1.0, 1000, 1e-6, false, true, msg)) ; + t = LAGraph_WallClockTime() - t ; + printf (" Time for LAGr_KatzCentrality: %g sec\n", t) ; + printf (" Iterations for LAGr_KatzCentrality: %" PRId64 "\n", niters) ; + + // Compare with reference values. + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + TEST_CHECK (cnvals == n) ; + + double err = difference (centrality, ldbc_directed_example_katz, 10) ; + printf (" ldbc-directed-example: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; + +#endif +} + +void test_katz_rand (void) +{ +#if LAGRAPH_SUITESPARSE + LAGraph_Graph G = NULL ; + OK (LAGraph_Init (msg)) ; + GrB_Matrix A = NULL ; + GrB_Vector centrality = NULL ; + int64_t niters = 0 ; + + // Create rand graph (directed), treated as unweighted + snprintf (filename, LEN, LG_DATA_DIR "%s", "rand.mtx") ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + // Check that AT is cached + int result = LAGraph_Cached_AT (G, msg) ; + TEST_CHECK (result == GrB_SUCCESS || result == LAGRAPH_CACHE_NOT_NEEDED) ; + + // Print graph stats + uint64_t n, nedges ; + OK (GrB_Matrix_nrows(&n, G->A)) ; + OK (GrB_Matrix_nvals(&nedges, G->A)) ; + printf ("\n\nrand graph (%" PRIu64 " nodes, %" PRIu64 " edges):\n", n, nedges) ; + + double alpha = 0.1 ; + + // normalize results + double t = LAGraph_WallClockTime() ; + OK (LAGr_KatzCentrality (¢rality, &niters, G, alpha, 1.0, 1000, 1e-6, true, false, msg)) ; + t = LAGraph_WallClockTime() - t ; + printf (" Time for LAGr_KatzCentrality: %g sec\n", t) ; + printf (" Iterations for LAGr_KatzCentrality: %" PRId64 "\n", niters) ; + + // Compare with reference values. + GrB_Index cn = 0, cnvals = 0 ; + OK (GrB_Vector_size (&cn, centrality)) ; + OK (GrB_Vector_nvals (&cnvals, centrality)) ; + TEST_CHECK (cn == n) ; + TEST_CHECK (cnvals == n) ; + + double err = difference (centrality, rand_katz, 20) ; + printf (" rand: err: %e\n", err) ; + TEST_CHECK (err < 1e-4) ; + + OK (GrB_free (¢rality)) ; + OK (LAGraph_Delete (&G, msg)) ; + OK (LAGraph_Finalize (msg)) ; + +#endif +} + + +//------------------------------------------------------------------------------ +// list of tests +//------------------------------------------------------------------------------ + +TEST_LIST = { + {"test_katz_diamonds", test_katz_diamonds}, + {"test_katz_karate", test_katz_karate}, + {"test_katz_ldbc_directed_example", test_katz_ldbc_directed_example}, + {"test_katz_rand", test_katz_rand}, + {NULL, NULL} +} ; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 2e2740940d..5251699387 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1317,6 +1317,54 @@ int LAGr_EdgeBetweennessCentrality char *msg ) ; +//------------------------------------------------------------------------------ +// Closeness centrality +//------------------------------------------------------------------------------ + +typedef enum +{ + CC_BFS, // unweighted BFS (unit edge weights) + CC_SSSP, // delta-stepping SSSP (non-negative weights) + CC_BELLMAN_FORD, // Bellman-Ford (general weights, negative-cycle check) + CC_FLOYD_WARSHALL // Floyd-Warshall (FW) APSP + // CC_DEFAULT, select algorithm automatically +} cc_algo_t ; + +LAGRAPHX_PUBLIC +int LAGr_ClosenessCentrality +( + // output: + GrB_Vector *centrality, + // input: + LAGraph_Graph G, + GrB_Vector sources, // nodes to score; NULL or empty => all nodes + bool use_weights, // if true, use edge weights in shortest paths + cc_algo_t algorithm, // shortest-path algorithm to use + GrB_Scalar Delta, // delta for SSSP; if NULL, derived from G->emin + char *msg +) ; + +//------------------------------------------------------------------------------ +// Katz centrality +//------------------------------------------------------------------------------ + +LAGRAPHX_PUBLIC +int LAGr_KatzCentrality +( + // output: + GrB_Vector *centrality, + int64_t *iters, + // input: + LAGraph_Graph G, + double alpha, + double beta, + int64_t max_iter, + double tol, + bool normalize, + bool use_weights, + char *msg +) ; + //------------------------------------------------------------------------------ // harmonic centrality (approximate via HLL sketches) //------------------------------------------------------------------------------ diff --git a/papers/BHALLA-PRIMARY-2026.pdf b/papers/BHALLA-PRIMARY-2026.pdf new file mode 100644 index 0000000000..9ce36e3a5a Binary files /dev/null and b/papers/BHALLA-PRIMARY-2026.pdf differ