Mailman 3 July 2011 - NumPy-Discussion

segfault on complex array on solaris x86
by John Hunter 17 Aug '11

17 Aug '11

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

numpy.interp running time
by Timo Kluck 16 Aug '11

16 Aug '11

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

inverting and calculating eigenvalues for many small matrices
by Daniel Wheeler 16 Aug '11

16 Aug '11

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

Error building numpy (1.5.1 and 1.6.1rc3) with python2.7 debug
by Sandro Tosi 06 Aug '11

06 Aug '11

Hello, while preparing a test upload for 1.6.1rc3 in Debian, I noticed that it gets an error when building blas with python 2.7 in the debug flavor, the build log is at [1]. It's also been confirmed it fails also with 1.5.1 [2] [1] http://people.debian.org/~morph/python-numpy_1.6.1~rc3-1_amd64.build [2] http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=634012 I think it might be a toolchain change in Debian (since 1.5.1 was built successfully and now it fails), but could you please give me a hand in debugging the issue? Thanks in advance, -- Sandro Tosi (aka morph, morpheus, matrixhasu) My website: http://matrixhasu.altervista.org/ Me at Debian: http://wiki.debian.org/SandroTosi

3 4

[ANN] IPython 0.11 is officially out
by Fernando Perez 03 Aug '11

03 Aug '11

Hi all, on behalf of the IPython development team, I'm thrilled to announce, after more than two years of development work, the official release of IPython 0.11. This release brings a long list of improvements and new features (along with hopefully few new bugs). We have completely refactored IPython, making it a much more friendly project to participate in by having better separated and organized internals. We hope you will not only use the new tools and libraries, but also join us with new ideas and development. After this very long development effort, we hope to make a few stabilization releases at a quicker pace, where we iron out the kinks in the new APIs and complete some remaining internal cleanup work. We will then make a (long awaited) IPython 1.0 release with these stable APIs. *Downloads* Download links and instructions are at: http://ipython.org/download.html And IPython is also on PyPI: http://pypi.python.org/pypi/ipython Those contain a built version of the HTML docs; if you want pure source downloads with no docs, those are available on github: Tarball: https://github.com/ipython/ipython/tarball/rel-0.11 Zipball: https://github.com/ipython/ipython/zipball/rel-0.11 * Features * Here is a quick listing of the major new features: - Standalone Qt console - High-level parallel computing with ZeroMQ - New model for GUI/plotting support in the terminal - A two-process architecture - Fully refactored internal project structure - Vim integration - Integration into Microsoft Visual Studio - Improved unicode support - Python 3 support - New profile model - SQLite storage for history - New configuration system - Pasting of code with prompts And many more... We closed over 500 tickets, merged over 200 pull requests, and more than 60 people contributed over 2200 commits for the final release. Please see our release notes for the full details on everything about this release: https://github.com/ipython/ipython/zipball/rel-0.11 As usual, if you find any problem, please file a ticket --or even better, a pull request fixing it-- on our github issues site (https://github.com/ipython/ipython/issues/). Many thanks to all who contributed! Fernando, on behalf of the IPython development team. http://ipython.org

1 1

numpy.sqrt behaving differently on MacOS Lion
by Ilan Schnell 02 Aug '11

02 Aug '11

MacOS Lion: >>> numpy.sqrt([complex(numpy.nan, numpy.inf)]) array([ nan+infj]) Other all system: array([ inf+infj]) This causes a few numpy tests to fail on Lion. The numpy was not compiled using the new LLVM based gcc, it is the same numpy binary I used on other MacOS systems, which was compiled using gcc-4.0.1. However on Lion it is linked to Lions LLVM based gcc runtime, which apparently has some different behavior when it comes to such strange complex values. - Ilan

4 8

recommendation for saving data
by Brian Blais 01 Aug '11

01 Aug '11

Hello, I was wondering if there are any recommendations for formats for saving scientific data. I am running a simulation, which has many somewhat-indepedent parts which have their own internal state and parameters. I've been using pickle (gzipped) to save the entire object (which contains subobjects, etc...), but it is getting too unwieldy and I think it is time to look for a more robust solution. Ideally I'd like to have something where I can call a save method on the simulation object, and it will call the save methods on all the children, on down the line all saving into one file. It'd also be nice if it were cross-platform, and I could depend on the files being readable into the future for a while. Are there any good standards for this? What do you use for saving scientific data? thank you, Brian Blais -- Brian Blais bblais(a)bryant.edu http://web.bryant.edu/~bblais http://bblais.blogspot.com/

4 3

Rationale for returning type-wrapped min() / max() scalars? (was: Problem with ufunc of a numpy.ndarray derived class)
by Hans Meine 31 Jul '11

31 Jul '11

Hi, I have the same problem as Martin DRUON, who wrote 10 days ago: > I have a problem with the ufunc return type of a numpy.ndarray derived > class. In fact, I subclass a numpy.ndarray using the tutorial : > http://docs.scipy.org/doc/numpy/user/basics.subclassing.html > > But, for example, if I execute the "max" ufunc from my subclass, the return > type differs from the return type of the numpy ufunc. BTW: http://projects.scipy.org/numpy/ticket/1904 is meant to describe this, although the example code misses the actual calls to a1.min() and a2.min() in the assertion: # --------------------------------------------------- import numpy class Test(numpy.ndarray): pass a1 = numpy.ndarray((1,)) a2 = Test((1,)) assert type(a1.min()) == type(a2.min()), \ "%s != %s" % (type(a1.min()), type(a2.min())) # --------------------------------------------------- This code fails with 1.6.0, while it worked in 1.3.0. I tend to think that this is a bug (after all, a1.min() does not return ndarray, but an array scalar), but maybe there is also a good reason for this (for us, unexpected) behavor change and a nice solution? Have a nice day, Hans

7 13

Quaternion dtype for NumPy - initial implementation available
by Martin Ling 31 Jul '11

31 Jul '11

Hi all, I have just pushed a package to GitHub which adds a quaternion dtype to NumPy: https://github.com/martinling/numpy_quaternion Some backstory: on Wednesday I gave a talk at SciPy 2011 about an inertial sensing simulation package I have been working on (http://www.imusim.org/). One component I suggested might be reusable from that code was the quaternion math implementation, written in Cython. One of its features is a wrapper class for Nx4 NumPy arrays that supports efficient operations using arrays of quaternion values. Travis Oliphant suggested that a quaternion dtype would be a better solution, and got me talking to Mark Weibe about this. With Mark's help I completed this initial version at yesterday's sprint session. Incidentally, how to do something like this isn't well documented and I would have had little hope without both Mark's in-person help and his previous code (for adding a half-precision float dtype) to refer to. I don't know what the consensus is about whether people writing custom dtypes is a desirable thing, but if it is then the process needs to be made a lot easier. That said, the fact this is doable without patching the numpy core at all is really, really nice. Example usage: >>> import numpy as np >>> import quaternion >>> np.quaternion(1,0,0,0) quaternion(1, 0, 0, 0) >>> q1 = np.quaternion(1,2,3,4) >>> q2 = np.quaternion(5,6,7,8) >>> q1 * q2 quaternion(-60, 12, 30, 24) >>> a = np.array([q1, q2]) >>> a array([quaternion(1, 2, 3, 4), quaternion(5, 6, 7, 8)], dtype=quaternion) >>> exp(a) array([quaternion(1.69392, -0.78956, -1.18434, -1.57912), quaternion(138.909, -25.6861, -29.9671, -34.2481)], dtype=quaternion) The following ufuncs are implemented: add, subtract, multiply, divide, log, exp, power, negative, conjugate, copysign, equal, not_equal, less, less_equal, isnan, isinf, isfinite, absolute Quaternion components are stored as doubles. The package could be extended to support e.g. qfloat, qdouble, qlongdouble Comparison operations follow the same lexicographic ordering as tuples. The unary tests isnan, isinf and isfinite return true if they would return true for any individual component. Real types may be cast to quaternions, giving quaternions with zero for all three imaginary components. Complex types may also be cast to quaternions, with their single imaginary component becoming the first imaginary component of the quaternion. Quaternions may not be cast to real or complex types. Comments very welcome. This is my first attempt at NumPy hacking :-) Martin

12 22

Error when building numpy with Py3k
by Dirk Ullrich 31 Jul '11

31 Jul '11

Hi, trying to build current Git master of numpy with Py3k (Pythin 3.2.1, to be precise) yields to an error: building 'numpy.lib._compiled_base' extension compiling C sources C compiler: gcc -pthread -DDYNAMIC_ANNOTATIONS_ENABLED=1 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -march=x86-64 -mtune=generic -O2 -pipe -fPIC creating build/temp.linux-x86_64-3.2/numpy/lib creating build/temp.linux-x86_64-3.2/numpy/lib/src compile options: '-Inumpy/core/include -Ibuild/src.linux-x86_64-3.2/numpy/core/include/numpy -Inumpy/core/src/private -Inumpy/core/src -Inumpy/core -Inumpy/core/src/npymath -Inumpy/core/src/multiarray -Inumpy/core/src/umath -Inumpy/core/src/npysort -Inumpy/core/include -I/usr/include/python3.2mu -Ibuild/src.linux-x86_64-3.2/numpy/core/src/multiarray -Ibuild/src.linux-x86_64-3.2/numpy/core/src/umath -c' gcc: numpy/lib/src/_compiled_base.c numpy/lib/src/_compiled_base.c: In function 'pack_or_unpack_bits': numpy/lib/src/_compiled_base.c:1317:51: error: 'PyArrayObject' has no member named 'ob_type' numpy/lib/src/_compiled_base.c:1357:43: error: 'PyArrayObject' has no member named 'ob_type' numpy/lib/src/_compiled_base.c: In function 'pack_or_unpack_bits': numpy/lib/src/_compiled_base.c:1317:51: error: 'PyArrayObject' has no member named 'ob_type' numpy/lib/src/_compiled_base.c:1357:43: error: 'PyArrayObject' has no member named 'ob_type' error: Command "gcc -pthread -DDYNAMIC_ANNOTATIONS_ENABLED=1 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -march=x86-64 -mtune=generic -O2 -pipe -fPIC -Inumpy/core/include -Ibuild/src.linux-x86_64-3.2/numpy/core/include/numpy -Inumpy/core/src/private -Inumpy/core/src -Inumpy/core -Inumpy/core/src/npymath -Inumpy/core/src/multiarray -Inumpy/core/src/umath -Inumpy/core/src/npysort -Inumpy/core/include -I/usr/include/python3.2mu -Ibuild/src.linux-x86_64-3.2/numpy/core/src/multiarray -Ibuild/src.linux-x86_64-3.2/numpy/core/src/umath -c numpy/lib/src/_compiled_base.c -o build/temp.linux-x86_64-3.2/numpy/lib/src/_compiled_base.o" failed with exit status 1 Dirk

1 0