Mailman 3 November 2010 - NumPy-Discussion

broadcasting with numpy.interp
by greg whittier 02 Dec '10

02 Dec '10

Hi all, I'd like to be able to speed up the following code. def replace_dead(cube, dead): # cube.shape == (320, 640, 1200) # dead.shape == (320, 640) # cube[i,j,:] are bad points to be replaced via interpolation if dead[i,j] == True bands = np.arange(0, cube.shape[0]) for line in range(cube.shape[1]): dead_bands = bands[dead[:, line] == True] good_bands = bands[dead[:, line] == False] for sample in range(cube.shape[2]): # interp returns fp[0] for x < xp[0] and fp[-1] for x > xp[-1] cube[dead_bands, line, sample] = \ np.interp(dead_bands, good_bands, cube[good_bands, line, sample]) Removing that last loop via some sort of broadcasting seems like it should be possible, but that doesn't seem to work with interp. While interp allows the x-coordinates of interpolation points to be a multi-dimensional array, it expects the x- and y-coordinates of the data points to be 1-d arrays. Any suggestions for speeding this up? Thanks, Greg

2 3

A faster median (Wirth's method)
by Sturla Molden 01 Dec '10

01 Dec '10

We recently has a discussion regarding an optimization of NumPy's median to average O(n) complexity. After some searching, I found out there is a selection algorithm competitive in speed with Hoare's quick select. It has the advantage of being a lot simpler to implement. In plain Python: import numpy as np def wirthselect(array, k): """ Niklaus Wirth's selection algortithm """ a = np.ascontiguousarray(array) if (a is array): a = a.copy() l = 0 m = a.shape[0] - 1 while l < m: x = a[k] i = l j = m while 1: while a[i] < x: i += 1 while x < a[j]: j -= 1 if i <= j: tmp = a[i] a[i] = a[j] a[j] = tmp i += 1 j -= 1 if i > j: break if j < k: l = i if k < i: m = j return a Now, the median can be obtained in average O(n) time as: def median(x): """ median in average O(n) time """ n = x.shape[0] k = n >> 1 s = wirthselect(x, k) if n & 1: return s[k] else: return 0.5*(s[k]+s[:k].max()) The beauty of this is that Wirth select is extremely easy to migrate to Cython: import numpy ctypedef numpy.double_t T # or whatever def wirthselect(numpy.ndarray[T, ndim=1] array, int k): cdef int i, j, l, m cdef T x, tmp cdef T *a _array = np.ascontiguousarray(array) if (_array is array): _array = _array.copy() a = <T *> _array.data l = 0 m = <int> a.shape[0] - 1 with nogil: while l < m: x = a[k] i = l j = m while 1: while a[i] < x: i += 1 while x < a[j]: j -= 1 if i <= j: tmp = a[i] a[i] = a[j] a[j] = tmp i += 1 j -= 1 if i > j: break if j < k: l = i if k < i: m = j return _array For example, we could have a small script that generates withselect for all NumPy dtypes (T as template), and use a dict as jump table. Chad, you can continue to write quick select using NumPy's C quick sort in numpy/core/src/_sortmodule.c.src. When you are done, it might be about 10% faster than this. :-) Reference: http://ndevilla.free.fr/median/median.pdf Best regards, Sturla Molden

12 28

numpy speed question
by Jean-Luc Menut 01 Dec '10

01 Dec '10

Hello all, I have a little question about the speed of numpy vs IDL 7.0. I did a very simple little check by computing just a cosine in a loop. I was quite surprised to see an order of magnitude of difference between numpy and IDL, I would have thought that for such a basic function, the speed would be approximatively the same. I suppose that some of the difference may come from the default data type of 64bits in numpy and 32 bits in IDL. Is there a way to change the numpy default data type (without recompiling) ? And I'm not an expert at all, maybe there is a better explanation, like a better use of the several CPU core by IDL ? I'm working with windows 7 64 bits on a core i7. any hint is welcome. Thanks. Here the IDL code : Julian1 = SYSTIME( /JULIAN , /UTC ) for j=0,9999 do begin for i=0,999 do begin a=cos(2*!pi*i/100.) endfor endfor Julian2 = SYSTIME( /JULIAN , /UTC ) print, (Julian2-Julian1)*86400.0 print,cpt end result: % Compiled module: $MAIN$. 2.9999837 The python code: from numpy import * from time import time time1 = time() for j in range(10000): for i in range(1000): a=cos(2*pi*i/100.) time2 = time() print time2-time1 result: In [2]: run python_test_speed.py 24.1809999943

9 11

Warning: invalid value encountered in subtract
by Keith Goodman 30 Nov '10

30 Nov '10

After upgrading from numpy 1.4.1 to 1.5.1 I get warnings like "Warning: invalid value encountered in subtract" when I run unit tests (or timeit) using "python -c 'blah'" but not from an interactive session. How can I tell the warnings to go away?

4 5

NumPy 1.5.1 on RedHat 5.5
by David Brodbeck 30 Nov '10

30 Nov '10

I'm trying to install NumPy 1.5.1 on RedHat 5.5 and I'm having trouble getting it to find the ATLAS libraries. I did a lot of Googling but didn't find anything that helped...also looked through the install instructions, but they focus mainly on Ubuntu. The problem I'm having is it's looking for the libraries in the right location, but not finding them. e.g.: atlas_blas_info: libraries f77blas,cblas,atlas not found in /opt/python-2.5/lib libraries f77blas,cblas,atlas not found in /usr/local/lib64 libraries f77blas,cblas,atlas not found in /usr/local/lib libraries f77blas,cblas,atlas not found in /usr/lib64/atlas libraries f77blas,cblas,atlas not found in /usr/lib64/sse2 libraries f77blas,cblas,atlas not found in /usr/lib64 libraries f77blas,cblas,atlas not found in /usr/lib/sse2 libraries f77blas,cblas,atlas not found in /usr/lib NOT AVAILABLE ...and yet... brodbd@patas:~/numpy-1.5.1$ locate libf77blas /usr/lib64/atlas/libf77blas.so.3 /usr/lib64/atlas/libf77blas.so.3.0 brodbd@patas:~/numpy-1.5.1$ locate libcblas /usr/lib64/atlas/libcblas.so.3 /usr/lib64/atlas/libcblas.so.3.0 brodbd@patas:~/numpy-1.5.1$ locate libatlas /usr/lib64/atlas/libatlas.so.3 /usr/lib64/atlas/libatlas.so.3.0 So the libraries are there, and they're where NumPy is looking for them, but it's still not finding them? Clearly I'm missing something, but I'm not sure what. -- David Brodbeck System Administrator, Linguistics University of Washington

2 4

subdivide array
by Gerrit Holl 30 Nov '10

30 Nov '10

On 30 November 2010 17:58, Pierre GM <pgmdevlist(a)gmail.com> wrote: > On Nov 30, 2010, at 5:40 PM, John wrote: >> I have an array of data for a global grid at 1 degree resolution. It's >> filled with 1s and 0s, and it is just a land sea mask (not only, but >> as an example). I want to be able to regrid the data to higher or >> lower resolutions (i.e. 0.5 or 2 degrees). But if I try to use any >> standard interp functions, such as mpl_toolkits.basemap.interp it >> fails -- I assume due to the data being binary. >> >> I guess there may be a fairly easy routine to do this?? Does someone >> have an example? > > Just a random idea: have you tried to convert your input data to float? Hopefully you could get some values between 0 and 1 for your interpolated values, that you'll have to transform back to integers following a scheme of your choosing... I would argue that some float between 0 and 1 is an excellent representation when regridding a binary land-sea-mask onto a higher resolution. After all, this information is not there. Why should a land-sea mask be binary anyway? As if a grid-cell can only be fully ocean or fully land... BTW, I just realised Pythons convention that negative indices count from the end of the array is perfect when using a 180x360 land-sea-mask, as lon[-30] and lon[330] mean and should mean the same :) Gerrit. -- Gerrit Holl PhD student at Department of Space Science, Luleå University of Technology, Kiruna, Sweden http://www.sat.ltu.se/members/gerrit/

1 0

subdivide array
by John 30 Nov '10

30 Nov '10

Hello, I have an array of data for a global grid at 1 degree resolution. It's filled with 1s and 0s, and it is just a land sea mask (not only, but as an example). I want to be able to regrid the data to higher or lower resolutions (i.e. 0.5 or 2 degrees). But if I try to use any standard interp functions, such as mpl_toolkits.basemap.interp it fails -- I assume due to the data being binary. I guess there may be a fairly easy routine to do this?? Does someone have an example? Thanks!

3 2

Type of init_numpy()?
by Bruce Sherwood 28 Nov '10

28 Nov '10

For Python 2.x, init_numpy() was void on all platforms. For Python 3.1, I find experimentally that init_numpy() is int on Windows (but still void on Mac, and I think also void on Ubuntu Linux). Is this a bug? Bruce Sherwood

2 2

[ANN] Nanny, faster NaN functions
by Keith Goodman 28 Nov '10

28 Nov '10

===== Nanny ===== Nanny uses the magic of Cython to give you a faster, drop-in replacement for the NaN functions in NumPy and SciPy. For example:: >> import nanny as ny >> import numpy as np >> arr = np.random.rand(100, 100) >> timeit np.nansum(arr) 10000 loops, best of 3: 67.5 us per loop >> timeit ny.nansum(arr) 100000 loops, best of 3: 18.2 us per loop Let's not forget to add some NaNs:: >> arr[arr > 0.5] = np.nan >> timeit np.nansum(arr) 1000 loops, best of 3: 411 us per loop >> timeit ny.nansum(arr) 10000 loops, best of 3: 65 us per loop Nanny uses a separate Cython function for each combination of ndim, dtype, and axis. You can get rid of a lot of overhead (useful in an inner loop, e.g.) by directly importing the function that matches your problem:: >> arr = np.random.rand(10, 10) >> from nansum import nansum_2d_float64_axis1 >> timeit np.nansum(arr, axis=1) 10000 loops, best of 3: 25.5 us per loop >> timeit ny.nansum(arr, axis=1) 100000 loops, best of 3: 5.15 us per loop >> timeit nansum_2d_float64_axis1(arr) 1000000 loops, best of 3: 1.75 us per loop I put together Nanny as a way to learn Cython. It currently only supports: - functions: nansum - Operating systems: 64-bit (accumulator for int32 is hard coded to int64) - dtype: int32, int64, float64 - ndim: 1, 2, and 3 If there is interest in the project, I could continue adding the remaining NaN functions from NumPy and SciPy: nanmin, nanmax, nanmean, nanmedian (using a partial sort), nanstd. But why stop there? How about nancumsum or nanprod? Or anynan, which could short-circuit once a NaN is found? Feedback on the code or the direction of the project are welcomed. So is coding help---without that I doubt the package will ever be completed. Once nansum is complete, many of the remaining functions will be copy, paste, touch up operations. Remember, Nanny quickly protects your precious data from the corrupting influence of Mr. Nan. License ======= Nanny is distributed under a Simplified BSD license. Parts of NumPy, which has a BSD licenses, are included in Nanny. See the LICENSE file, which is distributed with Nanny, for details. Installation ============ You can grab Nanny at http://github.com/kwgoodman/nanny. nansum of ints is only supported by 64-bit operating systems at the moment. **GNU/Linux, Mac OS X, et al.** To install Nanny:: $ python setup.py build $ sudo python setup.py install Or, if you wish to specify where Nanny is installed, for example inside ``/usr/local``:: $ python setup.py build $ sudo python setup.py install --prefix=/usr/local **Windows** In order to compile the C code in Nanny you need a Windows version of the gcc compiler. MinGW (Minimalist GNU for Windows) contains gcc and has been used to successfully compile Nanny on Windows. Install MinGW and add it to your system path. Then install Nanny with the commands:: python setup.py build --compiler=mingw32 python setup.py install **Post install** After you have installed Nanny, run the suite of unit tests:: >>> import nanny >>> nanny.test() <snip> Ran 1 tests in 0.008s OK <nose.result.TextTestResult run=1 errors=0 failures=0>

13 63

merge_arrays is very slow; alternatives?
by Gerrit Holl 26 Nov '10

26 Nov '10

Hi, upon profiling my code, I found that numpy.lib.recfunctions.merge_arrays is extremely slow; it does some 7000 rows/second. This is not acceptable for me. I have two large record arrays, or arrays with a complicated dtype. All I want to do is to merge them into one. I don't think that should have to be a very slow operation, I don't need to copy anything, I just want to view the two record arrays as one. How can I do this in a faster way? In [45]: cProfile.runctx("numpy.lib.recfunctions.merge_arrays([metarows, targetrows2], flatten=True)", globals(), locals()) 225381902 function calls (150254635 primitive calls) in 166.620 CPU seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.031 0.031 166.620 166.620 <string>:1(<module>) 68/1 0.000 0.000 0.000 0.000 _internal.py:82(_array_descr) 2 0.000 0.000 0.000 0.000 numeric.py:286(asanyarray) 2 0.000 0.000 0.000 0.000 recfunctions.py:135(flatten_descr) 1 0.000 0.000 0.001 0.001 recfunctions.py:161(zip_descr) 149165600/74038400 117.195 0.000 139.701 0.000 recfunctions.py:235(_izip_fields_flat) 1088801 12.146 0.000 151.847 0.000 recfunctions.py:263(izip_records) 3 0.000 0.000 0.000 0.000 recfunctions.py:277(sentinel) 1 4.599 4.599 166.589 166.589 recfunctions.py:328(merge_arrays) 3 0.000 0.000 0.000 0.000 recfunctions.py:406(<genexpr>) 75127201 22.506 0.000 22.506 0.000 {isinstance} 69 0.000 0.000 0.000 0.000 {len} 1 0.000 0.000 0.000 0.000 {map} 1 0.000 0.000 0.000 0.000 {max} 2 0.000 0.000 0.000 0.000 {method '__array__' of 'numpy.ndarray' objects} 136 0.000 0.000 0.000 0.000 {method 'append' of 'list' objects} 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 2 0.000 0.000 0.000 0.000 {method 'extend' of 'list' objects} 2 0.000 0.000 0.000 0.000 {method 'pop' of 'list' objects} 2 0.000 0.000 0.000 0.000 {method 'ravel' of 'numpy.ndarray' objects} 2 0.000 0.000 0.000 0.000 {numpy.core.multiarray.array} 1 10.142 10.142 10.142 10.142 {numpy.core.multiarray.fromiter} Gerrit. -- Gerrit Holl PhD student at Department of Space Science, Luleå University of Technology, Kiruna, Sweden http://www.sat.ltu.se/members/gerrit/

2 2