Mailman 3 August 2011 - NumPy-Discussion

bug with assignment into an indexed array?
by Benjamin Root Aug. 18, 2011

Aug. 18, 2011

Came across this today when trying to determine what was wrong with my code: import numpy as np matched_to = np.array([-1] * 5) in_ellipse = np.array([False, True, True, True, False]) match = np.array([False, True, True]) matched_to[in_ellipse][match] = 3 I would expect matched_to to now be "array([-1, -1, 3, 3, -1])", but instead, it is still all -1. It would seem that unless the view was created by a slice, then the assignment into the indexed view would not work as expected. This works: >>> matched_to[:3][match] = 3 but not: >>> matched_to[np.array([0, 1, 2])][match] = 3 Note that the following does work: >>> matched_to[np.array([0, 1, 2])] = 3 Is this a bug, or was I wrong to expect this to work this way? Thanks, Ben Root

3 7

Best way to construct/slice 3-dimensional ndarray from multiple 2d ndarrays?
by Keith Hughitt Aug. 18, 2011

Aug. 18, 2011

Hi all, I have a method which builds a single 3d ndarray from several equal-dimension 2d ndarrays, and another method which extracts the original 2d ndarrays back out from the 3d one. The way I'm doing this right now is pretty simple, e.g.: cube = np.asarray([arr1, arr2,...]) ... x = cube[0] I believe the way this is currently handled, is to use new memory locations first for the 3d array, and then later for the 2d slices. Does anyone know if there is a better way to handle this? Ideally, I would like to reuse the same memory locations instead of copying it anew each time. Also, when subclassing ndarray and calling obj = data.view(cls) for an ndarray "data", does this copy the data into the new object by value or reference? The method which extracts the 2d slice actually returns a subclass of ndarray created using the extracted data, so this is why I ask. Any insight or suggestions would be appreciated. Thanks! Keith

3 4

segfault on complex array on solaris x86
by John Hunter Aug. 18, 2011

Aug. 18, 2011

johnh@udesktop253:~> gcc --version gcc (GCC) 3.4.3 (csl-sol210-3_4-branch+sol_rpath) Copyright (C) 2004 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. johnh@udesktop253:~> uname -a SunOS udesktop253 5.10 Generic_142910-17 i86pc i386 i86pc johnh@udesktop253:~> cat test.py import numpy as np print np.__version__ fs = 1000 t = np.linspace(0, 0.3, 301) A = np.array([2, 8]).reshape(-1, 1) f = np.array([150, 140]).reshape(-1, 1) xn = (A * np.exp(2j * np.pi * f * t)).sum(axis=0) johnh@udesktop253:~> python test.py 2.0.0.dev-9451260 Segmentation Fault (core dumped) johnh@udesktop253:~> johnh@udesktop253:~> sudo pstack /var/core/core.python.957 core '/var/core/core.python.957' of 9397: python test.py febf1928 cexp (0, 0, 0, 0, 8060ab0, 84321ac) + 1b0 fe9657e0 npy_cexp (80458e0, 0, 0, 0, 0, 84e2530) + 30 fe95064f nc_exp (8045920, 84e72a0, 8045978, 8045920, 10, 10) + 3f fe937d5b PyUFunc_D_D (84e2530, 84e20f4, 84e25b0, fe950610, 1, 0) + 5b fe95e818 PyUFunc_GenericFunction (81e96e0, 807deac, 0, 80460b8, 2, 2) + 448 fe95fb10 ufunc_generic_call (81e96e0, 807deac, 0, fe98a820) + 70 feeb2d78 PyObject_Call (81e96e0, 807deac, 0, 80a24ec, 8061c08, 0) + 28 fef11900 PyEval_EvalFrame (80a2394, 81645a0, 8079824, 8079824) + 146c fef17708 PyEval_EvalCodeEx (81645a0, 8079824, 8079824, 0, 0, 0) + 620 fef178af PyEval_EvalCode (81645a0, 8079824, 8079824, 8061488, fef3d9ee, 0) + 2f fef3d095 PyRun_FileExFlags (feb91c98, 804687b, 101, 8079824, 8079824, 1) + 75 fef3d9ee PyRun_SimpleFileExFlags (feb91c98, 804687b, 1, 80465a8, fef454a1, 804687b) + 172 fef3e4fd PyRun_AnyFileExFlags (feb91c98, 804687b, 1, 80465a8) + 61 fef454a1 Py_Main (1, 80466b8, feb1cf35, fea935a1, 29, feb96750) + 9d9 08050862 main (2, 80466b8, 80466c4) + 22 08050758 _start (2, 8046874, 804687b, 0, 8046883, 80468ad) + 60

2 3

Re: [Numpy-discussion] f2py - undefined symbol: _intel_fast_memset [SEC=UNCLASSIFIED]
by Jin Lee Aug. 17, 2011

Aug. 17, 2011

Hello, This is my very first attempt at using f2py but I have come across a problem. If anyone can assist me I would appreciate it very much. I have a very simple test Fortran source, sub.f90 which is: subroutine sub1(x,y) implicit none integer, intent(in) :: x integer, intent(out) :: y ! start y = x end subroutine sub1 I then used f2py to produce an object file, sub.so: f2py -c -m sub sub.f90 --fcompiler='gfortran' After starting a Python interactive session I tried to import the Fortran-derived Python module but I get an error message: >>> import sub Traceback (most recent call last): File "<stdin>", line 1, in <module> ImportError: ./sub.so: undefined symbol: _intel_fast_memset Can anyone suggest what this error message means and how I can overcome it, please? Regards, Jin

2 2

Trim a numpy array in numpy.
by Hongchun Jin Aug. 17, 2011

Aug. 17, 2011

*Hi there, * * * *I have a question regarding how to trim a string array in numpy. * * * *>>> import numpy as np* *>>> x = np.array(['aaa.hdf', 'bbb.hdf', 'ccc.hdf', 'ddd.hdf'])* * * *I expect to trim a certain part of each element in the array, for example '.hdf', giving me ['aaa', 'bbb', 'ccc', 'ddd']. Of course, I can do a loop thing. However, in my actual dataset, I have more than one million elements in such an array. So I am wondering is there a faster and better way to do it, like STRMID function in IDL? I try to google it, but it turns out that I can not find any discussion about it. Thanks. * * Hongchun*

3 4

[f2py] How to specify compile options in setup.py
by Jose Gomez-Dans Aug. 17, 2011

Aug. 17, 2011

Hi, Up to now, I have managed to build Fortran extensions with f2py by ussing the following command: $ python setup.py config_fc --fcompiler=gnu95 --f77flags='-fmy_flags' --f90flags='-fmy_flags' build I think that these options should be able to go in a setup.py file, and use the f2py_options file. One way of doing this is to extend sys.argv with the required command line options: import sys sys.argv.extend ( ['config_fc', '--fcompiler=gnu95', '--f77flags="-fmy_flags"', "--f90flags='-fmy_flags"] ) This works well if all the extensions require the same flags. In my case, however, One of the extensions requires a different set of flags (in particular, it requires that flag -fdefault-real-8 isn't set, which is required by the extensions). I tried setting the f2py_options in the add_extension method call: config.add_extension( 'my_extension', sources = my_sources, f2py_options=['f77flags="-ffixed-line-length-0" -fdefault-real-8', 'f90flags="-fdefault-real-8"'] ) This compiles the extensions (using the two dashes in front of the f2py option eg --f77flags results in an unrecognised option), but the f2p_options goes unheeded. Here's the relevant bit of the output from python setup.py build: compiling Fortran sources Fortran f77 compiler: /usr/bin/gfortran -ffixed-line-length-0 -fPIC -O3 -march=native Fortran f90 compiler: /usr/bin/gfortran -ffixed-line-length-0 -fPIC -O3 -march=native Fortran fix compiler: /usr/bin/gfortran -Wall -ffixed-form -fno-second-underscore -ffixed-line-length-0 -fPIC -O3 -march=native compile options: '-Ibuild/src.linux-i686-2.7 -I/usr/lib/pymodules/python2.7/numpy/core/include -I/usr/include/python2.7 -c' extra options: '-Jbuild/temp.linux-i686-2.7/my_dir -Ibuild/temp.linux-i686-2.7/my_dir' How can I disable (or enable) one option for compiling one particular extension? Thanks! Jose

2 1

Segfault for np.lookfor
by Matthew Brett Aug. 17, 2011

Aug. 17, 2011

Hi, On current trunk, all tests pass but running the (forgive my language) doctests, I found this: In [1]: import numpy as np In [2]: np.__version__ Out[2]: '2.0.0.dev-730b861' In [3]: np.lookfor('cos') Segmentation fault on: Linux angela 2.6.38-10-generic #46-Ubuntu SMP Tue Jun 28 15:07:17 UTC 2011 x86_64 x86_64 x86_64 GNU/Linux Ubuntu Natty Python 2.7.1+ Best, Matthew

2 5

numpy.interp running time
by Timo Kluck Aug. 17, 2011

Aug. 17, 2011

Dear numpy developers, The current implementation of numpy.interp(x,xp,fp) comes down to: first calculating all the slopes of the linear interpolant (these are len(xp)-1), then use a binary search to find where x is in xp (running time log(len(xp)). So we obtain a running time of O( len(xp) + len(x)*log(len(xp) ) We could improve this to just O( len(x)*log(len(xp) ) by not caching the slopes. The point is, of course, that this is slightly slower in the common use case where x is is refinement of xp, and where you will have to compute all the slopes anyway. In my personal use case, however, I needed the value of the interp(x0,xp,fp) in order to calculate the next point x1 where I wanted to calculate interp(x1,xp,fp). The current implementation gave a severe running time penalty. I have looked at the source and I could easily produce a patch for this. Would you be interested in it? Cheers, Timo Kluck

2 5

Questionable reduceat behavior
by Wes McKinney Aug. 17, 2011

Aug. 17, 2011

I'm a little perplexed why reduceat was made to behave like this: In [26]: arr = np.ones((10, 4), dtype=bool) In [27]: arr Out[27]: array([[ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True], [ True, True, True, True]], dtype=bool) In [30]: np.add.reduceat(arr, [0, 3, 3, 7, 9], axis=0) Out[30]: array([[3, 3, 3, 3], [1, 1, 1, 1], [4, 4, 4, 4], [2, 2, 2, 2], [1, 1, 1, 1]]) this does not seem intuitively correct. Since we have: In [33]: arr[3:3].sum(0) Out[33]: array([0, 0, 0, 0]) I would expect array([[3, 3, 3, 3], [0, 0, 0, 0], [4, 4, 4, 4], [2, 2, 2, 2], [1, 1, 1, 1]]) Obviously I can RTFM and see why it does this ("if ``indices[i] >= indices[i + 1]``, the i-th generalized "row" is simply ``a[indices[i]]``"), but it doesn't make much sense to me, and I need work around it. Suggestions?

2 3

inverting and calculating eigenvalues for many small matrices
by Daniel Wheeler Aug. 17, 2011

Aug. 17, 2011

Hi, I am trying to find the eigenvalues and eigenvectors as well as the inverse for a large number of small matrices. The matrix size (MxM) will typically range from 2x2 to 8x8 at most. The number of matrices (N) can be from 100 up to a million or more. My current solution is to define "eig" and "inv" to be, def inv(A): """ Inverts N MxM matrices, A.shape = (M, M, N), inv(A).shape = (M, M, N). """ return np.array(map(np.linalg.inv, A.transpose(2, 0, 1))).transpose(1, 2, 0) def eig(A): """ Calculate the eigenvalues and eigenvectors of N MxM matrices, A.shape = (M, M, N), eig(A)[0].shape = (M, N), eig(A)[1].shape = (M, M, N) """ tmp = zip(*map(np.linalg.eig, A.transpose(2, 0, 1))) return (np.array(tmp[0]).swapaxes(0,1), np.array(tmp[1]).transpose(1,2,0)) The above uses "map" to fake a vector solution, but this is heinously slow. Are there any better ways to do this without resorting to cython or weave (would it even be faster (or possible) to use "np.linalg.eig" and "np.linalg.inv" within cython)? I could write specialized versions when M=2 or M=3, which could be fully vectorized, but I'd rather have a general solution. Are there better algorithms than the ones used in "np.linalg.inv" and "np.linalg.eig" when M < 9, say, that I could hand code using numpy in a fully vectorized way? My end goal is to implement a Riemann flux calculation as part of a finite volume solver. This requires calculating "Abar = R E inv(R)" where R are the right eigenvectors of a matrix A and E is the matrix with the absolute value of the eigenvalues of A along the diagonal (both R and E are sorted in ascending order of their corresponding eigenvalues). As well as using "eig" and "inv" defined above, matrix multiplication and sorted eigenvalues and eigenvectors are also required. Fortunately, these can be vectorized as follows, def sum(a, axis=0): """ Faster than using np.sum. """ return np.tensordot(np.ones(a.shape[axis], 'l'), a, (0, axis)) def mul(A, B): """ Matrix multiply N MxM matrices, A.shape = B.shape = (M, M, N), mul(A, B).shape = (M, M, N) """ return sum(A.swapaxes(0,1)[:, :, np.newaxis] * B[:, np.newaxis], 0) def sortedeig(A): """ Calculates the sorted eigenvalues and eigenvectors of N MxM matrices, A.shape = (M, M, N), sortedeig(A)[0].shape = (M, N), sortedeig(A)[1].shape = (M, M, N). """ N = A.shape[-1] eigenvalues, R = eig(A) order = eigenvalues.argsort(0).swapaxes(0, 1) Nlist = [[i] for i in xrange(N)] return (eigenvalues[order, Nlist].swapaxes(0, 1), R[:, order, Nlist].swapaxes(1, 2)) The above two functions take very little time compared with "eig" and "inv". Using "sum" helps too. Given A, calculating "Abar = R E inv(R)" can then simply be written eigenvalues, R = sortedeig(A) E = abs(eigenvalues) * numerix.identity(eigenvalues.shape[0])[..., np.newaxis] Abar = mul(mul(R, E), inv(R)) Maybe there is a better way to calculate "Abar" rather than explicitly calculating the eigenvalues and eigenvectors for every matrix. Any help is much appreciated. -- Daniel Wheeler

6 12