Projecting Coordinates in Numpy array

So I have come across a rather large bottle neck in my software. I have a set of co-ordinates in cords where each row corresponds to X,Y,Z co-ordinates. Each co-ordinate in cords has a defined area in atom_proj . The atoms variable corresponds to the cords variable and provides the key to the atom_proj .

I project the co-ordinates onto the grid array then rotate and repeat until the number of rotations is satisfied. I only project the X and Z co-ordinates ignoring the Y.

I have simplified version of my code below. The code runs relatively quick for small co-ordinate sets and number of rotations. But can take a long time if both co-ordinate set and rotation list is large. The number of co-ordinates can vary from a few hundred to tens of thousands. I project the area on the grid over a number or rotations to produce a heat map. An example of the heat map for a co-ordinate set is also shown below.

Question:

(i) - How can I decrease the projection time of the co-ordinates onto the matrix

(ii) - Is there a more pythonic way of applying the co-ordinate area to the grid rather than array splicing?

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

Heat map: 在这里输入图像描述


Something like this should work for the inner loop

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

A couple notes:

  • There's still a loop, but it is through the length- 2 array of atom types, not the length- N array of individual atoms.
  • I've left out the for rot in [] loop, since it wasn't doing anything here, but it should fit back in just fine.
  • This works by putting a one at the central location of each atom, in a stack of grids. Then, for each atom type, those ones are all added. Then, as @Joe suggested, the atom projection is convolved with those centers.
  • For testing, my atoms is 1d, yours is 2d. Not sure if this was on purpose or not.
  • Below, I've also added a third version, which is your algorithm but with variables that I was able to understand, it's called OP_simplified
  • Here's the full suite:

    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
    

    It works!

    In [957]: np.allclose(OP(cords, atoms, atom_proj), project_atom(cords, atoms, atom_proj))
    Out[957]: True
    

    And timing:

    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/75386.html

    上一篇: 如何从侧面折叠bootstrap导航栏

    下一篇: 在Numpy数组中投影坐标