pathrex-sys 0.1.0

//------------------------------------------------------------------------------
// LAGraph_FastGraphletTransform: fast graphlet transform
//------------------------------------------------------------------------------

// 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 Tanner Hoke, Texas A&M University

//------------------------------------------------------------------------------

// LAGraph_FastGraphletTransform: computes the Fast Graphlet Transform of
// an undirected graph.  No self edges are allowed on the input graph.

// TODO: rename this, and add to src (but it uses GxB)

// https://arxiv.org/pdf/2007.11111.pdf

//------------------------------------------------------------------------------
#define F_UNARY(f)  ((void (*)(void *, const void *)) f)

#define LG_FREE_WORK                \
{                                   \
    GrB_free (&C_3) ;               \
    GrB_free (&d_0) ;               \
    GrB_free (&d_1) ;               \
    GrB_free (&d_2) ;               \
    GrB_free (&d_3) ;               \
    GrB_free (&d_4) ;               \
    GrB_free (&d_5) ;               \
    GrB_free (&d_6) ;               \
    GrB_free (&d_7) ;               \
    GrB_free (&d_8) ;               \
    GrB_free (&d_9) ;               \
    GrB_free (&d_10) ;              \
    GrB_free (&d_11) ;              \
    GrB_free (&d_12) ;              \
    GrB_free (&d_13) ;              \
    GrB_free (&d_14) ;              \
    GrB_free (&d_15) ;              \
    GrB_free (&d_2) ;               \
    GrB_free (&v) ;                 \
    GrB_free (&p_1_minus_one) ;     \
    GrB_free (&p_1_minus_two) ;     \
    GrB_free (&two_c_3) ;           \
    GrB_free (&p_1_p_1_had) ;       \
    GrB_free (&C_42) ;              \
    GrB_free (&P_2) ;               \
    GrB_free (&D_1) ;               \
    GrB_free (&D_4c) ;              \
    GrB_free (&D_43) ;              \
    GrB_free (&U_inv) ;             \
    GrB_free (&F_raw) ;             \
    GrB_free (&C_4) ;               \
    GrB_free (&Sub_one_mult) ;      \
    GrB_free (&T) ;                 \
    if (A_Tiles != NULL)                                                \
    {                                                                   \
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&A_Tiles [i]) ;    \
    }                                                                   \
    if (D_Tiles != NULL)                                                \
    {                                                                   \
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&D_Tiles [i]) ;    \
    }                                                                   \
    if (C_Tiles != NULL)                                                \
    {                                                                   \
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&C_Tiles [i]) ;    \
    }                                                                   \
    LAGraph_Free ((void **) &neighbors, msg) ;      \
    LAGraph_Free ((void **) &k4cmn, msg) ;          \
    LAGraph_Free ((void **) &f15, msg) ;            \
    LAGraph_Free ((void **) &I, msg) ;              \
    LAGraph_Free ((void **) &isNeighbor, msg) ;     \
    LAGraph_Free ((void **) &J, msg) ;              \
    LAGraph_Free ((void **) &vals, msg) ;           \
    LAGraph_Free ((void **) &A_Tiles, msg) ;        \
    LAGraph_Free ((void **) &D_Tiles, msg) ;        \
    LAGraph_Free ((void **) &C_Tiles, msg) ;        \
    LAGraph_Free ((void **) &Tile_nrows, msg) ;     \
}

#define LG_FREE_ALL                 \
{                                   \
    LG_FREE_WORK ;                  \
}

#define F_UNARY(f)  ((void (*)(void *, const void *)) f)

#include "LG_internal.h"
#include "LAGraphX.h"

void sub_one_mult (int64_t *z, const int64_t *x) { (*z) = (*x) * ((*x)-1) ; }

int LAGraph_FastGraphletTransform
(
    // outputs:
    GrB_Matrix *F_net,  // 16-by-n matrix of graphlet counts
    // inputs:
    LAGraph_Graph G,
    bool compute_d_15,  // probably this makes most sense
    char *msg
)
{
#if LAGRAPH_SUITESPARSE

    LG_CLEAR_MSG ;
    GrB_Index const U_inv_I[] = {0, 1, 2, 2, 3, 3, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 12, 12, 12, 12, 13, 13, 14, 14, 15} ;
    GrB_Index const U_inv_J[] = {0, 1, 2, 4, 3, 4, 4, 5, 9, 10, 12, 13, 14, 15, 6, 10, 11, 12, 13, 14, 15, 7, 9, 10, 13, 14, 15, 8, 11, 14, 15, 9, 13, 15, 10, 13, 14, 15, 11, 14, 15, 12, 13, 14, 15, 13, 15, 14, 15, 15} ;
    int64_t const U_inv_X[] = {1, 1, 1, -2, 1, -1, 1, 1, -2, -1, -2, 4, 2, -6, 1, -1, -2, -2, 2, 4, -6, 1, -1, -1, 2, 1, -3, 1, -1, 1, -1, 1, -2, 3, 1, -2, -2, 6, 1, -2, 3, 1, -1, -1, 3, 1, -3, 1, -3, 1} ;
    GrB_Index const U_inv_nvals = 50;
    GrB_UnaryOp Sub_one_mult = NULL ;
    int tile_cnt = 0 ;
    GrB_Matrix *A_Tiles = NULL ;
    GrB_Matrix *D_Tiles = NULL ;
    GrB_Matrix *C_Tiles = NULL ;
    GrB_Index *Tile_nrows = NULL ;
    GrB_Matrix T = NULL ;

    GrB_Index *neighbors = NULL ;
    GrB_Index *k4cmn = NULL ;
    int64_t *f15 = NULL ;
    GrB_Index *I = NULL ;
    int *isNeighbor = NULL ;
    GrB_Index *J = NULL ;
    int64_t *vals = NULL ;

    GrB_Matrix C_3 = NULL,
	       A = NULL,
	       C_42 = NULL,
	       P_2 = NULL,
	       D_1 = NULL,
	       D_4c = NULL,
	       D_43 = NULL,
	       U_inv = NULL,
	       F_raw = NULL,
               C_4 = NULL ;

    GrB_Vector d_0 = NULL,
	       d_1 = NULL,
	       d_2 = NULL,
	       d_3 = NULL,
	       d_4 = NULL,
	       d_5 = NULL,
	       d_6 = NULL,
	       d_7 = NULL,
	       d_8 = NULL,
	       d_9 = NULL,
	       d_10 = NULL,
	       d_11 = NULL,
	       d_12 = NULL,
               d_13 = NULL,
               d_14 = NULL,
               d_15 = NULL;

    GrB_Vector v = NULL,
               two_c_3 = NULL,
	       p_1_minus_one = NULL,
	       p_1_minus_two = NULL,
	       p_1_p_1_had = NULL ;

    GrB_Index nvals ;
    int64_t ntri ;

    A = G->A ;

    GrB_Index n ;
    GRB_TRY (GrB_Matrix_nrows (&n, A)) ;

    //--------------------------------------------------------------------------
    // compute d_0 = e
    //--------------------------------------------------------------------------

    // d_0 = e
    GRB_TRY (GrB_Vector_new (&d_0, GrB_INT64, n)) ;
    GRB_TRY (GrB_assign (d_0, NULL, NULL, 1, GrB_ALL, n, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_1 = Ae (in_degree)
    //--------------------------------------------------------------------------

//  GRB_TRY (GrB_Vector_new (&d_1, GrB_INT64, n)) ;

    // d_1 = Ae (in_degree)
    LG_TRY (LAGraph_Cached_OutDegree (G, msg)) ;

    GRB_TRY (GrB_Vector_dup (&d_1, G->out_degree)) ;

    //--------------------------------------------------------------------------
    // compute d_2 = p_2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_2, GrB_INT64, n)) ;

    // TODO: use LAGraph_plus_second_int64
    // d_2 = p_2 = A*p_1 - c_2 = A*d_1 - d_1
    GRB_TRY (GrB_mxv (d_2, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_1, NULL)) ;
    GRB_TRY (GrB_eWiseMult (d_2, NULL, NULL, GrB_MINUS_INT64, d_2, d_1, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_3 = hadamard(p_1, p_1 - 1) / 2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_3, GrB_INT64, n)) ;

    GRB_TRY (GrB_UnaryOp_new (&Sub_one_mult, F_UNARY (sub_one_mult), GrB_INT64, GrB_INT64)) ;

    GRB_TRY (GrB_apply (d_3, NULL, NULL, Sub_one_mult, d_1, NULL)) ;
    GRB_TRY (GrB_apply (d_3, NULL, NULL, GrB_DIV_INT64, d_3, (int64_t) 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_4 = C_3e/2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&C_3, GrB_INT64, n, n)) ;
    GRB_TRY (GrB_Vector_new (&d_4, GrB_INT64, n)) ;

    // TODO: use LAGraph_plus_first_int64
    // C_3 = hadamard(A, A^2)
    GRB_TRY (GrB_mxm (C_3, A, NULL, GxB_PLUS_FIRST_INT64, A, A, GrB_DESC_ST1)) ;

    // d_4 = c_3 = C_3e/2
    GRB_TRY (GrB_reduce (d_4, NULL, NULL, GrB_PLUS_MONOID_INT64, C_3, NULL)) ;
    GRB_TRY (GrB_apply (d_4, NULL, NULL, GrB_DIV_INT64, d_4, (int64_t) 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_5 = p_3 = A*d_2 - hadamard(p_1, p_1 - 1) - 2c_3
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&v, GrB_INT64, n)) ;
    GRB_TRY (GrB_Vector_new (&two_c_3, GrB_INT64, n)) ;
    GRB_TRY (GrB_Vector_new (&d_5, GrB_INT64, n)) ;

    // v = hadamard(p_1, p_1 - 1)
    GRB_TRY (GrB_apply (v, NULL, NULL, Sub_one_mult, d_1, NULL)) ;

    // two_c_3 = 2 * c_3 = 2 * d_4
    GRB_TRY (GrB_apply (two_c_3, NULL, NULL, GrB_TIMES_INT64, 2, d_4, NULL)) ;

    // d_5 = A * d_2
    // TODO: use LAGraph_plus_second_int64
    GRB_TRY (GrB_mxv (d_5, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_2, NULL)) ;

    // d_5 -= hadamard(p_1, p_1 - 1)
    GRB_TRY (GrB_eWiseAdd (d_5, NULL, NULL, GrB_MINUS_INT64, d_5, v, NULL)) ;

    // d_5 -= two_c_3
    GRB_TRY (GrB_eWiseAdd (d_5, NULL, NULL, GrB_MINUS_INT64, d_5, two_c_3, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_6 = hadamard(d_2, p_1-1) - 2c_3
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&p_1_minus_one, GrB_INT64, n)) ;
    GRB_TRY (GrB_Vector_new (&d_6, GrB_INT64, n)) ;

    // p_1_minus_one = p_1 - 1
    GRB_TRY (GrB_apply (p_1_minus_one, NULL, NULL, GrB_MINUS_INT64, d_1, (int64_t) 1, NULL)) ;

    // d_6 = hadamard(d_2, p_1-1)
    GRB_TRY (GrB_eWiseMult (d_6, NULL, NULL, GrB_TIMES_INT64, d_2, p_1_minus_one, NULL)) ;

    // d_6 -= 2c_3
    GRB_TRY (GrB_eWiseAdd (d_6, NULL, NULL, GrB_MINUS_INT64, d_6, two_c_3, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_7 = A*hadamard(p_1-1, p_1-2) / 2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&p_1_minus_two, GrB_INT64, n)) ;
    GRB_TRY (GrB_Vector_new (&p_1_p_1_had, GrB_INT64, n)) ;
    GRB_TRY (GrB_Vector_new (&d_7, GrB_INT64, n)) ;

    GRB_TRY (GrB_apply (p_1_minus_two, NULL, NULL, GrB_MINUS_INT64, d_1, (int64_t) 2, NULL)) ;
    GRB_TRY (GrB_eWiseMult (p_1_p_1_had, NULL, NULL, GrB_TIMES_INT64, p_1_minus_one, p_1_minus_two, NULL)) ;

    // TODO: use LAGraph_plus_second_int64
    GRB_TRY (GrB_mxv (d_7, NULL, NULL, GxB_PLUS_SECOND_INT64, A, p_1_p_1_had, NULL)) ;
    GRB_TRY (GrB_apply (d_7, NULL, NULL, GrB_DIV_INT64, d_7, (int64_t) 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_8 = hadamard(p_1, p_1_p_1_had) / 6
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_8, GrB_INT64, n)) ;

    GRB_TRY (GrB_eWiseMult (d_8, NULL, NULL, GrB_TIMES_INT64, d_1, p_1_p_1_had, NULL)) ;
    GRB_TRY (GrB_apply (d_8, NULL, NULL, GrB_DIV_INT64, d_8, (int64_t) 6, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_9 = A*c_3 - 2*c_3
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_9, GrB_INT64, n)) ;

    // TODO: use LAGraph_plus_second_int64
    GRB_TRY (GrB_mxv (d_9, NULL, NULL, GxB_PLUS_SECOND_INT64, A, d_4, NULL)) ;
    GRB_TRY (GrB_eWiseAdd (d_9, NULL, NULL, GrB_MINUS_INT64, d_9, two_c_3, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_10 = C_3 * (p_1 - 2)
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_10, GrB_INT64, n)) ;

    // TODO: use GrB_PLUS_TIMES_INT64
    GRB_TRY (GrB_mxv (d_10, NULL, NULL, GxB_PLUS_TIMES_INT64, C_3, p_1_minus_two, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_11 = hadamard(p_1 - 2, c_3)
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Vector_new (&d_11, GrB_INT64, n)) ;

    GRB_TRY (GrB_eWiseMult (d_11, NULL, NULL, GrB_TIMES_INT64, p_1_minus_two, d_4, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_12 = c_4 = C_{4,2}e/2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&C_4, GrB_INT64, n, 1)) ;
    GRB_TRY (GrB_Matrix_new (&D_1, GrB_INT64, n, n)) ;

    GRB_TRY (GrB_Vector_new (&d_12, GrB_INT64, n)) ;

    // D_1 = diag(d_1)
    // TODO: use GrB_Matrix_diag
    GRB_TRY (GxB_Matrix_diag (D_1, d_1, (int64_t) 0, NULL)) ;

    // TODO: true GxB starts here.  For vanilla method, put this
    // in a single call to GrB_mxm to compute C_4

    GRB_TRY (GrB_Matrix_nvals (&nvals, A));

    const GrB_Index entries_per_tile = 1000;
    GrB_Index ntiles = (nvals + entries_per_tile - 1) / entries_per_tile ;

    LG_TRY (LAGraph_Calloc ((void **) &A_Tiles, ntiles, sizeof (GrB_Matrix),
        msg)) ;
    LG_TRY (LAGraph_Calloc ((void **) &D_Tiles, ntiles, sizeof (GrB_Matrix),
        msg)) ;
    LG_TRY (LAGraph_Calloc ((void **) &C_Tiles, ntiles, sizeof (GrB_Matrix),
        msg)) ;
    LG_TRY (LAGraph_Calloc ((void **) &Tile_nrows, ntiles, sizeof (GrB_Index),
        msg)) ;

    GrB_Index Tile_ncols [1] = {n} ;

    int64_t tot_deg = 0 ;
    GrB_Index last_row = -1 ;
    for (GrB_Index i = 0; i < n; ++i)
    {
        int64_t deg = 0 ;
        GRB_TRY (GrB_Vector_extractElement (&deg, d_1, i)) ;
        if (i == n - 1 || (tot_deg / entries_per_tile != (tot_deg + deg) / entries_per_tile))
        {
            Tile_nrows [tile_cnt++] = i - last_row ;
            last_row = i ;
        }
        tot_deg += deg ;
    }

    GRB_TRY (GxB_Matrix_split (A_Tiles, tile_cnt, 1, Tile_nrows, Tile_ncols, A, NULL)) ;
    GRB_TRY (GxB_Matrix_split (D_Tiles, tile_cnt, 1, Tile_nrows, Tile_ncols, D_1, NULL)) ;

#define TRY(method)                         \
    {                                       \
        GrB_Info info2 = method ;           \
        if (info2 != GrB_SUCCESS)           \
        {                                   \
            GrB_free (&A_i) ;               \
            GrB_free (&C_Tiles [i_tile]) ;  \
            GrB_free (&e) ;                 \
            info1 = info2 ;                 \
            continue ;                      \
        }                                   \
    }

    int save_nthreads_outer, save_nthreads_inner ;
    LG_TRY (LAGraph_GetNumThreads (&save_nthreads_outer, &save_nthreads_inner, msg)) ;
    LG_TRY (LAGraph_SetNumThreads (1, 1, msg)) ;

    int i_tile;
    GrB_Info info1 = GrB_SUCCESS ;
    #pragma omp parallel for num_threads(omp_get_max_threads()) schedule(dynamic,1)
    for (i_tile = 0; i_tile < tile_cnt; ++i_tile)
    {
        GrB_Matrix A_i = NULL, e = NULL ;

        TRY (GrB_Matrix_new (&e, GrB_INT64, n, 1)) ;
        TRY (GrB_assign (e, NULL, NULL, (int64_t) 1, GrB_ALL, n, GrB_ALL, 1, NULL)) ;

        TRY (GrB_Matrix_new (&A_i, GrB_INT64, Tile_nrows [i_tile], n)) ;
        TRY (GrB_Matrix_new (&C_Tiles [i_tile], GrB_INT64, Tile_nrows [i_tile], 1)) ;

        TRY (GrB_mxm (A_i, NULL, NULL, GxB_PLUS_PAIR_INT64, A_Tiles [i_tile], A, NULL)) ;
        TRY (GrB_eWiseAdd (A_i, NULL, NULL, GrB_MINUS_INT64, A_i, D_Tiles [i_tile], NULL)) ;
        TRY (GrB_apply (A_i, NULL, NULL, Sub_one_mult, A_i, NULL)) ;

        // multiply A_i by it on the right
        TRY (GrB_mxm (C_Tiles [i_tile], NULL, NULL, GxB_PLUS_FIRST_INT64, A_i, e, NULL)) ;

        GrB_free (&A_i) ;
        GrB_free (&e) ;
    }

    GRB_TRY (info1) ;

    LG_TRY (LAGraph_SetNumThreads (save_nthreads_outer, save_nthreads_inner, msg)) ;

    GRB_TRY (GxB_Matrix_concat (C_4, C_Tiles, tile_cnt, 1, NULL)) ;

    if (A_Tiles != NULL)
    {
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&A_Tiles [i]) ;
    }
    if (D_Tiles != NULL)
    {
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&D_Tiles [i]) ;
    }
    if (C_Tiles != NULL)
    {
        for (int i = 0; i < tile_cnt; ++i) GrB_free (&C_Tiles [i]) ;
    }
    LAGraph_Free ((void **) &Tile_nrows, msg) ;
    LAGraph_Free ((void **) &A_Tiles, msg) ;
    LAGraph_Free ((void **) &D_Tiles, msg) ;
    LAGraph_Free ((void **) &C_Tiles, msg) ;

    // d_12 = sum (C_4)
    GRB_TRY (GrB_reduce (d_12, NULL, NULL, GrB_PLUS_MONOID_INT64, C_4, NULL)) ;
    GRB_TRY (GrB_apply (d_12, NULL, NULL, GrB_DIV_INT64, d_12, 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_13 = D_{4,c}e/2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&D_4c, GrB_INT64, n, n)) ;
    GRB_TRY (GrB_Vector_new (&d_13, GrB_INT64, n)) ;

    GRB_TRY (GrB_eWiseMult (D_4c, NULL, NULL, GrB_MINUS_INT64, C_3, A, NULL)) ; // can be mult because we mask with A next
    GRB_TRY (GrB_mxm (D_4c, A, NULL, GxB_PLUS_SECOND_INT64, A, D_4c, GrB_DESC_S)) ;

    // d_13  = D_{4,c}*e/2
    GRB_TRY (GrB_reduce (d_13, NULL, NULL, GrB_PLUS_INT64, D_4c, NULL)) ;
    GRB_TRY (GrB_apply (d_13, NULL, NULL, GrB_DIV_INT64, d_13, (int64_t) 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_14 = D_{4,3}e/2 = hadamard(A, C_42)e/2
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&D_43, GrB_INT64, n, n)) ;
    GRB_TRY (GrB_Vector_new (&d_14, GrB_INT64, n)) ;
    GRB_TRY (GrB_Matrix_new (&C_42, GrB_INT64, n, n)) ;
    GRB_TRY (GrB_Matrix_new (&P_2, GrB_INT64, n, n)) ;

    // P_2 = A*A - diag(d_1)
    GRB_TRY (GrB_eWiseAdd (P_2, A, NULL, GrB_MINUS_INT64, C_3, D_1, NULL)) ;

    // C_42 = hadamard(P_2, P_2 - 1)
    GRB_TRY (GrB_apply (C_42, A, NULL, Sub_one_mult, P_2, NULL)) ;

    GRB_TRY (GrB_eWiseMult (D_43, NULL, NULL, GrB_TIMES_INT64, A, C_42, NULL)) ;

    // d_14  = D_{4,3}*e/2
    GRB_TRY (GrB_reduce (d_14, NULL, NULL, GrB_PLUS_INT64, D_43, NULL)) ;
    GRB_TRY (GrB_apply (d_14, NULL, NULL, GrB_DIV_INT64, d_14, (int64_t) 2, NULL)) ;

    //--------------------------------------------------------------------------
    // compute d_15 = Te/6
    //--------------------------------------------------------------------------

    // TODO: computing d_15 requires GxB iterators.  For vanilla GrB,
    // export the CSR matrix and use it directly.

    if (compute_d_15)
    {
        LG_TRY (LAGraph_KTruss (&T, G, 4, msg)) ;
        GRB_TRY (GrB_Vector_new (&d_15, GrB_INT64, n)) ;

        int nthreads = 1 ;
        // todo: parallelize this...
//#pragma omp parallel for num_threads(nthreads)
        //for (int tid = 0 ; tid < nthreads ; tid++)
        {

            // allocate workspcae
            LG_TRY (LAGraph_Malloc ((void **) &neighbors, n,
                sizeof (GrB_Index), msg)) ;
            LG_TRY (LAGraph_Malloc ((void **) &k4cmn, n,
                sizeof (GrB_Index), msg)) ;
            LG_TRY (LAGraph_Malloc ((void **) &f15, n,
                sizeof (int64_t), msg)) ;
            LG_TRY (LAGraph_Malloc ((void **) &I, n,
                sizeof (GrB_Index), msg)) ;
            LG_TRY (LAGraph_Malloc ((void **) &isNeighbor, n,
                sizeof (int), msg)) ;

            for (int i = 0; i < n; ++i) {
                neighbors [i] = k4cmn [i] = f15 [i] = isNeighbor [i] = 0 ;
                I [i] = i ;
            }

            // thread tid operates on T(row1:row2-1,:)
            GrB_Index row1 = 0;//tid * (n / nthreads) ;
            GrB_Index row2 = n;//(tid == nthreads - 1) ? n : ((tid+1) * (n / nthreads)) ;

            GxB_Iterator riterator ;
            GxB_Iterator_new (&riterator) ;
            GRB_TRY (GxB_rowIterator_attach (riterator, T, NULL)) ;

            GxB_Iterator iterator ;
            GxB_Iterator_new (&iterator) ;
            GRB_TRY (GxB_rowIterator_attach (iterator, T, NULL)) ;

            // seek to T(row1,:)
            GrB_Info info = GxB_rowIterator_seekRow (iterator, row1) ;
            while (info != GxB_EXHAUSTED)
            {
                // iterate over entries in T(i,:)
                GrB_Index idx2 = GxB_rowIterator_getRowIndex (iterator) ;
                if (idx2 >= row2) break ;
                int neighbor_cnt = 0 ;
                while (info == GrB_SUCCESS)
                {
                    // working with edge (idx2, j)
                    GrB_Index j = GxB_rowIterator_getColIndex (iterator) ;

                    if (j > idx2) {
                        neighbors [neighbor_cnt++] = j ;
                        isNeighbor [j] = 1 ;
                    }

                    info = GxB_rowIterator_nextCol (iterator) ;
                }

                for (int neighbor_id = 0 ; neighbor_id < neighbor_cnt ; ++neighbor_id) {
                    GrB_Index j = neighbors [neighbor_id] ;
                    int cmn_cnt = 0 ;
                    info = GxB_rowIterator_seekRow(riterator, j) ;

                    while (info == GrB_SUCCESS) { // iterate over neighbors of j
                        GrB_Index k = GxB_rowIterator_getColIndex (riterator) ;
                        if (k > j && isNeighbor [k]) {
                            k4cmn [cmn_cnt++] = k ;
                            isNeighbor [k] = -1 ;
                        }
                        info = GxB_rowIterator_nextCol (riterator) ;
                    }
                    // check every combination
                    for (int k_1 = 0 ; k_1 < cmn_cnt ; k_1++) {
                        GrB_Index k = k4cmn [k_1] ;
                        info = GxB_rowIterator_seekRow(riterator, k) ;

                        while (info == GrB_SUCCESS) { // iterate over neighbors of k
                            GrB_Index l = GxB_rowIterator_getColIndex (riterator) ;
                            if (l > k && isNeighbor [l] == -1) {
                                f15[idx2]++ ;
                                f15[j]++ ;
                                f15[k]++ ;
                                f15[l]++ ;
                            }
                            info = GxB_rowIterator_nextCol (riterator) ;
                        }
                    }
                    for (int k_1 = 0 ; k_1 < cmn_cnt ; k_1++) {
                        isNeighbor[k4cmn[k_1]] = 1 ;
                    }
                }

                for (int neighbor_id = 0 ; neighbor_id < neighbor_cnt ; ++neighbor_id) {
                    GrB_Index j = neighbors [neighbor_id] ;
                    isNeighbor [j] = 0 ;
                }

                // move to the next row, T(i+1,:)
                info = GxB_rowIterator_nextRow (iterator) ;
            }
            GrB_free (&iterator) ;
            GrB_free (&riterator) ;
            GrB_free (&T) ;
            GRB_TRY (GrB_Vector_build (d_15, I, f15, n, NULL)) ;

            // free workspace
            LAGraph_Free ((void **) &neighbors, msg) ;
            LAGraph_Free ((void **) &k4cmn, msg) ;
            LAGraph_Free ((void **) &f15, msg) ;
            LAGraph_Free ((void **) &I, msg) ;
            LAGraph_Free ((void **) &isNeighbor, msg) ;
        }
    }

    //--------------------------------------------------------------------------
    // construct raw frequencies matrix F_raw
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&F_raw, GrB_INT64, 16, n)) ;

    GrB_Vector d[16] = {d_0, d_1, d_2, d_3, d_4, d_5, d_6, d_7, d_8, d_9, d_10, d_11, d_12, d_13, d_14, d_15} ;

    for (int i = 0; i < 15 + (compute_d_15 ? 1 : 0); ++i)
    {
        GRB_TRY (GrB_Vector_nvals (&nvals, d[i]));

        // allocate workspace
        LG_TRY (LAGraph_Malloc ((void **) &J, nvals, sizeof (GrB_Index), msg)) ;
        LG_TRY (LAGraph_Malloc ((void **) &vals, nvals, sizeof (int64_t),
            msg)) ;

        GRB_TRY (GrB_Vector_extractTuples (J, vals, &nvals, d[i])) ;
        for (int j = 0; j < nvals; ++j) {
            GRB_TRY (GrB_Matrix_setElement (F_raw, vals[j], i, J[j])) ;
        }

        // free workspace
        LAGraph_Free ((void **) &J, msg) ;
        LAGraph_Free ((void **) &vals, msg) ;
    }

    //--------------------------------------------------------------------------
    // construct U_inv
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (&U_inv, GrB_INT64, 16, 16)) ;
    GRB_TRY (GrB_Matrix_build (U_inv, U_inv_I, U_inv_J, U_inv_X, U_inv_nvals, GrB_PLUS_INT64)) ;

    //--------------------------------------------------------------------------
    // construct net frequencies matrix F_net
    //--------------------------------------------------------------------------

    GRB_TRY (GrB_Matrix_new (F_net, GrB_INT64, 16, n)) ;
    GRB_TRY (GrB_mxm (*F_net, NULL, NULL, GxB_PLUS_TIMES_INT64, U_inv, F_raw, NULL)) ;
    GrB_Vector f_net = NULL ;
    GRB_TRY (GrB_Vector_new (&f_net, GrB_INT64, 16)) ;
    GRB_TRY (GrB_reduce (f_net, NULL, NULL, GrB_PLUS_INT64, *F_net, NULL)) ;
    GRB_TRY (GrB_free (&f_net)) ;

    //--------------------------------------------------------------------------
    // free work
    //--------------------------------------------------------------------------

    LG_FREE_WORK ;
    return (GrB_SUCCESS) ;
#else
    return (GrB_NOT_IMPLEMENTED) ;
#endif
}