Mailman 3 October 2006 - NumPy-Discussion

PyArray_ContiguousFromObject & PyDECREF
by Michiel Jan Laurens de Hoon 16 May '11

16 May '11

Dear pythoneers, In Chapter 12 ("Writing a C extension") in the Numerical Python manual, in the first example PyArray_ContiguousFromObject is used together with PyDECREF, while in the last example PyArray_ContiguousFromObject is used but the arrays are not PyDECREF'ed before the function returns. From Chapter 13 ("C API Reference") it seems that a PyDECREF is needed in both cases. Or am I missing something here? Is PyDECREF not needed for some reason in the last example? Thanks in advance, Michiel de Hoon U Tokyo First example: ============== static PyObject * trace(PyObject *self, PyObject *args) { PyObject *input; PyArrayObject *array; double sum; int i, n; if (!PyArg_ParseTuple(args, "O", &input)) return NULL; array = (PyArrayObject *) PyArray_ContiguousFromObject(input, PyArray_DOUBLE, 2, 2); if (array == NULL) return NULL; n = array->dimensions[0]; if (n > array->dimensions[1]) n = array->dimensions[1]; sum = 0.; for (i = 0; i < n; i++) sum += *(double *)(array->data + i*array->strides[0] + i*array->strides[1]); Py_DECREF(array); return PyFloat_FromDouble(sum); } Last example: ============= static PyObject * matrix_vector(PyObject *self, PyObject *args) { PyObject *input1, *input2; PyArrayObject *matrix, *vector, *result; int dimensions[1]; double factor[1]; double real_zero[1] = {0.}; long int_one[1] = {1}; long dim0[1], dim1[1]; extern dgemv_(char *trans, long *m, long *n, double *alpha, double *a, long *lda, double *x, long *incx, double *beta, double *Y, long *incy); if (!PyArg_ParseTuple(args, "dOO", factor, &input1, &input2)) return NULL; matrix = (PyArrayObject *) PyArray_ContiguousFromObject(input1, PyArray_DOUBLE, 2, 2); if (matrix == NULL) return NULL; vector = (PyArrayObject *) PyArray_ContiguousFromObject(input2, PyArray_DOUBLE, 1, 1); if (vector == NULL) return NULL; if (matrix->dimensions[1] != vector->dimensions[0]) { PyErr_SetString(PyExc_ValueError, "array dimensions are not compatible"); return NULL; } dimensions[0] = matrix->dimensions[0]; result = (PyArrayObject *)PyArray_FromDims(1, dimensions, PyArray_DOUBLE); if (result == NULL) return NULL; dim0[0] = (long)matrix->dimensions[0]; dim1[0] = (long)matrix->dimensions[1]; dgemv_("T", dim1, dim0, factor, (double *)matrix->data, dim1, (double *)vector->data, int_one, real_zero, (double *)result->data, int_one); return PyArray_Return(result); } -- Michiel de Hoon, Assistant Professor University of Tokyo, Institute of Medical Science Human Genome Center 4-6-1 Shirokane-dai, Minato-ku Tokyo 108-8639 Japan http://bonsai.ims.u-tokyo.ac.jp/~mdehoon

3 2

Re: [Numpy-discussion] Floating point exception with numpy and embedded python interpreter
by Andrew Straw 20 Sep '07

20 Sep '07

Arkaitz Bitorika wrote: > > On 19 Apr 2006, at 18:37, Andrew Straw wrote: > >> >>> I've been embedding Python for ages on the same program and other >>> modules work fine, only numpy fails. >> >> >> Most other modules don't use the SSE units, so wouldn't get hit by such >> a bug. > > > Is there a way of not using those units from numpy, to check if > that's what's going on? I think that numpy only accesses the SSE units through ATLAS or other external library. So, build numpy without ATLAS. But I'm not 100% sure anymore if there aren't any optimizations that directly use SSE if it's available. > Or alternatively, how would I check if my program is messing with the > SSE bits? Hmm, I think that's a bit hairy. I'd suggest simply asking the C++ library's mailing list if they alter the error bits on the control registers of the SSE unit. (Out of curiousity, what library is it?) If you want hairy, though, I think you'd have to check from C with the appropriate calls -- I'd start with the source code in that bug report. It looks like they're inlining an assembly statement to query a SSE control register.

5 7

Numpy x Matlab: some synthetic benchmarks
by Paulo Jose da Silva e Silva 11 Feb '07

11 Feb '07

Hello, Travis asked me to benchmark numpy versus matlab in some basic linear algebra operations. Here are the resuts for matrices/vectors of dimensions 5, 50 and 500: Operation x'*y x*y' A*x A*B A'*x Half 2in2 Dimension 5 Array 0.94 0.7 0.22 0.28 1.12 0.98 1.1 Matrix 7.06 1.57 0.66 0.79 1.6 3.11 4.56 Matlab 1.88 0.44 0.41 0.35 0.37 1.2 0.98 Dimension 50 Array 9.74 3.09 0.56 18.12 13.93 4.2 4.33 Matrix 81.99 3.81 1.04 19.13 14.58 6.3 7.88 Matlab 16.98 1.94 1.07 17.86 0.73 1.57 1.77 Dimension 500 Array 1.2 8.97 2.03 166.59 20.34 3.99 4.31 Matrix 17.95 9.09 2.07 166.62 20.67 4.11 4.45 Matlab 2.09 6.07 2.17 169.45 2.1 2.56 3.06 Obs: The operation Half is actually A*x using only the lower half of the matrix and vector. The operation 2in2 is A*x using only the even indexes. Of course there are many repetitions of the same operation: 100000 for dim 5 and 50 and 1000 for dim 500. The inner product is number of repetitions is multiplied by dimension (it is very fast). The software is numpy svn version 1926 Matlab 6.5.0.180913a Release 13 (Jun 18 2002) Both softwares are using the *same* BLAS and LAPACK (ATLAS for sse). As you can see, numpy array looks very competitive. The matrix class in numpy has too much overhead for small dimension though. This overhead is very small for medium size arrays. Looking at the results above (specially the small dimensions ones, for higher dimensions the main computations are being performed by the same BLAS) I believe we can say: 1) Numpy array is faster on usual operations but outerproduct (I believe the reason is that the dot function uses the regular matrix multiplication to compute outer-products, instead of using a special function. This can "easily" changes). In particular numpy was faster in matrix times vector operations, which is the most usual in numerical linear algebra. 2) Any operation that involves transpose suffers a very big penalty in numpy. Compare A'*x and A*x, it is 10 times slower. In contrast Matlab deals with transpose quite well. Travis is already aware of this and it can be probably solved. 3) When using subarrays, numpy is a slower. The difference seems acceptable. Travis, can this be improved? Best, Paulo Obs: Latter on (in a couple of days) I may present less synthetic benchmarks (a QR factorization and a Modified Cholesky).

13 27

A reimplementation of MaskedArray
by Pierre GM 23 Nov '06

23 Nov '06

Folks, I just posted on the scipy/developers zone wiki (http://projects.scipy.org/scipy/numpy/wiki/MaskedArray) a reimplementation of the masked_array mopdule, motivated by some problems I ran into while subclassing MaskedArray. The main differences with the initial numpy.core.ma package are that MaskedArray is now a subclass of ndarray and that the _data section can now be any subclass of ndarray (well, it should work in most cases, some tweaking might required here and there). Apart from a couple of issues listed below, the behavior of the new MaskedArray class reproduces the old one. It is quite likely to be significantly slower, though: I was more interested in a clear organization than in performance, so I tended to use wrappers liberally. I'm sure we can improve that rather easily. The new module, along with a test suite and some utilities, are available here: http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/maskedarr… http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/masked_te… http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/test_mask… Please note that it's still a work in progress (even if it seems to work quite OK when I use it). Suggestions, comments, improvements and general feedback are more than welcome ! ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

8 27

A way to compute a 2d histogram
by Xavier Gnata 14 Nov '06

14 Nov '06

Hi, I have a set of 3 1D large arrays. The first 2 one stand for the coordinates of particules and the last one for their masses. I would like to be able to plot this data ie to compute a 2D histogram summing the masses in each bin. I cannot find a way to do that without any loop on the indices resulting too a very slow function. I'm looking for an elegant way to do that with numpy (or scipy??) function. For instance, scipy.histogram2d cannot do the job because it only counts the number of samples in each bin. There is no way to deal with weights. Xavier. -- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata(a)obs.univ-lyon1.fr ############################################ ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

3 9

Helper function to "unroll" a array
by Gael Varoquaux 12 Nov '06

12 Nov '06

Hi, There is an operation I do a lot, I would call it "unrolling" a array. The best way to describe it is probably to give the code: def unroll(M): """ Flattens the array M and returns a 2D array with the first columns being the indices of M, and the last column the flatten M. """ return hstack((indices(M.shape).reshape(-1,M.ndim),M.reshape(-1,1))) Example: >>> M array([[ 0.73530097, 0.3553424 , 0.3719772 ], [ 0.83353373, 0.74622133, 0.14748905], [ 0.72023762, 0.32306969, 0.19142366]]) >>> unroll(M) array([[ 0. , 0. , 0.73530097], [ 0. , 1. , 0.3553424 ], [ 1. , 1. , 0.3719772 ], [ 2. , 2. , 0.83353373], [ 2. , 0. , 0.74622133], [ 1. , 2. , 0.14748905], [ 0. , 1. , 0.72023762], [ 2. , 0. , 0.32306969], [ 1. , 2. , 0.19142366]]) The docstring sucks. The function is trivial (when you know numpy a bit). Maybe this function already exists in numpy, if so I couldn't find it. Elsewhere I propose it for inclusion. Cheers, Gaël ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

2 3

Strange and hard to reproduce crash
by Fernando Perez 06 Nov '06

06 Nov '06

Hi all, two colleagues have been seeing occasional crashes from very long-running code which uses numpy. We've now gotten a backtrace from one such crash, unfortunately it uses a build from a few days ago: In [3]: numpy.__version__ Out[3]: '1.0b5.dev3097' In [4]: scipy.__version__ Out[4]: '0.5.0.2180' Because it takes so long to get the code to crash (several days of 100%CPU usage), I can't make a new one right now, but I'll be happy to restart the same run with a current SVN build if necessary, and post the results in a few days. In the meantime, here's a gdb backtrace we were able to get by setting MALLOC_CHECK_ to 2 and running the python process from within gdb: Program received signal SIGABRT, Aborted. [Switching to Thread 1073880896 (LWP 26280)] 0x40000402 in __kernel_vsyscall () (gdb) bt #0 0x40000402 in __kernel_vsyscall () #1 0x0042c7d5 in raise () from /lib/tls/libc.so.6 #2 0x0042e149 in abort () from /lib/tls/libc.so.6 #3 0x0046b665 in free_check () from /lib/tls/libc.so.6 #4 0x00466e65 in free () from /lib/tls/libc.so.6 #5 0x005a4ab7 in PyObject_Free () from /usr/lib/libpython2.3.so.1.0 #6 0x403f6336 in arraydescr_dealloc (self=0x40424020) at arrayobject.c:10455 #7 0x403fab3e in PyArray_FromArray (arr=0xe081cb0, newtype=0x40424020, flags=0) at arrayobject.c:7725 #8 0x403facc3 in PyArray_FromAny (op=0xe081cb0, newtype=0x0, min_depth=0, max_depth=0, flags=0, context=0x0) at arrayobject.c:8178 #9 0x4043bc45 in PyUFunc_GenericFunction (self=0x943a660, args=0xa9dbf2c, mps=0xbfc83730) at ufuncobject.c:906 #10 0x40440a04 in ufunc_generic_call (self=0x943a660, args=0xa9dbf2c) at ufuncobject.c:2742 #11 0x0057d607 in PyObject_Call () from /usr/lib/libpython2.3.so.1.0 #12 0x0057d6d4 in PyObject_CallFunction () from /usr/lib/libpython2.3.so.1.0 #13 0x403eabb6 in PyArray_GenericBinaryFunction (m1=Variable "m1" is not available. ) at arrayobject.c:3296 #14 0x0057b7e1 in PyNumber_Check () from /usr/lib/libpython2.3.so.1.0 #15 0x0057c1e0 in PyNumber_Multiply () from /usr/lib/libpython2.3.so.1.0 #16 0x005d16a3 in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #17 0x005d509e in PyEval_EvalCodeEx () from /usr/lib/libpython2.3.so.1.0 #18 0x005d3d8f in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #19 0x005d509e in PyEval_EvalCodeEx () from /usr/lib/libpython2.3.so.1.0 #20 0x00590e2e in PyFunction_SetClosure () from /usr/lib/libpython2.3.so.1.0 #21 0x0057d607 in PyObject_Call () from /usr/lib/libpython2.3.so.1.0 #22 0x00584d98 in PyMethod_New () from /usr/lib/libpython2.3.so.1.0 #23 0x0057d607 in PyObject_Call () from /usr/lib/libpython2.3.so.1.0 #24 0x005b584c in _PyObject_SlotCompare () from /usr/lib/libpython2.3.so.1.0 #25 0x005aec2c in PyType_IsSubtype () from /usr/lib/libpython2.3.so.1.0 #26 0x0057d607 in PyObject_Call () from /usr/lib/libpython2.3.so.1.0 #27 0x005d2b7f in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #28 0x005d509e in PyEval_EvalCodeEx () from /usr/lib/libpython2.3.so.1.0 #29 0x005d3d8f in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #30 0x005d509e in PyEval_EvalCodeEx () from /usr/lib/libpython2.3.so.1.0 #31 0x005d3d8f in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #32 0x005d497b in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #33 0x005d497b in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #34 0x005d497b in _PyEval_SliceIndex () from /usr/lib/libpython2.3.so.1.0 #35 0x005d509e in PyEval_EvalCodeEx () from /usr/lib/libpython2.3.so.1.0 #36 0x005d5362 in PyEval_EvalCode () from /usr/lib/libpython2.3.so.1.0 #37 0x005ee817 in PyErr_Display () from /usr/lib/libpython2.3.so.1.0 #38 0x005ef942 in PyRun_SimpleFileExFlags () from /usr/lib/libpython2.3.so.1.0 #39 0x005f0994 in PyRun_AnyFileExFlags () from /usr/lib/libpython2.3.so.1.0 #40 0x005f568e in Py_Main () from /usr/lib/libpython2.3.so.1.0 #41 0x080485b2 in main () # End of BT. This code is running on a Fedora Core 3 box, with python 2.3.4 and numpy/scipy compiled using gcc 3.4.4. I realize that it's extremely difficult to help with so little information, but unfortunately we have no small test that can reproduce the problem. Only our large research codes, when running for multiple days on a single run, cause this. Even very intensive uses of the same code but which last only a few hours never show this. This code is a long-runing iterative algorithm, so it's basically applying the same (complex) loop over and over until convergence, using numpy and scipy pretty extensively throughout. If super Travis (or anyone else) can have a Eureka moment from the above backtrace, that would be fantastic. If there's any other information you think I may be able to provide, I'll be happy to do my best. Cheers, f ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

10 22

Re: Need more comments from scientific community on python-dev
by Josh Marshall 04 Nov '06

04 Nov '06

On 11/31/06, Fernando Perez <fperez.net(a)gmail.com> wrote: > Fernando Perez wrote: > ps - one more thing. This guy: > > http://blog.vrplumber.com/ > > has been rewriting the OpenGL bindings using ctypes, and I've seen > posts from him about numpy (in his blog). He might be able to > contribute something... I've been working on the OSX port and extensions for OpenGL-ctypes. (now released as PyOpenGL-3.00a, please test!). What Mike has done is define setuptools plugins to interface to different kinds of array data. The default type is ctypes arrays, or numpy arrays if numpy is installed. The data types handled are ctypes sized arrays, ctypes pointers, strings (read-only), and Python lists. (and old Numeric/ numarray via a non-default build.) Summaries of development and usage, resp, are found: http://pyopengl.sourceforge.net/ctypes/development.html http://pyopengl.sourceforge.net/ctypes/using.html Notes on array handling from the above: > Perhaps the most complex mechanisms in OpenGL-ctypes are those > which implement the array-based operations which allow for using > low-level blocks of formatted data to communicate with the OpenGL > implementation. OpenGL-ctypes preferred basic array implementation > is the (new) numpy reimplementation of the original Numeric Python. > > The array handling functionality provided within OpenGL-ctypes is > localised to the OpenGL.arrays sub-package. Within the package, > there are two major classes, one (the FormatHandler) which > implements an interface to a way of storing data in Python, and > another (the ArrayDatatype) which models an OpenGL array format. > The ArrayDatatype classes use FormatHandlers to manipulate array- > compatible objects for use in the system. and on FormatHandlers, for each of the datatypes I mentioned earlier: > Each format handler is responsible for implementing an API that > ArrayDatatypes can use to work with the Python data-format. Data- > formats can support a subset of the API, they only need to support > those aspects of the data-format which make sense. Now, I haven't spent much time looking at these parts of OpenGL- ctypes, as they have just worked for me. I would think that it would be trivial to write a FormatHandler which uses the ndarray interface and data type description to use any object implenting it as an input for OpenGL. This would include things such as PIL images. Mike, can you give us your opinion on how a standardised data type descriptor would be helpful for PyOpenGL? The PEP and some information about it can be found here: http://www.scipy.org/ArrayInterfacePEP Cheers, Josh ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

3 2

Need more comments from scientific community on python-dev
by Travis Oliphant 04 Nov '06

04 Nov '06

I'm recruiting more comments on python-dev regarding my two proposals for improving Python's native ability to share ndarray-like information. There is a dearth of scientific-computing and number-crunching-aware people on python-dev. The result is that I sound like a lone voice arguing for something that nobody cares about. When, I don't think that is true. Please, please. If you want Python to grow support for the array interface (or something like it), then please speak up on python-dev. Even something as simple as I really see the need for a way to exchange data-format information between two objects sharing the buffer protocol can be helpful. You can post through the gmane newsgroup interface: gmane.comp.python.devel Find any of the posts on the PEP's I've introduced. Thanks for your help. -Travis ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

11 21

Reading records from file and sorting
by George Sakkis 03 Nov '06

03 Nov '06

Is there a more elegant and/or faster way to read some records from a file and then sort them by different fields ? What I have now is too specific and error-prone in general: import numpy as N records = N.fromfile(a_file, dtype=N.dtype('i2,i4')) records_by_f0 = records.take(records.getfield('i2').argsort()) records_by_f1 = records.take(records.getfield('i4',2).argsort()) If there's a better way, I'd like to see it; bonus points for in-place sorting. Thanks, George ------------------------------------------------------------------------- Using Tomcat but need to do more? Need to support web services, security? Get stuff done quickly with pre-integrated technology to make your job easier Download IBM WebSphere Application Server v.1.0.1 based on Apache Geronimo http://sel.as-us.falkag.net/sel?cmd=lnk&kid=120709&bid=263057&dat=121642

6 9

NumPy-Discussion October 2006