Mailman 3 April 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

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

Guide to Numpy book
by Gennan Chen 08 May '06

08 May '06

Hi! What is the newest version of Guide to numpy? The recent one I got is dated at Jan 9 2005 on the cover. Gen-Nan Chen, PhD Chief Scientist Research and Development Group CorTechs Labs Inc (www.cortechs.net) 1020 Prospect St., #304, La Jolla, CA, 92037 Tel: 1-858-459-9700 ext 16 Fax: 1-858-459-9705 Email: gnchen(a)cortechs.net

5 4

numexpr thoughts
by Tim Hochberg 04 May '06

04 May '06

1. David already mentioned opcodes to copy from an integer array into a float register. Another idea would be to have an opcode to copy from a strided array into a register. This would save a copy for noncontiguous matrices and opens up the door to doing broadcasting. I think we'd have to tweak the main interpreter loop a bit, but I think it's doable. Depending how crunched we feel for opcode space, the casting array and this array could overlap (OP_CAST_AND_COPY for instance). 2. Something I noticed while writing tests is that the rearangement of operations for 'a/const' to '(1/const)*a' changes the results slightly (try const=3). I don't think I care -- I just thought I'd point it out. 3. It may simplify things to use copy_c to eliminate a bunch of the extra bytecodes for operating on functions of more than one argument. I need to check the timings on this -- I don't know how much of a speed hit using copy_c will cause. 4. At some point I plan to add complex operations. My original thought was to just duplicate the current, floating-point interpreter for complex numbers. However, if we have opcodes for casting, it might instead make sense to have one interpreter that works with both. This would be more efficient for mixed expressions since we wouldnt' have to case everything complex right away. It looks like there's enough opcode space, although I worry a bit whether making that main switch loop huge is going to slow things down. Probably not, but it's something to keep an eye on. 5. Same thoughts as 4., but for integer operations. Potentially the interpreter could work with longs, doubles and cdoubles, downcasting as appropriate on the way out. It's looking like numexpr could come pretty close to evaluating most numpy expressions if we put enough work into it. None of the required changes look like they should slow down the original, simple, fast case. However we should keep an eye on that as we go along. It's all very cool; I'd like to thank David for providing such a fun toy. -tim

4 12

Scalar math module is ready for testing
by Travis Oliphant 04 May '06

04 May '06

The scalar math module is complete and ready to be tested. It should speed up code that relies heavily on scalar arithmetic by by-passing the ufunc machinery. It needs lots of testing to be sure that it is doing the "right" thing. To enable scalarmath you need to import numpy.core.scalarmath You cannot disable it once it's enabled except by restarting Python. If we need that feature we can add it. The array scalars respond to the error modes of ufuncs. There is an experimental function called alter_scalars that replaces the Python int, float, and complex number tables with the array scalar equivalents. Thus, to amaze (or seriously annoy) your Python friends you can do import numpy.core.scalarmath as ncs ncs.alter_scalars(int) 1 / 0 This will return 0 unless you change the error modes... ncs.retore_scalars(int) Will put things back the way Guido intended.... Please try it out and send us error reports. Many thanks to Sasha for his help in getting all the code so it at least compiles and loads. All bugs should be blamed on me, though... Best, -Travis

7 8

basearray lives!
by Tim Hochberg 30 Apr '06

30 Apr '06

After a fashion anyway. I implemented the simplest thing that could possibly work and I've left out some stuff that even I think we need (docstring, repr and str). Still it exists, ndarray inherits from it and some stuff seems to work automagically. >>> import numpy as n >>> ba = n.basearray([3,3], int, n.arange(9)) >>> ba <numpy.basearray object at 0x00B29690> >>> a = asarray(ba) >>> a array([[0, 1, 2], [3, 4, 5], [6, 7, 8]]) >>> a + ba array([[ 0, 2, 4], [ 6, 8, 10], [12, 14, 16]]) >>> isinstance(a, n.basearray) True >>> type(ba) <type 'numpy.basearray'> >>> type(a) <type 'numpy.ndarray'> >>> len(dir(ba)) 19 >>> len(dir(a)) 156 Travis: should I go ahead and check this into the trunk? It shouldn't interfear with anything. The only change to ndarray is the tp_base, which sets up the inheritance. -tim

3 2

Creating a descr with aligned=1 using the C API
by Albert Strasheim 30 Apr '06

30 Apr '06

Hello all I was wondering what the best way would be to create the following descr using the C API: descr = dtype({'names' : ['index', 'value'], 'formats' : [intc, 'f8']}, align=1) One could use PyArray_DescrConverter in multiarraymodule.c, but there doesn't seem to be a way to specify aligned=1 and one would have to build the dict object before being able to pass it on for conversion. Unless there's another easy way I'm missing, the API could possibly do with a function like PyArray_DescrFromCommaString(const char*, int align) which calls _convert_from_commastring. By the way, what is the general format of these commastrings? Comments appreciated. Regards, Albert

2 1

[Numeric] "put" into object array corrupts memory
by Sasha 30 Apr '06

30 Apr '06

I know that Numeric is no longer maintained, but since this bug cost me two sleepless nights, I think it is appropriate to announce the bug and the fix to the list. ---------- Forwarded message ---------- From: SourceForge.net <noreply(a)sourceforge.net> Date: Apr 30, 2006 12:58 PM Subject: [ numpy-Bugs-1479376 ] [Numeric] "put" into object array corrupts memory To: noreply(a)sourceforge.net Bugs item #1479376, was opened at 2006-04-30 12:46 Message generated for change (Comment added) made by belopolsky You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101369&aid=1479376&group_… Please note that this message will contain a full copy of the comment thread, including the initial issue submission, for this request, not just the latest update. Category: Fatal Error Group: Normal bug Status: Open Priority: 5 Submitted By: Alexander Belopolsky (belopolsky) Assigned to: Nobody/Anonymous (nobody) Summary: [Numeric] "put" into object array corrupts memory Initial Comment: This is one of those bugs that are easier to fix than to reproduce: $ cat test-put.py class A(object): def __del__(self): print "deleting %r" % self a = A() from Numeric import * x = array([None], 'O') y = array([a], 'O') put(x,[0],y) del a,y print "exiting" $ python test-put.py deleting <__main__.A object at 0xf7e4d24c> exiting Fatal Python error: deletion of interned string failed Aborted (core dumped) Numeric version: 24.2 ---------------------------------------------------------------------- >Comment By: Alexander Belopolsky (belopolsky) Date: 2006-04-30 12:58 Message: Logged In: YES user_id=835142 Attached patch fixes the bug. ---------------------------------------------------------------------- You can respond by visiting: https://sourceforge.net/tracker/?func=detail&atid=101369&aid=1479376&group_…

1 0