Mailman 3 July 2009 - NumPy-Discussion

saving complex vector
by Neal Becker 15 Jul '09

15 Jul '09

Simple question. I want to save a complex vector as text in format real_0 imag_0\n real_1 imag_1\n ... I thought to use something like: np.savetxt ('gen_qpsk.txt', (mod_out.real, mod_out.imag), fmt='%g %g\n') I need a way to reinterpret the complex data as an array with 2 columns to make this work. Ideas?

4 3

mean of two or more arrays while ignoring a specific value
by Greg Fiske 15 Jul '09

15 Jul '09

Dear list, I'm learning to work with numpy arrays. Can somebody explain how to get the average of two separate arrays while ignoring a user defined value in one array? For example: >>>a = numpy.array([1,5,4,99]) >>>b = numpy.array([3,7,2,8]) Ignoring the value 99, the result should be an array like c= ([2,6,3,8]) Thanks for any advice, Greg Greg Fiske Research Associate/GIS Manager The Woods Hole Research Center 149 Woods Hole Road Falmouth MA, 02540 508-540-9900x139 gfiske(a)whrc.org http://whrc.org

2 1

mean of two or more arrays while ignoring a specific value
by Greg Fiske 14 Jul '09

14 Jul '09

Dear list, I'm learning to work with numpy arrays. Can somebody explain how to get the average of two separate arrays while ignoring a user defined value in one array? For example: >>>a = numpy.array([1,5,4,99]) >>>b = numpy.array([3,7,2,8]) Ignoring the value 99, the result should be an array like c= ([2,6,3,8]) Thanks for any advice, Greg Greg Fiske Research Associate/GIS Manager The Woods Hole Research Center 149 Woods Hole Road Falmouth MA, 02540 508-540-9900x139 gfiske(a)whrc.org http://whrc.org

5 5

Integer Overflow?
by Dave 14 Jul '09

14 Jul '09

I got stung when taking an ordinary python integer to the power of a numpy integer - the result wasn't what I was expecting (see below)! Taking a wild guess I expect this is due to integer overflow (since it doesn't show up with int64). When working with an int32 type one has to be aware of such issues and I'm guessing the casting behaviour is a design decision, nevertheless I thought I'd post to make others aware of such considerations. HTH, Dave Python 2.5.4 (r254:67916, Dec 23 2008, 15:10:54) [MSC v.1310 32 bit (Intel)] Type "copyright", "credits" or "license" for more information. IPython 0.10.bzr.r1163 -- An enhanced Interactive Python. ? -> Introduction and overview of IPython's features. %quickref -> Quick reference. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [2]: import numpy as np In [3]: np.__version__ Out[3]: '1.4.0.dev7056' In [4]: 10**-10 Out[4]: 1e-010 In [5]: 10**np.float32(-10) Out[5]: 1e-010 In [6]: 10**np.int32(-10) Out[6]: 7.0918695992859933e-010 ???????????????????????????????????? In [7]: 10**np.int32(-1) Out[7]: 0.10000000000000001 In [8]: 10**np.int32(-2) Out[8]: 0.01 In [9]: 10**np.int32(-3) Out[9]: 0.001 In [10]: 10**np.int32(-4) Out[10]: 0.0001 In [11]: 10**np.int32(-5) Out[11]: 1.0000000000000001e-005 In [12]: 10**np.int32(-6) Out[12]: 9.9999999999999995e-007 In [13]: 10**np.int32(-7) Out[13]: 9.9999999999999995e-008 In [14]: 10**np.int32(-8) Out[14]: 1e-008 In [15]: 10**np.int32(-9) Out[15]: 1.0000000000000001e-009

2 1

Using interpolate with zero-rank array raises error
by Tony Yu 13 Jul '09

13 Jul '09

(Sorry if this is a duplicate; I think sent this from the wrong email the first time) When using interpolate with a zero-rank array, I get "ValueError: object of too small depth for desired array". The following code reproduces this issue >>> import numpy as np >>> x0 = np.array(0.1) >>> x = np.linspace(0, 1) >>> y = np.linspace(0, 1) >>> np.interp(x0, x, y) I hesitate to call this behavior a bug (when I've done this in the past, I find out I'm just doing something wrong), but I find the error unnecessary (and the error output a bit obscure). Below is a (poorly tested) fix, which seems to work for me. (Sorry I'm not on svn, so it isn't a proper diff) Cheers, -Tony Python 2.5.1 numpy 1.3.0 numpy/lib/function_base.py: line 1144 ============================= if isinstance(x, (float, int, number)): return compiled_interp([x], xp, fp, left, right).item() + elif isinstance(x, np.ndarray) and x.ndim == 0: + return compiled_interp(x[np.newaxis], xp, fp, left, right)[0] else: return compiled_interp(x, xp, fp, left, right) =============================

1 0

Optimizing reduction loops (sum(), prod(), et al.)
by Pauli Virtanen 13 Jul '09

13 Jul '09

Hi, Ticket #1143 points out that Numpy's reduction operations are not always cache friendly. I worked a bit on tuning them. Just to tickle some interest, a "pathological" case before optimization: In [1]: import numpy as np In [2]: x = np.zeros((80000, 256)) In [3]: %timeit x.sum(axis=0) 10 loops, best of 3: 850 ms per loop After optimization: In [1]: import numpy as np In [2]: x = np.zeros((80000, 256)) In [3]: %timeit x.sum(axis=0) 10 loops, best of 3: 78.5 ms per loop For comparison, a reduction operation on a contiguous array of the same size: In [4]: x = np.zeros((256, 80000)) In [5]: %timeit x.sum(axis=1) 10 loops, best of 3: 88.9 ms per loop Funnily enough, it's actually slower than the reduction over the axis with the larger stride. The improvement factor depends on the CPU and its cache size. The code is here: http://github.com/pv/numpy-work/ticket/1143-speedup-reduce A benchmark script numpy/core/src/bench-reduce.py is included. Just do ./bench-reduce.py ../../../build/lib.linux-i686-2.6 | tee output.txt to run it against the system Numpy, and ./bench-reduce.py plot < output.txt to plot the results. Here's one result set: http://www.iki.fi/pav/tmp/xeon.png http://www.iki.fi/pav/tmp/xeon.txt The y-axis is proportional to the number of elemental (addition) operations per second, x-axis the total size of the array. Vertical lines show how the result changed by the cache optimization. Factor of 2 improvement seems common in the test cases (typically, arrays that are "skinny" in some dimension). There are also some factors of 10, especially for arrays larger than the 6144 kb cache of this particular CPU. Also, for this CPU, there are no significant regressions. On an older CPU (slower, smaller cache), the situation is slightly different: http://www.iki.fi/pav/tmp/athlon.png http://www.iki.fi/pav/tmp/athlon.txt On average, it's still an improvement in many cases. However, now there are more regressions. The significant ones (factor of 1/2) are N-D arrays where the reduction runs over an axis with a small number of elements. The optimization is in evaluating the reduction operation with a blocked loop (pseudocode) for block_idx in loop_over_blocks_except_axis(arr): result[block_idx] = arr[block_idx,0] for j in xrange(1, N): result[block_idx] += arr[block_idx,j] } ufunc loop instead of the usual ufunc reduction loop for idx in loop_over_array_except_axis(arr): result[idx] = arr[idx,0] } for j in xrange(1, N): } ufunc loop result[idx] += arr[idx,j] } The blocks are chosen to minimize the stride of the inmost ufunc loop, maximizing reuse of elements fitting in a single cache line. This has to be balanced against avoiding inner loops over very small dimensions, which incur some overhead. For N-D arrays, I also attempt to lump several dimensions into the ufunc loop, if the data can be iterated over with a single stride. The decision heuristics are in ufunc_object.c:_construct_reduce, and the new outer loop NOBUFFER_UFUNCREDUCE_TRANSPOSE in ufunc_object.c:PyUFunc_Reduce. Unfortunately, improving the performance using the above scheme comes at the cost of some slightly murky heuristics. I didn't manage to come up with an optimal decision rule, so they are partly empirical. There is one parameter tuning the cross-over between minimizing stride and avoiding small dimensions. (This is more or less straightforward.) Another empirical decision is required in choosing whether to use the usual reduction loop, which is better in some cases, or the blocked loop. How to make this latter choice is not so clear to me. The present approach seems to mostly work, as far as the benchmark produces improvements and produces no severe regressions. However, I think the regressions for the older CPU point out that something still needs to be tuned. The main suspect is the latter decision rule, but I'm not sure what the correct solution here is. Comments? I'm not a cache optimization expert, so bright ideas about better decision rules and further optimizations would be appreciated. Also, I think this addition could be nice to have included in Numpy. The main value IMHO is in boosting the most cache-inefficient corner cases. I think many users do not realize that np.zeros((256, 80000)).sum(axis=1) is currently more efficient than np.zeros((80000, 256)).sum(axis=0). And as shown above, there is no real reason to have this performance difference... Cheers, Pauli

10 14

numpy for 2.6 on mac os x
by Tommy Grav 12 Jul '09

12 Jul '09

The current dmg on the numpy download pages is buildt against 2.5. Is there any plans to make one for 2.6 or do I have to compile from the source? Cheers Tommy

3 2

Small fix to nosetester
by Gael Varoquaux 11 Jul '09

11 Jul '09

Right now the logics to detect package name in nosetester is fairly simple, and fails for non-standard instal path. I have opened up a ticket for that, with a fix: http://projects.scipy.org/numpy/ticket/1168 Gaël

1 0

speed up np.diag
by Citi, Luca 11 Jul '09

11 Jul '09

Hello, I happened to have a look at the code for np.diag and found it more complicated that necessary. I think it can be rewritten more cleanly and efficiently. Appended you can find both versions. The speed improvement is significant: In [145]: x = S.rand(1000,1300) In [146]: assert all(diag(x,-13) == diag2(x,-13)) In [147]: %timeit diag(x,-13) 10000 loops, best of 3: 37.4 µs per loop In [148]: %timeit diag2(x,-13) 100000 loops, best of 3: 14.1 µs per loop In [149]: x = x.T In [150]: assert all(diag(x,-13) == diag2(x,-13)) In [151]: %timeit diag(x,-13) 10000 loops, best of 3: 81.1 µs per loop In [152]: %timeit diag2(x,-13) 100000 loops, best of 3: 14.8 µs per loop It takes slightly more than one third for the C-contiguous input and slightly more than one sixth for the F-contiguous input. Best, Luca ## ORIGINAL def diag(v, k=0): v = asarray(v) s = v.shape if len(s)==1: n = s[0]+abs(k) res = zeros((n,n), v.dtype) if (k>=0): i = arange(0,n-k) fi = i+k+i*n else: i = arange(0,n+k) fi = i+(i-k)*n res.flat[fi] = v return res elif len(s)==2: N1,N2 = s if k >= 0: M = min(N1,N2-k) i = arange(0,M) fi = i+k+i*N2 else: M = min(N1+k,N2) i = arange(0,M) fi = i + (i-k)*N2 return v.flat[fi] else: raise ValueError, "Input must be 1- or 2-d." ## SUGGESTED def diag(v, k=0): v = asarray(v) s = v.shape if len(s) == 1: n = s[0]+abs(k) res = zeros((n,n), v.dtype) i = k if k >= 0 else (-k) * s[1] res.flat[i::n+1] = v return res elif len(s) == 2: if v.flags.f_contiguous: v, k, s = v.T, -k, s[::-1] i = k if k >= 0 else (-k) * s[1] return v.flat[i::s[1]+1] else: raise ValueError, "Input must be 1- or 2-d."

3 4

argmin and argmax
by f0X_in_s0X 11 Jul '09

11 Jul '09

Do argmin and argmax skip nan values? I cant find it anywhere in the numpy documentation. In any case, if I want to find max and min values in an array that may contain nan values also, what would the be the most efficient way to do that? Any pointers?? -- View this message in context: http://www.nabble.com/argmin-and-argmax-tp24428341p24428341.html Sent from the Numpy-discussion mailing list archive at Nabble.com.

2 2