高效的债券清单到三角转换？ (python)

Question

我需要将 'bonds' 的数组（点索引对）转换为 'triangles' 的数组（点索引的三元组，表示三角剖分）。我的方法对于大 (N~100000+) 个点来说太慢了，因为我需要对 ~1,000,000 个点的网格执行多次此操作。

我得到了一个数组 '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]

在这里，我使用字典，这意味着我们受益于各种哈希 table 即使有大量节点也能加速。我还通过确保它们都是 (a,b) where a<b.

来减少要检查的节点数

值得注意的是，一般来说，对于许多涉及缺少 numpy（和相关库 - scitools 等）的大型数组的问题，我发现它通常更容易（并且通常稍微更干净，具体取决于你对某些东西的喜爱程度）更晦涩的 numpy 语法）将繁重的工作卸载到 C - 通过诸如 ctypes (https://docs.python.org/2/library/ctypes.html) 的库。这很容易上手，所以如果您对编码的内容很灵活并且需要进行额外的计算，那么值得一看。

高效的债券清单到三角转换？ (python)

Efficient bond list to triangulation conversion? (python)

python

arrays

triangulation