在Numpy数组中投影坐标
所以我在我的软件中遇到了一个相当大的瓶颈。 我有一组坐标的cords
,其中每一行对应于X,Y,Z
坐标。 每个坐标的cords
在限定区域atom_proj
。 atoms
变量对应于cords
变量,并将关键字提供给atom_proj
。
我将坐标投影到grid
阵列上,然后旋转并重复,直到满足旋转次数。 我只投影X和Z坐标忽略Y.
我在下面简化了我的代码版本。 代码运行相对较快,适用于小坐标集和旋转次数。 但是如果坐标集和旋转列表都很大,可能需要很长时间。 坐标的数量可以从几百到几万不等。 我通过一些数字或旋转来投影grid
上的区域以产生热图。 下面还显示了一个坐标集的热图的例子。
题:
(i) - 如何减少坐标到矩阵上的投影时间
(ii) - 将坐标区域应用到grid
而不是阵列拼接有更多的pythonic方法吗?
import numpy as np
cords = np.array([[5,4,5],[5,4,3],[6,4,6]])
atoms = np.array([['C'],['H'],['C']])
atom_proj = {'H':np.array([[0,0,0,0,0],[0,0,1,0,0],[0,1,1,1,0],[0,0,1,0,0],[0,0,0,0,0]]),'C':np.array([[0,0,0,0,0,0,0],[0,0,0,0,0,0,0],[0,0,1,1,1,0,0],[0,0,1,1,1,0,0],[0,0,1,1,1,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,0,0]])}
grid = np.zeros((10,10))
for rot in xrange(1,10):
# This for loop would contain a list of list of rotations to apply which are calculated before hand.
# apply rotation
for values in zip(cords, atoms):
atom_shape = np.shape(atom_proj[values[1][0]])
rad = (atom_shape[0]-1)/2
grid[values[0][2]-rad:values[0][2]+rad+1,values[0][0]-rad:values[0][0]+rad+1] += atom_proj[values[1][0]]
print grid
热图:
像这样的东西应该适用于内部循环
extruded = np.zeros((N, 10,10))
extruded[range(N), cords[:,2], cords[:,0]] = 1
grid = np.zeros((10,10))
for atom, proj in atom_proj.iteritems():
centers = extruded[atoms==atom].sum(0)
projected = nd.convolve(centers, proj)
grid += projected
一对夫妇笔记:
2
的原子类型阵列,而不是单个原子的长度为N
阵列。 for rot in []
循环中忽略了for rot in []
,因为它在这里没有做任何事情,但它应该适得其反。 atoms
是1d,你的是2d。 不确定这是否有意或无意。 OP_simplified
这是完整的套件:
import numpy as np
import scipy.ndimage as nd
N = 1000
cords = np.random.randint(3, 7, (N, 3)) #np.array([[5,4,5],[5,4,3],[6,4,6]])
atoms = np.random.choice(list('HC'), N) #np.array([['C'],['H'],['C']])
atom_proj = {'H': np.array([[0,0,0,0,0],
[0,0,1,0,0],
[0,1,1,1,0],
[0,0,1,0,0],
[0,0,0,0,0]]),
'C': np.array([[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0],
[0,0,1,1,1,0,0],
[0,0,1,1,1,0,0],
[0,0,1,1,1,0,0],
[0,0,0,0,0,0,0],
[0,0,0,0,0,0,0]])}
def project_atom(cords, atoms, atom_proj):
extruded = np.zeros((N, 10,10))
extruded[range(N), cords[:,2], cords[:,0]] = 1
grid = np.zeros((10,10))
for atom, proj in atom_proj.iteritems():
grid += nd.convolve(extruded[atoms.squeeze()==atom].sum(0), proj, mode='constant')
return grid
def OP_simplified(cords, atoms, atom_proj):
rads = {atom: (proj.shape[0] - 1)/2 for atom, proj in atom_proj.iteritems()}
grid = np.zeros((10,10))
for (x,y,z), atom in zip(cords, atoms):
rad = rads[atom]
grid[z-rad:z+rad+1, x-rad:x+rad+1] += atom_proj[atom]
return grid
def OP(cords, atoms, atom_proj):
grid = np.zeros((10,10))
for values in zip(cords, atoms):
atom_shape = np.shape(atom_proj[values[1][0]])
rad = (atom_shape[0]-1)/2
grid[values[0][2]-rad:values[0][2]+rad+1,values[0][0]-rad:values[0][0]+rad+1] += atom_proj[values[1][0]]
return grid
有用!
In [957]: np.allclose(OP(cords, atoms, atom_proj), project_atom(cords, atoms, atom_proj))
Out[957]: True
和时机:
In [907]: N = 1000
In [910]: timeit OP(cords, atoms, atom_proj)
10 loops, best of 3: 30.7 ms per loop
In [911]: timeit project_atom(cords, atoms, atom_proj)
100 loops, best of 3: 2.97 ms per loop
In [913]: N = 10000
In [916]: timeit project_atom(cords, atoms, atom_proj)
10 loops, best of 3: 33.3 ms per loop
In [917]: timeit OP(cords, atoms, atom_proj)
1 loops, best of 3: 314 ms per loop
链接地址: http://www.djcxy.com/p/75385.html