Optimize speed of for loop using numpy
Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this:
tmp1=zeros((eta,xi),float)
tmp2=zeros((eta,xi),float)
tmp1=tmp1+10000
tmp2=tmp2+10000
for i in range(xi):
for j in range(eta):
for k in range(s):
if z_r[k,j,i] < depth:
if tmp1[j,i]==10000:
if (depth - z_r[k,j,i]) <= (depth - z_r[0,j,i]):
tmp1[j,i]=k
else:
if (depth - z_r[k,j,i]) <= (depth - z_r[int(tmp1[j,i]),j,i]):
tmp1[j,i]=k
elif z_r[k,j,i] >= depth:
if tmp2[j,i]==10000:
if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[s-1,j,i]) :
tmp2[j,i]=k
else:
if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[int(tmp2[j,i]),j,i]) :
tmp2[j,i]=k
Not very impressive. My problem is that I can not find any numpy functions that actually can do the tests I do for each k value in the arrays. I need to identify the position (i,j) of k-values that meet the specified requirement. There are way too many if-else tests here as well. The scripts takes about 10 seconds to run (for 238MB input file), but these arrays are read from netCDF files and can be much larger and easily grow to enormous dimensions. It is therefore crucial for me to speed things up. Hope some of you can help. I really appreciate all feedback on this. I am just a rooky to numpy.
Cheers and thanks for all help, Trond
On Mon, Feb 25, 2008 at 8:08 PM, Trond Kristiansen trond@unc.edu wrote:
Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this:
Can you try to repost this in such a way that the indentation is preserved? Feel free to attach it as a text file. Also, can you describe at a higher level what it is you are trying to accomplish and what the arrays mean?
Hi again.
I have attached the function that the FOR loop is part of as a python file. What I am trying to do is to create a set of functions that will read the output files (NetCDF) from running the ROMS model (ocean model). The output file is organized in xi (x-direction), eta (y-direction), and s (z-direction) where the s-values are vertical layers and not depth. This particular function (z_slice) will find the closest upper and lower s-layer for a given depth in meters (e.g. -100). Then values from the two selcted layers will be interpolated to create a new layer at the selected depth (-100). The problem is that the s-layers follow the bathymetry and a particular s-layer will therefore sometimes be above and sometimes below the selected depth that we want to interpolate to. That's why I need a quick script that searches all of the layers and find the upper and lower layers for a given depth value (which is negative). The z_r is a 3D array (s,eta,xi) that is created using the netcdf module.
The main goal of these set of functions is to move away from using matlab, but also to speed things up. The sliced data array will be plotted using GMT or pyNGL.
Thanks for helping me. Cheers, Trond
On 2/25/08 9:15 PM, "Robert Kern" robert.kern@gmail.com wrote:
On Mon, Feb 25, 2008 at 8:08 PM, Trond Kristiansen trond@unc.edu wrote:
Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this:
Can you try to repost this in such a way that the indentation is preserved? Feel free to attach it as a text file. Also, can you describe at a higher level what it is you are trying to accomplish and what the arrays mean?
I would definitely suggest using scipy's weave.inline for this. It seems like this particular function can be translated into C code really easily, which would give you a HUGE speed up. Look at some of the examples in scipy/weave/examples to see how to do this. The numpy book also has a section on it.
One of the reasons I've left matlab and never looked back is how easy it is to interweave bits of compiled C code for loops like this.
--Hoyt
On Mon, Feb 25, 2008 at 6:32 PM, Trond Kristiansen trond@unc.edu wrote:
Hi again.
I have attached the function that the FOR loop is part of as a python file. What I am trying to do is to create a set of functions that will read the output files (NetCDF) from running the ROMS model (ocean model). The output file is organized in xi (x-direction), eta (y-direction), and s (z-direction) where the s-values are vertical layers and not depth. This particular function (z_slice) will find the closest upper and lower s-layer for a given depth in meters (e.g. -100). Then values from the two selcted layers will be interpolated to create a new layer at the selected depth (-100). The problem is that the s-layers follow the bathymetry and a particular s-layer will therefore sometimes be above and sometimes below the selected depth that we want to interpolate to. That's why I need a quick script that searches all of the layers and find the upper and lower layers for a given depth value (which is negative). The z_r is a 3D array (s,eta,xi) that is created using the netcdf module.
The main goal of these set of functions is to move away from using matlab, but also to speed things up. The sliced data array will be plotted using GMT or pyNGL.
Thanks for helping me. Cheers, Trond
On 2/25/08 9:15 PM, "Robert Kern" robert.kern@gmail.com wrote:
On Mon, Feb 25, 2008 at 8:08 PM, Trond Kristiansen trond@unc.edu wrote:
Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this:
Can you try to repost this in such a way that the indentation is preserved? Feel free to attach it as a text file. Also, can you describe at a higher level what it is you are trying to accomplish and what the arrays mean?
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Mon, Feb 25, 2008 at 8:32 PM, Trond Kristiansen trond@unc.edu wrote:
Hi again.
I have attached the function that the FOR loop is part of as a python file. What I am trying to do is to create a set of functions that will read the output files (NetCDF) from running the ROMS model (ocean model). The output file is organized in xi (x-direction), eta (y-direction), and s (z-direction) where the s-values are vertical layers and not depth. This particular function (z_slice) will find the closest upper and lower s-layer for a given depth in meters (e.g. -100). Then values from the two selcted layers will be interpolated to create a new layer at the selected depth (-100). The problem is that the s-layers follow the bathymetry and a particular s-layer will therefore sometimes be above and sometimes below the selected depth that we want to interpolate to. That's why I need a quick script that searches all of the layers and find the upper and lower layers for a given depth value (which is negative). The z_r is a 3D array (s,eta,xi) that is created using the netcdf module.
Ah, that makes things clearer. You should be able to remove the innermost k-loop by using searchsorted().
On 25/02/2008, Trond Kristiansen trond@unc.edu wrote:
I have attached the function that the FOR loop is part of as a python file. What I am trying to do is to create a set of functions that will read the output files (NetCDF) from running the ROMS model (ocean model). The output file is organized in xi (x-direction), eta (y-direction), and s (z-direction) where the s-values are vertical layers and not depth. This particular function (z_slice) will find the closest upper and lower s-layer for a given depth in meters (e.g. -100). Then values from the two selcted layers will be interpolated to create a new layer at the selected depth (-100). The problem is that the s-layers follow the bathymetry and a particular s-layer will therefore sometimes be above and sometimes below the selected depth that we want to interpolate to. That's why I need a quick script that searches all of the layers and find the upper and lower layers for a given depth value (which is negative). The z_r is a 3D array (s,eta,xi) that is created using the netcdf module.
If I understand what you're doing correctly, you have a 3D array of depths, indexed by two unproblematic coordinates (eta and xi), and by a third coordinate whose values are peculiar. For each pair (eta, xi), you want to find the value of the third coordinate that occurs at a depth closest to some given depth.
Are you looking to do this for many values of depth? It seems like what you're trying to do is (approximately) invert a function - you have z as a function of s, and you want s as a function of z. Is the mapping between depths and s values monotonic?
First of all, numpy.searchsorted() is the right tool for finding where depth occurs in a list of z values. It gives you the position of the next lower value; it's pretty straightforward to write down a linear interpolation (you should be able to do it without any for loop at all) and more sophisticated interpolation schemes may be straightforward as well.
If you want a smoother interpolation scheme, you may want to look at scipy.interpolate. It is not always as vectorial as you might wish, but it implements splines (optionally with smoothing) quite efficiently. I believe scipy.ndimage may well have some suitable interpolation code as well.
In my opinion, the value of numpy/scipy comes from the tools that allow you to express major parts of your algorithm as a line or two of code. But it is often difficult to take advantage of this by "trying to accelerate for loops". I find it really pays to step back and look at my algorithm and think how to express it as uniform operations on arrays. (Of course it helps to have a rough idea of what operations are available; the numpy documentation is sometimes sadly lacking in terms of letting you know what tools are there.)
You might find useful the new http://www.scipy.org/Numpy_Functions_by_Category
Good luck, Anne
Trond,
See if the attached file contains something close to what you need. It has no loops at all; I have not timed it, but it should be quite quick. I have given it only a cursory check, so I don't guarantee it works correctly.
Depending on how your particular NetCDF interface works, you might need to copy arrays from NetCDF to ensure they are genuine ndarray objects.
For plotting, you might want to try matplotlib. I think you will find it easier to use than GMT, especially if you are accustomed to matlab.
Eric
Trond Kristiansen wrote:
Hi again.
I have attached the function that the FOR loop is part of as a python file. What I am trying to do is to create a set of functions that will read the output files (NetCDF) from running the ROMS model (ocean model). The output file is organized in xi (x-direction), eta (y-direction), and s (z-direction) where the s-values are vertical layers and not depth. This particular function (z_slice) will find the closest upper and lower s-layer for a given depth in meters (e.g. -100). Then values from the two selcted layers will be interpolated to create a new layer at the selected depth (-100). The problem is that the s-layers follow the bathymetry and a particular s-layer will therefore sometimes be above and sometimes below the selected depth that we want to interpolate to. That's why I need a quick script that searches all of the layers and find the upper and lower layers for a given depth value (which is negative). The z_r is a 3D array (s,eta,xi) that is created using the netcdf module.
The main goal of these set of functions is to move away from using matlab, but also to speed things up. The sliced data array will be plotted using GMT or pyNGL.
Thanks for helping me. Cheers, Trond
Hey all.
I would just like to thank you all for extremely good feedback on my problem with optimizing loops. Thank you all for being so helpful. Cheers, Trond
On Mon, Feb 25, 2008 at 7:08 PM, Trond Kristiansen trond@unc.edu wrote:
Hi all. This is my first email to the discussion group. I have spent two days trying to get a particular loop to speed up, and the best result I got was this:
tmp1=zeros((eta,xi),float)
tmp2=zeros((eta,xi),float)
tmp1=tmp1+10000
tmp2=tmp2+10000
for i in range(xi):
for j in range(eta):
for k in range(s):
if z_r[k,j,i] < depth:
if tmp1[j,i]==10000:
if (depth - z_r[k,j,i]) <= (depth - z_r[0,j,i]):
tmp1[j,i]=k
else:
if (depth - z_r[k,j,i]) <= (depth - z_r[int(tmp1[j,i]),j,i]):
tmp1[j,i]=k
elif z_r[k,j,i] >= depth:
if tmp2[j,i]==10000:
if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[s-1,j,i]) :
tmp2[j,i]=k
else:
if abs(depth - z_r[k,j,i]) <= abs(depth - z_r[int(tmp2[j,i]),j,i]) :
tmp2[j,i]=k
Not very impressive. My problem is that I can not find any numpy functions that actually can do the tests I do for each k value in the arrays. I need to identify the position (i,j) of k-values that meet the specified requirement. There are way too many if-else tests here as well. The scripts takes about 10 seconds to run (for 238MB input file), but these arrays are read from netCDF files and can be much larger and easily grow to enormous dimensions. It is therefore crucial for me to speed things up. Hope some of you can help. I really appreciate all feedback on this. I am just a rooky to numpy.
Python if statements are certain death.
Chuck
participants (6)
-
Anne Archibald -
Charles R Harris -
Eric Firing -
Hoyt Koepke -
Robert Kern -
Trond Kristiansen