Mailman 3 October 2006 - NumPy-Discussion

PyArray_ContiguousFromObject & PyDECREF
by Michiel Jan Laurens de Hoon May 17, 2011

May 17, 2011

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 … [View More]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 [View Less]

3 2

Re: [Numpy-discussion] Floating point exception with numpy and embedded python interpreter
by Andrew Straw Sept. 20, 2007

Sept. 20, 2007

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 … [View More]

5 7

Numpy x Matlab: some synthetic benchmarks
by Paulo Jose da Silva e Silva Feb. 11, 2007

Feb. 11, 2007

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.… [View More]

13 27

A reimplementation of MaskedArray
by Pierre GM Nov. 23, 2006

Nov. 23, 2006

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 … [View More]

8 27

A way to compute a 2d histogram
by Xavier Gnata Nov. 14, 2006

Nov. 14, 2006

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 … [View More]

3 9

Helper function to "unroll" a array
by Gael Varoquaux Nov. 12, 2006

Nov. 12, 2006

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.… [View More]

2 3

Strange and hard to reproduce crash
by Fernando Perez Nov. 6, 2006

Nov. 6, 2006

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 … [View More]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 [View Less]

10 22

Re: Need more comments from scientific community on python-dev
by Josh Marshall Nov. 4, 2006

Nov. 4, 2006

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 … [View More]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 [View Less]

3 2

Need more comments from scientific community on python-dev
by Travis Oliphant Nov. 4, 2006

Nov. 4, 2006

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 … [View More]

11 21

Reading records from file and sorting
by George Sakkis Nov. 4, 2006

Nov. 4, 2006

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 -------------------------… [View More]

6 9