Mailman 3 March 2007 - 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 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

rant against from numpy import * / from pylab import *
by Sebastian Haase Aug. 9, 2007

Aug. 9, 2007

Hi! I use the wxPython PyShell. I like especially the feature that when typing a module and then the dot "." I get a popup list of all available functions (names) inside that module. Secondly, I think it really makes code clearer when one can see where a function comes from. I have a default import numpy as N executed before my shell even starts. In fact I have a bunch of my "standard" modules imported as <some single capital letter>. This - I think - is a good compromise to the commonly used "extra typing" and "unreadable" argument. a = sin(b) * arange(10,50, .1) * cos(d) vs. a = N.sin(b) * N.arange(10,50, .1) * N.cos(d) I would like to hear some comments by others. On a different note: I just started using pylab, so I did added an automatic "from matplotlib import pylab as P" -- but now P contains everything that I already have in N. It makes it really hard to *find* (as in *see* n the popup-list) the pylab-only functions. -- what can I do about this ? Thanks, Sebastian

2 3

numexpr efficency depends on the size of the computing kernel
by Francesc Altet April 26, 2007

April 26, 2007

Hi, Now that I'm commanding my old AMD Duron machine, I've made some benchmarks just to prove that the numexpr computing is not influenced by the size of the CPU cache, but I failed miserably (and Tim was right: there is a dependency of the numexpr efficency on CPU cache size). Provided that the pytables instance of the computing kernel of numexpr is quite larger (it supports more datatypes) than the original, comparing the performance of both versions can be a good way to check the influence of CPU cache on the computing efficency. The attached benchmark is a small modification of the timing.py that comes with the numexpr package (the modification was needed to allow the numexpr version of pytables to run all the cases). Basically, the expressions tested operations with arrays of 1 million of elements, with a mix of contiguous and strided arrays (no unaligned arrays are present here). See the code in benchmark for the details. The speed-ups of numexpr over plain numpy on a AMD Duron machine (64 + 64 KB L1 cache, 64 KB L2 cache) are: For the original numexpr package: 2.14, 2.21, 2.21 (these represent averages for 3 complete runs) For the modified pytables version (enlarged computing kernel): 1.32, 1.34, 1.37 So, with a CPU with a very small cache, the original numexpr kernel is 1.6x faster than the pytables one. However, using a AMD Opteron which has a much bigger L2 cache (64 + 64 KB L1 cache, 1 MB L2 cache), the speed-ups are quite similar: For the original numexpr package: 3.10, 3.35, 3.35 For the modified pytables version (enlarged computing kernel): 3.37, 3.50, 3.45 So, there is effectively a dependency on the CPU cache size. It would be nice to run the benchmark with other CPUs with a L2 cache in the range between 64 KB and 1 MB so as to find the point where the performance starts to be similar (this should be a good guess on the size of the computing kernel). Meanwhile, the lesson learned is that Tim worries were correct: one should be very careful on adding more opcodes (at least, until CPUs with a very small L2 cache are in use). With this, perhaps we will have to reduce the opcodes in the numexpr version for pytables to a bare minimum :-/ Cheers, -- Francesc Altet | Be careful about using the following code -- Carabos Coop. V. | I've only proven that it works, www.carabos.com | I haven't tested it. -- Donald Knuth

1 1

Re: [Numpy-discussion] nd_image.affine_transform edge effects
by James Turner April 13, 2007

April 13, 2007

Hi Stefan, Thanks for the suggestions! > Is this related to > > http://projects.scipy.org/scipy/scipy/ticket/213 > > in any way? As far as I can see, the problems look different, but thanks for the example of how to document this. I did confirm that your example exhibits the same behaviour under numarray, in case that is useful. > Code snippets to illustrate the problem would be welcome. OK. I have had a go at producing a code snippet. I apologize that this is based on numarray rather than numpy, since I'm using STScI Python, but I think it should be very easy to convert if you have numpy instead. What I am doing is to transform overlapping input images onto a common, larger grid and co-adding them. Although I'm using affine_ transform on 3D data from FITS images, the issue can be illustrated using a simple 1D translation of a single 2D test array. The input values are just [4., 3., 2., 1.] in each row. With a translation of -0.1, the values should therefore be something like [X, 3.1, 2.1, 1.1, X, X], where the Xs represent points outside the original data range. What I actually get, however, is roughly [X, 3.1, 2.1, 1.0, 1.9, X]. The 5th value of 1.9 contaminates the co-added data in the final output array. Now I'm looking at this element-by-element, I suppose the bad value of 1.9 is just a result of extrapolating in order to preserve the original number of data points, isn't it? Sorry I wasn't clear on that in my original post -- but surely a blank value (as specified by cval) would be better? I suppose I could work around this by blanking out the extrapolated column after doing the affine_transform. I could calculate which is the column to blank out based on the sense of the offset and the input array dimensions. It seems pretty messy and inefficient though. Another idea is to split the translation into integer and fractional parts, keep the input and output array dimensions the same initially and then copy the output into a larger array with integer offsets. That is messy to keep track of though. Maybe a parameter could instead be added to affine_transform that tells it to shrink the number of elements instead of extrapolating? I'd be a bit out of my depth trying to implement that though, even if the authors agree... (maybe in a few months though). Can anyone comment on whether this problem should be considered a bug, or whether it's intended behaviour that I should work around? The code snippet follows below. Thanks for your patience with someone who isn't accustomed to posting questions like this routinely :-). James. ----- import numarray as N import numarray.nd_image as ndi # Create a 2D test pattern: I = N.zeros((2,4),N.Float32) I[:,:] = N.arange(4.0, 0.0, -1.0) # Transformation parameters for a simple translation in 1D: trmatrix = N.array([[1,0],[0,1]]) troffset = (0.0, -0.1) # Apply the offset to the test pattern: I_off1 = ndi.affine_transform(I, trmatrix, troffset, order=3, mode='constant', cval=-1.0, output_shape=(2,6)) I_off2 = ndi.affine_transform(I, trmatrix, troffset, order=3, mode='constant', cval=-1.0, output_shape=(2,6), prefilter=False) # Compare the data before and after interpolation: print I print I_off1 print I_off2

6 40

Dynamic module not initialized properly
by Mark Janikas April 2, 2007

April 2, 2007

Hello all, I am having an issue importing numpy on subsequent (I.e. not on first load) attempts in our software. The majority of the code is written in C, C++ and I am a python developer and do not have direct access to a lot of it. This is a bit of a difficult question to ask all of you because I cant provide you a direct example. All I can do is point to a numpy thread that discusses the issue: http://groups.google.com/group/Numpy-discussion/browse_thread/thread/321 77a82deab05ae/d8eecaf494ba5ad5?lnk=st&q=dynamic+module+not+initialized+p roperly+numpy&rnum=1&hl=en#d8eecaf494ba5ad5 ERROR: exceptions.SystemError: dynamic module not initialized properly What is really odd about my specific issue is that if I don't change anything in the source code.... Then the error doesn't pop up. Furthermore, the error doesn't show on some attempts even after I make a change!!!! Not sure whether there is anything I can do from the scripting side (some alternative form of reload?)... or if I have to forward it along to the C developers. You have my appreciation ahead of time. Mark Janikas Product Engineer ESRI, Geoprocessing 380 New York St. Redlands, CA 92373 909-793-2853 (2563) mjanikas(a)esri.com

2 2

organizational question
by Alan G Isaac April 1, 2007

April 1, 2007

Is either NumPy or SciPy substantially supported by an identifiable and actual non-profit organization? I ask because I think both fit under http://www.mellon.org/grant_programs/programs/copy_of_research item 4. Here is the announcement: http://matc.mellon.org/ Note that universities are among the nominees: http://matc.mellon.org/2007_nominations Cheers, Alan Isaac

2 1

comments to ticket #454 !? : Importing numpy prevents decrementing references for local variables
by Sebastian Haase March 31, 2007

March 31, 2007

Another rather old ticket is this one: http://projects.scipy.org/scipy/numpy/ticket/454> Any comments on this !? Thanks, -Sebastian On 2/16/07, NumPy <numpy-tickets(a)scipy.net> wrote: > #454: Importing numpy prevents decrementing references for local variables > ------------------------+--------------------------------------------------- > Reporter: goddard | Owner: somebody > Type: defect | Status: new > Priority: normal | Milestone: > Component: numpy.core | Version: 1.0.1 > Severity: major | Keywords: > ------------------------+--------------------------------------------------- > If the first import of numpy is within a function it prevents that > function's > local variables and all locals in higher call frames that led to the > import > from being deleted. This can prevent freeing up large numpy arrays. > > The problem is that numpy __init__.py imports numpy _import_tools.py > and creates a PackageLoader object that saves the parent frame > > self.parent_frame = frame = sys._getframe(1) > > Saving the parent frame results in saving the entire call stack even after > all those functions have returned. > > I've attached Python code illustrating the problem.

1 0

ticket #450
by Sebastian Haase March 31, 2007

March 31, 2007

Hi! Could someone comment on ticket #450 [Numpy-tickets] [NumPy] #450: Make a.min() not copy data It seems odd that calculating the min of an array causes a memory error. The demonstrated data array is a 3D cube of uint8 and about 1GB in size. even creation of an extra (2d) array when reducing the 1st dimension, would only need about 1MB -- that can not really explain the memory error. Thanks, Sebastian

1 0

matrix indexing question
by Alan G Isaac March 31, 2007

March 31, 2007

One thing keeps bugging me when I use numpy.matrix. All this is fine:: >>> x=N.mat('1 1;1 0') >>> x matrix([[1, 1], [1, 0]]) >>> x[1,:] matrix([[1, 0]]) But it seems to me that I should be able to extract a matrix row as an array. So this :: >>> x[1] matrix([[1, 0]]) feels wrong. (Similarly when iterating across rows.) Of course I realize that I can just :: >>> x.A[1] array([1, 0]) but since the above keeps feeling wrong I felt I should raise this as a possible design issue, better discussed early than latter. Cheers, Alan Isaac

17 73