From d4c4b1920656d3fec7ee541ebf25e6215151bb3a Mon Sep 17 00:00:00 2001 From: GomezGab Date: Thu, 19 Mar 2026 18:01:58 -0500 Subject: [PATCH 1/5] add harmonic centrality method, demo, and test --- Contributors.txt | 2 +- .../algorithm/LAGr_HarmonicCentrality.c | 571 ++++++++++++++++++ experimental/algorithm/LAGr_SwapEdges.c | 2 +- .../algorithm/LAGraph_RichClubCoefficient.c | 2 +- experimental/algorithm/LAGraph_SwapEdges.c | 2 +- experimental/algorithm/LAGraph_msf.c | 2 +- experimental/algorithm/LG_CC_FastSV7_FA.c | 4 +- .../benchmark/harmonicCentrality_demo.c | 129 ++++ experimental/test/LG_check_RCC.c | 2 +- experimental/test/test_HarmonicCentrality.c | 236 ++++++++ experimental/utility/LAGraph_FastAssign.c | 2 +- include/LAGraphX.h | 35 ++ 12 files changed, 980 insertions(+), 9 deletions(-) create mode 100644 experimental/algorithm/LAGr_HarmonicCentrality.c create mode 100644 experimental/benchmark/harmonicCentrality_demo.c create mode 100644 experimental/test/test_HarmonicCentrality.c diff --git a/Contributors.txt b/Contributors.txt index 70e737fd95..94bf15e27a 100644 --- a/Contributors.txt +++ b/Contributors.txt @@ -15,7 +15,7 @@ Contributors (in alphabetical order): Florentin Dorre Technische Univeritat Dresden, Neo4j Marton Elekes Budapest University of Technology and Economics, Hungary Alexandra Goff Texas A&M University - Gabriel Gomez Texas A&M University + Gabriel A. Gomez Texas A&M University, FalkorDB Semyon Grigoriev St. Petersburg State University Balint Hegyi Budapest University of Technology and Economics, Hungary Tanner Hoke Texas A&M University diff --git a/experimental/algorithm/LAGr_HarmonicCentrality.c b/experimental/algorithm/LAGr_HarmonicCentrality.c new file mode 100644 index 0000000000..89ccf80090 --- /dev/null +++ b/experimental/algorithm/LAGr_HarmonicCentrality.c @@ -0,0 +1,571 @@ +//------------------------------------------------------------------------------ +// LAGr_HarmonicCentrality: Estimate the harmonic centrality for all nodes +//------------------------------------------------------------------------------ + +// 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 Gabriel A. Gomez, FalkorDB + +//------------------------------------------------------------------------------ + +// TODO: ready for src +// FIXME: add HLL licence + +// TODO: References: + +#include "GraphBLAS.h" +#include "LAGraphX.h" +#include "LG_internal.h" + +#define CENTRALITY_MAX_ITER 100 + +#define HLL_P 10 // Precision +#define HLL_REGISTERS (1 << HLL_P) // number of registers + +#define HC_JIT_DEFINITIONS \ +"#define HLL_P 10 \n"\ + "#define HLL_REGISTERS (1 << HLLP) \n" + +#define HC_MAKE_STR(x) #x +#define HC_TO_STR(x) HC_MAKE_STR(x) +// register a function (or type) and the jit string for its definition +#define JIT_STR(f, var) static char *var = #f; f + +JIT_STR( +typedef struct { + uint8_t registers[(1 << 10)]; +} HLL; +, HLL_jit) + +static __inline uint8_t _hll_rank(uint32_t hash, uint8_t bits) { + uint8_t i; + + for (i = 1; i <= 32 - bits; i++) { + if (hash & 1) + break; + + hash >>= 1; + } + + return i; +} + +static __inline void _hll_add_hash(HLL *hll, uint32_t hash) { + uint32_t index = hash >> (32 - HLL_P); + uint8_t rank = _hll_rank(hash, HLL_P); + + if (rank > hll->registers[index]) { + hll->registers[index] = rank; + } +} + +double _hll_count(const HLL *hll) { + uint32_t i; + + double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)HLL_REGISTERS); + + alpha_mm *= ((double)HLL_REGISTERS * (double)HLL_REGISTERS); + + double sum = 0; + for (uint32_t i = 0; i < HLL_REGISTERS; i++) { + sum += 1.0 / (1 << hll->registers[i]); + } + + double estimate = alpha_mm / sum; + + if (estimate <= 5.0 / 2.0 * (double)HLL_REGISTERS) { + int zeros = 0; + + for (i = 0; i < HLL_REGISTERS; i++) + zeros += (hll->registers[i] == 0); + + if (zeros) + estimate = (double)HLL_REGISTERS * log((double)HLL_REGISTERS / zeros); + + } else if (estimate > (1.0 / 30.0) * 4294967296.0) { + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + } + + return estimate; +} + +//------------------------------------------------------------------------------ +// GraphBLAS Ops +//------------------------------------------------------------------------------ +#define GOLDEN_GAMMA 0x9E3779B97F4A7C15LL + +// TODO: does this hash have good properties? Good spread is much more important +// than speed here. +uint64_t LG_HC_hash(GrB_Index i) { + uint64_t result = (i + GOLDEN_GAMMA); + result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9LL; + result = (result ^ (result >> 27)) * 0x94D049BB133111EBLL; + result = (result ^ (result >> 31)); + return result; +} + +// no need to register in JIT, not proformance critical +// load a new HLL set with *x hashes in the set. +void fdb_hll_init(HLL *z, const uint64_t *x, GrB_Index i, GrB_Index j, + bool theta) { + int64_t weight = *x; + memset(z->registers, 0, HLL_REGISTERS); + // build a hash chain seeded from the node index i: + // each iteration hashes the previous result, yielding 'weight' + // distinct pseudo-random values for the HLL sketch + uint32_t hash = LG_HC_hash(i); + for (int64_t h = 0; h < weight; h++) { + _hll_add_hash(z, hash); + hash = LG_HC_hash(hash); + } +} + +JIT_STR( +void fdb_hll_merge(HLL *z, const HLL *x, const HLL *y) { + for (uint32_t i = 0; i < (1 << 10); i++) { + z->registers[i] = y->registers[i] > x->registers[i] ? y->registers[i] + : x->registers[i]; + } +}, +FDB_HLL_MERGE_STR) + +JIT_STR( + void fdb_hll_delta(double *z, const HLL *x, const HLL *y) { + *z = 0; + bool diff = 0 != memcmp(x->registers, y->registers, (1 << 10)); + if (diff) { + uint32_t i; + + double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(1 << 10)); + + alpha_mm *= ((double)(1 << 10) * (double)(1 << 10)); + + double sum = 0; + for (uint32_t i = 0; i < (1 << 10); i++) { + sum += 1.0 / (1 << x->registers[i]); + } + + double estimate = alpha_mm / sum; + + if (estimate <= 5.0 / 2.0 * (double)(1 << 10)) { + int zeros = 0; + + for (i = 0; i < (1 << 10); i++) + zeros += (x->registers[i] == 0); + + if (zeros) + estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); + + } else if (estimate > (1.0 / 30.0) * 4294967296.0) { + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + } + + *z = estimate; + + sum = 0; + for (uint32_t i = 0; i < (1 << 10); i++) { + sum += 1.0 / (1 << y->registers[i]); + } + + estimate = alpha_mm / sum; + + if (estimate <= 5.0 / 2.0 * (double)(1 << 10)) { + int zeros = 0; + + for (i = 0; i < (1 << 10); i++) + zeros += (y->registers[i] == 0); + + if (zeros) + estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); + + } else if (estimate > (1.0 / 30.0) * 4294967296.0) { + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + } + *z -= estimate; + } +}, +FDB_HLL_DELTA_STR) + +JIT_STR( + void fdb_hll_second(HLL *z, + bool *x, // unused + const HLL *y) { + memcpy(z->registers, y->registers, 1 << 10); +}, +FDB_HLL_SECOND_STR) + +int64_t print_hll(char *string, // value is printed to the string + size_t string_size, // size of the string array + const void *value, // HLL value to print + int verbose // if >0, print verbosely; else tersely + ) { + const HLL *hll = (const HLL *)value; + + if (verbose > 0) { + return snprintf(string, string_size, "HLL{bits=%u, size=%u, count=%.2f}", + HLL_P, HLL_REGISTERS, _hll_count(hll)); + } + + return snprintf(string, string_size, "HLL{count=%.2f}", _hll_count(hll)); +} + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ +GrB_free(&_A); \ +GrB_free(&score_cont); \ +GrB_free(&hll_t); \ +GrB_free(&new_sets); \ +GrB_free(&old_sets); \ +GrB_free(&old_cont); \ +GrB_free(&flat_scores); \ +GrB_free(&flat_weight); \ +GrB_free(&delta_vec); \ +GrB_free(&desc); \ +GrB_free(&init_hlls); \ +GrB_free(&shallow_second); \ +GrB_free(&merge_hll_biop); \ +GrB_free(&merge_hll); \ +GrB_free(&merge_second); \ +GrB_free(&delta_hll); + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ +LG_FREE_WORK ; \ +GrB_free(scores) ; + +// compute harmonic closeness centrality estimates using HLL BFS propagation +// +// each node maintains an HLL sketch of "reachable nodes seen so far" +// at BFS level d, sketches are propagated along edges +// the harmonic contribution delta/d is accumulated into each node's score, +// where delta is the number of new nodes discovered at that level +int LAGr_HarmonicCentrality( + // outputs: + GrB_Vector *scores, // FP64 scores by original node ID + GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- + // able nodes from the given node. + // inputs: + const GrB_Matrix A, // adjacency matrix + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) { + GrB_Matrix _A = NULL; + GxB_Container score_cont = NULL; + GrB_Index nrows = 0; + GrB_Index nvals = 0; + + GrB_Type hll_t = NULL; + GrB_Vector new_sets = NULL; + GrB_Vector old_sets = NULL; + GxB_Container old_cont = NULL; + GrB_Vector flat_scores = NULL; + GrB_Vector flat_weight = NULL; + GrB_Vector delta_vec = NULL; + + GrB_Descriptor desc = NULL; + GrB_IndexUnaryOp init_hlls = NULL; + GrB_BinaryOp shallow_second = NULL; + GrB_BinaryOp merge_hll_biop = NULL; + GrB_Monoid merge_hll = NULL; + GrB_Semiring merge_second = NULL; + GrB_BinaryOp delta_hll = NULL; + + // TODO: + LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); + LG_ASSERT(A != NULL, GrB_NULL_POINTER); + LG_ASSERT(scores != NULL, GrB_NULL_POINTER); + LG_ASSERT(node_weights != NULL, GrB_NULL_POINTER); + + GRB_TRY(GrB_Vector_size(&nrows, node_weights)); + GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)); + GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)); + + if (nvals == 0) { + return GrB_SUCCESS; + } + + // double check weight type and maximum weight requirements + GrB_Type weight_t = NULL; + GRB_TRY(GxB_Vector_type(&weight_t, node_weights)); + LG_ASSERT_MSG(weight_t == GrB_BOOL || weight_t == GrB_INT64, + GrB_DOMAIN_MISMATCH, "weight must be integer or uint64"); + + if (weight_t == GrB_INT64) { + int64_t max_w = 0, min_w = 0; + GRB_TRY(GrB_Vector_reduce_INT64(&max_w, NULL, GrB_MAX_MONOID_INT64, + node_weights, NULL)); + GRB_TRY(GrB_Vector_reduce_INT64(&min_w, NULL, GrB_MIN_MONOID_INT64, + node_weights, NULL)); + LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, + "Negative node weights not supported"); + // TODO: is this cap reasonable? + LG_ASSERT_MSG(max_w < 1000000, GrB_NOT_IMPLEMENTED, + "Node weights over 1000000 not supported"); + } + + //-------------------------------------------------------------------------- + // build compact submatrix _A + //-------------------------------------------------------------------------- + + GRB_TRY(GrB_Descriptor_new(&desc)); + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)); + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)); + GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)); + GRB_TRY(GxB_Matrix_extract_Vector(_A, NULL, NULL, A, node_weights, + node_weights, desc)); + GRB_TRY(GrB_free(&desc)); + + //-------------------------------------------------------------------------- + // create scores vector (0.0 at each participating node) + //-------------------------------------------------------------------------- + + GRB_TRY(GrB_Vector_assign_FP64(*scores, node_weights, NULL, 0.0, GrB_ALL, 0, + GrB_DESC_S)); + GRB_TRY(GrB_set(*scores, GxB_SPARSE | GxB_FULL, GxB_SPARSITY_CONTROL)); + GRB_TRY(GxB_Container_new(&score_cont)); + GRB_TRY(GxB_unload_Vector_into_Container(*scores, score_cont, NULL)); + + //-------------------------------------------------------------------------- + // Initialize HLL + //-------------------------------------------------------------------------- + + GRB_TRY(GxB_Type_new(&hll_t, sizeof(HLL), "HLL", HLL_jit)); + GRB_TRY(GrB_Vector_new(&new_sets, hll_t, nvals)); + GRB_TRY(GrB_Vector_new(&old_sets, hll_t, nvals)); + GRB_TRY(GxB_Container_new(&old_cont)); + GRB_TRY(GrB_Vector_new(&flat_scores, GrB_FP64, nvals)); + GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)); + GRB_TRY( + GxB_Vector_extractTuples_Vector(NULL, flat_weight, node_weights, NULL)); + + + // init op: weight (INT64) at row index i → HLL seeded with 'weight' hashes + GRB_TRY(GrB_IndexUnaryOp_new(&init_hlls, + (GxB_index_unary_function)fdb_hll_init, hll_t, + GrB_INT64, GrB_BOOL)); + + // merge binary op (HLL, HLL) -> HLL in-place: z == x required + GRB_TRY(GxB_BinaryOp_new(&merge_hll_biop, (GxB_binary_function)fdb_hll_merge, + hll_t, hll_t, hll_t, "fdb_hll_merge", + FDB_HLL_MERGE_STR)); + + // second op + GRB_TRY(GxB_BinaryOp_new(&shallow_second, (GxB_binary_function)fdb_hll_second, + hll_t, GrB_BOOL, hll_t, "fdb_hll_second", + FDB_HLL_SECOND_STR)); + + // merge monoid - identity is an empty (all-zero) HLL sketch + HLL hll_zero = {0}; + GRB_TRY(GrB_Monoid_new_UDT(&merge_hll, merge_hll_biop, &hll_zero)); + + // semiring: add = merge monoid, multiply = copy (pass-through second operand) + GRB_TRY(GrB_Semiring_new(&merge_second, merge_hll, shallow_second)); + + // delta op: (HLL_old, HLL_new) → FP64 cardinality change + GRB_TRY(GxB_BinaryOp_new(&delta_hll, (GxB_binary_function)fdb_hll_delta, + GrB_FP64, hll_t, hll_t, "fdb_hll_delta", + FDB_HLL_DELTA_STR)); + + // delta output: one FP64 entry per participating node + GRB_TRY(GrB_Vector_new(&delta_vec, GrB_FP64, nvals)); + + //-------------------------------------------------------------------------- + // Load vector values + //-------------------------------------------------------------------------- + // initialize HLLs + GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(new_sets, NULL, NULL, init_hlls, + flat_weight, false, NULL)); + GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(old_sets, NULL, NULL, init_hlls, + flat_weight, false, NULL)); + GRB_TRY(GrB_set(old_sets, GxB_BITMAP, GxB_SPARSITY_CONTROL)); + + GRB_TRY(GrB_free(&flat_weight)); + GRB_TRY( + GrB_Vector_assign_FP64(flat_scores, NULL, NULL, 0.0, GrB_ALL, 0, NULL)); + + //-------------------------------------------------------------------------- + // HLL BFS propagation + //-------------------------------------------------------------------------- + + int64_t changes = 0; + for (int d = 1; d <= CENTRALITY_MAX_ITER; d++) { + changes = 0; + + // foward bfs + // merge each neighbor's pre-round set into this node's set + // target kernel for inplace adjustments + GRB_TRY(GrB_mxv(new_sets, NULL, merge_hll_biop, merge_second, _A, old_sets, + NULL)); + + GRB_TRY(GxB_unload_Vector_into_Container(old_sets, old_cont, NULL)); + // find the delta between last round and this one + + GRB_TRY(GrB_eWiseMult(delta_vec, NULL, NULL, delta_hll, new_sets, + old_cont->x, NULL)); + GRB_TRY( + GrB_Vector_assign(old_cont->x, NULL, NULL, new_sets, GrB_ALL, 0, NULL)); + + // old_set bitmap is the set of nodes with non-zero deltas + GRB_TRY( + GrB_apply(old_cont->b, NULL, NULL, GrB_IDENTITY_BOOL, delta_vec, NULL)); + GRB_TRY(GrB_Vector_reduce_INT64(&changes, NULL, GrB_PLUS_MONOID_INT64, + old_cont->b, NULL)); + old_cont->nvals = changes; + GRB_TRY(GxB_load_Vector_from_Container(old_sets, old_cont, NULL)); + + // use the deltas to update the score + GRB_TRY(GrB_apply(flat_scores, NULL, GrB_PLUS_FP64, GrB_DIV_FP64, delta_vec, + (double)d, NULL)); + + // stop when no HLL cardinality changed this round + if (changes == 0) { + break; + } + } + + //-------------------------------------------------------------------------- + // write flat_scores into score_cont->x, clear iso, reload scores vector + //-------------------------------------------------------------------------- + GRB_TRY (GrB_free (&score_cont->x)); + score_cont->iso = false; + score_cont->x = flat_scores; + flat_scores = NULL ; + GRB_TRY (GxB_load_Vector_from_Container (*scores, score_cont, NULL)); + + //-------------------------------------------------------------------------- + // cleanup + //-------------------------------------------------------------------------- + + // free operators (semiring before monoid) + + LG_FREE_WORK; + return GrB_SUCCESS; +} + +// TODO: add a descriptor? let nodes be values or indecies + +#undef LG_FREE_WORK +#define LG_FREE_WORK \ + GrB_free(&_A); \ + GrB_free(&flat_weight); \ + GrB_free(&level); \ + GrB_free(&it); \ + GrB_free(&score); \ + GrB_free(&desc); \ + LAGraph_Delete(&G_compact, NULL); + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + LG_FREE_WORK ; \ + GrB_free(scores) ; + +// This is an exact (but much slower) calculation of harmonic centrality +int LAGr_HarmonicCentrality_exact( + // outputs: + GrB_Vector *scores, // FP64 scores by original node ID + GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- + // able nodes from the given node. + // inputs: + const GrB_Matrix A, // adjacency matrix + const GrB_Vector nodes, // nodes to calculate centrality of + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) { + GrB_Matrix _A = NULL; + GrB_Index nrows = 0; + GrB_Index nvals = 0; + + GrB_Vector level = NULL; + GrB_Vector flat_weight = NULL; + GxB_Iterator it = NULL; + GrB_Scalar score = NULL; + + GrB_Descriptor desc = NULL; + LAGraph_Graph G_compact = NULL; + + // TODO: + LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); + LG_ASSERT(A != NULL, GrB_NULL_POINTER); + LG_ASSERT(scores != NULL, GrB_NULL_POINTER); + LG_ASSERT(node_weights != NULL, GrB_NULL_POINTER); + + GRB_TRY(GrB_Vector_size(&nrows, node_weights)); + GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)); + GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)); + + if (nvals == 0) { + return GrB_SUCCESS; + } + + // double check weight type and maximum weight requirements + GrB_Type weight_t = NULL; + GRB_TRY(GxB_Vector_type(&weight_t, node_weights)); + LG_ASSERT_MSG(weight_t == GrB_BOOL || weight_t == GrB_INT64, + GrB_DOMAIN_MISMATCH, "weight must be integer or uint64"); + + if (weight_t == GrB_INT64) { + int64_t min_w = 0; + GRB_TRY(GrB_Vector_reduce_INT64(&min_w, NULL, GrB_MIN_MONOID_INT64, + node_weights, NULL)); + LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, + "Negative node weights not supported"); + } + + //-------------------------------------------------------------------------- + // build compact submatrix _A + //-------------------------------------------------------------------------- + + GRB_TRY(GrB_Descriptor_new(&desc)); + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)); + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)); + GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)); + GRB_TRY(GxB_Matrix_extract_Vector(_A, NULL, NULL, A, node_weights, + node_weights, desc)); + GRB_TRY(GrB_free(&desc)); + LG_TRY (LAGraph_New (&G_compact, &_A, LAGraph_ADJACENCY_DIRECTED, msg)); + + GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)); + GRB_TRY( GxB_Vector_extractTuples_Vector( + NULL, flat_weight, node_weights, NULL)); + //-------------------------------------------------------------------------- + // calculate centrality via BFS + //-------------------------------------------------------------------------- + + GRB_TRY (GxB_Scalar_new (&score, GrB_FP64)); + GRB_TRY(GxB_Iterator_new(&it)); + GRB_TRY(GxB_Vector_Iterator_attach(it, nodes, NULL)); + GrB_Info info = GxB_Vector_Iterator_seek(it, 0); + while (info != GxB_EXHAUSTED) { + GrB_Index nodeIdx = GxB_Vector_Iterator_getIndex(it); + + LG_TRY (LAGr_BreadthFirstSearch(&level, NULL, G_compact, nodeIdx, msg)); + + // remove 0 to prevent div by 0 + GRB_TRY (GrB_Vector_removeElement (level, nodeIdx)); + + // dot product (plus div) + GRB_TRY(GrB_vxm((GrB_Vector) score, NULL, NULL, GxB_PLUS_DIV_FP64, + flat_weight, (GrB_Matrix) level, NULL)); + GrB_free (&level); + + GRB_TRY (GrB_Vector_setElement_Scalar(*scores, score, nodeIdx)); + + info = GxB_Vector_Iterator_next(it); + } + + //-------------------------------------------------------------------------- + // cleanup + //-------------------------------------------------------------------------- + + // free operators (semiring before monoid) + + LG_FREE_WORK; + return GrB_SUCCESS; +} diff --git a/experimental/algorithm/LAGr_SwapEdges.c b/experimental/algorithm/LAGr_SwapEdges.c index a72ffc636e..398957a6e3 100755 --- a/experimental/algorithm/LAGr_SwapEdges.c +++ b/experimental/algorithm/LAGr_SwapEdges.c @@ -11,7 +11,7 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 -// Contributed by Gabriel Gomez, Texas A&M University +// Contributed by Gabriel A. Gomez, Texas A&M University //------------------------------------------------------------------------------ diff --git a/experimental/algorithm/LAGraph_RichClubCoefficient.c b/experimental/algorithm/LAGraph_RichClubCoefficient.c index d20ed9dac1..2b6f67539f 100755 --- a/experimental/algorithm/LAGraph_RichClubCoefficient.c +++ b/experimental/algorithm/LAGraph_RichClubCoefficient.c @@ -11,7 +11,7 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 -// Contributed by Gabriel Gomez, Texas A&M University +// Contributed by Gabriel A. Gomez, Texas A&M University //------------------------------------------------------------------------------ diff --git a/experimental/algorithm/LAGraph_SwapEdges.c b/experimental/algorithm/LAGraph_SwapEdges.c index 81d03b5fad..76eb03aed9 100644 --- a/experimental/algorithm/LAGraph_SwapEdges.c +++ b/experimental/algorithm/LAGraph_SwapEdges.c @@ -11,7 +11,7 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 -// Contributed by Gabriel Gomez, Texas A&M University +// Contributed by Gabriel A. Gomez, Texas A&M University //------------------------------------------------------------------------------ diff --git a/experimental/algorithm/LAGraph_msf.c b/experimental/algorithm/LAGraph_msf.c index 193c6fe322..c8279c54a9 100755 --- a/experimental/algorithm/LAGraph_msf.c +++ b/experimental/algorithm/LAGraph_msf.c @@ -12,7 +12,7 @@ // DM22-0790 // Contributed by Yongzhe Zhang (zyz915@gmail.com) -// Revised by Gabriel Gomez and Tim Davis +// Revised by Gabriel A. Gomez and Tim Davis //------------------------------------------------------------------------------ diff --git a/experimental/algorithm/LG_CC_FastSV7_FA.c b/experimental/algorithm/LG_CC_FastSV7_FA.c index 3060159227..8b478bca73 100644 --- a/experimental/algorithm/LG_CC_FastSV7_FA.c +++ b/experimental/algorithm/LG_CC_FastSV7_FA.c @@ -34,8 +34,8 @@ // Changed to use GxB load/unload. Converted to use the LAGraph_Graph object. // Exploiting iso status for the temporary matrices Parent and T. -// Modified by Gabriel Gomez, Texas A&M University: moved Parent matrix trick -// out to LAGraph_FastAssign. +// Modified by Gabriel A. Gomez, Texas A&M University: moved Parent matrix +// trick out to LAGraph_FastAssign. // The input graph G must be undirected, or directed and with an adjacency // matrix that has a symmetric structure. Self-edges (diagonal entries) are diff --git a/experimental/benchmark/harmonicCentrality_demo.c b/experimental/benchmark/harmonicCentrality_demo.c new file mode 100644 index 0000000000..3f3845b5ff --- /dev/null +++ b/experimental/benchmark/harmonicCentrality_demo.c @@ -0,0 +1,129 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/benchmark/harmonicCentrality_demo.c: benchmark for +// harmonic centrality +//------------------------------------------------------------------------------ + +// 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 Gabriel A. Gomez, FalkorDB + +//------------------------------------------------------------------------------ + +// Usage: harmonicCentrality_demo matrixmarketfile.mtx +// harmonicCentrality_demo matrixmarketfile.mtx exact + +#include "../../src/benchmark/LAGraph_demo.h" +#include "LAGraphX.h" + +#define LG_FREE_ALL \ +{ \ + GrB_free (&A) ; \ + GrB_free (&scores_approx) ; \ + GrB_free (&scores_exact) ; \ + GrB_free (&node_weights) ; \ +} + +int main (int argc, char **argv) +{ + + //-------------------------------------------------------------------------- + // initialize LAGraph and GraphBLAS + //-------------------------------------------------------------------------- + + char msg [LAGRAPH_MSG_LEN] ; + + GrB_Matrix A = NULL ; + GrB_Vector scores_approx = NULL ; + GrB_Vector scores_exact = NULL ; + GrB_Vector node_weights = NULL ; + + // start GraphBLAS and LAGraph + bool burble = false ; + demo_init (burble) ; + + //-------------------------------------------------------------------------- + // read in the graph + //-------------------------------------------------------------------------- + + if (argc < 2) + { + printf ("Usage: %s [exact]\n", argv [0]) ; + return (GrB_INVALID_VALUE) ; + } + + bool run_exact = (argc >= 3 && strcmp (argv [2], "exact") == 0) ; + + char *matrix_name = argv [1] ; + FILE *f = fopen (matrix_name, "r") ; + if (f == NULL) + { + printf ("Error: unable to open file %s\n", matrix_name) ; + return (GrB_INVALID_VALUE) ; + } + + double t = LAGraph_WallClockTime ( ) ; + LAGRAPH_TRY (LAGraph_MMRead (&A, f, msg)) ; + fclose (f) ; + t = LAGraph_WallClockTime ( ) - t ; + printf ("Time to read the graph: %g sec\n", t) ; + + GrB_Index n, nvals ; + GRB_TRY (GrB_Matrix_nrows (&n, A)) ; + GRB_TRY (GrB_Matrix_nvals (&nvals, A)) ; + printf ("Graph: %s (%" PRIu64 " nodes, %" PRIu64 " edges)\n", + matrix_name, (uint64_t) n, (uint64_t) nvals) ; + + //-------------------------------------------------------------------------- + // create boolean node_weights (all nodes participate with weight = 1) + //-------------------------------------------------------------------------- + + GRB_TRY (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; + GRB_TRY (GrB_Vector_assign_BOOL ( + node_weights, NULL, NULL, true, GrB_ALL, n, NULL)) ; + + //-------------------------------------------------------------------------- + // compute approximate harmonic centrality (HLL) + //-------------------------------------------------------------------------- + + t = LAGraph_WallClockTime ( ) ; + LAGRAPH_TRY (LAGr_HarmonicCentrality ( + &scores_approx, NULL, A, node_weights, msg)) ; + t = LAGraph_WallClockTime ( ) - t ; + printf ("Time for LAGr_HarmonicCentrality (approx): %g sec\n", t) ; + + // print a summary of the scores + LAGraph_PrintLevel pr = (n <= 40) ? LAGraph_COMPLETE : LAGraph_SHORT ; + printf ("\napproximate scores:\n") ; + LAGRAPH_TRY (LAGraph_Vector_Print (scores_approx, pr, stdout, msg)) ; + + //-------------------------------------------------------------------------- + // optionally compute exact harmonic centrality (BFS) + //-------------------------------------------------------------------------- + + if (run_exact) + { + t = LAGraph_WallClockTime ( ) ; + LAGRAPH_TRY (LAGr_HarmonicCentrality_exact ( + &scores_exact, NULL, A, node_weights, node_weights, msg)) ; + t = LAGraph_WallClockTime ( ) - t ; + printf ("\nTime for LAGr_HarmonicCentrality_exact: %g sec\n", t) ; + + printf ("\nexact scores:\n") ; + LAGRAPH_TRY (LAGraph_Vector_Print (scores_exact, pr, stdout, msg)) ; + } + + //-------------------------------------------------------------------------- + // free everything and finish + //-------------------------------------------------------------------------- + + LG_FREE_ALL ; + LAGRAPH_TRY (LAGraph_Finalize (msg)) ; + return (GrB_SUCCESS) ; +} diff --git a/experimental/test/LG_check_RCC.c b/experimental/test/LG_check_RCC.c index 3a9da8ad81..32c2aad7d3 100644 --- a/experimental/test/LG_check_RCC.c +++ b/experimental/test/LG_check_RCC.c @@ -11,7 +11,7 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 -// Contributed by Gabriel Gomez, Texas A&M University +// Contributed by Gabriel A. Gomez, Texas A&M University //------------------------------------------------------------------------------ diff --git a/experimental/test/test_HarmonicCentrality.c b/experimental/test/test_HarmonicCentrality.c new file mode 100644 index 0000000000..782bc60d96 --- /dev/null +++ b/experimental/test/test_HarmonicCentrality.c @@ -0,0 +1,236 @@ +//------------------------------------------------------------------------------ +// LAGraph/experimental/test/test_HarmonicCentrality.c: test cases for +// harmonic centrality (approximate HLL and exact BFS) +//------------------------------------------------------------------------------ + +// 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 Gabriel A. Gomez, FalkorDB + +//------------------------------------------------------------------------------ + +// NOTE: these tests require SuiteSparse:GraphBLAS + +#include +#include + +#include +#include + +char msg [LAGRAPH_MSG_LEN] ; +GrB_Matrix A = NULL ; +GrB_Vector scores_approx = NULL ; +GrB_Vector scores_exact = NULL ; +GrB_Vector node_weights = NULL ; +#define LEN 512 +char filename [LEN+1] ; + +const char *files [ ] = +{ + "karate.mtx", + "A.mtx", + "jagmesh7.mtx", + "ldbc-directed-example.mtx", + "ldbc-undirected-example.mtx", + "", +} ; + +//------------------------------------------------------------------------------ +// max_relative_error: max percentage difference between two FP64 vectors +//------------------------------------------------------------------------------ +// For each entry i where exact[i] != 0, compute |approx[i] - exact[i]| / +// |exact[i]|. Returns the maximum of these ratios. + +double max_relative_error +( + GrB_Vector approx, + GrB_Vector exact +) +{ + GrB_Index n ; + GrB_Vector_size (&n, exact) ; + + double max_err = 0 ; + for (GrB_Index i = 0 ; i < n ; i++) + { + double a = 0, e = 0 ; + GrB_Info info_a = GrB_Vector_extractElement_FP64 (&a, approx, i) ; + GrB_Info info_e = GrB_Vector_extractElement_FP64 (&e, exact, i) ; + + // skip entries not present in both vectors + if (info_a != GrB_SUCCESS || info_e != GrB_SUCCESS) continue ; + + if (fabs (e) > 0) + { + double rel = fabs (a - e) / fabs (e) ; + if (rel > max_err) max_err = rel ; + } + } + return max_err ; +} + +//------------------------------------------------------------------------------ +// test_HarmonicCentrality: compare approximate vs exact on small graphs +//------------------------------------------------------------------------------ + +void test_HarmonicCentrality (void) +{ + #if LAGRAPH_SUITESPARSE + LAGraph_Init (msg) ; + + for (int k = 0 ; ; k++) + { + const char *aname = files [k] ; + if (strlen (aname) == 0) break ; + printf ("\n================================== %s:\n", aname) ; + TEST_CASE (aname) ; + snprintf (filename, LEN, LG_DATA_DIR "%s", aname) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + TEST_MSG ("Loading of adjacency matrix failed") ; + + // get matrix dimensions + GrB_Index n ; + OK (GrB_Matrix_nrows (&n, A)) ; + + // create boolean node_weights (all nodes, weight = true) + OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; + OK (GrB_Vector_assign_BOOL ( + node_weights, NULL, NULL, true, GrB_ALL, n, NULL)) ; + + // compute approximate harmonic centrality + OK (LAGr_HarmonicCentrality ( + &scores_approx, NULL, A, node_weights, msg)) ; + printf ("%s", msg); + + // compute exact harmonic centrality + OK (LAGr_HarmonicCentrality_exact ( + &scores_exact, NULL, A, node_weights, node_weights, msg)) ; + + // print results for small graphs + LAGraph_PrintLevel pr = + (n <= 40) ? LAGraph_COMPLETE : LAGraph_SHORT ; + printf ("\napproximate scores:\n") ; + OK (LAGraph_Vector_Print (scores_approx, pr, stdout, msg)) ; + printf ("\nexact scores:\n") ; + OK (LAGraph_Vector_Print (scores_exact, pr, stdout, msg)) ; + + // compare: max relative (percentage) error + double rel_err = max_relative_error (scores_approx, scores_exact) ; + printf ("max relative error: %.2f%%\n", rel_err * 100) ; + + // allow up to 50% relative error for HLL on small graphs + TEST_CHECK (rel_err < 0.5) ; + TEST_MSG ("Relative error too large: %.2f%%", rel_err * 100) ; + + // cleanup for this iteration + OK (GrB_free (&A)) ; + OK (GrB_free (&scores_approx)) ; + OK (GrB_free (&scores_exact)) ; + OK (GrB_free (&node_weights)) ; + } + + LAGraph_Finalize (msg) ; + #endif +} + +//------------------------------------------------------------------------------ +// test_HarmonicCentrality_empty: test with empty node set +//------------------------------------------------------------------------------ + +void test_HarmonicCentrality_empty (void) +{ + #if LAGRAPH_SUITESPARSE + LAGraph_Init (msg) ; + + 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)) ; + + GrB_Index n ; + OK (GrB_Matrix_nrows (&n, A)) ; + + // empty node_weights — no participating nodes + OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; + + OK (LAGr_HarmonicCentrality (&scores_approx, NULL, A, node_weights, msg)) ; + + // scores should exist but have no entries + GrB_Index nvals ; + OK (GrB_Vector_nvals (&nvals, scores_approx)) ; + TEST_CHECK (nvals == 0) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&scores_approx)) ; + OK (GrB_free (&node_weights)) ; + + LAGraph_Finalize (msg) ; + #endif +} + +//------------------------------------------------------------------------------ +// test_errors: test error handling +//------------------------------------------------------------------------------ + +void test_errors (void) +{ + #if LAGRAPH_SUITESPARSE + LAGraph_Init (msg) ; + + 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)) ; + + GrB_Index n ; + OK (GrB_Matrix_nrows (&n, A)) ; + OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; + OK (GrB_Vector_assign_BOOL ( + node_weights, NULL, NULL, true, GrB_ALL, n, NULL)) ; + + int result ; + + // scores is NULL + result = LAGr_HarmonicCentrality (NULL, NULL, A, node_weights, msg) ; + printf ("\nresult: %d %s\n", result, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // A is NULL + result = LAGr_HarmonicCentrality ( + &scores_approx, NULL, NULL, node_weights, msg) ; + printf ("\nresult: %d %s\n", result, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + // node_weights is NULL + result = LAGr_HarmonicCentrality ( + &scores_approx, NULL, A, NULL, msg) ; + printf ("\nresult: %d %s\n", result, msg) ; + TEST_CHECK (result == GrB_NULL_POINTER) ; + + OK (GrB_free (&A)) ; + OK (GrB_free (&node_weights)) ; + + LAGraph_Finalize (msg) ; + #endif +} + +//**************************************************************************** + +TEST_LIST = { + {"HarmonicCentrality", test_HarmonicCentrality}, + {"HarmonicCentrality_empty", test_HarmonicCentrality_empty}, + {"HarmonicCentrality_errors", test_errors}, + {NULL, NULL} +} ; diff --git a/experimental/utility/LAGraph_FastAssign.c b/experimental/utility/LAGraph_FastAssign.c index 8374b35b97..25717aff73 100755 --- a/experimental/utility/LAGraph_FastAssign.c +++ b/experimental/utility/LAGraph_FastAssign.c @@ -13,7 +13,7 @@ // funding and support from the U.S. Government (see Acknowledgments.txt file). // DM22-0790 -// Contributed by Gabriel Gomez, Texas A&M University +// Contributed by Gabriel A. Gomez, Texas A&M University //------------------------------------------------------------------------------ diff --git a/include/LAGraphX.h b/include/LAGraphX.h index fa5058fe5e..6d9be99c38 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1317,6 +1317,41 @@ int LAGr_EdgeBetweennessCentrality char *msg ) ; +//------------------------------------------------------------------------------ +// harmonic centrality (approximate via HLL sketches) +//------------------------------------------------------------------------------ + +LAGRAPHX_PUBLIC +int LAGr_HarmonicCentrality +( + // outputs: + GrB_Vector *scores, // FP64 harmonic centrality scores + GrB_Vector *reachable_nodes, // [optional] estimated reachable node count + // (pass NULL, not yet implemented) + // inputs: + const GrB_Matrix A, // adjacency matrix + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) ; + +//------------------------------------------------------------------------------ +// harmonic centrality (exact via BFS) +//------------------------------------------------------------------------------ + +LAGRAPHX_PUBLIC +int LAGr_HarmonicCentrality_exact +( + // outputs: + GrB_Vector *scores, // FP64 harmonic centrality scores + GrB_Vector *reachable_nodes, // [optional] estimated reachable node count + // (pass NULL, not yet implemented) + // inputs: + const GrB_Matrix A, // adjacency matrix + const GrB_Vector nodes, // nodes to calculate centrality of + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) ; + //------------------------------------------------------------------------------ // graph clustering with quality metrics //------------------------------------------------------------------------------ From d07663829a773476fc221ab0d23ae99d02c04f76 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Fri, 20 Mar 2026 01:08:35 -0500 Subject: [PATCH 2/5] +// Permission is hereby granted, free of charge, to any person obtaining a copy +// of this software and associated documentation files (the "Software"), to deal +// in the Software without restriction, including without limitation the rights +// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the Software is +// furnished to do so, subject to the following conditions: +// The above copyright notice and this permission notice shall be included +// in all copies or substantial portions of the Software. +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR +// OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, +// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. #include "GraphBLAS.h" #include "LAGraphX.h" @@ -26,15 +48,11 @@ #define CENTRALITY_MAX_ITER 100 +// WARNING: these macros are not recognized by the JIT, so have to be changed +// manually in the functions below #define HLL_P 10 // Precision #define HLL_REGISTERS (1 << HLL_P) // number of registers -#define HC_JIT_DEFINITIONS \ -"#define HLL_P 10 \n"\ - "#define HLL_REGISTERS (1 << HLLP) \n" - -#define HC_MAKE_STR(x) #x -#define HC_TO_STR(x) HC_MAKE_STR(x) // register a function (or type) and the jit string for its definition #define JIT_STR(f, var) static char *var = #f; f @@ -90,7 +108,7 @@ double _hll_count(const HLL *hll) { estimate = (double)HLL_REGISTERS * log((double)HLL_REGISTERS / zeros); } else if (estimate > (1.0 / 30.0) * 4294967296.0) { - estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; } return estimate; @@ -107,13 +125,13 @@ uint64_t LG_HC_hash(GrB_Index i) { uint64_t result = (i + GOLDEN_GAMMA); result = (result ^ (result >> 30)) * 0xBF58476D1CE4E5B9LL; result = (result ^ (result >> 27)) * 0x94D049BB133111EBLL; - result = (result ^ (result >> 31)); + result = (result ^ (result >> 31)) ; return result; } // no need to register in JIT, not proformance critical // load a new HLL set with *x hashes in the set. -void fdb_hll_init(HLL *z, const uint64_t *x, GrB_Index i, GrB_Index j, +void lg_hll_init(HLL *z, const uint64_t *x, GrB_Index i, GrB_Index j, bool theta) { int64_t weight = *x; memset(z->registers, 0, HLL_REGISTERS); @@ -128,24 +146,24 @@ void fdb_hll_init(HLL *z, const uint64_t *x, GrB_Index i, GrB_Index j, } JIT_STR( -void fdb_hll_merge(HLL *z, const HLL *x, const HLL *y) { +void lg_hll_merge(HLL *z, const HLL *x, const HLL *y) { for (uint32_t i = 0; i < (1 << 10); i++) { z->registers[i] = y->registers[i] > x->registers[i] ? y->registers[i] : x->registers[i]; } }, -FDB_HLL_MERGE_STR) +LG_HLL_MERGE_STR) JIT_STR( - void fdb_hll_delta(double *z, const HLL *x, const HLL *y) { + void lg_hll_delta(double *z, const HLL *x, const HLL *y) { *z = 0; - bool diff = 0 != memcmp(x->registers, y->registers, (1 << 10)); + bool diff = 0 != memcmp(x->registers, y->registers, (1 << 10)) ; if (diff) { uint32_t i; - double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(1 << 10)); + double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(1 << 10)) ; - alpha_mm *= ((double)(1 << 10) * (double)(1 << 10)); + alpha_mm *= ((double)(1 << 10) * (double)(1 << 10)) ; double sum = 0; for (uint32_t i = 0; i < (1 << 10); i++) { @@ -164,7 +182,7 @@ JIT_STR( estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); } else if (estimate > (1.0 / 30.0) * 4294967296.0) { - estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; } *z = estimate; @@ -186,35 +204,36 @@ JIT_STR( estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); } else if (estimate > (1.0 / 30.0) * 4294967296.0) { - estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)); + estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; } *z -= estimate; } }, -FDB_HLL_DELTA_STR) +LG_HLL_DELTA_STR) JIT_STR( - void fdb_hll_second(HLL *z, + void lg_hll_second(HLL *z, bool *x, // unused const HLL *y) { memcpy(z->registers, y->registers, 1 << 10); }, -FDB_HLL_SECOND_STR) - -int64_t print_hll(char *string, // value is printed to the string - size_t string_size, // size of the string array - const void *value, // HLL value to print - int verbose // if >0, print verbosely; else tersely - ) { - const HLL *hll = (const HLL *)value; - - if (verbose > 0) { - return snprintf(string, string_size, "HLL{bits=%u, size=%u, count=%.2f}", - HLL_P, HLL_REGISTERS, _hll_count(hll)); - } - - return snprintf(string, string_size, "HLL{count=%.2f}", _hll_count(hll)); -} +LG_HLL_SECOND_STR) + +// int64_t print_hll ( +// char *string, // value is printed to the string +// size_t string_size, // size of the string array +// const void *value, // HLL value to print +// int verbose // if >0, print verbosely; else tersely +// ) { +// const HLL *hll = (const HLL *)value; +// +// if (verbose > 0) { +// return snprintf(string, string_size, "HLL{bits=%u, size=%u, count=%.2f}", +// HLL_P, HLL_REGISTERS, _hll_count(hll)) ; +// } +// +// return snprintf(string, string_size, "HLL{count=%.2f}", _hll_count(hll)) ; +// } #undef LG_FREE_WORK #define LG_FREE_WORK \ @@ -246,13 +265,14 @@ GrB_free(scores) ; // at BFS level d, sketches are propagated along edges // the harmonic contribution delta/d is accumulated into each node's score, // where delta is the number of new nodes discovered at that level +#if LG_SUITESPARSE_GRAPHBLAS_V10 int LAGr_HarmonicCentrality( // outputs: GrB_Vector *scores, // FP64 scores by original node ID GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- // able nodes from the given node. // inputs: - const GrB_Matrix A, // adjacency matrix + const LAGraph_Graph G, // input graph const GrB_Vector node_weights, // participating nodes and their weights char *msg ) { @@ -279,13 +299,15 @@ int LAGr_HarmonicCentrality( // TODO: LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); - LG_ASSERT(A != NULL, GrB_NULL_POINTER); + + LG_ASSERT(G != NULL, GrB_NULL_POINTER); + LG_ASSERT(G->A != NULL, GrB_NULL_POINTER); LG_ASSERT(scores != NULL, GrB_NULL_POINTER); LG_ASSERT(node_weights != NULL, GrB_NULL_POINTER); - GRB_TRY(GrB_Vector_size(&nrows, node_weights)); - GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)); - GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)); + GRB_TRY(GrB_Vector_size(&nrows, node_weights)) ; + GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)) ; + GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)) ; if (nvals == 0) { return GrB_SUCCESS; @@ -293,16 +315,16 @@ int LAGr_HarmonicCentrality( // double check weight type and maximum weight requirements GrB_Type weight_t = NULL; - GRB_TRY(GxB_Vector_type(&weight_t, node_weights)); + GRB_TRY(GxB_Vector_type(&weight_t, node_weights)) ; LG_ASSERT_MSG(weight_t == GrB_BOOL || weight_t == GrB_INT64, GrB_DOMAIN_MISMATCH, "weight must be integer or uint64"); if (weight_t == GrB_INT64) { int64_t max_w = 0, min_w = 0; GRB_TRY(GrB_Vector_reduce_INT64(&max_w, NULL, GrB_MAX_MONOID_INT64, - node_weights, NULL)); + node_weights, NULL)) ; GRB_TRY(GrB_Vector_reduce_INT64(&min_w, NULL, GrB_MIN_MONOID_INT64, - node_weights, NULL)); + node_weights, NULL)) ; LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, "Negative node weights not supported"); // TODO: is this cap reasonable? @@ -314,81 +336,81 @@ int LAGr_HarmonicCentrality( // build compact submatrix _A //-------------------------------------------------------------------------- - GRB_TRY(GrB_Descriptor_new(&desc)); - GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)); - GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)); - GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)); - GRB_TRY(GxB_Matrix_extract_Vector(_A, NULL, NULL, A, node_weights, - node_weights, desc)); - GRB_TRY(GrB_free(&desc)); + GRB_TRY (GrB_Descriptor_new (&desc)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)) ; + GRB_TRY (GrB_set (desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + GRB_TRY (GrB_Matrix_new (&_A, GrB_BOOL, nvals, nvals)) ; + GRB_TRY (GxB_Matrix_extract_Vector ( + _A, NULL, NULL, G->A, node_weights, node_weights, desc)) ; + GRB_TRY(GrB_free(&desc)) ; //-------------------------------------------------------------------------- // create scores vector (0.0 at each participating node) //-------------------------------------------------------------------------- - GRB_TRY(GrB_Vector_assign_FP64(*scores, node_weights, NULL, 0.0, GrB_ALL, 0, - GrB_DESC_S)); - GRB_TRY(GrB_set(*scores, GxB_SPARSE | GxB_FULL, GxB_SPARSITY_CONTROL)); - GRB_TRY(GxB_Container_new(&score_cont)); - GRB_TRY(GxB_unload_Vector_into_Container(*scores, score_cont, NULL)); + GRB_TRY (GrB_Vector_assign_FP64 ( + *scores, node_weights, NULL, 0.0, GrB_ALL, 0, GrB_DESC_S)) ; + GRB_TRY (GrB_set (*scores, GxB_SPARSE | GxB_FULL, GxB_SPARSITY_CONTROL)) ; + GRB_TRY (GxB_Container_new (&score_cont)) ; + GRB_TRY (GxB_unload_Vector_into_Container (*scores, score_cont, NULL)) ; //-------------------------------------------------------------------------- // Initialize HLL //-------------------------------------------------------------------------- - GRB_TRY(GxB_Type_new(&hll_t, sizeof(HLL), "HLL", HLL_jit)); - GRB_TRY(GrB_Vector_new(&new_sets, hll_t, nvals)); - GRB_TRY(GrB_Vector_new(&old_sets, hll_t, nvals)); - GRB_TRY(GxB_Container_new(&old_cont)); - GRB_TRY(GrB_Vector_new(&flat_scores, GrB_FP64, nvals)); - GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)); - GRB_TRY( - GxB_Vector_extractTuples_Vector(NULL, flat_weight, node_weights, NULL)); + GRB_TRY (GxB_Type_new (&hll_t, sizeof(HLL), "HLL", HLL_jit)) ; + GRB_TRY (GrB_Vector_new (&new_sets, hll_t, nvals)) ; + GRB_TRY (GrB_Vector_new (&old_sets, hll_t, nvals)) ; + GRB_TRY (GxB_Container_new (&old_cont)) ; + GRB_TRY (GrB_Vector_new (&flat_scores, GrB_FP64, nvals)) ; + GRB_TRY (GrB_Vector_new (&flat_weight, GrB_INT64, nvals)) ; + GRB_TRY (GxB_Vector_extractTuples_Vector ( + NULL, flat_weight, node_weights, NULL)) ; // init op: weight (INT64) at row index i → HLL seeded with 'weight' hashes - GRB_TRY(GrB_IndexUnaryOp_new(&init_hlls, - (GxB_index_unary_function)fdb_hll_init, hll_t, - GrB_INT64, GrB_BOOL)); + GRB_TRY(GrB_IndexUnaryOp_new( + &init_hlls, (GxB_index_unary_function)lg_hll_init, hll_t, + GrB_INT64, GrB_BOOL)) ; // merge binary op (HLL, HLL) -> HLL in-place: z == x required - GRB_TRY(GxB_BinaryOp_new(&merge_hll_biop, (GxB_binary_function)fdb_hll_merge, - hll_t, hll_t, hll_t, "fdb_hll_merge", - FDB_HLL_MERGE_STR)); + GRB_TRY(GxB_BinaryOp_new( + &merge_hll_biop, (GxB_binary_function)lg_hll_merge, + hll_t, hll_t, hll_t, "lg_hll_merge", LG_HLL_MERGE_STR)) ; // second op - GRB_TRY(GxB_BinaryOp_new(&shallow_second, (GxB_binary_function)fdb_hll_second, - hll_t, GrB_BOOL, hll_t, "fdb_hll_second", - FDB_HLL_SECOND_STR)); + GRB_TRY(GxB_BinaryOp_new( + &shallow_second, (GxB_binary_function)lg_hll_second, + hll_t, GrB_BOOL, hll_t, "lg_hll_second", LG_HLL_SECOND_STR)) ; // merge monoid - identity is an empty (all-zero) HLL sketch HLL hll_zero = {0}; - GRB_TRY(GrB_Monoid_new_UDT(&merge_hll, merge_hll_biop, &hll_zero)); + GRB_TRY(GrB_Monoid_new_UDT(&merge_hll, merge_hll_biop, &hll_zero)) ; // semiring: add = merge monoid, multiply = copy (pass-through second operand) - GRB_TRY(GrB_Semiring_new(&merge_second, merge_hll, shallow_second)); + GRB_TRY(GrB_Semiring_new(&merge_second, merge_hll, shallow_second)) ; // delta op: (HLL_old, HLL_new) → FP64 cardinality change - GRB_TRY(GxB_BinaryOp_new(&delta_hll, (GxB_binary_function)fdb_hll_delta, - GrB_FP64, hll_t, hll_t, "fdb_hll_delta", - FDB_HLL_DELTA_STR)); + GRB_TRY(GxB_BinaryOp_new( + &delta_hll, (GxB_binary_function)lg_hll_delta, + GrB_FP64, hll_t, hll_t, "lg_hll_delta", LG_HLL_DELTA_STR)) ; // delta output: one FP64 entry per participating node - GRB_TRY(GrB_Vector_new(&delta_vec, GrB_FP64, nvals)); + GRB_TRY(GrB_Vector_new(&delta_vec, GrB_FP64, nvals)) ; //-------------------------------------------------------------------------- // Load vector values //-------------------------------------------------------------------------- // initialize HLLs GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(new_sets, NULL, NULL, init_hlls, - flat_weight, false, NULL)); + flat_weight, false, NULL)) ; GRB_TRY(GrB_Vector_apply_IndexOp_BOOL(old_sets, NULL, NULL, init_hlls, - flat_weight, false, NULL)); - GRB_TRY(GrB_set(old_sets, GxB_BITMAP, GxB_SPARSITY_CONTROL)); + flat_weight, false, NULL)) ; + GRB_TRY(GrB_set(old_sets, GxB_BITMAP, GxB_SPARSITY_CONTROL)) ; - GRB_TRY(GrB_free(&flat_weight)); - GRB_TRY( - GrB_Vector_assign_FP64(flat_scores, NULL, NULL, 0.0, GrB_ALL, 0, NULL)); + GRB_TRY(GrB_free(&flat_weight)) ; + GRB_TRY (GrB_Vector_assign_FP64 ( + flat_scores, NULL, NULL, 0.0, GrB_ALL, 0, NULL)) ; //-------------------------------------------------------------------------- // HLL BFS propagation @@ -401,28 +423,28 @@ int LAGr_HarmonicCentrality( // foward bfs // merge each neighbor's pre-round set into this node's set // target kernel for inplace adjustments - GRB_TRY(GrB_mxv(new_sets, NULL, merge_hll_biop, merge_second, _A, old_sets, - NULL)); + GRB_TRY(GrB_mxv( + new_sets, NULL, merge_hll_biop, merge_second, _A, old_sets, NULL)) ; - GRB_TRY(GxB_unload_Vector_into_Container(old_sets, old_cont, NULL)); + GRB_TRY(GxB_unload_Vector_into_Container(old_sets, old_cont, NULL)) ; // find the delta between last round and this one - GRB_TRY(GrB_eWiseMult(delta_vec, NULL, NULL, delta_hll, new_sets, - old_cont->x, NULL)); - GRB_TRY( - GrB_Vector_assign(old_cont->x, NULL, NULL, new_sets, GrB_ALL, 0, NULL)); + GRB_TRY(GrB_eWiseMult( + delta_vec, NULL, NULL, delta_hll, new_sets, old_cont->x, NULL)) ; + GRB_TRY(GrB_Vector_assign( + old_cont->x, NULL, NULL, new_sets, GrB_ALL, 0, NULL)) ; // old_set bitmap is the set of nodes with non-zero deltas - GRB_TRY( - GrB_apply(old_cont->b, NULL, NULL, GrB_IDENTITY_BOOL, delta_vec, NULL)); - GRB_TRY(GrB_Vector_reduce_INT64(&changes, NULL, GrB_PLUS_MONOID_INT64, - old_cont->b, NULL)); + GRB_TRY(GrB_apply( + old_cont->b, NULL, NULL, GrB_IDENTITY_BOOL, delta_vec, NULL)) ; + GRB_TRY(GrB_Vector_reduce_INT64( + &changes, NULL, GrB_PLUS_MONOID_INT64, old_cont->b, NULL)) ; old_cont->nvals = changes; - GRB_TRY(GxB_load_Vector_from_Container(old_sets, old_cont, NULL)); + GRB_TRY(GxB_load_Vector_from_Container(old_sets, old_cont, NULL)) ; // use the deltas to update the score - GRB_TRY(GrB_apply(flat_scores, NULL, GrB_PLUS_FP64, GrB_DIV_FP64, delta_vec, - (double)d, NULL)); + GRB_TRY(GrB_apply(flat_scores, NULL, GrB_PLUS_FP64, + GrB_DIV_FP64, delta_vec, (double)d, NULL)) ; // stop when no HLL cardinality changed this round if (changes == 0) { @@ -433,11 +455,11 @@ int LAGr_HarmonicCentrality( //-------------------------------------------------------------------------- // write flat_scores into score_cont->x, clear iso, reload scores vector //-------------------------------------------------------------------------- - GRB_TRY (GrB_free (&score_cont->x)); + GRB_TRY (GrB_free (&score_cont->x)) ; score_cont->iso = false; score_cont->x = flat_scores; flat_scores = NULL ; - GRB_TRY (GxB_load_Vector_from_Container (*scores, score_cont, NULL)); + GRB_TRY (GxB_load_Vector_from_Container (*scores, score_cont, NULL)) ; //-------------------------------------------------------------------------- // cleanup @@ -448,8 +470,20 @@ int LAGr_HarmonicCentrality( LG_FREE_WORK; return GrB_SUCCESS; } - -// TODO: add a descriptor? let nodes be values or indecies +#else +int LAGr_HarmonicCentrality( + // outputs: + GrB_Vector *scores, // FP64 scores by original node ID + GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- + // able nodes from the given node. + // inputs: + const LAGraph_Graph G, // input graph + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) { + return (GrB_NOT_IMPLEMENTED) ; +} +#endif #undef LG_FREE_WORK #define LG_FREE_WORK \ @@ -467,13 +501,14 @@ int LAGr_HarmonicCentrality( GrB_free(scores) ; // This is an exact (but much slower) calculation of harmonic centrality +#if LG_SUITESPARSE_GRAPHBLAS_V10 int LAGr_HarmonicCentrality_exact( // outputs: GrB_Vector *scores, // FP64 scores by original node ID GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- // able nodes from the given node. // inputs: - const GrB_Matrix A, // adjacency matrix + const LAGraph_Graph G, // input graph const GrB_Vector nodes, // nodes to calculate centrality of const GrB_Vector node_weights, // participating nodes and their weights char *msg @@ -492,13 +527,14 @@ int LAGr_HarmonicCentrality_exact( // TODO: LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); - LG_ASSERT(A != NULL, GrB_NULL_POINTER); + LG_ASSERT(G != NULL, GrB_NULL_POINTER); + LG_ASSERT(G->A != NULL, GrB_NULL_POINTER); LG_ASSERT(scores != NULL, GrB_NULL_POINTER); LG_ASSERT(node_weights != NULL, GrB_NULL_POINTER); - GRB_TRY(GrB_Vector_size(&nrows, node_weights)); - GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)); - GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)); + GRB_TRY(GrB_Vector_size(&nrows, node_weights)) ; + GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)) ; + GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)) ; if (nvals == 0) { return GrB_SUCCESS; @@ -506,14 +542,14 @@ int LAGr_HarmonicCentrality_exact( // double check weight type and maximum weight requirements GrB_Type weight_t = NULL; - GRB_TRY(GxB_Vector_type(&weight_t, node_weights)); + GRB_TRY(GxB_Vector_type(&weight_t, node_weights)) ; LG_ASSERT_MSG(weight_t == GrB_BOOL || weight_t == GrB_INT64, - GrB_DOMAIN_MISMATCH, "weight must be integer or uint64"); + GrB_DOMAIN_MISMATCH, "weight must be integer or uint64"); if (weight_t == GrB_INT64) { int64_t min_w = 0; - GRB_TRY(GrB_Vector_reduce_INT64(&min_w, NULL, GrB_MIN_MONOID_INT64, - node_weights, NULL)); + GRB_TRY(GrB_Vector_reduce_INT64( + &min_w, NULL, GrB_MIN_MONOID_INT64, node_weights, NULL)) ; LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, "Negative node weights not supported"); } @@ -522,40 +558,40 @@ int LAGr_HarmonicCentrality_exact( // build compact submatrix _A //-------------------------------------------------------------------------- - GRB_TRY(GrB_Descriptor_new(&desc)); - GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)); - GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)); - GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)); - GRB_TRY(GxB_Matrix_extract_Vector(_A, NULL, NULL, A, node_weights, - node_weights, desc)); - GRB_TRY(GrB_free(&desc)); - LG_TRY (LAGraph_New (&G_compact, &_A, LAGraph_ADJACENCY_DIRECTED, msg)); + GRB_TRY(GrB_Descriptor_new(&desc)) ; + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_COLINDEX_LIST)) ; + GRB_TRY(GrB_set(desc, GxB_USE_INDICES, GxB_ROWINDEX_LIST)) ; + GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)) ; + GRB_TRY(GxB_Matrix_extract_Vector( + _A, NULL, NULL, G->A, node_weights, node_weights, desc)) ; + GRB_TRY(GrB_free(&desc)) ; + LG_TRY (LAGraph_New (&G_compact, &_A, LAGraph_ADJACENCY_DIRECTED, msg)) ; - GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)); + GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)) ; GRB_TRY( GxB_Vector_extractTuples_Vector( - NULL, flat_weight, node_weights, NULL)); + NULL, flat_weight, node_weights, NULL)) ; //-------------------------------------------------------------------------- // calculate centrality via BFS //-------------------------------------------------------------------------- - GRB_TRY (GxB_Scalar_new (&score, GrB_FP64)); - GRB_TRY(GxB_Iterator_new(&it)); - GRB_TRY(GxB_Vector_Iterator_attach(it, nodes, NULL)); + GRB_TRY (GxB_Scalar_new (&score, GrB_FP64)) ; + GRB_TRY(GxB_Iterator_new(&it)) ; + GRB_TRY(GxB_Vector_Iterator_attach(it, nodes, NULL)) ; GrB_Info info = GxB_Vector_Iterator_seek(it, 0); while (info != GxB_EXHAUSTED) { GrB_Index nodeIdx = GxB_Vector_Iterator_getIndex(it); - LG_TRY (LAGr_BreadthFirstSearch(&level, NULL, G_compact, nodeIdx, msg)); + LG_TRY (LAGr_BreadthFirstSearch(&level, NULL, G_compact, nodeIdx, msg)) ; // remove 0 to prevent div by 0 - GRB_TRY (GrB_Vector_removeElement (level, nodeIdx)); + GRB_TRY (GrB_Vector_removeElement (level, nodeIdx)) ; // dot product (plus div) GRB_TRY(GrB_vxm((GrB_Vector) score, NULL, NULL, GxB_PLUS_DIV_FP64, - flat_weight, (GrB_Matrix) level, NULL)); + flat_weight, (GrB_Matrix) level, NULL)) ; GrB_free (&level); - GRB_TRY (GrB_Vector_setElement_Scalar(*scores, score, nodeIdx)); + GRB_TRY (GrB_Vector_setElement_Scalar(*scores, score, nodeIdx)) ; info = GxB_Vector_Iterator_next(it); } @@ -569,3 +605,18 @@ int LAGr_HarmonicCentrality_exact( LG_FREE_WORK; return GrB_SUCCESS; } +#else +int LAGr_HarmonicCentrality_exact( + // outputs: + GrB_Vector *scores, // FP64 scores by original node ID + GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- + // able nodes from the given node. + // inputs: + const LAGraph_Graph G, // input graph + const GrB_Vector nodes, // nodes to calculate centrality of + const GrB_Vector node_weights, // participating nodes and their weights + char *msg +) { + return (GrB_NOT_IMPLEMENTED) ; +} +#endif diff --git a/experimental/benchmark/harmonicCentrality_demo.c b/experimental/benchmark/harmonicCentrality_demo.c index 3f3845b5ff..d060dcf3b1 100644 --- a/experimental/benchmark/harmonicCentrality_demo.c +++ b/experimental/benchmark/harmonicCentrality_demo.c @@ -40,6 +40,7 @@ int main (int argc, char **argv) char msg [LAGRAPH_MSG_LEN] ; GrB_Matrix A = NULL ; + LAGraph_Graph G = NULL ; GrB_Vector scores_approx = NULL ; GrB_Vector scores_exact = NULL ; GrB_Vector node_weights = NULL ; @@ -80,6 +81,9 @@ int main (int argc, char **argv) printf ("Graph: %s (%" PRIu64 " nodes, %" PRIu64 " edges)\n", matrix_name, (uint64_t) n, (uint64_t) nvals) ; + // construct a graph + LAGRAPH_TRY (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + //-------------------------------------------------------------------------- // create boolean node_weights (all nodes participate with weight = 1) //-------------------------------------------------------------------------- @@ -94,7 +98,7 @@ int main (int argc, char **argv) t = LAGraph_WallClockTime ( ) ; LAGRAPH_TRY (LAGr_HarmonicCentrality ( - &scores_approx, NULL, A, node_weights, msg)) ; + &scores_approx, NULL, G, node_weights, msg)) ; t = LAGraph_WallClockTime ( ) - t ; printf ("Time for LAGr_HarmonicCentrality (approx): %g sec\n", t) ; @@ -111,7 +115,7 @@ int main (int argc, char **argv) { t = LAGraph_WallClockTime ( ) ; LAGRAPH_TRY (LAGr_HarmonicCentrality_exact ( - &scores_exact, NULL, A, node_weights, node_weights, msg)) ; + &scores_exact, NULL, G, node_weights, node_weights, msg)) ; t = LAGraph_WallClockTime ( ) - t ; printf ("\nTime for LAGr_HarmonicCentrality_exact: %g sec\n", t) ; @@ -124,6 +128,7 @@ int main (int argc, char **argv) //-------------------------------------------------------------------------- LG_FREE_ALL ; + LAGraph_Delete (&G, msg) ; LAGRAPH_TRY (LAGraph_Finalize (msg)) ; return (GrB_SUCCESS) ; } diff --git a/experimental/test/test_HarmonicCentrality.c b/experimental/test/test_HarmonicCentrality.c index 782bc60d96..6f1fba4a35 100644 --- a/experimental/test/test_HarmonicCentrality.c +++ b/experimental/test/test_HarmonicCentrality.c @@ -26,6 +26,7 @@ char msg [LAGRAPH_MSG_LEN] ; GrB_Matrix A = NULL ; +LAGraph_Graph G = NULL ; GrB_Vector scores_approx = NULL ; GrB_Vector scores_exact = NULL ; GrB_Vector node_weights = NULL ; @@ -80,9 +81,9 @@ double max_relative_error // test_HarmonicCentrality: compare approximate vs exact on small graphs //------------------------------------------------------------------------------ +#if LG_SUITESPARSE_GRAPHBLAS_V10 void test_HarmonicCentrality (void) { - #if LAGRAPH_SUITESPARSE LAGraph_Init (msg) ; for (int k = 0 ; ; k++) @@ -102,6 +103,10 @@ void test_HarmonicCentrality (void) GrB_Index n ; OK (GrB_Matrix_nrows (&n, A)) ; + // construct a graph + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; // A has been moved into G->A + // create boolean node_weights (all nodes, weight = true) OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; OK (GrB_Vector_assign_BOOL ( @@ -109,12 +114,11 @@ void test_HarmonicCentrality (void) // compute approximate harmonic centrality OK (LAGr_HarmonicCentrality ( - &scores_approx, NULL, A, node_weights, msg)) ; - printf ("%s", msg); + &scores_approx, NULL, G, node_weights, msg)) ; // compute exact harmonic centrality OK (LAGr_HarmonicCentrality_exact ( - &scores_exact, NULL, A, node_weights, node_weights, msg)) ; + &scores_exact, NULL, G, node_weights, node_weights, msg)) ; // print results for small graphs LAGraph_PrintLevel pr = @@ -133,23 +137,24 @@ void test_HarmonicCentrality (void) TEST_MSG ("Relative error too large: %.2f%%", rel_err * 100) ; // cleanup for this iteration - OK (GrB_free (&A)) ; + OK (LAGraph_Delete (&G, msg)) ; OK (GrB_free (&scores_approx)) ; OK (GrB_free (&scores_exact)) ; OK (GrB_free (&node_weights)) ; } LAGraph_Finalize (msg) ; - #endif } +#endif + //------------------------------------------------------------------------------ // test_HarmonicCentrality_empty: test with empty node set //------------------------------------------------------------------------------ +#if LG_SUITESPARSE_GRAPHBLAS_V10 void test_HarmonicCentrality_empty (void) { - #if LAGRAPH_SUITESPARSE LAGraph_Init (msg) ; snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; @@ -161,23 +166,28 @@ void test_HarmonicCentrality_empty (void) GrB_Index n ; OK (GrB_Matrix_nrows (&n, A)) ; + // construct a graph + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; // A has been moved into G->A + // empty node_weights — no participating nodes OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; - OK (LAGr_HarmonicCentrality (&scores_approx, NULL, A, node_weights, msg)) ; + OK (LAGr_HarmonicCentrality (&scores_approx, NULL, G, node_weights, msg)) ; // scores should exist but have no entries GrB_Index nvals ; OK (GrB_Vector_nvals (&nvals, scores_approx)) ; TEST_CHECK (nvals == 0) ; - OK (GrB_free (&A)) ; + OK (LAGraph_Delete (&G, msg)) ; OK (GrB_free (&scores_approx)) ; OK (GrB_free (&node_weights)) ; LAGraph_Finalize (msg) ; - #endif } +#endif + //------------------------------------------------------------------------------ // test_errors: test error handling @@ -185,7 +195,6 @@ void test_HarmonicCentrality_empty (void) void test_errors (void) { - #if LAGRAPH_SUITESPARSE LAGraph_Init (msg) ; snprintf (filename, LEN, LG_DATA_DIR "%s", "karate.mtx") ; @@ -196,18 +205,24 @@ void test_errors (void) GrB_Index n ; OK (GrB_Matrix_nrows (&n, A)) ; + + // construct a graph + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; // A has been moved into G->A + OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; OK (GrB_Vector_assign_BOOL ( node_weights, NULL, NULL, true, GrB_ALL, n, NULL)) ; int result ; + #if LG_SUITESPARSE_GRAPHBLAS_V10 // scores is NULL - result = LAGr_HarmonicCentrality (NULL, NULL, A, node_weights, msg) ; + result = LAGr_HarmonicCentrality (NULL, NULL, G, node_weights, msg) ; printf ("\nresult: %d %s\n", result, msg) ; TEST_CHECK (result == GrB_NULL_POINTER) ; - // A is NULL + // G is NULL result = LAGr_HarmonicCentrality ( &scores_approx, NULL, NULL, node_weights, msg) ; printf ("\nresult: %d %s\n", result, msg) ; @@ -215,22 +230,35 @@ void test_errors (void) // node_weights is NULL result = LAGr_HarmonicCentrality ( - &scores_approx, NULL, A, NULL, msg) ; + &scores_approx, NULL, G, NULL, msg) ; printf ("\nresult: %d %s\n", result, msg) ; TEST_CHECK (result == GrB_NULL_POINTER) ; + #else + // Below V10, functions should return GrB_NOT_IMPLEMENTED + result = LAGr_HarmonicCentrality ( + &scores_approx, NULL, G, node_weights, msg) ; + printf ("\nresult: %d (expected GrB_NOT_IMPLEMENTED)\n", result) ; + TEST_CHECK (result == GrB_NOT_IMPLEMENTED) ; + + result = LAGr_HarmonicCentrality_exact ( + &scores_approx, NULL, G, node_weights, node_weights, msg) ; + printf ("\nresult: %d (expected GrB_NOT_IMPLEMENTED)\n", result) ; + TEST_CHECK (result == GrB_NOT_IMPLEMENTED) ; + #endif - OK (GrB_free (&A)) ; + OK (LAGraph_Delete (&G, msg)) ; OK (GrB_free (&node_weights)) ; LAGraph_Finalize (msg) ; - #endif } //**************************************************************************** TEST_LIST = { + #if LG_SUITESPARSE_GRAPHBLAS_V10 {"HarmonicCentrality", test_HarmonicCentrality}, {"HarmonicCentrality_empty", test_HarmonicCentrality_empty}, + #endif {"HarmonicCentrality_errors", test_errors}, {NULL, NULL} } ; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index 6d9be99c38..8e72e4b16c 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1329,7 +1329,7 @@ int LAGr_HarmonicCentrality GrB_Vector *reachable_nodes, // [optional] estimated reachable node count // (pass NULL, not yet implemented) // inputs: - const GrB_Matrix A, // adjacency matrix + const LAGraph_Graph G, // input graph const GrB_Vector node_weights, // participating nodes and their weights char *msg ) ; @@ -1346,7 +1346,7 @@ int LAGr_HarmonicCentrality_exact GrB_Vector *reachable_nodes, // [optional] estimated reachable node count // (pass NULL, not yet implemented) // inputs: - const GrB_Matrix A, // adjacency matrix + const LAGraph_Graph G, // input graph const GrB_Vector nodes, // nodes to calculate centrality of const GrB_Vector node_weights, // participating nodes and their weights char *msg From 39868a214630e553f2c2d839ed53ae7dd0062a59 Mon Sep 17 00:00:00 2001 From: GomezGab Date: Tue, 7 Apr 2026 01:08:36 -0500 Subject: [PATCH 3/5] added return for the number of reachable nodes --- .../algorithm/LAGr_HarmonicCentrality.c | 97 +++++++++++-------- experimental/test/test_HarmonicCentrality.c | 80 +++++++++++++++ 2 files changed, 136 insertions(+), 41 deletions(-) diff --git a/experimental/algorithm/LAGr_HarmonicCentrality.c b/experimental/algorithm/LAGr_HarmonicCentrality.c index ca97ecfb8e..c8d661ad97 100644 --- a/experimental/algorithm/LAGr_HarmonicCentrality.c +++ b/experimental/algorithm/LAGr_HarmonicCentrality.c @@ -17,7 +17,7 @@ // The LAGr_HarmonicCentrality algorithm calculates an estimate of the "harmonic // centrality" of all nodes in a graph. Harmonic Centrality is defined on a non- // weighted graph as follows: -// HC(u) = \sum_{v\in neighbors(u)} d(u,v) +// HC(u) = \sum_{v \ne u} 1/d(u,v) // where d(u,v) is the shortest path distance from u to v. // // HyperLogLog allows us to estimate the cardinality of the subsets which each @@ -47,16 +47,10 @@ #include "LG_internal.h" #define CENTRALITY_MAX_ITER 100 - -// WARNING: these macros are not recognized by the JIT, so have to be changed -// manually in the functions below #define HLL_P 10 // Precision #define HLL_REGISTERS (1 << HLL_P) // number of registers -// register a function (or type) and the jit string for its definition -#define JIT_STR(f, var) static char *var = #f; f - -JIT_STR( +LG_JIT_STRING( typedef struct { uint8_t registers[(1 << 10)]; } HLL; @@ -84,8 +78,9 @@ static __inline void _hll_add_hash(HLL *hll, uint32_t hash) { } } -double _hll_count(const HLL *hll) { - uint32_t i; +LG_JIT_STRING ( +void lg_hll_count(double *z, const HLL *x) { + const HLL *hll = x; double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)HLL_REGISTERS); @@ -101,7 +96,7 @@ double _hll_count(const HLL *hll) { if (estimate <= 5.0 / 2.0 * (double)HLL_REGISTERS) { int zeros = 0; - for (i = 0; i < HLL_REGISTERS; i++) + for (uint32_t i = 0; i < HLL_REGISTERS; i++) zeros += (hll->registers[i] == 0); if (zeros) @@ -111,8 +106,8 @@ double _hll_count(const HLL *hll) { estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; } - return estimate; -} + *z = estimate; +}, LG_HLL_COUNT) //------------------------------------------------------------------------------ // GraphBLAS Ops @@ -145,41 +140,41 @@ void lg_hll_init(HLL *z, const uint64_t *x, GrB_Index i, GrB_Index j, } } -JIT_STR( +LG_JIT_STRING( void lg_hll_merge(HLL *z, const HLL *x, const HLL *y) { - for (uint32_t i = 0; i < (1 << 10); i++) { - z->registers[i] = y->registers[i] > x->registers[i] ? y->registers[i] - : x->registers[i]; + for (uint32_t i = 0; i < HLL_REGISTERS; i++) { + z->registers[i] = y->registers[i] > x->registers[i] ? + y->registers[i] : x->registers[i]; } }, LG_HLL_MERGE_STR) -JIT_STR( +LG_JIT_STRING( void lg_hll_delta(double *z, const HLL *x, const HLL *y) { *z = 0; - bool diff = 0 != memcmp(x->registers, y->registers, (1 << 10)) ; + bool diff = 0 != memcmp(x->registers, y->registers, HLL_REGISTERS) ; if (diff) { uint32_t i; double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(1 << 10)) ; - alpha_mm *= ((double)(1 << 10) * (double)(1 << 10)) ; + alpha_mm *= ((double)(HLL_REGISTERS) * (double)(HLL_REGISTERS)) ; double sum = 0; - for (uint32_t i = 0; i < (1 << 10); i++) { + for (uint32_t i = 0; i < (HLL_REGISTERS); i++) { sum += 1.0 / (1 << x->registers[i]); } double estimate = alpha_mm / sum; - if (estimate <= 5.0 / 2.0 * (double)(1 << 10)) { + if (estimate <= 5.0 / 2.0 * (double)(HLL_REGISTERS)) { int zeros = 0; - for (i = 0; i < (1 << 10); i++) + for (i = 0; i < (HLL_REGISTERS); i++) zeros += (x->registers[i] == 0); if (zeros) - estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); + estimate = (double)(HLL_REGISTERS) * log((double)(HLL_REGISTERS) / zeros); } else if (estimate > (1.0 / 30.0) * 4294967296.0) { estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; @@ -188,20 +183,20 @@ JIT_STR( *z = estimate; sum = 0; - for (uint32_t i = 0; i < (1 << 10); i++) { + for (uint32_t i = 0; i < (HLL_REGISTERS); i++) { sum += 1.0 / (1 << y->registers[i]); } estimate = alpha_mm / sum; - if (estimate <= 5.0 / 2.0 * (double)(1 << 10)) { + if (estimate <= 5.0 / 2.0 * (double)(HLL_REGISTERS)) { int zeros = 0; - for (i = 0; i < (1 << 10); i++) + for (i = 0; i < (HLL_REGISTERS); i++) zeros += (y->registers[i] == 0); if (zeros) - estimate = (double)(1 << 10) * log((double)(1 << 10) / zeros); + estimate = (double)(HLL_REGISTERS) * log((double)(HLL_REGISTERS) / zeros); } else if (estimate > (1.0 / 30.0) * 4294967296.0) { estimate = -4294967296.0 * log(1.0 - (estimate / 4294967296.0)) ; @@ -211,11 +206,9 @@ JIT_STR( }, LG_HLL_DELTA_STR) -JIT_STR( - void lg_hll_second(HLL *z, - bool *x, // unused - const HLL *y) { - memcpy(z->registers, y->registers, 1 << 10); +LG_JIT_STRING( + void lg_hll_second(HLL *z, bool *x, const HLL *y) { + memcpy(z->registers, y->registers, sizeof(z->registers)); }, LG_HLL_SECOND_STR) @@ -252,12 +245,14 @@ GrB_free(&shallow_second); \ GrB_free(&merge_hll_biop); \ GrB_free(&merge_hll); \ GrB_free(&merge_second); \ -GrB_free(&delta_hll); +GrB_free(&delta_hll); \ +GrB_free(&count_hll); #undef LG_FREE_ALL #define LG_FREE_ALL \ LG_FREE_WORK ; \ -GrB_free(scores) ; +GrB_free(scores) ; \ +GrB_free(reachable_nodes) ; // compute harmonic closeness centrality estimates using HLL BFS propagation // @@ -291,15 +286,13 @@ int LAGr_HarmonicCentrality( GrB_Descriptor desc = NULL; GrB_IndexUnaryOp init_hlls = NULL; + GrB_UnaryOp count_hll = NULL; GrB_BinaryOp shallow_second = NULL; GrB_BinaryOp merge_hll_biop = NULL; GrB_Monoid merge_hll = NULL; GrB_Semiring merge_second = NULL; GrB_BinaryOp delta_hll = NULL; - // TODO: - LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); - LG_ASSERT(G != NULL, GrB_NULL_POINTER); LG_ASSERT(G->A != NULL, GrB_NULL_POINTER); LG_ASSERT(scores != NULL, GrB_NULL_POINTER); @@ -327,6 +320,7 @@ int LAGr_HarmonicCentrality( node_weights, NULL)) ; LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, "Negative node weights not supported"); + // TODO: is this cap reasonable? LG_ASSERT_MSG(max_w < 1000000, GrB_NOT_IMPLEMENTED, "Node weights over 1000000 not supported"); @@ -367,6 +361,9 @@ int LAGr_HarmonicCentrality( GRB_TRY (GxB_Vector_extractTuples_Vector ( NULL, flat_weight, node_weights, NULL)) ; + // count op + GRB_TRY (GxB_UnaryOp_new (&count_hll, (GxB_unary_function) lg_hll_count, + GrB_FP64, hll_t, "lg_hll_count", LG_HLL_COUNT)) ; // init op: weight (INT64) at row index i → HLL seeded with 'weight' hashes GRB_TRY(GrB_IndexUnaryOp_new( @@ -459,6 +456,17 @@ int LAGr_HarmonicCentrality( score_cont->iso = false; score_cont->x = flat_scores; flat_scores = NULL ; + + if (reachable_nodes) { + // make reachable_nodes vector with same sparsity pattern as scores + GRB_TRY (GrB_Vector_new (reachable_nodes, GrB_FP64, nrows)) ; + GRB_TRY (GrB_apply (delta_vec, NULL, NULL, count_hll, new_sets, NULL)) ; + GrB_Vector I_vec = (score_cont->format == GxB_FULL) ? + NULL : score_cont->i; + GRB_TRY (GxB_Vector_assign_Vector ( + *reachable_nodes, NULL, NULL, delta_vec, I_vec, NULL)) ; + } + GRB_TRY (GxB_load_Vector_from_Container (*scores, score_cont, NULL)) ; //-------------------------------------------------------------------------- @@ -470,6 +478,7 @@ int LAGr_HarmonicCentrality( LG_FREE_WORK; return GrB_SUCCESS; } + #else int LAGr_HarmonicCentrality( // outputs: @@ -493,6 +502,7 @@ int LAGr_HarmonicCentrality( GrB_free(&it); \ GrB_free(&score); \ GrB_free(&desc); \ + GrB_free(&node_compact); \ LAGraph_Delete(&G_compact, NULL); #undef LG_FREE_ALL @@ -519,6 +529,7 @@ int LAGr_HarmonicCentrality_exact( GrB_Vector level = NULL; GrB_Vector flat_weight = NULL; + GrB_Vector node_compact = NULL; GxB_Iterator it = NULL; GrB_Scalar score = NULL; @@ -564,7 +575,6 @@ int LAGr_HarmonicCentrality_exact( GRB_TRY(GrB_Matrix_new(&_A, GrB_BOOL, nvals, nvals)) ; GRB_TRY(GxB_Matrix_extract_Vector( _A, NULL, NULL, G->A, node_weights, node_weights, desc)) ; - GRB_TRY(GrB_free(&desc)) ; LG_TRY (LAGraph_New (&G_compact, &_A, LAGraph_ADJACENCY_DIRECTED, msg)) ; GRB_TRY(GrB_Vector_new(&flat_weight, GrB_INT64, nvals)) ; @@ -575,8 +585,13 @@ int LAGr_HarmonicCentrality_exact( //-------------------------------------------------------------------------- GRB_TRY (GxB_Scalar_new (&score, GrB_FP64)) ; - GRB_TRY(GxB_Iterator_new(&it)) ; - GRB_TRY(GxB_Vector_Iterator_attach(it, nodes, NULL)) ; + GRB_TRY (GxB_Iterator_new(&it)) ; + GRB_TRY (GrB_Vector_new (&node_compact, GrB_INT64, nvals)); + GRB_TRY (GxB_Vector_extract_Vector( + node_compact, NULL, NULL, nodes, node_weights, desc)) ; + GRB_TRY(GrB_free(&desc)) ; + + GRB_TRY(GxB_Vector_Iterator_attach(it, node_compact, NULL)) ; GrB_Info info = GxB_Vector_Iterator_seek(it, 0); while (info != GxB_EXHAUSTED) { GrB_Index nodeIdx = GxB_Vector_Iterator_getIndex(it); diff --git a/experimental/test/test_HarmonicCentrality.c b/experimental/test/test_HarmonicCentrality.c index 6f1fba4a35..04f6440080 100644 --- a/experimental/test/test_HarmonicCentrality.c +++ b/experimental/test/test_HarmonicCentrality.c @@ -18,6 +18,7 @@ // NOTE: these tests require SuiteSparse:GraphBLAS +#include "LG_internal.h" #include #include @@ -30,6 +31,7 @@ LAGraph_Graph G = NULL ; GrB_Vector scores_approx = NULL ; GrB_Vector scores_exact = NULL ; GrB_Vector node_weights = NULL ; +GrB_Vector reachable_approx = NULL ; #define LEN 512 char filename [LEN+1] ; @@ -189,6 +191,83 @@ void test_HarmonicCentrality_empty (void) #endif +//------------------------------------------------------------------------------ +// test_HarmonicCentrality_reachable: compare reachable_nodes to BFS nvals +//------------------------------------------------------------------------------ + +#if LG_SUITESPARSE_GRAPHBLAS_V10 +void test_HarmonicCentrality_reachable (void) +{ + LAGraph_Init (msg) ; + + for (int k = 0 ; ; k++) + { + const char *aname = files [k] ; + if (strlen (aname) == 0) break ; + printf ("\n================================== %s:\n", aname) ; + TEST_CASE (aname) ; + snprintf (filename, LEN, LG_DATA_DIR "%s", aname) ; + FILE *f = fopen (filename, "r") ; + TEST_CHECK (f != NULL) ; + OK (LAGraph_MMRead (&A, f, msg)) ; + OK (fclose (f)) ; + + GrB_Index n ; + OK (GrB_Matrix_nrows (&n, A)) ; + + OK (LAGraph_New (&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)) ; + TEST_CHECK (A == NULL) ; + + OK (GrB_Vector_new (&node_weights, GrB_BOOL, n)) ; + OK (GrB_Vector_assign_BOOL ( + node_weights, NULL, NULL, true, GrB_ALL, n, NULL)) ; + + // compute approximate harmonic centrality with reachable_nodes output + GrB_Info info = LAGr_HarmonicCentrality ( + &scores_approx, &reachable_approx, G, node_weights, msg); + OK (info) ; + + // for each node, run BFS and compare nvals(level) to reachable_approx + double max_err = 0 ; + for (GrB_Index i = 0 ; i < n ; i++) + { + double approx_count = 0 ; + GrB_Info info = GrB_Vector_extractElement_FP64 ( + &approx_count, reachable_approx, i) ; + if (info != GrB_SUCCESS) continue ; + + GrB_Vector level = NULL ; + OK (LAGr_BreadthFirstSearch (&level, NULL, G, i, msg)) ; + + GrB_Index bfs_nvals = 0 ; + OK (GrB_Vector_nvals (&bfs_nvals, level)) ; + OK (GrB_free (&level)) ; + + // bfs_nvals includes the source node (level 0) + if (bfs_nvals > 0) + { + double rel = fabs (approx_count - (double) bfs_nvals) + / (double) bfs_nvals ; + if (rel > max_err) max_err = rel ; + } + } + + printf ("max relative error (reachable): %.2f%%\n", max_err * 100) ; + TEST_CHECK (max_err < 0.5) ; + TEST_MSG ("Reachable node count error too large: %.2f%%", + max_err * 100) ; + + OK (LAGraph_Delete (&G, msg)) ; + OK (GrB_free (&scores_approx)) ; + OK (GrB_free (&reachable_approx)) ; + OK (GrB_free (&node_weights)) ; + } + + LAGraph_Finalize (msg) ; +} +#endif + + //------------------------------------------------------------------------------ // test_errors: test error handling //------------------------------------------------------------------------------ @@ -258,6 +337,7 @@ TEST_LIST = { #if LG_SUITESPARSE_GRAPHBLAS_V10 {"HarmonicCentrality", test_HarmonicCentrality}, {"HarmonicCentrality_empty", test_HarmonicCentrality_empty}, + {"HarmonicCentrality_reachable", test_HarmonicCentrality_reachable}, #endif {"HarmonicCentrality_errors", test_errors}, {NULL, NULL} From 837ac7ee46b1a4fc8c4ac77d776d3ee6b905570b Mon Sep 17 00:00:00 2001 From: GomezGab Date: Tue, 7 Apr 2026 11:05:48 -0500 Subject: [PATCH 4/5] correct assertions and types --- experimental/algorithm/LAGr_HarmonicCentrality.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/experimental/algorithm/LAGr_HarmonicCentrality.c b/experimental/algorithm/LAGr_HarmonicCentrality.c index c8d661ad97..e987a35860 100644 --- a/experimental/algorithm/LAGr_HarmonicCentrality.c +++ b/experimental/algorithm/LAGr_HarmonicCentrality.c @@ -318,11 +318,11 @@ int LAGr_HarmonicCentrality( node_weights, NULL)) ; GRB_TRY(GrB_Vector_reduce_INT64(&min_w, NULL, GrB_MIN_MONOID_INT64, node_weights, NULL)) ; - LG_ASSERT_MSG(min_w > 0, GrB_INVALID_VALUE, + LG_ASSERT_MSG(min_w >= 0, GrB_INVALID_VALUE, "Negative node weights not supported"); // TODO: is this cap reasonable? - LG_ASSERT_MSG(max_w < 1000000, GrB_NOT_IMPLEMENTED, + LG_ASSERT_MSG(max_w <= 1000000, GrB_NOT_IMPLEMENTED, "Node weights over 1000000 not supported"); } @@ -459,7 +459,7 @@ int LAGr_HarmonicCentrality( if (reachable_nodes) { // make reachable_nodes vector with same sparsity pattern as scores - GRB_TRY (GrB_Vector_new (reachable_nodes, GrB_FP64, nrows)) ; + GRB_TRY (GrB_Vector_new (reachable_nodes, GrB_UINT64, nrows)) ; GRB_TRY (GrB_apply (delta_vec, NULL, NULL, count_hll, new_sets, NULL)) ; GrB_Vector I_vec = (score_cont->format == GxB_FULL) ? NULL : score_cont->i; From c3d6fbc48fc2e68d682298b4151e2118764b57f5 Mon Sep 17 00:00:00 2001 From: Gabriel Gomez Date: Thu, 9 Apr 2026 12:14:31 -0500 Subject: [PATCH 5/5] Comment variables for HLL --- .../algorithm/LAGr_HarmonicCentrality.c | 89 ++++++++++++++----- 1 file changed, 67 insertions(+), 22 deletions(-) diff --git a/experimental/algorithm/LAGr_HarmonicCentrality.c b/experimental/algorithm/LAGr_HarmonicCentrality.c index e987a35860..c7a759a715 100644 --- a/experimental/algorithm/LAGr_HarmonicCentrality.c +++ b/experimental/algorithm/LAGr_HarmonicCentrality.c @@ -52,7 +52,7 @@ LG_JIT_STRING( typedef struct { - uint8_t registers[(1 << 10)]; + uint8_t registers[HLL_REGISTERS]; } HLL; , HLL_jit) @@ -156,7 +156,7 @@ LG_JIT_STRING( if (diff) { uint32_t i; - double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(1 << 10)) ; + double alpha_mm = 0.7213 / (1.0 + 1.079 / (double)(HLL_REGISTERS)) ; alpha_mm *= ((double)(HLL_REGISTERS) * (double)(HLL_REGISTERS)) ; @@ -241,7 +241,7 @@ GrB_free(&flat_weight); \ GrB_free(&delta_vec); \ GrB_free(&desc); \ GrB_free(&init_hlls); \ -GrB_free(&shallow_second); \ +GrB_free(&hll_second); \ GrB_free(&merge_hll_biop); \ GrB_free(&merge_hll); \ GrB_free(&merge_second); \ @@ -271,26 +271,51 @@ int LAGr_HarmonicCentrality( const GrB_Vector node_weights, // participating nodes and their weights char *msg ) { + // compact G->A: without the nodes with no specified weights GrB_Matrix _A = NULL; - GxB_Container score_cont = NULL; + GrB_Index nrows = 0; GrB_Index nvals = 0; + // HC score of the nodes + GxB_Container score_cont = NULL; + GrB_Vector flat_scores = NULL; + + // HyperLogLog GraphBLAS type GrB_Type hll_t = NULL; + + // Vectors holding the hll sets. After step i, holds the ball of radius i + // around each node GrB_Vector new_sets = NULL; GrB_Vector old_sets = NULL; GxB_Container old_cont = NULL; - GrB_Vector flat_scores = NULL; + + // compact wieghts GrB_Vector flat_weight = NULL; + + // |new_sets| - |old_sets| GrB_Vector delta_vec = NULL; GrB_Descriptor desc = NULL; + + // Operators + // Add the given node with (weight) hashes into its HLL set GrB_IndexUnaryOp init_hlls = NULL; + + // Query set cardinality GrB_UnaryOp count_hll = NULL; - GrB_BinaryOp shallow_second = NULL; + + // z = y + GrB_BinaryOp hll_second = NULL; + + // z = x union y GrB_BinaryOp merge_hll_biop = NULL; GrB_Monoid merge_hll = NULL; + + // [merge_hll.hll_second] GrB_Semiring merge_second = NULL; + + // z = |y| - |x| GrB_BinaryOp delta_hll = NULL; LG_ASSERT(G != NULL, GrB_NULL_POINTER); @@ -361,33 +386,26 @@ int LAGr_HarmonicCentrality( GRB_TRY (GxB_Vector_extractTuples_Vector ( NULL, flat_weight, node_weights, NULL)) ; - // count op GRB_TRY (GxB_UnaryOp_new (&count_hll, (GxB_unary_function) lg_hll_count, GrB_FP64, hll_t, "lg_hll_count", LG_HLL_COUNT)) ; - // init op: weight (INT64) at row index i → HLL seeded with 'weight' hashes - GRB_TRY(GrB_IndexUnaryOp_new( + GRB_TRY (GrB_IndexUnaryOp_new ( &init_hlls, (GxB_index_unary_function)lg_hll_init, hll_t, GrB_INT64, GrB_BOOL)) ; - // merge binary op (HLL, HLL) -> HLL in-place: z == x required GRB_TRY(GxB_BinaryOp_new( &merge_hll_biop, (GxB_binary_function)lg_hll_merge, hll_t, hll_t, hll_t, "lg_hll_merge", LG_HLL_MERGE_STR)) ; - // second op GRB_TRY(GxB_BinaryOp_new( - &shallow_second, (GxB_binary_function)lg_hll_second, + &hll_second, (GxB_binary_function)lg_hll_second, hll_t, GrB_BOOL, hll_t, "lg_hll_second", LG_HLL_SECOND_STR)) ; - // merge monoid - identity is an empty (all-zero) HLL sketch HLL hll_zero = {0}; GRB_TRY(GrB_Monoid_new_UDT(&merge_hll, merge_hll_biop, &hll_zero)) ; - // semiring: add = merge monoid, multiply = copy (pass-through second operand) - GRB_TRY(GrB_Semiring_new(&merge_second, merge_hll, shallow_second)) ; + GRB_TRY(GrB_Semiring_new(&merge_second, merge_hll, hll_second)) ; - // delta op: (HLL_old, HLL_new) → FP64 cardinality change GRB_TRY(GxB_BinaryOp_new( &delta_hll, (GxB_binary_function)lg_hll_delta, GrB_FP64, hll_t, hll_t, "lg_hll_delta", LG_HLL_DELTA_STR)) ; @@ -503,16 +521,19 @@ int LAGr_HarmonicCentrality( GrB_free(&score); \ GrB_free(&desc); \ GrB_free(&node_compact); \ + GrB_free(&score_compact); \ + GrB_free(&reach_compact); \ LAGraph_Delete(&G_compact, NULL); #undef LG_FREE_ALL #define LG_FREE_ALL \ LG_FREE_WORK ; \ - GrB_free(scores) ; + GrB_free(scores) ; \ + GrB_free(reachable_nodes) ; // This is an exact (but much slower) calculation of harmonic centrality #if LG_SUITESPARSE_GRAPHBLAS_V10 -int LAGr_HarmonicCentrality_exact( +int LAGr_HarmonicCentrality_exact ( // outputs: GrB_Vector *scores, // FP64 scores by original node ID GrB_Vector *reachable_nodes, // [optional] estimate the number of reach- @@ -523,21 +544,29 @@ int LAGr_HarmonicCentrality_exact( const GrB_Vector node_weights, // participating nodes and their weights char *msg ) { + // compact G->A GrB_Matrix _A = NULL; + + // Dimensions of node_weights vector GrB_Index nrows = 0; GrB_Index nvals = 0; + // BFS result GrB_Vector level = NULL; + + // compacted vectors GrB_Vector flat_weight = NULL; GrB_Vector node_compact = NULL; + GrB_Vector score_compact = NULL; + GrB_Vector reach_compact = NULL; + GxB_Iterator it = NULL; GrB_Scalar score = NULL; + GrB_Index n_reachable = 0; GrB_Descriptor desc = NULL; LAGraph_Graph G_compact = NULL; - // TODO: - LG_ASSERT(reachable_nodes == NULL, GrB_NOT_IMPLEMENTED); LG_ASSERT(G != NULL, GrB_NULL_POINTER); LG_ASSERT(G->A != NULL, GrB_NULL_POINTER); LG_ASSERT(scores != NULL, GrB_NULL_POINTER); @@ -546,6 +575,10 @@ int LAGr_HarmonicCentrality_exact( GRB_TRY(GrB_Vector_size(&nrows, node_weights)) ; GRB_TRY(GrB_Vector_nvals(&nvals, node_weights)) ; GRB_TRY(GrB_Vector_new(scores, GrB_FP64, nrows)) ; + if (reachable_nodes != NULL) + { + GRB_TRY(GrB_Vector_new(reachable_nodes, GrB_INT64, nrows)) ; + } if (nvals == 0) { return GrB_SUCCESS; @@ -587,6 +620,8 @@ int LAGr_HarmonicCentrality_exact( GRB_TRY (GxB_Scalar_new (&score, GrB_FP64)) ; GRB_TRY (GxB_Iterator_new(&it)) ; GRB_TRY (GrB_Vector_new (&node_compact, GrB_INT64, nvals)); + GRB_TRY (GrB_Vector_new (&score_compact, GrB_FP64, nvals)); + GRB_TRY (GrB_Vector_new (&reach_compact, GrB_INT64, nvals)); GRB_TRY (GxB_Vector_extract_Vector( node_compact, NULL, NULL, nodes, node_weights, desc)) ; GRB_TRY(GrB_free(&desc)) ; @@ -602,15 +637,25 @@ int LAGr_HarmonicCentrality_exact( GRB_TRY (GrB_Vector_removeElement (level, nodeIdx)) ; // dot product (plus div) - GRB_TRY(GrB_vxm((GrB_Vector) score, NULL, NULL, GxB_PLUS_DIV_FP64, + GRB_TRY (GrB_vxm ((GrB_Vector) score, NULL, NULL, GxB_PLUS_DIV_FP64, flat_weight, (GrB_Matrix) level, NULL)) ; + GRB_TRY (GrB_Vector_nvals (&n_reachable, level)); GrB_free (&level); - GRB_TRY (GrB_Vector_setElement_Scalar(*scores, score, nodeIdx)) ; + GRB_TRY (GrB_Vector_setElement_INT64 (reach_compact, n_reachable, nodeIdx)); + GRB_TRY (GrB_Vector_setElement_Scalar (score_compact, score, nodeIdx)) ; info = GxB_Vector_Iterator_next(it); } + GRB_TRY (GxB_Vector_assign_Vector ( + *scores, NULL, NULL, score_compact, node_weights, desc)) ; + if (reachable_nodes != NULL) + { + GRB_TRY (GxB_Vector_assign_Vector ( + *reachable_nodes, NULL, NULL, reach_compact, node_weights, desc)) ; + } + //-------------------------------------------------------------------------- // cleanup //--------------------------------------------------------------------------