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

Scipy - Numpy incompatiblility when calling upon Old Numeric
by frbeaxs Jan. 3, 2007

Jan. 3, 2007

I am using Python 2.4 with Numpy 0.9.8. Matpylib graphs function under these two version except when old numeric must be called to utilized the spline function in contouring the graphs leading to the error: Traceback (most recent call last): File "C:\File\XXX", line 3, in ? from scipy.sandbox.delaunay import * File "C:\Python24\lib\site-packages\scipy\__init__.py", line 32, in ? from numpy import oldnumeric File "C:\Python24\lib\site-packages\numpy\oldnumeric\__init__.py", line 3, in ? from compat import * File "C:\Python24\lib\site-packages\numpy\oldnumeric\compat.py", line 133, in ? from numpy import deprecate ImportError: cannot import name deprecate Scipy 0.5.0 seems to be compatible only with Numpy 101b, not version 0.9.8 and I have to switch version of Numpy to should I need to utilize Scipy, however certain matpylib graphs will not function. Upgrading Scipy to 0.5.2 and Numpy 101 has no negative effects on the programs I wrote or the matplotlib sample programs, however the following error message replaces the Import error above Matrix = matrix NameError:name matrix is not defined. Downgrading NUMPY to ver 9.6 produces even more problems. The best match seems to be Mumpy 0.9.8 used with SciPy 0.50 seems to be the best combination however it prohibits calling upon Oldnemeric which seems to be necessary to use the spline function. Does anyone have any idea on how to get around this? frbeaxs

2 1

newbie: attempt at data frame
by Vincent Nijs Jan. 3, 2007

Jan. 3, 2007

I just started working on a time-series module/class in scipy/numpy and it seemed useful to have some of the R data-frame functionality (i.e., select columns of data based on variable names). I tried rec-arrays but couldn't get them to work the way I wanted. I also looked at the Dataframe class by Andrew Straw but at over 400 lines of code that seemed pretty complicated, to me at least. I searched the mailing-list archives and found a discussion on 'Table like array' (see exert below). To get the minimal functionality discussed, I wrote a simple class (see attached) to try and implement X.get('a','c') where 'a' and 'c' are variables names linked to columns of data in X. I added some test code so that if you run the code in the attachment you will see that is seems to work. However, since this is my first class I'd appreciate your input on the approach I used and any suggestions on how to improve the class (or use something else). I'd like to read the data and variable names directly from a single csv file. I tried this through the python csv module but it would read all data as strings and I couldn't figure out how to easily separate the variable names and the data. Thanks, Vincent > [Numpy-discussion] Re: [SciPy-user] Table like array > Paul Barrett pebarrett at gmail.com > Wed Mar 1 06:45:02 CST 2006 > > On 3/1/06, Travis Oliphant <oliphant.travis at ieee.org> wrote: >> >> >> How many people would like to see x['f1','f2','f5'] return a new array >> with a new data-type descriptor constructed from the provided fields? >> >> > > I'm surprised that it's not already available. > > -- Paul

3 7

numpy install on mac os x 10.4
by belinda thom Dec. 31, 2006

Dec. 31, 2006

Hi, I just used easy_install to get the latest version of numpy. Is this the "preferred" method for installing? I'm on a G5 w/mac OS X 10.4 and MacPython 2.4. I'm wondering why the os-related file names that easy_install creates have macosx-10.3 in them (as opposed to 10.4), e.g. creating /Library/Frameworks/Python.framework/Versions/2.4/lib/ python2.4/site-packages/numpy-1.0.1-py2.4-macosx-10.3-fat.egg Is this something I should be concerned about? Thanks, --b

5 10

PyArray_FromDims segfault??
by Xavier Gnata Dec. 31, 2006

Dec. 31, 2006

Hello, I would like to use PyObject_CallObject to call the imshow function from matplolib from a C code. So, the first step is to create a 2D PyArrayObject from a C array. I have read the Numpy book which is great but now I'm puzzled: My goal is to convert a C array of doubles into a 2D numpy array object. PyArray_SimpleNewFromData(nd, dims, typenum, data) seems to be the function I need but I did not manage to write/find any compiling *and* not segfaulting code performing this convertion :( Ok, maybe PyArray_SimpleNewFromData is to complex to begin with. Let's try with PyArray_FromDims... segfault also :( My code looks so simpleKK I cannot see what I'm doing the wrong way :( #include <python2.4/Python.h> #include <numpy/arrayobject.h> int main (int argc, char *argv[]) { PyArrayObject *array; int i, length = 100; double* a; double x; Py_Initialize (); array = (PyArrayObject *) PyArray_FromDims(1, &length, PyArray_DOUBLE); a = new double[100]; for (i = 0; i < length; i++) { x = (double) i/(length-1); a[i] = x; } Py_Finalize (); return 0; } g++ imshow.cpp -o imshow -lpython2.4 -I/usr/lib/python2.4/site-packages/numpy/core/include/ ./imshow -> Segmentation fault I really need the Py_Initialize/Py_Finalize calls because it is 'embedding python' (and not extending :)). To sum up, any code taking an double* and 2 int standind for the dimX and the dimY and returning the corresponding PyObject would be very much appreciated :) I don't care if it has to copy the data. The simpler the better. In [2]: numpy.__version__ Out[2]: '1.0.2.dev3491' Xavier. ps : I really want to use only the C API first at all to learn it and because it really looks like as the simple way to fit my needs. Ok PyObject_CallObject needs quite a lot of additionnal code but it is always almost the same :) -- ############################################ Xavier Gnata CRAL - Observatoire de Lyon 9, avenue Charles André 69561 Saint Genis Laval cedex Phone: +33 4 78 86 85 28 Fax: +33 4 78 86 83 86 E-mail: gnata(a)obs.univ-lyon1.fr ############################################

2 1

Re: [Numpy-discussion] Advice please on efficient subtotal function
by Andreas Eisele Dec. 30, 2006

Dec. 30, 2006

Hi Stephen, > > I'm looking for efficient ways to subtotal a 1-d array onto a 2-D > grid. This > is more easily explained in code that words, thus: > > for n in xrange(len(data)): > totals[ i[n], j[n] ] += data[n] > > data comes from a series of PyTables files with ~200m rows. Each row > has ~20 > cols, and I use the first three columns (which are 1-3 char strings) > to form > the indexing functions i[] and j[], then want to calc averages of the > remaining 17 numerical cols. > > I have tried various indirect ways of doing this with searchsorted and > bincount, but intuitively they feel overly complex solutions to what > is > essentially a very simple problem. > > My work involved comparing the subtotals for various different > segmentation > strategies (the i[] and j[] indexing functions). Efficient solutions > are > important because I need to make many passes through the 200m rows of > data. > Memory usage is the easiest thing for me to adjust by changing how > many rows > of data to read in for each pass and then reusing the same array data > buffers. > It looks as if the values in your i and j columns come from a limited range, so you may consider encoding pairs of (i,j) values into one int using a suitable encoding function (e.g. ij = i+K*j if both i and j are non-negative and K=max(i)+1). You could then use bincount(ij, data) to get the sums per encoded (i,j) pair. This should be efficient, and the complexity is only in the encoding/decoding steps. Best regards, Andreas ---- Dr. Andreas Eisele Senior Researcher DFKI GmbH, Language Technology Lab eisele(a)dfki.de Stuhlsatzenhausweg 3 tel: +49-681-302-5285 D-66123 Saarbrücken fax: +49-681-302-5338

1 0

test issue
by belinda thom Dec. 30, 2006

Dec. 30, 2006

Hello, I've been going thru Dave Kuhlman's "SciPy Course Outline" (http:// www.rexx.com/~dkuhlman/scipy_course_01.html) and found out about test functions -- very cool. Except that on my end, not all tests pass (appended below). Is this a problem for other people? Is it something I should worry about? Here's my setup: Mac G5 w/OS X 10.4.8, using MacPython 2.4, numpy.__version__ is 1.0, matplotlib.__version__ 0.87.7 and Numeric.__version__ 24.2 Thanks, --b ========== In [94]: import numpy In [95]: numpy.test() Found 13 tests for numpy.core.umath Found 9 tests for numpy.lib.arraysetops Found 3 tests for numpy.fft.helper Found 1 tests for numpy.lib.ufunclike Found 4 tests for numpy.ctypeslib Found 2 tests for numpy.lib.polynomial Found 8 tests for numpy.core.records Found 26 tests for numpy.core.numeric Found 5 tests for numpy.distutils.misc_util Found 3 tests for numpy.lib.getlimits Found 31 tests for numpy.core.numerictypes Found 4 tests for numpy.core.scalarmath Found 12 tests for numpy.lib.twodim_base Found 47 tests for numpy.lib.shape_base Found 4 tests for numpy.lib.index_tricks Found 32 tests for numpy.linalg.linalg Found 42 tests for numpy.lib.type_check Found 184 tests for numpy.core.multiarray Found 36 tests for numpy.core.ma Found 10 tests for numpy.core.defmatrix Found 41 tests for numpy.lib.function_base Found 0 tests for __main__ ........................................................................ ........................................................................ ........................................................................ ........................................................................ ........................................................................ F....................................................................... ........................................................................ ............. ====================================================================== FAIL: Ticket #112 ---------------------------------------------------------------------- Traceback (most recent call last): File "/Library/Frameworks/Python.framework/Versions/2.4/lib/ python2.4/site-packages/numpy/core/tests/test_regression.py", line 220, in check_longfloat_repr assert(str(a)[1:9] == str(a[0])[:8]) AssertionError ---------------------------------------------------------------------- Ran 517 tests in 1.241s FAILED (failures=1) <unittest.TextTestRunner object at 0x6194f30>

2 4

Advice please on efficient subtotal function
by Stephen Simmons Dec. 30, 2006

Dec. 30, 2006

Hi, I'm looking for efficient ways to subtotal a 1-d array onto a 2-D grid. This is more easily explained in code that words, thus: for n in xrange(len(data)): totals[ i[n], j[n] ] += data[n] data comes from a series of PyTables files with ~200m rows. Each row has ~20 cols, and I use the first three columns (which are 1-3 char strings) to form the indexing functions i[] and j[], then want to calc averages of the remaining 17 numerical cols. I have tried various indirect ways of doing this with searchsorted and bincount, but intuitively they feel overly complex solutions to what is essentially a very simple problem. My work involved comparing the subtotals for various different segmentation strategies (the i[] and j[] indexing functions). Efficient solutions are important because I need to make many passes through the 200m rows of data. Memory usage is the easiest thing for me to adjust by changing how many rows of data to read in for each pass and then reusing the same array data buffers. Thanks in advance for any suggestions! Stephen

3 3