<p>This looks like the perfect work for cython. It it's great opp optimizing loops.<br>
Another option is the new<br>
Numba, an automatic compiler.</p>
<p>David.</p>
<div class="gmail_quote">El 26/12/2012 10:09, "Eric Emsellem" <<a href="mailto:eric.emsellem@eso.org">eric.emsellem@eso.org</a>> escribió:<br type="attribution"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Hi!<br>
<br>
I am looking for an efficient way of doing some simple binning of points<br>
and then applying some functions to points within each bin.<br>
<br>
I have tried several ways, including crude looping over the indices, or<br>
using digitize (see below) but I cannot manage to get it as efficient as<br>
I need it to be. I have a comparison with a similar (although complex)<br>
code in idl, and thought I would ask the forum. In idl there is a way to<br>
"invert" an histogram and get a reverse set of indices (via the<br>
histogram function) which seems to do the trick (or maybe it is faster<br>
for another reason).<br>
<br>
Below I provide a (dummy) example of what I wish to achieve. Any hint on<br>
how to do this EFFICIENTLY using numpy is most welcome. I need to speed<br>
things up quite a bit (at the moment, the version I have, see below, is<br>
10 times slower than the more complex idl routine..., I must be doing<br>
something wrong!).<br>
<br>
thanks!!<br>
Eric<br>
========================================================<br>
# I have a random set of data points in 2D with coordinates x,y :<br>
import numpy as np<br>
x = np.random.random(1000)<br>
y = np.random.random(1000)<br>
<br>
# I have now a 2D grid given by let's say 10x10 grid points:<br>
nx = 11<br>
ny = 21<br>
lx = linspace(0,1,nx)<br>
ly = linspace(0,1,ny)<br>
gx, gy = np.meshgrid(lx, ly)<br>
<br>
# So my set of 2D bins are (not needed in the solution I present but<br>
just for clarity)<br>
bins = np.dstack((gx.ravel(), gy.ravel()))[0]<br>
<br>
# Now I want to have the list of points in each bin and<br>
# if the number of points in that bin is larger than 10, apply (dummy)<br>
function func1 (see below)<br>
# If less than 10, apply (dummy) function func2 so (dum?)<br>
# if 0, do nothing<br>
# for two dummy functions like (for example):<br>
def func1(x) : return x.mean()<br>
<br>
def func2(x) : return x.std()<br>
<br>
# One solution would be to use digitize in 1D and histogram in 2D (don't<br>
need gx, gy for this one):<br>
<br>
h = histogram2d(x, y, bins=[lx, ly])[0]<br>
<br>
digitX = np.digitize(x, lx)<br>
digitY = np.digitize(y, ly)<br>
<br>
# create the output array, with -999 values to make sure I see which<br>
ones are not filled in<br>
result = np.zeros_like(h) - 999<br>
<br>
for i in range(nx-1) :<br>
     for j in range(ny-1) :<br>
         selectionofpoints = (digitX == i+1) & (digitY == j+1)<br>
         if h[i,j] > 10 : result[i,j] = func1(x[selectionofpoints])<br>
         elif h[i,j] > 0 : result[i,j] = func2(x[selectionofpoints])<br>
_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
</blockquote></div>