Note that there are some explanatory texts on larger screens.

plurals
  1. POSpeed up array query in Numpy/Python
    text
    copied!<p>I have an array of points (called <em>points</em>), consisting of ~30000 x,y, and z values. I also have a separate array of points (called <em>vertices</em>), about ~40000 x,y, and z values. The latter array indexes the lower-front-left corners of some cubes of side length <em>size.</em> I would like to find out which points reside in which cubes, and how many points reside in each cube. I wrote a loop to do this, which works like this:</p> <pre><code>for i in xrange(len(vertices)): cube=((vertices[i,0]&lt;= points[:,0]) &amp; (points[:,0]&lt;(vertices[i,0]+size)) &amp; (vertices[i,1]&lt;= points[:,1]) &amp; (points[:,1] &lt; (vertices[i,1]+size)) &amp; (vertices[i,2]&lt;= points[:,2]) &amp; (points[:,2] &lt; (vertices[i,2]+size)) ) numpoints[i]=len(points[cube]) </code></pre> <p>(The loop is order the individual cubes, and the "cube" creates a boolean array of indices.) I then store points[cube] somewhere, but this isn't what's slowing me down; it's the creating of "cube=".</p> <p>I would like to speed this loop up (it takes tens of seconds to complete on a macbook pro). I tried rewriting the "cube=" part in C, as follows:</p> <pre><code>for i in xrange(len(vertices)): cube=zeros(pp, dtype=bool) code=""" for (int j=0; j&lt;pp; ++j){ if (vertices(i,0)&lt;= points(j,0)) if (points(j,0) &lt; (vertices(i,0)+size)) if (vertices(i,1)&lt;= points(j,1)) if (points(j,1) &lt; (vertices(i,1)+size)) if (vertices(i,2)&lt;= points(j,2)) if (points(j,2) &lt; (vertices(i,2)+size)) cube(j)=1; } return_val = 1;""" weave.inline(code, ['vertices', 'points','size','pp','cube', 'i']) numpoints[i]=len(points[cube]) </code></pre> <p>This sped it up by a bit more than a factor of two. Rewriting <em>both</em> loops in C actually made it only slightly faster than the original numpy-only version, because of the frequent references to the array objects necessary to keep track of which points are in which cubes. I suspect it's possible to do this much more rapidly, and that I'm missing something. Can anyone suggest how to speed this up?? I'm new to numpy/python, and thanks in advance.</p>
 

Querying!

 
Guidance

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload