Mailman 3 June 2006 - NumPy-Discussion

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

17 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

updated Ubuntu Dapper packages for numpy, matplotlib, and scipy online
by Andrew Straw 18 Jul '06

18 Jul '06

I have updated the apt repository I maintain for Ubuntu's Dapper, which now includes: numpy matplotlib scipy Each package is from a recent SVN checkout and should thus be regarded as "bleeding edge". The repository has a new URL: http://debs.astraw.com/dapper/ I intend to keep this repository online for an extended duration. If you want to put this repository in your sources list, you need to add the following lines to /etc/apt/sources.list:: deb http://debs.astraw.com/ dapper/ deb-src http://debs.astraw.com/ dapper/ I have not yet investigated the use of ATLAS in building or using the numpy binaries, and if performance is critical for you, please evaluate speed before using it. I intend to visit this issue, but I cannot say when. The Debian source packages were generated using stdeb, [ http://stdeb.python-hosting.com/ ] a Python to Debian source package conversion utility I wrote. stdeb does not build packages that follow the Debian Python Policy, so the packages here may be slighly unusual compared to Python packages in the official Debian or Ubuntu repositiories. For example, example scripts do not get installed, and no documentation is installed. Future releases of stdeb may resolve these issues. As always, feedback is very appreciated. Cheers! Andrew

3 5

unique() should return a sorted array
by David Huard 17 Jul '06

17 Jul '06

Hi, Numpy's unique(x) returns an array x with repetitions removed. However, since it returns asarray(dict.keys()), the resulting array is not sorted, worse, the original order may not be conserved. I think that unique() should return a sorted array, like its matlab homonym. Regards, David Huard

10 17

logspace behaviour/documentation
by Arnd Baecker 10 Jul '06

10 Jul '06

Hi, I am wondering a bit about the the behaviour of logspace: Definition: numpy.logspace(start, stop, num=50, endpoint=True, base=10.0) Reading this I would assume that numpy.logspace(10**-12, 0.0, 100) gives 100 values, from start=10**-12 to stop=0.0, equispaced on a logarithmic scale. But this is not the case. Instead one has to do: numpy.logspace(-12, 0.0, 100) Docstring: Evenly spaced numbers on a logarithmic scale. Computes int(num) evenly spaced exponents from start to stop. If endpoint=True, then last exponent is stop. Returns base**exponents. My impression is that only the very last line is clearly saying what logspace does. And of course the code itself: y = linspace(start,stop,num=num,endpoint=endpoint) return _nx.power(base,y) Possible solutions (see below): a) modify logspace so that numpy.logspace(10**-12, 0.0, 100) works b) keep the current behaviour and improve the doc-string I would be interested in opinions on this. Best, Arnd Possible solution for (a) (no error checking yet): def logspace_modified(start, stop, num=50, endpoint=True): """Evenly spaced numbers on a logarithmic scale. Computes `num` evenly spaced numbers on a logarithmic scale from `start` to `stop`. If endpoint=True, then last exponent is `stop`. """ lstart = log(start) lstop = log(stop) y = linspace(lstart, lstop, num=num, endpoint=endpoint) return exp(y) Possible improvent of the doc-string (b) - due to Lars Bittrich: def logspace(start,stop,num=50,endpoint=True,base=10.0): """Evenly spaced numbers on a logarithmic scale. Return 'int(num)' evenly spaced samples on a logarithmic scale from 'base'**'start' to 'base'**'stop'. If 'endpoint' is True, the last sample is 'base'**'stop'."""

2 2

recarray field names
by Erin Sheldon 08 Jul '06

08 Jul '06

Hi all- I'm new to numpy and I'm slowly weening myself off of IDL. So far the experience has been very positive. I use recarrays a lot. I often read these recarrays from fits files using pyfits (modified to work with numpy). I find myself doing the following more than I would like: if 'TAG' in rec._coldefs.names: .... It seems messy to be accessing this "hidden" attribute in this way. Is there a plan to add methods to more transparently do such things? e.g. def fieldnames(): return _coldefs.names def field_exists(fieldname): return fieldname.upper() in _coldefs.names def field_index(fieldname): if field_exists(fieldname): return _coldefs.names.index(fieldname.upper()) else: return -1 # or None maybe Thanks, Erin

5 13

MA bug or feature?
by Michael Sorich 04 Jul '06

04 Jul '06

When transposing a masked array of dtype '<f8' I noticed that an ndarray of dtype '|O4' was returned. I found this a little strange as the masked array did not contain any masked values. Upon closer examination (see script below) it seems that if the mask is a boolean array with all false values and the fill_value is a string this will occur. However if the fill_value is a number or the mask is simply False, the dtype stays as '<f8'. I was going to submit this as a bug but wanted to check that this was not a deliberate feature. I am using the numpy version that comes with a recent version of python enthought edition. import numpy print numpy.version.version ma1 = numpy.ma.array(((1.,2,3),(4,5,6)), mask=((0,0,0),(0,0,0))) print ma1.filled('NA').dtype print ma1.filled(-999).dtype ma2 = numpy.ma.array(((1.,2,3),(4,5,6)), mask=False) print ma2.filled('NA').dtype print ma2.filled(-999).dtype --------------- output: 0.9.9.2538 '|O4' '<f8' '<f8' '<f8'

3 8

iterate along a ray: linear algebra?
by stephen emslie 04 Jul '06

04 Jul '06

I am in the process of implementing an image processing algorithm that requires following rays extending outwards from a starting point and calculating the intensity derivative at each point. The idea is to find the point where the difference in intensity goes beyond a particular threshold. Specifically I'm examining an image of an eye to find the pupil, and the edge of the pupil is a sharp change in intensity. How does one iterate along a line in a 2d matrix, and is there a better way to do this? Is this a problem that linear algebra can help with? Thanks Stephen Emslie

3 3

Time for beta1 of NumPy 1.0
by Travis Oliphant 03 Jul '06

03 Jul '06

I think it's time for the first beta-release of NumPy 1.0 I'd like to put it out within 2 weeks. Please make any comments or voice major concerns so that the 1.0 release series can be as stable as possible. -Travis

38 86