diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index ff5b7ff06e..c3523ac8d0 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -19,7 +19,7 @@ jobs: # if there are multiple items in this list, only use should # deployit=true for just one of them. - {grb_version: 9.0.0, deployit: false} - - {grb_version: 10.1.0, deployit: true} + - {grb_version: 10.2.0, deployit: true} steps: @@ -69,7 +69,7 @@ jobs: matrix: config: - {grb_version: 9.0.0} - - {grb_version: 10.1.0} + - {grb_version: 10.2.0} steps: @@ -99,5 +99,5 @@ jobs: cd build cmake .. JOBS=2 make - ctest . + ctest . --rerun-failed --output-on-failure diff --git a/data/50node.mtx b/data/50node.mtx new file mode 100644 index 0000000000..b4a6e2afed --- /dev/null +++ b/data/50node.mtx @@ -0,0 +1,336 @@ +%%MatrixMarket matrix coordinate pattern symmetric +%%GraphBLAS type double +50 50 283 +1 2 +1 3 +1 4 +1 5 +1 6 +1 20 +1 19 +1 22 +1 29 +1 41 +1 11 +2 7 +2 30 +2 37 +2 28 +2 13 +2 38 +2 43 +2 16 +2 33 +2 5 +2 18 +2 24 +2 12 +2 42 +2 35 +2 39 +3 37 +3 28 +3 19 +3 42 +3 29 +3 20 +3 12 +3 13 +3 47 +3 17 +3 6 +3 4 +3 46 +3 14 +4 17 +4 40 +4 5 +4 19 +4 6 +4 22 +4 41 +4 10 +4 27 +4 49 +4 29 +4 32 +5 20 +5 29 +5 11 +5 32 +5 15 +5 19 +6 17 +6 8 +6 19 +6 22 +6 50 +6 13 +6 32 +6 40 +6 41 +6 29 +7 8 +7 9 +7 10 +7 11 +7 12 +7 37 +7 33 +7 15 +7 38 +8 9 +8 38 +8 14 +8 42 +8 30 +8 13 +8 15 +8 28 +8 44 +8 21 +8 24 +8 43 +9 13 +9 14 +9 15 +9 16 +9 31 +9 33 +9 38 +9 43 +9 42 +9 11 +10 28 +10 33 +10 27 +10 23 +10 37 +10 39 +10 42 +10 14 +11 17 +11 20 +11 40 +11 21 +11 32 +12 16 +12 31 +12 37 +12 13 +12 33 +12 24 +12 27 +12 30 +12 36 +12 25 +12 43 +13 30 +13 27 +13 18 +13 39 +13 42 +13 40 +13 14 +13 22 +13 15 +13 24 +14 16 +14 44 +14 22 +14 18 +14 25 +14 45 +14 50 +14 34 +14 32 +14 36 +14 49 +14 43 +14 23 +14 28 +15 30 +15 23 +15 42 +16 27 +16 28 +16 44 +16 24 +16 42 +17 18 +17 19 +17 20 +17 21 +17 22 +17 23 +17 33 +18 28 +18 49 +18 44 +18 47 +18 34 +18 48 +19 46 +19 34 +19 21 +19 40 +19 32 +19 41 +19 29 +19 47 +19 38 +20 21 +20 29 +20 22 +21 40 +21 26 +21 27 +21 41 +22 48 +22 34 +22 32 +22 49 +22 39 +23 24 +23 25 +23 26 +23 44 +23 47 +23 35 +24 37 +24 33 +24 39 +24 47 +24 42 +24 27 +24 28 +24 35 +24 48 +25 33 +25 44 +25 47 +25 26 +25 50 +25 34 +25 36 +25 46 +26 45 +26 47 +26 36 +26 50 +26 41 +26 35 +27 33 +27 43 +27 30 +27 28 +27 39 +28 42 +28 38 +29 30 +29 32 +29 45 +30 31 +30 32 +30 33 +30 43 +30 42 +30 39 +31 34 +31 35 +31 36 +31 47 +31 50 +33 37 +33 43 +33 38 +33 39 +33 42 +33 40 +34 44 +34 50 +34 47 +34 48 +34 35 +35 47 +35 49 +35 46 +35 36 +35 45 +35 48 +35 39 +36 48 +36 47 +36 49 +37 38 +37 39 +37 42 +38 39 +38 42 +39 43 +40 41 +42 44 +42 50 +42 43 +44 45 +44 46 +44 50 +44 48 +45 46 +45 47 +46 47 +46 48 +47 49 +47 48 +47 50 +48 50 +1 1 +2 2 +3 3 +4 4 +5 5 +6 6 +7 7 +8 8 +9 9 +10 10 +11 11 +12 12 +13 13 +14 14 +15 15 +16 16 +17 17 +18 18 +19 19 +20 20 +21 21 +22 22 +23 23 +24 24 +25 25 +26 26 +27 27 +28 28 +29 29 +30 30 +31 31 +32 32 +33 33 +34 34 +35 35 +36 36 +37 37 +38 38 +39 39 +40 40 +41 41 +42 42 +43 43 +44 44 +45 45 +46 46 +47 47 +48 48 +49 49 +50 50 \ No newline at end of file diff --git a/data/comm0.mtx b/data/comm0.mtx new file mode 100644 index 0000000000..6640e67974 --- /dev/null +++ b/data/comm0.mtx @@ -0,0 +1,10 @@ +%%MatrixMarket matrix coordinate pattern symmetric +%%GraphBLAS type double_t +6 6 7 +1 2 +1 3 +2 3 +3 4 +4 5 +4 6 +5 6 \ No newline at end of file diff --git a/data/comm0_S.mtx b/data/comm0_S.mtx new file mode 100644 index 0000000000..c40906d850 --- /dev/null +++ b/data/comm0_S.mtx @@ -0,0 +1,9 @@ +%%MatrixMarket matrix coordinate pattern general +%%GraphBLAS type double_t +6 6 6 +1 1 +2 2 +3 3 +4 4 +5 5 +6 6 \ No newline at end of file diff --git a/data/comm0_Sa.mtx b/data/comm0_Sa.mtx new file mode 100644 index 0000000000..2d9a489192 --- /dev/null +++ b/data/comm0_Sa.mtx @@ -0,0 +1,9 @@ +%%MatrixMarket matrix coordinate pattern general +%%GraphBLAS type double_t +6 6 6 +1 3 +2 3 +3 3 +4 6 +5 6 +6 6 \ No newline at end of file diff --git a/experimental/algorithm/LAGr_ModularityMatrix.c b/experimental/algorithm/LAGr_ModularityMatrix.c new file mode 100644 index 0000000000..0482959717 --- /dev/null +++ b/experimental/algorithm/LAGr_ModularityMatrix.c @@ -0,0 +1,129 @@ +//------------------------------------------------------------------------------ +// LAGr_Modularity2.c: Calculates the modularity of a graph +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2024 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 Olumayowa Olowomeye, Texas A&M University + +//------------------------------------------------------------------------------ +// Q = (1/2m)*Trace(S^T*B*S) +// B the modularity matrix = A-(kk^t/2m) + +// Given a symmetric graph A with no-self edges, LAGr_Modularity calculates the +// Modularity of that graph and a given community matrix +// + +// Newman ME. Modularity and community structure in networks. +// Proc Natl Acad Sci U S A. 2006 Jun 6;103(23):8577-82. doi: 10.1073/pnas.0601602103. +// pub 2006 May 24. PMID: 16723398; PMCID: PMC1482622. + +// Current Test File: experimental/test/test_modularity.c + +#include "LG_internal.h" +#include +#include + +#define LG_FREE_MOD \ + { \ + GrB_free(&k); \ + GrB_free(&kk_); \ + GrB_free(&B); \ + GrB_free(&BS); \ + GrB_free(&S_BS); \ + GrB_free(&Diag); \ + } +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + LG_FREE_MOD; \ + } +#define DEBUG 0 +#if DEBUG +#define check() printf("here") +#define dbg(x) \ + if (DEBUG) \ + GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +#else +#define check() +#define dbg(x) +#define err(x, info) +#endif + +int LAGr_AdjModularity( + double *Q, + double gamma, + GrB_Matrix A, + GrB_Matrix S, + char *msg) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LG_CLEAR_MSG; + char MATRIX_TYPE[LAGRAPH_MSG_LEN]; + + // GrB_set(GrB_GLOBAL, false, GxB_BURBLE); + + GrB_Index n; + GrB_Matrix k = NULL; + GrB_Matrix kk_ = NULL; + GrB_Matrix B = NULL; + GrB_Matrix BS = NULL; + GrB_Matrix S_BS = NULL; + GrB_Matrix Diag = NULL; + GrB_Matrix mask = NULL; + + GRB_TRY(GrB_Matrix_nrows(&n, A)); + + GRB_TRY(GrB_Matrix_new(&B, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&k, GrB_FP64, n, 1)); + GRB_TRY(GrB_Matrix_new(&kk_, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&BS, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&S_BS, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&Diag, GrB_FP64, n, n)); + + double m = 0.0; + double Q_ = 0.0; + + GRB_TRY(GrB_Matrix_reduce_Monoid((GrB_Vector)k, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + + GRB_TRY(GrB_Matrix_reduce_FP64(&m, GrB_PLUS_FP64, GrB_PLUS_MONOID_FP64, A, NULL)); + m /= 2.0; + if (m == 0.0) + { + *Q = 0.0; + LG_FREE_ALL; + return GrB_SUCCESS; + } + GRB_TRY(GrB_mxm(kk_, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, k, k, GrB_DESC_T1)); + double inv_m = -gamma / (2.0 * m); + GRB_TRY(GrB_Matrix_apply_BinaryOp2nd_FP64(kk_, NULL, NULL, GrB_TIMES_FP64, kk_, inv_m, GrB_DESC_R)); + GRB_TRY(GrB_eWiseAdd(B, NULL, NULL, GrB_PLUS_FP64, A, kk_, NULL)); + GRB_TRY(GrB_mxm(BS, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, B, S, NULL)); + GRB_TRY(GrB_mxm(S_BS, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, S, BS, GrB_DESC_T0)); + + GRB_TRY(GrB_select(Diag, NULL, NULL, GrB_DIAG, S_BS, 0, NULL)); + + GRB_TRY(GrB_Matrix_reduce_FP64(&Q_, NULL, GrB_PLUS_MONOID_FP64, Diag, NULL)); + Q_ *= -inv_m; + *Q = Q_; + + LG_FREE_ALL; + return (GrB_SUCCESS); +#else + return (GrB_NOT_IMPLEMENTED); +#endif +} diff --git a/experimental/algorithm/LAGraph_IsolateSets.c b/experimental/algorithm/LAGraph_IsolateSets.c new file mode 100644 index 0000000000..27f858ac48 --- /dev/null +++ b/experimental/algorithm/LAGraph_IsolateSets.c @@ -0,0 +1,230 @@ +//------------------------------------------------------------------------------ +// LAGraph_IsolateSets.c: Returns Isolate Sets or a covering 2-hop independent set of a +// of the giving adjacency matrix +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2024 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 Olumayowa Olowomeye, Texas A&M University +// Test File: experimental/test/test_IsolateSets +//------------------------------------------------------------------------------ + +#include "LG_internal.h" +#include +#include +#include +#include +#define DEBUG 0 +#if DEBUG +#define check() printf("here") +#define dbg(x) \ + if (DEBUG) \ + GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +#else +#define check() +#define dbg(x) +#define err(x, info) +#endif +#undef LG_FREE_IS +#define LG_FREE_IS \ + { \ + GrB_free(&score); \ + GrB_free(&scoreA); \ + GrB_free(&neighbor_max); \ + GrB_free(&new_members); \ + GrB_free(&new_membersA); \ + GrB_free(&new_neighbors); \ + GrB_free(&candidates); \ + GrB_free(&empty); \ + GrB_free(&Seed); \ + GrB_free(°ree); \ + GrB_free(&iset); \ + } + +int LAGraph_IsolateSet( + // output + GrB_Vector *isolate_set, + // input + GrB_Matrix A, + GrB_Vector ignore_node, + uint64_t seed, + char *msg) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LG_CLEAR_MSG; + + char MATRIX_TYPE[LAGRAPH_MSG_LEN]; + // GrB_set(GrB_GLOBAL, DEBUG, GxB_BURBLE); + GrB_Vector iset = NULL; // independent set (output vector) + GrB_Vector score = NULL; + GrB_Vector scoreA = NULL; // random score for each node + // GrB_Vector scoreAA = NULL; + GrB_Vector neighbor_max = NULL; // value of max neighbor score + GrB_Vector new_members = NULL; // set of new members to add to iset + GrB_Vector new_membersA = NULL; + GrB_Vector new_neighbors = NULL; // new neighbors to new iset members + GrB_Vector candidates = NULL; // candidate nodes + GrB_Vector empty = NULL; // an empty vector + GrB_Vector Seed = NULL; // random number seed vector + GrB_Vector degree = NULL; // (float) G->out_degree + // GrB_Matrix A ; // G->A, the adjacency matrix + GrB_Index n; // # of nodes + // printf("in Isolate set algorithm"); + // LG_TRY (LAGraph_CheckGraph (G, msg)) ; + LG_ASSERT(isolate_set != NULL, GrB_NULL_POINTER); + // A = G->A; + dbg(A); + + GRB_TRY(GrB_Matrix_nrows(&n, A)); + GRB_TRY(GrB_Vector_new(&iset, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&neighbor_max, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(°ree, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&new_members, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&new_neighbors, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&new_membersA, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&candidates, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&empty, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&Seed, GrB_UINT64, n)); + GRB_TRY(GrB_Vector_new(&score, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&scoreA, GrB_FP64, n)); + + // rand + // seed = 6247; + // printf("%ld",seed); + GRB_TRY(GrB_Matrix_reduce_Monoid(degree, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + dbg(degree); + + GrB_Index ncandidates; + GRB_TRY(GrB_assign(candidates, ignore_node, NULL, (bool)true, GrB_ALL, n, GrB_DESC_C)); + GRB_TRY(GrB_assign(Seed, candidates, NULL, 1, GrB_ALL, n, GrB_DESC_S)); + dbg(candidates); + GRB_TRY(LAGraph_Random_Seed(Seed, seed, msg)); + dbg(Seed); + + GRB_TRY(GrB_Vector_nvals(&ncandidates, candidates)); + + GRB_TRY(GrB_assign(score, NULL, NULL, Seed, GrB_ALL, n, NULL)); + GRB_TRY(GrB_eWiseMult(score, NULL, NULL, GrB_DIV_FP64, score, degree, NULL)); + dbg(score); + + dbg(candidates); + GRB_TRY(GrB_vxm(scoreA, candidates, NULL, GrB_MAX_FIRST_SEMIRING_FP64, score, A, GrB_DESC_RS)); + dbg(scoreA); + GRB_TRY(GrB_vxm(neighbor_max, candidates, NULL, GrB_MAX_FIRST_SEMIRING_FP64, scoreA, A, GrB_DESC_RS)); + GRB_TRY(GrB_vxm(neighbor_max, candidates, NULL, GrB_MAX_FIRST_SEMIRING_FP64, neighbor_max, A, GrB_DESC_RS)); + dbg(neighbor_max); + dbg(score); + + GRB_TRY(GrB_eWiseAdd(new_members, NULL, NULL, GrB_GE_FP64, score, neighbor_max, NULL)); + dbg(new_members); + GRB_TRY(GrB_select(new_members, NULL, NULL, GrB_VALUEEQ_BOOL, + new_members, (bool)true, NULL)); + dbg(new_members); + GRB_TRY(GrB_assign(iset, new_members, NULL, true, GrB_ALL, n, NULL)); + (*isolate_set) = iset; + + // printf("done iset"); + iset = NULL; + dbg(*isolate_set); + LG_FREE_IS; + return (GrB_SUCCESS); +#else + return (GrB_NOT_IMPLEMENTED); +#endif +} + +// #undef LG_FREE_WORK +#undef LG_FREE_ISS +#define LG_FREE_ISS \ + { \ + GrB_free(&ignore_nodes); \ + } + +int LAGraph_IsolateSets( + GrB_Matrix *IsolateSets, // Output: k x n Boolean matrix + // LAGraph_Graph G, // Input: graph + GrB_Matrix A, + // GrB_Vector ignore_nodes, + uint64_t seed, // Input: RNG seed + char *msg // Error message buffer +) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LG_CLEAR_MSG; + // LG_TRY(LAGraph_CheckGraph(G, msg)); + LG_ASSERT(IsolateSets != NULL, GrB_NULL_POINTER); + LG_ASSERT(A != NULL, GrB_NULL_POINTER); + GrB_Index nvals; + GrB_Matrix_nvals(&nvals, A); + if (nvals == 0) + { + *IsolateSets = NULL; + return GrB_SUCCESS; + } + // GrB_Matrix A = G->A; + // dbg(A); + GrB_Index n; + // start NULL -- LAGraph_IsolateSet will allocate + GrB_Index k = 0; + GrB_Index vals_res = 0; + GrB_Vector ignore_nodes = NULL; + GrB_Matrix result = NULL; + GrB_Vector iset = NULL; + + GRB_TRY(GrB_Matrix_nrows(&n, A)); + GrB_Index max_k = n; // Max possible number of isolate sets is <= n + + GRB_TRY(GrB_Vector_new(&ignore_nodes, GrB_BOOL, n)); + GRB_TRY(GrB_Matrix_new(&result, GrB_BOOL, max_k, n)); + + GRB_TRY(GrB_assign(ignore_nodes, NULL, NULL, (bool)false, GrB_ALL, n, NULL)); + + while (true) + { + // LAGraph_IsolateSet allocates a new vector and stores it into 'iset' + GRB_TRY(LAGraph_IsolateSet(&iset, A, ignore_nodes, seed, msg)); + + // if 'iset' is NULL or empty, break + GRB_TRY(GrB_Vector_nvals(&vals_res, iset)); + if (vals_res == 0) + { + GrB_free(&iset); // free the empty iset returned + break; + } + + // mark ignored nodes (update ignore_nodes) BEFORE copying into matrix if desired + GRB_TRY(GrB_Vector_eWiseAdd_BinaryOp(ignore_nodes, NULL, NULL, GrB_LOR, ignore_nodes, iset, NULL)); + + // copy the iset vector into the result matrix row k + GRB_TRY(GxB_Row_assign_Vector(result, NULL, NULL, iset, k, NULL, NULL)); + + // free the iset after copying into the matrix to avoid leaking + GrB_free(&iset); + iset = NULL; + + k++; + } + + GRB_TRY(GrB_Matrix_resize(result, k, n)); + *IsolateSets = result; + + LG_FREE_ISS; + return (GrB_SUCCESS); +#else + return (GrB_NOT_IMPLEMENTED); +#endif +} diff --git a/experimental/algorithm/LAGraph_LouvainIS.c b/experimental/algorithm/LAGraph_LouvainIS.c new file mode 100644 index 0000000000..f5e0084eba --- /dev/null +++ b/experimental/algorithm/LAGraph_LouvainIS.c @@ -0,0 +1,544 @@ +//------------------------------------------------------------------------------ +// LAGraph_LouvainIS.c: Runs the first phase of the Louvain Algorithm with +// Isolate Sets on a given Graph +// +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2024 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 Olumayowa Olowomeye, Texas A&M University +//------------------------------------------------------------------------------ + +// Current Test File: experimental/test/test_louvain.c +#include "LG_internal.h" +#include +#include +#include +#include +#include + +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + +// #define TIMING +#define DEBUG 0 + +#if DEBUG +#define dbg(x) \ + if (DEBUG) GxB_print(x, 5) \ +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +#else +#define dbg(x) +#define err(x, info) +#endif + +typedef struct // Theta +{ + double *d; + double m; + uint64_t seed; +} Theta; +#define THETA_DEFN \ + "typedef struct /* Theta */" \ + "{" \ + " double *d;" \ + " double m;" \ + " uint64_t seed;" \ + "} Theta;" + +typedef struct // argmax_tup +{ + double score; /* change in modularity */ + int64_t comm; /* who */ + double tb; +} argmax_tup; + +#define AM_TUP \ + "typedef struct /* argmax_tup */ \n" \ + "{\n" \ + " double score; /* change in modularity */\n" \ + " int64_t comm; /* who */\n" \ + " double tb;\n" \ + "} argmax_tup;\n" + +void make_argmax_tup(argmax_tup *z, + const double *x, GrB_Index ix, GrB_Index jx, + const double *y, GrB_Index iy, GrB_Index jy, + const void *theta) +{ + Theta *_theta = (Theta *)theta; + double ki_to_C = *x; // Edge weight from node i to community C + double ki = _theta->d[ix]; // Degree of node i + double kC = *y; // Sum of degrees in community C + double m = _theta->m; // Total edge weight of graph (already divided by 2) + + // Compute modularity gain (delta Q) + z->score = (ki_to_C / m) - ((ki * kC) / (m * m)); + + // Generate tiebreaker (seed) + uint64_t seed = _theta->seed + (*y + ix + iy + jy); + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + + z->comm = (int64_t)jx; + z->tb = seed; +} + +#define MAKE_AM_TUP \ + "void make_argmax_tup(argmax_tup *z,\n" \ + " const double *x, GrB_Index ix, GrB_Index jx,\n" \ + " const double *y, GrB_Index iy, GrB_Index jy,\n" \ + " const void *theta)\n" \ + "{\n" \ + " Theta *_theta = (Theta *)theta;\n" \ + " double ki_to_C = *x;\n" \ + " double ki = _theta->d[ix];\n" \ + " double kC = *y;\n" \ + " double m = _theta->m;\n" \ + " z->score = (ki_to_C / m) - ((ki * kC) / (m * m));\n" \ + " uint64_t seed = _theta->seed + (*y + ix + iy + jy);\n" \ + " seed ^= seed << 13;\n" \ + " seed ^= seed >> 7;\n" \ + " seed ^= seed << 17;\n" \ + " z->comm = (int64_t)jx;\n" \ + " z->tb = seed;\n" \ + "}\n" + +void argmax_op(argmax_tup *z, argmax_tup *x, argmax_tup *y) +{ + if (x->score > y->score) + { + z->score = x->score; + z->comm = x->comm; + z->tb = x->tb; + } + else if (x->score == y->score) + { + if (x->tb > y->tb) + { + z->score = x->score; + z->comm = x->comm; + z->tb = x->tb; + } + else + { + z->score = y->score; + z->comm = y->comm; + z->tb = y->tb; + } + } + else + { + z->score = y->score; + z->comm = y->comm; + z->tb = y->tb; + } +} +#define AM_OP \ + "void argmax_op(argmax_tup *z, argmax_tup *x, argmax_tup *y)\n" \ + "{\n" \ + " if (x->score > y->score)\n" \ + " {\n" \ + " z->score = x->score;\n" \ + " z->comm = x->comm;\n" \ + " z->tb = x->tb;\n" \ + " }\n" \ + " else if (x->score == y->score)\n" \ + " {\n" \ + " if (x->tb > y->tb)\n" \ + " {\n" \ + " z->score = x->score;\n" \ + " z->comm = x->comm;\n" \ + " z->tb = x->tb;\n" \ + " }\n" \ + " else\n" \ + " {\n" \ + " z->score = y->score;\n" \ + " z->comm = y->comm;\n" \ + " z->tb = y->tb;\n" \ + " }\n" \ + " }\n" \ + " else\n" \ + " {\n" \ + " z->score = y->score;\n" \ + " z->comm = y->comm;\n" \ + " z->tb = y->tb;\n" \ + " }\n" \ + "}\n" +static GrB_Info build_argmax_operator( + GrB_Type Theta_UDT, GrB_Type Tuple, + double *d, double m, uint64_t seed, + GrB_BinaryOp *MAKEAMTUP_Bop, + GrB_BinaryOp *AM_Bop, + GrB_Monoid *AM_mon, + GrB_Semiring *AM_Semiring, + GxB_IndexBinaryOp *MAKEAMTUP_op, + GrB_Scalar *_0, + char *msg) +{ + Theta theta_scalar; + theta_scalar.d = d; + theta_scalar.seed = seed; + theta_scalar.m = m; + + GRB_TRY(GrB_Scalar_setElement_UDT(*_0, (void *)&theta_scalar)); + GRB_TRY(GxB_IndexBinaryOp_new(MAKEAMTUP_op, + (GxB_index_binary_function)make_argmax_tup, + Tuple, GrB_FP64, GrB_FP64, Theta_UDT, + "make_argmax_tup", MAKE_AM_TUP)); + GRB_TRY(GxB_BinaryOp_new_IndexOp(MAKEAMTUP_Bop, *MAKEAMTUP_op, *_0)); + + argmax_tup id; + memset(&id, 0, sizeof(argmax_tup)); + id.tb = 0; + id.comm = INT64_MAX; + id.score = (double)(-INFINITY); + + GRB_TRY(GxB_BinaryOp_new(AM_Bop, + (GxB_binary_function)argmax_op, Tuple, Tuple, Tuple, + "argmax_op", AM_OP)); + + GRB_TRY(GrB_Monoid_new_UDT(AM_mon, *AM_Bop, &id)); + + GRB_TRY(GrB_Semiring_new(AM_Semiring, *AM_mon, *MAKEAMTUP_Bop)); + + return GrB_SUCCESS; +} + +#define LG_FREE_EXTRACT \ + { \ + GrB_free(&extract_k_if_gain_op); \ + } +void extract_k_if_gain(void *out, const void *in) +{ + const argmax_tup *a = in; + int64_t *k_out = out; + if (a->score > 0.0) + { + *k_out = a->comm; + } + else + { + *k_out = a->comm; + } +} +#define EXTRACT_K_IF_GAIN_SRC \ + "void extract_k_if_gain(void *out, const void *in) {\n" \ + " const argmax_tup *a = in;\n" \ + " int64_t *k_out = out;\n" \ + " if (a->score > 0.0) {\n" \ + " *k_out = a->comm;\n" \ + " } else {\n" \ + " *k_out = a->comm;\n" \ + " }\n" \ + "}\n" + +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + /* vectors */ \ + GrB_free(&x); \ + GrB_free(&y); \ + GrB_free(&k); \ + GrB_free(&iset); \ + GrB_free(&Wy); \ + GrB_free(&k_values); \ + GrB_free(&Si_old); \ + GrB_free(&A); \ + GrB_free(&Miset); \ + GrB_free(&W); \ + GrB_free(&A_rows); \ + GxB_Container_free(&S_container); \ + GxB_Container_free(&k_container); \ + GrB_free(&ri); \ + GrB_free(&rv); \ + GrB_free(&AM_Semiring); \ + GrB_free(&AM_mon); \ + GrB_free(&AM_Bop); \ + GrB_free(&MAKEAMTUP_Bop); \ + GrB_free(&MAKEAMTUP_op); \ + GrB_free(&argmax_0); \ + GrB_free(&Theta_UDT); \ + GrB_Type_free(&Tuple); \ + if (A_iset != NULL) \ + { \ + for (int i = 0; i < loop; i++) \ + { \ + GrB_free(&A_iset[i]); \ + } \ + LAGraph_Free ((void **) &A_iset, NULL) ; \ + } \ + } + +#endif + +int LAGraph_LouvainIS( + // output + GrB_Matrix *S_result, + uint64_t seed, + // input + LAGraph_Graph G, + char *msg) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + char MATRIX_TYPE[LAGRAPH_MSG_LEN]; + // GrB_set(GrB_GLOBAL, false, GxB_BURBLE); + + GrB_Descriptor ri = NULL; + GrB_Descriptor rv = NULL; + + GrB_Vector iset = NULL; + GrB_Vector k = NULL; + GrB_Vector x = NULL; + GrB_Vector y = NULL; + GrB_Vector Wy = NULL; + GrB_Vector k_values = NULL; + GrB_Vector Si_old = NULL; + GrB_Vector gain_values = NULL; + GrB_Vector gain_mask = NULL; + + GrB_Matrix S = NULL; + GrB_Matrix A = NULL; + GrB_Matrix W = NULL; + GrB_Matrix Miset = NULL; + GrB_Matrix A_rows = NULL; + + GrB_Index n; + GrB_Index niset; + GrB_Index ncols; + + GxB_Container S_container = NULL; + GxB_Container k_container = NULL; + + GrB_Scalar argmax_0; + GrB_Type Tuple = NULL; + GxB_IndexBinaryOp MAKEAMTUP_op = NULL; + GrB_BinaryOp MAKEAMTUP_Bop = NULL, AM_Bop = NULL; + GrB_Monoid AM_mon = NULL; + GrB_Semiring AM_Semiring = NULL; + GrB_Type Theta_UDT = NULL; + + GrB_Matrix *A_iset = NULL ; + GrB_Index loop = 0 ; + + GRB_TRY(GxB_Type_new(&Theta_UDT, sizeof(Theta), "Theta", THETA_DEFN)); + GRB_TRY(GrB_Scalar_new(&argmax_0, Theta_UDT)); + + GRB_TRY(GxB_Type_new(&Tuple, sizeof(argmax_tup), "argmax_tup", AM_TUP)); + + GrB_UnaryOp extract_k_if_gain_op = NULL; + GRB_TRY(GxB_UnaryOp_new(&extract_k_if_gain_op, extract_k_if_gain, GrB_INT64, Tuple, "extract_k_if_gain", EXTRACT_K_IF_GAIN_SRC)); + + LG_ASSERT(S_result != NULL, GrB_NULL_POINTER); + GrB_Info info; + dbg(A); + GRB_TRY(GrB_Matrix_dup(&A, G->A)); + dbg(Tuple); + + GRB_TRY(GrB_Matrix_nrows(&n, A)); + GRB_TRY(GrB_Matrix_ncols(&ncols, A)); + // printf("n: %lu\n", n); + + GRB_TRY(GrB_Vector_new(&y, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&x, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&k, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&Wy, Tuple, n)); + GRB_TRY(GrB_Vector_new(&iset, GrB_BOOL, n)); + GRB_TRY(GrB_Vector_new(&gain_values, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&k_values, GrB_INT64, n)); + GRB_TRY(GrB_Vector_new(&gain_mask, GrB_BOOL, n)); + + GRB_TRY(GrB_Matrix_new(&W, GrB_FP64, n, n)); + GRB_TRY(GrB_Descriptor_new(&ri)); + GRB_TRY(GrB_set(ri, GxB_USE_INDICES, GxB_ROWINDEX_LIST)); + GRB_TRY(GrB_Descriptor_new(&rv)); + GRB_TRY(GrB_set(rv, GxB_USE_VALUES, GxB_ROWINDEX_LIST)); + + GRB_TRY(GxB_Container_new(&S_container)); + GRB_TRY(GxB_Container_new(&k_container)); + + double m = 0.0; + + void *f = NULL; + uint64_t f_size, f_nvals = 0, f_nheld = 0; + GrB_Type ftype = NULL; + int f_handling; + + void *c = NULL; + uint64_t c_size, c_nvals = 0, c_nheld = 0; + GrB_Type ctype = NULL; + int c_handling; + + // S <- I + + GRB_TRY(GrB_assign(x, NULL, NULL, true, GrB_ALL, n, NULL)); + GRB_TRY(GrB_Matrix_diag(&S, x, 0)); + GRB_TRY(GrB_set(S, GxB_SPARSE, GxB_SPARSITY_CONTROL)); + dbg(S); + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + GRB_TRY(GrB_Vector_dup(&Si_old, S_container->i)); + dbg(Si_old); + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + GRB_TRY(GrB_Matrix_reduce_Monoid(k, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + GRB_TRY(GrB_set(k, GxB_SPARSE, GxB_SPARSITY_CONTROL)); + dbg(k); + GRB_TRY(GrB_Vector_reduce_FP64(&m, NULL, GrB_PLUS_MONOID_FP64, k, NULL)); + m /= 2; + if (m == 0) + { + *S_result = NULL; + LG_FREE_ALL; + return 0; + } + GRB_TRY(GxB_unload_Vector_into_Container(k, k_container, NULL)); + GRB_TRY(GxB_Vector_unload(k_container->x, &f, &ftype, &f_nheld, &f_size, &f_handling, NULL)); + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + GRB_TRY(GxB_Vector_unload(S_container->i, &c, &ctype, &c_nheld, &c_size, &c_handling, NULL)); + + GRB_TRY(build_argmax_operator( + Theta_UDT, Tuple, f, m, seed, + &MAKEAMTUP_Bop, &AM_Bop, &AM_mon, &AM_Semiring, + &MAKEAMTUP_op, &argmax_0, msg)); + GRB_TRY(GxB_Vector_load(k_container->x, &f, ftype, f_nheld, f_size, f_handling, NULL)); + dbg(k_container->x); + GRB_TRY(GxB_load_Vector_from_Container(k, k_container, NULL)); + dbg(k); + GRB_TRY(GxB_Vector_load(S_container->i, &c, ctype, c_nheld, c_size, c_handling, NULL)); + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + dbg(S); + bool changed = true; + + int iter = 0; + double Q = 0; + double gamma = 1; + int max_iter = 20; + + #ifdef TIMING + double tsimple = LAGraph_WallClockTime(); + #endif + GRB_TRY(LAGraph_IsolateSets(&Miset, A, seed, msg)); + #ifdef TIMING + tsimple = LAGraph_WallClockTime() - tsimple; + printf("Isolate Set calc time: %10.10f\n", tsimple); + #endif + + GRB_TRY(GrB_Matrix_nrows(&loop, Miset)); + + // allocate an array of GrB_Matrix handles + // GrB_Matrix *A_iset = malloc(loop * sizeof(GrB_Matrix)); + LG_TRY (LAGraph_Calloc ((void **) &A_iset, loop, sizeof (GrB_Matrix), msg)) ; + + GrB_Index n_k_values; + #ifdef TIMING + double tsimple = LAGraph_WallClockTime(); + #endif + for (int i = 0; i < loop; i++) + { + // extract column into iset + GRB_TRY(GrB_Col_extract(iset, NULL, NULL, Miset, GrB_ALL, n, i, GrB_DESC_T0)); + // dbg(iset); + + GrB_Index niset; + GrB_Vector_nvals(&niset, iset); + + // temp matrix to hold extracted rows + GrB_Matrix A_rows = NULL; + GRB_TRY(GrB_Matrix_new(&A_rows, GrB_FP64, niset, ncols)); + + // fill A_rows from A + GRB_TRY(GxB_Matrix_extract_Vector(A_rows, NULL, NULL, A, iset, NULL, ri)); + // dbg(A_rows); + + // now create the ith matrix in A_iset + GRB_TRY(GrB_Matrix_new(&A_iset[i], GrB_FP64, n, n)); + + // assign A_rows into it + GRB_TRY(GxB_Matrix_assign_Vector(A_iset[i], NULL, NULL, A_rows, iset, NULL, ri)); + + GrB_free(&A_rows); + A_rows = NULL; + } + + #ifdef TIMING + tsimple = LAGraph_WallClockTime() - tsimple; + printf("A_ist array set calc time: %10.10f\n", tsimple); + double tsimple2 = LAGraph_WallClockTime(); + #endif + bool recompute = false; + + while (changed && iter < max_iter) + { + recompute = changed; + changed = false; + for (int i = 0; i < loop; i++) + { + double tsimple3 = LAGraph_WallClockTime(); + + dbg(S); + GRB_TRY(GrB_mxm(W, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, A_iset[i], S, NULL)); + dbg(W); + if(recompute){ + GRB_TRY(GrB_vxm(y, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, k, S, GrB_DESC_T0)); + dbg(y); + } + GRB_TRY(GrB_mxv(Wy, NULL, NULL, AM_Semiring, W, y, NULL)); + dbg(Wy); + + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + + GRB_TRY(GrB_Vector_apply(k_values, NULL, NULL, extract_k_if_gain_op, Wy, NULL)); + GRB_TRY(GrB_Vector_nvals(&n_k_values, k_values)); + if (n_k_values > 0) + { + GRB_TRY(GrB_assign(S_container->i, k_values, NULL, k_values, GrB_ALL, n, GrB_DESC_S)); + } + + dbg(S_container->i); + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + dbg(S); + } + + // break; + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + GRB_TRY(LAGraph_Vector_IsEqual(&changed, Si_old, S_container->i, msg)); + changed = !changed; + if (changed) + { // skip if they are the same + GRB_TRY(GrB_assign(Si_old, NULL, NULL, S_container->i, GrB_ALL, n, GrB_DESC_S)); + } + + + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + dbg(Si_old); + + iter++; + } + #ifdef TIMING + tsimple2 = LAGraph_WallClockTime() - tsimple2; + printf("Main calc time: %10.10f\n", tsimple2); + printf("Iterations: %d\n", iter); + #endif + + (*S_result) = S; + S = NULL; + + LG_FREE_EXTRACT; + LG_FREE_ALL; + return (GrB_SUCCESS) ; +#else + return (GrB_NOT_IMPLEMENTED); +#endif +} diff --git a/experimental/algorithm/LAGraph_LouvainSeq.c b/experimental/algorithm/LAGraph_LouvainSeq.c new file mode 100644 index 0000000000..36d6dcd589 --- /dev/null +++ b/experimental/algorithm/LAGraph_LouvainSeq.c @@ -0,0 +1,379 @@ +//------------------------------------------------------------------------------ +// LAGraph_LouvainSeq.c: Runs the first phase of the Louvain Algorithm on a given graph +//------------------------------------------------------------------------------ + +// LAGraph, (c) 2019-2024 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 Olumayowa Olowomeye, Texas A&M University +// Originally written by T. M. Low, D. G. Spampinato, S. McMillan, and M. Pelletier, “Linear Algebraic Louvain Method in Python,” in +// 2020 IEEE International Parallel and Distributed Processing Symposium Workshops +// (IPDPSW), New Orleans, LA, USA: IEEE, May 2020, pp. 223–226. doi: 10.1109/IPDPSW50202.2020.00050. +// Modified to fit GraphBLAS +//------------------------------------------------------------------------------ + +// Current Test File: experimental/test/test_louvain.c + +#include "LG_internal.h" +#include +#include +#include +#include +#undef LG_FREE_ALL +#define LG_FREE_ALL \ + { \ + GrB_free(&k); \ + GrB_free(&x); \ + GrB_free(&v); \ + GrB_free(&srxq); \ + GrB_free(&sr); \ + GrB_free(&q1); \ + GrB_free(&t); \ + GrB_free(&t_q); \ + GrB_free(&dS); \ + GrB_free(&dSTk); \ + GrB_free(&vtS); \ + GrB_free(&temp); \ + GrB_free(&y_rand); \ + GrB_free(&max_q1); \ + GrB_free(&za); \ + GrB_free(&z_dSTk); \ + GrB_free(&AS); \ + GrB_free(&StAS); \ + GrB_free(&z); \ + GrB_Scalar_free(&Theta); \ + GrB_Type_free(&Tuple); \ + GrB_Semiring_free(&IBop_MAX); \ + GrB_BinaryOp_free(&MAKEFP64_Bop); \ + GxB_IndexBinaryOp_free(&MAKEFP64_op); \ + GrB_BinaryOp_free(&MAXFP64_op); \ + GrB_Monoid_free(&MAXFP64_mon); \ + GxB_Container_free(&S_container); \ + } +#define DEBUG 0 +#if DEBUG +#define check() printf("here") +#define dbg(x) \ + if (DEBUG) \ + GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +#else +#define check() +#define dbg(x) +#define err(x, info) +#endif +typedef struct tuple_fp64 +{ + int64_t k; + double v; + uint64_t tb; +} tuple_fp64; +#define FP64_K "typedef struct tuple_fp64 { int64_t k ; double v ; uint64_t tb ;} tuple_fp64 ;" +void make_fp64(tuple_fp64 *z, + const double *x, GrB_Index ix, GrB_Index jx, + const uint64_t *y, GrB_Index iy, GrB_Index jy, + const void *theta) +{ + z->k = (int64_t)jx; + z->v = (*x); + uint64_t seed = (*y + ix + iy + jy); + seed ^= seed << 13; + seed ^= seed >> 7; + seed ^= seed << 17; + z->tb = seed; +} +void max_fp64(tuple_fp64 *z, const tuple_fp64 *x, const tuple_fp64 *y) +{ + + if (x->v > y->v) + { + z->k = x->k; + z->v = x->v; + } + else if (x->v == y->v) + { + if (x->tb > y->tb) + { + z->k = y->k; + z->v = y->v; + } + else + { + z->k = x->k; + z->v = x->v; + } + } + else + { + z->k = y->k; + z->v = y->v; + } +} +#define MAX_FP64 \ + "void max_fp64(tuple_fp64 *z, const tuple_fp64 *x, const tuple_fp64 *y){ \n" \ + " if (x->v > y->v) \n" \ + " { \n" \ + " z->k = x->k; \n" \ + " z->v = x->v; \n" \ + " }else if(x->v == y->v){ \n" \ + " if(x->tb > y->tb){z->k = y->k;z->v = y->v;} else {z->k = x->k;z->v = x->v;}\n" \ + " }else{ \n" \ + " z->k = y->k; \n" \ + " z->v = y->v; \n" \ + " } \n" \ + "}" +#define MAKE_FP64 \ + "void make_fp64(tuple_fp64 *z, \n" \ + " const double *x, GrB_Index ix, GrB_Index jx, \n" \ + " const uint64_t *y, GrB_Index iy, GrB_Index jy, \n" \ + " const void *theta) \n" \ + "{ \n" \ + " z->k = (int64_t)jx; \n" \ + " z->v = (*x); \n" \ +" uint64_t seed = (*y + ix + iy + jy); \n" \ + " seed ^= seed << 13 ; \n" \ + " seed ^= seed >> 7 ;\n" \ + " seed ^= seed << 17 ;\n" \ + " z->tb = seed; \n" \ + "}" + + +int LAGraph_LouvainSeq( + // output + GrB_Matrix *S_result, + // input + LAGraph_Graph G, + uint64_t seed, + char *msg) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + char MATRIX_TYPE[LAGRAPH_MSG_LEN]; + // GrB_set(GrB_GLOBAL, false, GxB_BURBLE); + + + GrB_Vector t_q = NULL, q1 = NULL, t = NULL, v = NULL; + GrB_Vector k = NULL; + GrB_Vector x = NULL; + GrB_Vector z = NULL; + GrB_Vector za = NULL; + GrB_Vector z_dSTk = NULL; + GrB_Index n, b; + GrB_Matrix S = NULL; + GrB_Matrix dS = NULL; + GrB_Vector sr = NULL; + GrB_Vector srxq = NULL; + GrB_Index vals_srxq; + // GrB_Matrix dS = NULL ; + GrB_Vector dSTk = NULL, vtS = NULL; + GrB_Vector temp = NULL; + GrB_Vector y_rand = NULL; + GrB_Vector max_q1 = NULL; + GrB_Matrix AS = NULL; + GrB_Matrix StAS = NULL; + GrB_Semiring IBop_MAX = NULL; + GxB_IndexBinaryOp MAKEFP64_op = NULL; + GrB_BinaryOp MAKEFP64_Bop = NULL, MAXFP64_op = NULL; + GrB_Scalar Theta = NULL; + GrB_Type Tuple = NULL; + GrB_Monoid MAXFP64_mon = NULL; + + GxB_Container S_container = NULL; + GrB_Index q1_size; + // GrB_Container + tuple_fp64 o; + double o1; + double k_i = 0; + LG_ASSERT(S_result != NULL, GrB_NULL_POINTER); + + GrB_Matrix A = G->A; + // GxB_print(A,5); + // printf("rand init"); + double test = 23; + // index bin op definitions + //------------------------------------------------------------------------------------------------------------------ + + GRB_TRY(GrB_Scalar_new(&Theta, GrB_BOOL)); + GRB_TRY(GrB_Scalar_setElement_BOOL(Theta, 0)); + GRB_TRY(GxB_Type_new(&Tuple, sizeof(tuple_fp64), "tuple_fp64", FP64_K)); + GRB_TRY(GxB_IndexBinaryOp_new(&MAKEFP64_op, (GxB_index_binary_function)make_fp64, Tuple, GrB_FP64, GrB_UINT64, GrB_BOOL, "make_fp64", MAKE_FP64)); + GRB_TRY(GxB_BinaryOp_new_IndexOp(&MAKEFP64_Bop, MAKEFP64_op, Theta)); + tuple_fp64 id; + memset(&id, 0, sizeof(tuple_fp64)); + id.k = INT64_MAX; + id.v = (double)(-INFINITY); + GRB_TRY(GxB_BinaryOp_new(&MAXFP64_op, (GxB_binary_function)max_fp64, Tuple, Tuple, Tuple, "max_fp64", MAX_FP64)); + GRB_TRY(GrB_Monoid_new_UDT(&MAXFP64_mon, MAXFP64_op, &id)); + GRB_TRY(GrB_Semiring_new(&IBop_MAX, MAXFP64_mon, MAKEFP64_Bop)); + //------------------------------------------------------------------------------------------------------------------ + + GRB_TRY(GrB_Matrix_nrows(&n, A)); + GRB_TRY(GrB_Matrix_ncols(&b, A)); + + // k = [+_j A(:,j)] + GRB_TRY(GrB_Vector_new(&k, GrB_FP64, n)); + GRB_TRY(GrB_Matrix_reduce_Monoid(k, NULL, NULL, GrB_PLUS_MONOID_FP64, A, NULL)); + // GxB_print(k,5); + + // m = .5[+_i k(i)] + double m; + GRB_TRY(GrB_Vector_reduce_FP64(&m, NULL, GrB_PLUS_MONOID_FP64, k, NULL)); + m *= .5; + if (m == 0) + { + *S_result = NULL; + LG_FREE_ALL; + return 0; + } + // printf("m= %f\n", m); + + // S <- I + GRB_TRY(GrB_Vector_new(&x, GrB_BOOL, n)); + GRB_TRY(GrB_assign(x, NULL, NULL, true, GrB_ALL, n, NULL)); + // GxB_print(i,5); + GRB_TRY(GrB_Matrix_diag(&S, x, 0)); + GRB_TRY(GrB_set(S, false, GxB_ISO)); + GrB_set(S, GxB_SPARSE, GxB_SPARSITY_CONTROL); + // GxB_print(S,5); + GRB_TRY(GrB_Vector_new(&v, GrB_FP64, n)); + // var used in for loop + GRB_TRY(GrB_Vector_new(&t_q, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&q1, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&z, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&za, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&dSTk, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&vtS, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&z_dSTk, GrB_FP64, n)); + // temp used to set dS to 0 matrix + GRB_TRY(GrB_Vector_new(&temp, GrB_FP64, n)); + GRB_TRY(GrB_assign(temp, NULL, NULL, 0, GrB_ALL, n, NULL)); + GRB_TRY(GxB_Container_new(&S_container)); + GRB_TRY(GrB_Vector_new(&max_q1, Tuple, 1)); + GrB_Vector_new(&srxq, GrB_FP64, n); + GRB_TRY(GrB_Matrix_diag(&dS, temp, 0)); + GRB_TRY(GrB_Vector_new(&y_rand, GrB_UINT64, n)); + GRB_TRY(GrB_Vector_new(&sr, GrB_FP64, n)); + GRB_TRY(GrB_Vector_new(&srxq, GrB_FP64, n)); + GRB_TRY(GrB_Matrix_new(&AS, GrB_FP64, n, n)); + GRB_TRY(GrB_Matrix_new(&StAS, GrB_FP64, n, n)); + GRB_TRY(GrB_assign(y_rand, NULL, NULL, 1, GrB_ALL, n, NULL)); + + bool changed = true; + int max_iter = 20; + int iter = 0; + int aggr_iter = 0; + // uint64_t seed = 1212312224; + + // GxB_print(y_rand,5); + GRB_TRY(GrB_mxv(z, NULL, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, S, k, NULL)); + // GxB_print(z,5); + while (aggr_iter < max_iter) + { + while (changed && iter < max_iter) + { + changed = false; + for (int i = 0; i < n; i++) + { // extract tuples + // printf("%d",i); + + // v = A(i,:) + GRB_TRY(GrB_Col_extract(v, NULL, NULL, A, GrB_ALL, b, i, GrB_DESC_T0)); + // GxB_print(v,5); + + // -- extract k_i + GRB_TRY(GrB_Vector_extractElement_FP64(&k_i, k, i)); + // GxB_print(S,5); + + // t_q =v any.pair S O(|v|) + GRB_TRY(GrB_vxm(t_q, NULL, NULL, GxB_ANY_PAIR_FP64, v, S, GrB_DESC_T0)); + // GxB_print(t_q,5); + + // sr = S(i,:) + GRB_TRY(GrB_Col_extract(sr, NULL, NULL, S, GrB_ALL, 1, i, GrB_DESC_T0)); + // GxB_print(sr,5); + + // S(i,:) = 0/false + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + GRB_TRY(GrB_Vector_setElement_BOOL(S_container->x, false, i)); + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + dbg(S); + + //////////////////////////////////////////////////////////// + //-------------q1 = a(kTS)+vTS----------- -----------// + //-------------q1 = z+vTS----------------------------// + // double alpha_p = 1; + double alpha = -k_i / m; + + GRB_TRY(GrB_Vector_apply_BinaryOp2nd_FP64(za, NULL, NULL, GrB_TIMES_FP64, z, alpha, GrB_DESC_T0)); + dbg(za); + + GRB_TRY(GrB_eWiseAdd(za, NULL, NULL, GrB_PLUS_FP64, za, v, NULL)); + // printf("z*alpha + v\n"); + dbg(za); + GRB_TRY(GrB_vxm(q1, t_q, NULL, GrB_PLUS_TIMES_SEMIRING_FP64, za, S, GrB_DESC_R)); + dbg(q1); + /////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////// + //-------------Index Binary OP Rand_argminmax -----------// + GRB_TRY(GrB_Vector_nvals(&q1_size, q1)); + // GxB_print(q1,5); + // printf("Size of q1: %ld\n",q1_size); + // GRB_TRY(GrB_Vector_setElement_UINT64(y_rand,seed,0)); + seed += 3; + GRB_TRY(GrB_assign(y_rand, t_q, NULL, seed, GrB_ALL, n, GrB_DESC_S)); + + // GxB_print(q1,5); + GRB_TRY(GrB_mxv(max_q1, NULL, NULL, IBop_MAX, (GrB_Matrix)q1, y_rand, GrB_DESC_T0)); + // GxB_print(q1,5); + + GRB_TRY(GrB_Vector_extractElement_UDT((void *)&o, max_q1, 0)); + // printf("choice:%ld tb: %ld\n",(long)o.k, (long)o.tb); + // dbg(S); + + GRB_TRY(GxB_unload_Matrix_into_Container(S, S_container, NULL)); + GRB_TRY(GrB_Vector_setElement(S_container->i, o.k, i)); + GRB_TRY(GrB_Vector_setElement_BOOL(S_container->x, true, i)); + GRB_TRY(GxB_load_Matrix_from_Container(S, S_container, NULL)); + // GxB_print(S,5); + + ////////////////////////////////////////////////////////// + + // GxB_print(sr,5); + GRB_TRY(GrB_Vector_eWiseMult_BinaryOp(srxq, NULL, NULL, GrB_TIMES_FP64, sr, q1, NULL)); + GRB_TRY(GrB_Vector_nvals(&vals_srxq, srxq)); + // GxB_print(srxq,5); + + if (vals_srxq == 0) + { + changed = true; + } + // if(i==3)break; + } + iter++; + // break; + + } + + aggr_iter++; + } + // printf("Iterations: %d\n", iter); + (*S_result) = S; + S = NULL; + LG_FREE_ALL; + return (GrB_SUCCESS) ; +#else + return (GrB_NOT_IMPLEMENTED); +#endif +} diff --git a/experimental/test/test_IsolateSets.c b/experimental/test/test_IsolateSets.c new file mode 100644 index 0000000000..9395769c3f --- /dev/null +++ b/experimental/test/test_IsolateSets.c @@ -0,0 +1,82 @@ +#include +#include + +#include +#include +#include "LG_Xtest.h" + +char msg[LAGRAPH_MSG_LEN]; +LAGraph_Graph G; +GrB_Matrix A; + +#define dbg(x) GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +#define LEN 512 +char filename[LEN + 1]; +typedef struct +{ + const char *matrix_file; // Adjeancy matrix or graph +} matrix_info; +const matrix_info files[] = { + {"empty.mtx"}, + {"comm0.mtx"}, + {"karate.mtx"}, + {""}}; + +void test_IsolateSets(void) +{ +#if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + printf("\n"); + // GrB_set(GrB_GLOBAL, 0, GxB_BURBLE); + GrB_Info info; + for (int k = 0;; k++) + { + if (strlen(files[k].matrix_file) == 0) + break; + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + OK(fclose(f)); + // GxB_print(A,5); + GrB_Matrix MIset = NULL; + + + double tsimple = LAGraph_WallClockTime(); + OK (LAGraph_IsolateSets(&MIset, A, 1231245, msg)); + tsimple = LAGraph_WallClockTime() - tsimple; + printf(" time: %f\n", tsimple); + + if (MIset == NULL) + { + TEST_CHECK(true); + } + else + { + GrB_Index nrows; + GrB_Matrix_nrows(&nrows, MIset); + printf("%lu\n", nrows); + GrB_Index nvals; + GrB_Matrix_nvals(&nvals, MIset); + printf("%lu\n", nvals); + TEST_CHECK(nrows <= nvals); + } + // GxB_print(MIset,5); + OK(LAGraph_Delete(&G, msg)); + GrB_Matrix_free(&A); + GrB_Matrix_free(&MIset); + } + LAGraph_Finalize(msg); +#endif +}; + +TEST_LIST = { + {"IsolateSets", test_IsolateSets}, + {NULL, NULL}}; diff --git a/experimental/test/test_louvain.c b/experimental/test/test_louvain.c new file mode 100644 index 0000000000..6fc4a1638b --- /dev/null +++ b/experimental/test/test_louvain.c @@ -0,0 +1,152 @@ +#include +#include + +#include +#include +#include "LG_Xtest.h" + +#define dbg(x) GxB_print(x, 5) +#define err(x, info) \ + if (!(info == GrB_SUCCESS || info == GrB_NO_VALUE)) \ + { \ + char **err; \ + GrB_error(err, x); \ + printf("\ninfo: %d error: %s\n", info, err); \ + } +char msg[LAGRAPH_MSG_LEN]; +LAGraph_Graph G = NULL; +GrB_Matrix A = NULL; +#define LEN 512 +char filename[LEN + 1]; +typedef struct +{ + const char *matrix_file; // Adjeancy matrix or graph + const double mod; +} matrix_info; + +const matrix_info files[] = { + {"empty.mtx",0}, + {"comm0.mtx", 0.357142857142857}, + {"karate.mtx", .42}, + {"", -1}}; + + +void test_LouvainSeq(void) +{ + #if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + // LG_SET_BURBLE (false) ; + + for (int k = 0;; k++) + { + uint64_t seed = 1224; + + const char *aname = files[k].matrix_file; + if (strlen(aname) == 0) + break; + printf("\n================================== %s:\n", aname); + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + fclose(f); + + OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); + TEST_CHECK(A == NULL); + + + OK(LAGraph_Cached_AT(G, msg)); + // check if the pattern is symmetric - if it isn't make it. + OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); + GrB_Matrix S = NULL; + printf("Sequntial Louvain"); + for (int jit = 0 ; jit <= 1 ; jit++) + { + OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; + + double tsimple = LAGraph_WallClockTime(); + OK(LAGraph_LouvainSeq(&S, G,seed, msg)); + + // OK(LAGraph_Louvain_res(&S,G,.3,msg)); + tsimple = LAGraph_WallClockTime() - tsimple; + double Q = 0.0; + OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); + printf("Q:%f\n", Q); + // printf("Number of Communities: %d",comms); + printf(" time: %f\n", tsimple); + } + + OK(LAGraph_Delete(&G, msg)); + } + LAGraph_Finalize(msg); + #endif +} +void test_LouvainIS(void) +{ + #if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + // LG_SET_BURBLE (true) ; + + for (int k = 0;; k++) + { + uint64_t seed = 1249141465; + const char *aname = files[k].matrix_file; + if (strlen(aname) == 0) + break; + printf("\n================================== %s:\n", aname); + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + fclose(f); + + OK(LAGraph_New(&G, &A, LAGraph_ADJACENCY_DIRECTED, msg)); + TEST_CHECK(A == NULL); + + + OK(LAGraph_Cached_AT(G, msg)); + // check if the pattern is symmetric - if it isn't make it. + OK(LAGraph_Cached_IsSymmetricStructure(G, msg)); + GrB_Matrix S = NULL; + printf("Isolate Sets Louvain\n"); + for (int jit = 0 ; jit <= 1 ; jit++) + { + OK (LG_SET_JIT (jit ? GxB_JIT_ON : GxB_JIT_OFF)) ; + + double tsimple = LAGraph_WallClockTime(); + OK(LAGraph_LouvainIS(&S,seed, G, msg)); + // GrB_Info info = (LAGraph_LouvainIS(&S,seed, G, msg)); + // printf ("info %d\nmsg: %s\n", info, msg) ; + // OK (info) ; + + // OK(LAGraph_Louvain_res(&S,G,.3,msg)); + tsimple = LAGraph_WallClockTime() - tsimple; + double Q = 0.0; + double tsimple2 = LAGraph_WallClockTime(); + OK(LAGr_AdjModularity(&Q, 1.0, G->A, S, msg)); + tsimple2 = LAGraph_WallClockTime() - tsimple2; + + printf("Q:%f time to calc Q: %f\n", Q,tsimple2); + // printf("Number of Communities: %d",comms); + printf(" time: %f\n", tsimple); + GrB_free(&S); + } + + OK(LAGraph_Delete(&G, msg)); + + } + + LAGraph_Finalize(msg); + #endif +} + +TEST_LIST = { + + {"LouvainSeq", test_LouvainSeq}, + {"LouvainIS", test_LouvainIS}, + + {NULL, NULL}}; + + + + diff --git a/experimental/test/test_modularity.c b/experimental/test/test_modularity.c new file mode 100644 index 0000000000..7a22a5e73f --- /dev/null +++ b/experimental/test/test_modularity.c @@ -0,0 +1,79 @@ + +#include +#include + +#include "LG_Xtest.h" +#include +#include +char msg[LAGRAPH_MSG_LEN]; +LAGraph_Graph G = NULL; +GrB_Matrix A = NULL; +GrB_Matrix S = NULL; +#define LEN 512 +char filename[LEN + 1]; +char filename2[LEN + 1]; +typedef struct +{ + const char *matrix_file; // Adjeancy matrix or graph + const char *comm_matrix; // comm matrix define as pattern as general + const double gamma; // resolution of communities 1 for now + const double mod; // expected modularity from Q = 1/2m (S^TBS) +} matrix_info; + +const matrix_info files[] = { + {"empty.mtx","empty.mtx",1,0}, + {"comm0.mtx", "comm0_S.mtx", 1, -0.17347}, + {"comm0.mtx", "comm0_Sa.mtx", 1, 0.35714}, + // {"com-Amazon.mtx", "comm0_Sa.mtx",1, -1}, + {"", "", -1, -1}}; +void test_modularity(void) +{ + #if LG_SUITESPARSE_GRAPHBLAS_V10_2 + LAGraph_Init(msg); + for (int k = 0;; k++) + { + if (strlen(files[k].matrix_file) == 0) + break; + snprintf(filename, LEN, LG_DATA_DIR "%s", files[k].matrix_file); + FILE *f = fopen(filename, "r"); + TEST_CHECK(f != NULL); + OK(LAGraph_MMRead(&A, f, msg)); + OK(fclose(f)); + printf("\nInput of Matrix:\n"); + GxB_print(A, 3); + + snprintf(filename2, LEN, LG_DATA_DIR "%s", files[k].comm_matrix); + FILE *t = fopen(filename2, "r"); + TEST_CHECK(t != NULL); + OK(LAGraph_MMRead(&S, t, msg)); + OK(fclose(t)); + printf("\nInput of Matrix S:\n"); + GxB_print(S, 3); + + double gamma = files[k].gamma; + double Q; + OK(LAGr_AdjModularity(&Q, gamma, A, S, msg)); + Q = floor(100000*Q)/100000; + TEST_CHECK(Q == files[k].mod); + printf("Q:%.15g\n", Q); + OK(GrB_free(&A)); + OK(GrB_free(&S)); + } + + OK(GrB_free(&A)); + OK(GrB_free(&S)); + + OK(LAGraph_Delete(&G, msg)); + + LAGraph_Finalize(msg); + #endif +} + +//---------------------------------------------------------------------------- +// the make program is created by acutest, and it runs a list of tests: +//---------------------------------------------------------------------------- + +TEST_LIST = + { + {"Modularity", test_modularity}, // just one test in this example + {NULL, NULL}}; diff --git a/include/LAGraphX.h b/include/LAGraphX.h index e3e593d0f7..6d4f3235c0 100644 --- a/include/LAGraphX.h +++ b/include/LAGraphX.h @@ -1530,6 +1530,20 @@ int LAGraph_MinCut( char *msg ); +//------------------------------------------------------------------------------ +// Louvain sub-algorithms +//------------------------------------------------------------------------------ + +LAGRAPHX_PUBLIC +int LAGr_AdjModularity( + //output + double *Q, + //input + double gamma, + GrB_Matrix A, + GrB_Matrix S, + char* msg +); LAGRAPHX_PUBLIC int LAGr_Jaccard( // output @@ -1540,6 +1554,46 @@ int LAGr_Jaccard( char *msg ); +LAGRAPHX_PUBLIC +int LAGraph_IsolateSet( + //output + GrB_Vector *isolate_set, + //input + GrB_Matrix A, + GrB_Vector ignore_node, + uint64_t seed, + char* msg +); + +LAGRAPHX_PUBLIC +int LAGraph_IsolateSets( + GrB_Matrix *IsolateSets, + // LAGraph_Graph G, + GrB_Matrix A, + // GrB_Vector ignore_nodes, + uint64_t seed, + char *msg +); + +LAGRAPHX_PUBLIC +int LAGraph_LouvainSeq( + // output + GrB_Matrix *S_result, + // input + LAGraph_Graph G, + uint64_t seed, + char *msg +); + +LAGRAPHX_PUBLIC +int LAGraph_LouvainIS( + // output + GrB_Matrix *S_result, + uint64_t seed, + // input + LAGraph_Graph G, + char *msg +); #if defined ( __cplusplus ) } #endif