I need to convert an array of 'bonds' (pairs of point indices) to an array of 'triangles' (triplets of point indices, representing a triangulation). My method is too slow for large (N~100000+) points, since I need to do this many times for meshes of ~1,000,000 pts.
I am given an array 'BL' in which each row contains indices of two points which have a connection. My approach is to make an array 'NL' in which the ith row contains indices of the neighbors of the ith point. Then for each row of NL (say the ith row), I make a boolean array that tells me which rows of BL contain indices in which both points are neighbors of the ith particle. If the kth row of BL satisfies this condition for the ith row of NL, then I have a triangle between the points [i, BL[k,0], BL[k,1]].
I am sure that this process can be made more efficient. Any suggestions? My function 'BL2TRI' is below:
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
To use this function, here is a short working example:
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)
Note that the result has repeated rows and is not sorted, but these are easy steps that I omit here.
Here's a method that's faster, according to my timing, and one which scales better according to some quick tests I've done. It's also a little cleaner:
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]
Here, I use dictionaries which means we benefit from all sorts of hash table speed up even with a large number of nodes. I've also cut down the number of nodes to check by making sure they are all (a,b) where a<b.
It's worth noting that in general, for many problems involving large arrays where numpy (and associated libraries - scitools etc) is lacking, I find it is frequently easier (and generally slightly cleaner, depending on how much you enjoy some of the more obscure numpy syntax) to offload the grunt work to C - through a library such as ctypes (https://docs.python.org/2/library/ctypes.html). This is pretty easy to get going so if you are flexible about what you code in and you need to do additional calculation then it's worth taking a look at.