有效的键列表到三角剖分转换？（Python）

Question

我需要将一个“bonds”（点索引对）数组转换为一个“triangles”（点索引的三元组，表示三角剖分）数组。 我的方法对于较大的（N~100000+）点来说太慢，因为对于~1000000点的网格，我需要多次这样做。你知道吗

我得到一个数组'BL'，其中每一行包含两个有连接点的索引。我的方法是制作一个数组'NL'，其中第i行包含第i个点的邻居的索引。然后，对于NL的每一行（比如第I行），我创建一个布尔数组，告诉我BL的哪些行包含索引，其中两个点都是第I个粒子的邻居。如果BL的第k行对于NL的第I行满足这个条件，那么在点[I，BL[k，0]，BL[k，1]]之间有一个三角形。你知道吗

我相信这一进程可以更有效率。有什么建议吗？ 我的函数“BL2TRI”如下：

def BL2TRI(BL,nn):
    '''Convert bond list (#bonds x 2) to Triangulation array (#tris x 3)

    Parameters
    ----------
    BL : array of dimension #bonds x 2
        Each row is a bond and contains indices of connected points
    nn : int
        max number of nearest neighbors expected in NL and KL

    Returns
    ----------
    TRI : array of dimension #tris x 3
        Each row contains indices of the 3 points lying at the vertices of the tri.
    '''
    print('BL2TRI: computing NL...')
    NL, KL = BL2NLandKL(BL,nn)
    print('BL2TRI: assembling TRI...')
    ind = 0 #index to keep track of where we are in TRI
    # make TRItmp longer than necessary (to fill it in without appending)
    TRItmp = np.zeros((10*len(NL),3),dtype='int')
    # add 1 to all BL values and to all NL values for which KL!=0 to avoid finding that zero is a neighbor of every particle with fewer than nn neighbors
    BLp = BL + np.ones(np.shape(BL))
    NLp = copy.deepcopy(NL)
    NLp[KL!=0] +=1

    for kk in range(len(NLp)):
        idx = np.logical_and( ismember(BLp[:,0], NLp[kk,:]), ismember(BLp[:,1], NLp[kk,:]) )
        TRIS = BL[idx,:]
        TRItmp[ind:ind+len(TRIS),:]  = np.hstack(( TRIS, kk*np.ones((len(TRIS),1)) ))
        ind = ind+len(TRIS)

    #trim off the extra zeros at the end
    TRI = TRItmp[0:ind,:]
    return TRI

要使用此函数，下面是一个简短的工作示例：

import numpy as np
import copy

def BL2NLandKL(BL,nn=6):
    '''Convert bond list (#bonds x 2) to neighbor list (#pts x max# neighbors) for lattice of bonded points. Also returns KL: ones where there is a bond and zero where there is not.
    '''
    NL = np.zeros((max(BL.ravel())+1,nn))
    KL = np.zeros((max(BL.ravel())+1,nn))
    for row in BL:
        col = np.where(KL[row[0],:]==0)[0][0]
        NL[row[0],col] = row[1]
        KL[row[0],col] = 1
        col = np.where(KL[row[1],:]==0)[0][0]
        NL[row[1],col] = row[0]
        KL[row[1],col] = 1        
    return NL, KL


def ismember(a, b):
    '''Return logical array (c) testing where elements of a are members of b. 
    The speed is O(len(a)+len(b)), so it's fast.
    '''
    bind = {}
    for i, elt in enumerate(b):
        if elt not in bind:
            bind[elt] = True
    return np.array([bind.get(itm, False) for itm in a])  

# make some points in 2D
pts = np.array([[-0.5,-0.5],[-0.5, 0.0],[-0.5,0.5],[0.0,-0.5],[0.0,0.0],[0.0,0.5],[0.5,-0.5],[0.5,0.0],[0.5,0.5]])
# Give the connectivity between the pts as BL
BL = np.array([[0, 1],[1, 2],[0, 3],[1, 3],[1, 4],[2, 4],[3, 4],[2, 5],[4, 5], [3, 6],[4, 6],[4, 7],[5, 7],[6, 7],[5, 8],[7, 8]], dtype='int32')
# Convert BL to triangulation (TRI)
TRI = BL2TRI(BL,8) 

请注意，结果有重复的行并且没有排序，但是这些是我在这里省略的简单步骤。你知道吗

Answer 1

根据我的时间安排，这里有一个更快的方法，根据我做的一些快速测试，它的伸缩性更好。也比较干净：

def BL2TRI(BL):
    d = {}
    tri = np.zeros((len(BL),3), dtype=np.int)
    c = 0
    for i in BL:
        if(i[0] > i[1]):
            t = i[0]
            i[0] = i[1]
            i[1] = t
        if(i[0] in d):
            d[i[0]].append(i[1])
        else:
            d[i[0]] = [i[1]]
    for key in d:
        for n in d[key]:
            for n2 in d[key]:
                if (n>n2) or n not in d:
                    continue
                if (n2 in d[n]):
                    tri[c,:] = [key,n,n2]
                    c += 1
    return tri[0:c]

在这里，我使用字典，这意味着我们受益于各种哈希表的加速，即使有大量的节点。我还减少了要检查的节点数，确保它们都是(a,b)其中a<b。你知道吗

值得注意的是，一般来说，对于许多涉及大型数组的问题，如果缺少numpy（以及相关的库-scitools等），我发现它通常更容易（而且通常稍微干净一些，根据您对一些更晦涩的numpy语法的喜爱程度，您可以通过诸如ctypes（https://docs.python.org/2/library/ctypes.html）之类的库将繁重的工作卸载到C中。这是很容易开始，所以如果你是灵活的代码，你需要做额外的计算，然后值得一看。你知道吗