triangles_in_sparsegraph.py triangles_to_edges.py 17 Nov 2020

triangles_in_sparsegraph.py

#!/usr/bin/env python3
""" triangles in a sparse graph """
# Keywords: sparse-graph, sparse-matrix, scipy, mesh, triangles
# w Polygon_mesh Mesh_generation ...

import numpy as np
from scipy import sparse
    # https://docs.scipy.org/doc/scipy/reference/sparse.csgraph.html

# some type names --
Node        = "id 0 .. Nnode"
Nodes       = "[node ...]"  # np.ndarray, think list / row / bitset of Nodes
Edges       = "[[i j: Node] ...]"  # ndarray, i < j ?
Triangles   = "[[i j k: Node] ...]"  # ndarray, i < j < k ?
Sparsegraph = "scipy.sparse matrix of 0 / 1"  # no edge weights, .data ignored
Sparserow   = "Sparsegraph[i]"  # ins, outs: Node -> Nodes

#...............................................................................
def triangles( G: Sparsegraph ) -> Triangles:
    """ -> [[i j k] ...] with edges i -- j -- k in G """
    G = G.tocsr()
    GT = G.T.tocsc()  # faster, 2*mem: outs() row, ins() col of G
    G2 = G * G  # G2 ik = sum Gij Gjk = nr paths i -- any j -- k
    GandG2 = G.multiply( G2 )  # pointwise, edges direct or pathlen 2
    I, K = GandG2.nonzero()
    tris = []
    for i, k in zip( I, K ):
        J = intersect( outs( G, i ), ins( GT, k ))  # i -- j -- k
        if len(J) > 0:
            tris.extend( [[i, j, k] for j in J] )
    if len(tris) == 0:
        print( "Warning: no triangles in G %s  %.0f nnz/row" % (
                G.shape, G.nnz / G.shape[0] ))
        return []
    return np.sort( tris, axis=1 )  # , GandG2, I, K


def intersect( I, J: Nodes ) -> Nodes:
    I = _toix( I )
    J = _toix( J )
    return np.intersect1d( I, J, assume_unique=True )
        # if sorted, O( min( len I, len J ))
        # historical footnote: fault list Bell Labs


def _ix( S: Sparserow, i: Node ) -> Nodes:
    """ -> [j ...] where S[i,j] != 0 """
    assert isinstance( S, (sparse.csr_matrix, sparse.csc_matrix) ), \
            [type(S).__name__, S.shape]
    return S.indices[ S.indptr[i] : S.indptr[i+1] ]
        # return S[i].nonzero()[1]  # .tocoo ?

outs = _ix  # outs( S, node n ): [j: n -> j]
ins = _ix   # ins( S, node n ): [i: i -> n]

def _toix( S: "Sparserow" ) -> Nodes:
    if isinstance( S, np.ndarray ):
        return S
    else:
        return _ix( S, 0 )  # csr / csc


def edges_to_sparsegraph( I, J: Nodes, n=0 ) -> Sparsegraph:
    """ I, J = edges.T of a mesh """
    if n == 0:
        n = max( I.max(), J.max() ) + 1
    S = sparse.lil_matrix( (n, n), dtype=np.int8 )
    S[ np.ix_( I, J )] = 1
    return S.tocsr()


#...............................................................................
if __name__ == "__main__":
    import sys
    from randomsparse import randomsparse  # ~ sparse.random

    np.set_printoptions( threshold=20, edgeitems=10, linewidth=120,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    print( 80 * "▄" )

    n = 2e5  # sparse n x n
    avrow = 200

        # to change these params, run this.py  a=1  b=None  'c = expr' ... in sh or ipython --
    for arg in sys.argv[1:]:
        exec( arg )

    print( "%s  n %g  avrow %d" % (
            sys.argv[0], n, avrow ))

        # randomsparse 0/1 matrix --
    n = int( n )
    S = randomsparse( n, n, density=avrow, seed=0,
                distrib = lambda n: np.ones( n, np.int ))
    S = sparse.triu( S, k=1 ).tocsr()  # edges i < j only
    print( "S randomsparse triu: %.0f mb  row sums %s ... " % (
            S.nnz * 24 / 1e6, S[:10].sum( axis=1 ).reshape(-1) ))

    #...........................................................................
    tri = triangles( S )
    # n 2e5  avrow 200: 480 mb  64244 triangles -> 189632 uniq edges  wall time 127 sec

triangles_to_edges.py

#!/usr/bin/env python3
""" triangles_to_edges: n x 3 -> uniq edges, tri->edgeno """

import numpy as np

#...............................................................................
def triangles_to_edges( tri: "n x 3", verbose=1 ) -> "uniq edges, tri->edgeno":
    # cf. create_edges in $fem/meshzoo-master/meshzoo/helpers.py
    assert tri.shape[1] == 3, tri.shape
    # [[a b c] ...] ->
    # [[a b]
    #  [b c]
    #  [c a] ...
    pairs = tri[:, [0,1, 1,2, 2,0]] .reshape( -1, 2  )

    uniqedges, tritoedges = np.unique( pairs, return_inverse=True, axis=0 )
    tritoedges = tritoedges.reshape( -1, 3 )
    if verbose:
        print( "triangles_to_edges: %d triangles -> %d uniq edges" % (
                len(tri), len(uniqedges) ))

    return uniqedges, tritoedges


#...............................................................................
if __name__ == "__main__":
    import sys
    import meshzoo

    np.set_printoptions( threshold=20, edgeitems=5, linewidth=120,
            formatter = dict( float = lambda x: "%.2g" % x ))  # float arrays %.2g
    print( 80 * "▄" )

    n = 2000
        # to change these params, run this.py  a=1  b=None  'c = expr' ... in sh or ipython --
    for arg in sys.argv[1:]:
        exec( arg )

    xy, tri = meshzoo.tube( n=n )  # >= 10
    print( "meshzoo.tube( n=%d ): %d triangles" % (
            n, len(tri) ))

    edges, tritoedges = triangles_to_edges( tri )
    # meshzoo.tube( n=2000 ): 1268000 triangles
    # triangles_to_edges: 1268000 triangles -> 3804000 uniq edges
    # Wall time:       6.80 s