Mailman 3 July 2006 - NumPy-Discussion

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

May 16, 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 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 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 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 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.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

Patch against Image.py in the PIL
by Travis Oliphant Aug. 20, 2006

Aug. 20, 2006

Attached is a patch that makes PIL Image objects both export and consume the array interface. -Travis Index: PIL/Image.py =================================================================== --- PIL/Image.py (revision 358) +++ PIL/Image.py (working copy) @@ -187,6 +187,42 @@ } +if sys.byteorder == 'little': + _ENDIAN = '<' +else: + _ENDIAN = '>' + +_MODE_CONV = { + + # official modes + "1": ('|b1', None), + "L": ('|u1', None), + "I": ('%si4' % _ENDIAN, None), + "F": ('%sf4' % _ENDIAN, None), + "P": ('|u1', None), + "RGB": ('|u1', 3), + "RGBX": ('|u1', 4), + "RGBA": ('|u1', 4), + "CMYK": ('|u1', 4), + "YCbCr": ('|u1', 4), + + # Experimental modes include I;16, I;16B, RGBa, BGR;15, + # and BGR;24. Use these modes only if you know exactly + # what you're doing... + +} + +def _conv_type_shape(im): + shape = im.size[::-1] + typ, extra = _MODE_CONV[im.mode] + if extra is None: + return shape, typ + else: + return shape+(extra,), typ + + + + MODES = _MODEINFO.keys() MODES.sort() @@ -491,6 +527,22 @@ return string.join(data, "") ## + # Returns the array_interface dictionary + # + # @return A dictionary with keys 'shape', 'typestr', 'data' + + def __get_array_interface__(self): + new = {} + shape, typestr = _conv_type_shape(self) + new['shape'] = shape + new['typestr'] = typestr + new['data'] = self.tostring() + return new + + __array_interface__ = property(__get_array_interface__, None, doc="array interface") + + + ## # Returns the image converted to an X11 bitmap. This method # only works for mode "1" images. # @@ -1749,7 +1801,61 @@ return apply(fromstring, (mode, size, data, decoder_name, args)) + ## +# (New in 1.1.6) Create an image memory from an object exporting +# the array interface (using the buffer protocol). +# +# If obj is not contiguous, then the tostring method is called +# and frombuffer is used +# +# @param obj Object with array interface +# @param mode Mode to use (will be determined from type if None) + +def fromarray(obj, mode=None): + arr = obj.__array_interface__ + shape = arr['shape'] + ndim = len(shape) + try: + strides = arr['strides'] + except KeyError: + strides = None + if mode is None: + typestr = arr['typestr'] + if not (typestr[0] == '|' or typestr[0] == _ENDIAN or + typestr[1:] not in ['u1', 'b1', 'i4', 'f4']): + raise TypeError, "cannot handle data-type" + typestr = typestr[:2] + if typestr == 'i4': + mode = 'I' + elif typestr == 'f4': + mode = 'F' + elif typestr == 'b1': + mode = '1' + elif ndim == 2: + mode = 'L' + elif ndim == 3: + mode = 'RGB' + elif ndim == 4: + mode = 'RGBA' + else: + raise TypeError, "Do not understand data." + ndmax = 4 + bad_dims=0 + if mode in ['1','L','I','P','F']: + ndmax = 2 + elif mode == 'RGB': + ndmax = 3 + if ndim > ndmax: + raise ValueError, "Too many dimensions." + + size = shape[:2][::-1] + if strides is not None: + obj = obj.tostring() + + return frombuffer(mode, size, obj) + +## # Opens and identifies the given image file. # <p> # This is a lazy operation; this function identifies the file, but the

2 1

trivial question: how to compare dtype - but ignoring byteorder ?
by Sebastian Haase Aug. 14, 2006

Aug. 14, 2006

Hi, if I have a numpy array 'a' and say: a.dtype == numpy.float32 Is the result independent of a's byteorder ? (That's what I would expect ! Just checking !) Thanks, Sebastian Haase

5 11

Re: [Numpy-discussion] Profiling line-by-line
by David Grant Aug. 13, 2006

Aug. 13, 2006

Does anyone know if this issue related to profiling with numpy is a python problem or a numpy problem? Dave On 7/20/06, David Grant <davidgrant(a)gmail.com> wrote: > > > > On 7/20/06, Arnd Baecker <arnd.baecker(a)web.de> wrote: > > > > > > More importantly note that profiling in connection > > with ufuncs seems problematic: > > > Yes, that seems to be my problem... I read the threads you provided links > to. Do you know why this is the case? > > I have tried hotshot2calltree by the way, and I didn't find out anything > new. > > -- > David Grant > -- David Grant

1 1

Moving docstrings from C to Python
by Stefan van der Walt Aug. 9, 2006

Aug. 9, 2006

Hi, Travis recently advertised the add_docstring and add_newdoc functions, which are used to add documentation to functions written in C. This brings the advantage that you don't have to recompile everytime you change a docstrings, and also that the docstrings are much simpler to format (no more "\n\" tokens all over the place). I took all the docstrings from numpy.core.src.multiarray and formatted them to be used with add_docstring. I have two questions: Would anyone mind if the change was made? If not, where should they go? (numpy/add_newdocs.py or numpy/core/something) Thanks for the feedback. Stéfan

3 3

comments on r_ and c_ ?
by Bill Baxter Aug. 6, 2006

Aug. 6, 2006

When you have a chance, could the powers that be make some comment on the r_ and c_ situation? From what I understand Travis added these fairly recently and they aren't officially documented anywhere. There was discussion previously about whether to rename them or not, but no real conclusion was reached on better names (I've kind of grown to like the x_ style naming in the mean time). I do think the subtle differences in behavior between these and existing methods (hstack,vstack,column_stack) are an issue, though, which I outlined in a previous email. Thanks, --Bill

3 3

fixing diag() for matrices
by Sven Schreiber Aug. 6, 2006

Aug. 6, 2006

Hi, there was a thread about this before, diag() is currently only partly useful if you work with numpy-matrices, because the 1d->2d direction doesn't work, as there are no 1d-numpy-matrices. This is unfortunate because a numpy-matrix with shape (n,1) or (1,m) should be naturally treated as a vector, imho. So it would be nice if this could be fixed. It's probably not the most efficient solution, but what I want for numpy-matrix input x is to get: mat(diag(x.A.squeeze)) where diag is the current implementation. This means that if x is not a vector ("truly 2d"), then nothing is changed. But if one of the dimensions of x is ==1, then it's turned into a 1d-array, and diag works as it should. Does that sound reasonable? Thanks, Sven

6 13

numpy vs numarray
by Louis Cordier Aug. 1, 2006

Aug. 1, 2006

Hmmm, I think people are spreading "fud" (lower case)... http://media.libsyn.com/media/awaretek/Python411_060530_Numeric.mp3 -- Louis Cordier <lcordier(a)point45.com> cell: +27721472305 Point45 Entertainment (Pty) Ltd. http://www.point45.org

6 8