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数组中投影坐标