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

Different results from repeated calculation
by Keith Goodman 16 Feb '07

16 Feb '07

I get slightly different results when I repeat a calculation. In a long simulation the differences snowball and swamp the effects I am trying to measure. In the attach script there are three tests. In test1, I construct matrices x and y and then repeatedly calculate z = calc(x,y). The result z is the same every time. So this test passes. In test2, I construct matrices x and y each time before calculating z = calc(x,y). Sometimes z is slightly different---of the order of 1e-21 to 1e-18. But the x's test to be equal and so do the y's. This test fails. (It doesn't fail on my friend's windows box; I'm running linux.) test3 is the same as test2 but I calculate z like this: z = calc(100*x ,y) / (100 * 100). This test passes. What is going on? Here is some sample output: >> import repeat >> repeat.all() 4 z different 8.47032947254e-22 5 z different 8.47032947254e-22 7 z different 8.47032947254e-22 9 z different 8.47032947254e-22 10 z different 8.47032947254e-22 16 z different 8.47032947254e-22 24 z different 8.47032947254e-22 25 z different 8.47032947254e-22 26 z different 8.47032947254e-22 27 z different 8.47032947254e-22 30 z different 8.47032947254e-22 32 z different 8.47032947254e-22 34 z different 8.47032947254e-22 35 z different 8.47032947254e-22 36 z different 8.47032947254e-22 39 z different 8.47032947254e-22 40 z different 8.47032947254e-22 41 z different 8.47032947254e-22 45 z different 8.47032947254e-22 46 z different 8.47032947254e-22 50 z different 8.47032947254e-22 52 z different 8.47032947254e-22 53 z different 8.47032947254e-22 55 z different 8.47032947254e-22 56 z different 8.47032947254e-22 58 z different 8.47032947254e-22 66 z different 8.47032947254e-22 67 z different 8.47032947254e-22 71 z different 8.47032947254e-22 73 z different 8.47032947254e-22 74 z different 8.47032947254e-22 83 z different 8.47032947254e-22 86 z different 8.47032947254e-22 87 z different 8.47032947254e-22 88 z different 8.47032947254e-22 89 z different 8.47032947254e-22 90 z different 8.47032947254e-22 92 z different 8.47032947254e-22 test1: 0 differences test2: 38 differences test3: 0 differences Repeated runs tend to give me the same number of differences in test2 for several runs. Then I get a new number of differences which last for severals runs...

11 54

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

Latest Array-Interface PEP
by Travis Oliphant 08 Feb '07

08 Feb '07

I'm attaching my latest extended buffer-protocol PEP that is trying to get the array interface into Python. Basically, it is a translation of the numpy header files into something as simple as possible that can still be used to describe a complicated block of memory to another user. My purpose is to get feedback and criticisms from this community before display before the larger Python community. -Travis PEP: <unassigned> Title: Extending the buffer protocol to include the array interface Version: $Revision: $ Last-Modified: $Date: $ Author: Travis Oliphant <oliphant(a)ee.byu.edu> Status: Draft Type: Standards Track Created: 28-Aug-2006 Python-Version: 2.6 Abstract This PEP proposes extending the tp_as_buffer structure to include function pointers that incorporate information about the intended shape and data-format of the provided buffer. In essence this will place an array interface directly into Python. Rationale Several extensions to Python utilize the buffer protocol to share the location of a data-buffer that is really an N-dimensional array. However, there is no standard way to exchange the additional N-dimensional array information so that the data-buffer is interpreted correctly. The NumPy project introduced an array interface (http://numpy.scipy.org/array_interface.shtml) through a set of attributes on the object itself. While this approach works, it requires attribute lookups which can be expensive when sharing many small arrays. One of the key reasons that users often request to place something like NumPy into the standard library is so that it can be used as standard for other packages that deal with arrays. This PEP provides a mechanism for extending the buffer protocol (which already allows data sharing) to add the additional information needed to understand the data. This should be of benefit to all third-party modules that want to share memory through the buffer protocol such as GUI toolkits, PIL, PyGame, CVXOPT, PyVoxel, PyMedia, audio libraries, video libraries etc. Proposal Add bf_getarrview and bf_relarrview function pointers to the buffer protocol to allow objects to share a view on a memory pointer including information about accessing it as an N-dimensional array. Add the TP_HAS_ARRAY_BUFFER flag to types that define this extended buffer protocol. Also a few additionsl C-API calls should perhaps be added to Python to facilitate creating new PyArrViewObjects. Specification: static PyObject* bf_getarrayview (PyObject *obj) This function must return a new reference to a PyArrViewObject which contains the details of the array information exposed by the object. If failure occurs, then NULL is returned and an exception set. static int bf_relarrayview(PyObject *obj) If not NULL then this will be called when the object returned by bf_getarrview is destroyed so that the underlying object can be aware when acquired "views" are released. The object that defines bf_getarrview should not re-allocate memory (re-size itself) while views are extant. A 0 is returned on success and a -1 and an error condition set on failure. The PyArrayViewObject has the structure typedef struct { PyObject_HEAD void *data; /* pointer to the beginning of data */ int nd; /* the number of dimensions */ Py_ssize_t *shape; /* c-array of size nd giving shape */ Py_ssize_t *strides; /* SEE BELOW */ PyObject *base; /* the object this is a "view" of */ PyObject *format; /* SEE BELOW */ int flags; /* SEE BELOW */ } PyArrayViewObject; strides -- a c-array of size nd providing the striding information which is the number of bytes to skip to get to the next element in that dimension. format -- a Python data-format object (PyDataFormatObject) which contains information about how each item in the array should be interpreted. flags -- an integer of flags. PYARR_WRITEABLE is the only flag that must be set appropriately by types. Other flags: PYARR_ALIGNED, PYARR_C_CONTIGUOUS, PYARR_F_CONTIGUOUS, and PYARR_NOTSWAPPED can all be determined from the rest of the PyArrayViewObject using the UpdateFlags C-API. The PyDataFormatObject has the structure typedef struct { PyObject_HEAD PySimpleformat primitive; /* basic primitive type */ int flags; /* byte-order, isaligned */ int itemsize; /* SEE BELOW */ int alignment; /* SEE BELOW */ PyObject *extended; /* SEE BELOW */ } PyDataFormatObject; enum Pysimpleformat {PY_BIT='1', PY_BOOL='?', PY_BYTE='b', PY_SHORT='h', PY_INT='i', PY_LONG='l', PY_LONGLONG='q', PY_UBYTE='B', PY_USHORT='H', PY_UINT='I', PY_ULONG='L', PY_ULONGLONG='Q', PY_FLOAT='f', PY_DOUBLE='d', PY_LONGDOUBLE='g', PY_CFLOAT='F', PY_CDOUBLE='D', PY_CLONGDOUBLE='G', PY_OBJECT='O', PY_CHAR='c', PY_UCS2='u', PY_UCS4='w', PY_FUNCPTR='X', PY_VOIDPTR='V'}; Each of these simple formats has a special character code which can be used to identify this primitive in a nested python list. flags -- flags for the data-format object. Specified masks are PY_NATIVEORDER PY_BIGENDIAN PY_LITTLEENDIAN PY_IGNORE itemsize -- the total size represented by this data-format in bytes unless the primitive is PY_BIT in which case it is the size in bits. For data-formats that are simple 1-d arrays of the underlying primitive, this total size can represent more than one primitive (with extended still NULL). alignment -- For the primitive types this is offsetof(struct {char c; type v;},v) extended -- NULL if this is a primitive data-type or no additional information is available. If primitive is PY_FUNCPTR, then this can be a tuple with >=1 element: (args, {dim0, dim1, dim2, ...}). args -- A list (of at least length 2) of data-format objects specifying the input argument formats with the last argument specifying the output argument data-format (use None for void inputs and/or outputs). For other primitives, this can be a tuple with >=2 elements: (names, fields, {dim0, dim1, dim2, ...}) Use None for both names and fields if they should be ignored. names -- An ordered list of string or unicode objects giving the names of the fields for a structure data-format. fields -- a Python dictionary with ordered-keys given by the list in names. Each entry in the dictionary is a 3-tuple containing (data-format-object, offset, meta-information) where meta-information is Py_None if there is no meta-information. Offset is given in bytes from the start of the record or in bits if PY_BIT is the primitive. Any additional entries in the extended tuple (dim0, dim1, etc.) are interpreted as integers which specify that this data-format is an array of the given shape of the fundamental data-format specified by the remainder of the DataFormat Object. The dimensions are specified so that the last-index is always assumed to vary the fastest (C-order). The constructor of a PyArrViewObject allocates the memory for shape and strides and the destructor frees that memory. The constructor of a PyDataFormatObject allocates the objects it needs for fields, names, and shape. C-API void PyArrayView_UpdateFlags(PyObject *view, int flags) /* update the flags on the array view object provided */ PyDataFormatObject *Py_NewSimpleFormat(Pysimpleformat primitive) /* return a new primitive data-format object */ PyDataFormatObject *Py_DataFormatFromCType(PyObject *ctype) /* return a new data-format object from a ctype */ int Py_GetPrimitiveSize(Pysimpleformat primitive) /* return the size (in bytes) of the provided primitive */ PyDataFormatObject *Py_AlignDataFormat(PyObject *format) /* take a data-format object and construct an aligned data-format object where all fields are aligned on appropriate boundaries for the compiler */ Discussion The information provided in the array view object is patterned after the way a multi-dimensional array is defined in NumPy -- including the data-format object which allows a variety of descriptions of memory depending on the need. Reference Implementation Supplied when the PEP is accepted. Copyright This document is placed in the public domain.

15 50

Requests for NumPy Ports?
by Travis Oliphant 03 Feb '07

03 Feb '07

I'm trying to help out the conversion to NumPy by offering patches to various third-party packages that have used Numeric in the past. Does anybody here have requests for which packages they would like to see converted to use NumPy? Ultimately, we are at the mercy of the package-owners to complete the conversion by applying the patches. But, I think it helps when there are patches already available for trying out. The packages I've worked on so far are here: http://www.scipy.org/Porting_to_NumPy Thanks, -Travis

7 9

Release 0.6.1 of pyaudio, renamed pyaudiolab
by David Cournapeau 02 Feb '07

02 Feb '07

Hi There, A few months ago, I posted the first release of pyaudio, a python module to give numpy/scipy environment audio file IO capabilities (ala matlab wavread and co). I recently took time to update it significantly, and as several people showed interest in pyaudio recently, I thought this may interest some people here. The main improvements since last public annoucement: - renamed to pyaudiolab to avoid nameclash with other package pyaudio. - ability to seek into audio files and sync files for flushing IO buffers - matlab-like API: wavread and wavwrite - improvements of the API: more similar to numpy conventions, ability to retrieve most meta-data of a file. - some bug fixing (one bug was quite severe, leading to data corruption on ubuntu with python2.5, due to some weirdness in ubuntu python2.5 packaging) - a real documentation If some people manage to use it on something else than linux, I would happy to hear (particularly Mac OS X, which I cannot test myself). Cheers, David ====== pyaudiolab ====== * WHAT FOR ?: The Goal of pyaudiolab is to give to a numpy/scipy environment some basic audio IO facilities (ala sound, wavread, wavwrite of matlab). With pyaudiolab, you should be able to read and write most common audio files from and to numpy arrays. The underlying IO operations are done using libsndfile from Erik Castro Lopo (http://www.mega-nerd.com/libsndfile/), so he is the one who did the hard work. As libsndfile has support for a vast number of audio files (including wav, aiff, but also htk, ircam, and flac, an open source lossless codec), pyaudiolab enables the import from and export to a fairly large number of file formats. * WHERE TO GET ?: http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/pyaudiolab/index.… bzr archive: http://www.ar.media.kyoto-u.ac.jp/members/david/archives/pyaudiolab/pyaudio… # HOW TO USE ?: http://www.ar.media.kyoto-u.ac.jp/members/david/softwares/pyaudiolab/index.…

2 2

Cutting 1.0.2 release
by Travis Oliphant 01 Feb '07

01 Feb '07

I think it's time for the 1.0.2 release of NumPy. What outstanding issues need to be resolved before we do it? Hopefully, we can do it by the end of the week. -Travis

7 9

query to array
by BBands 31 Jan '07

31 Jan '07

Good afternoon, The following works: import pyodbc import numpy as np connection = pyodbc.connect('DSN=DSNname') cursor = connection.cursor() symbol = 'ibm' request = "select to_days(Date), Close from price where symbol = '" + symbol + "' and date > '2006-01-01'" for row in cursor.execute(request): price.append([row[0], row[1]]) priceArray = np.array(price) Is there a better way to do this? Can fromiter be used? Could this query result be put into a wider array from the start such that columns 2... could be filled with calcs later? Thanks in advance, jab -- John Bollinger, CFA, CMT www.BollingerBands.com If you advance far enough, you arrive at the beginning.

3 4

numpy.scipy.org refers to sourceforge mailing list
by Sebastian Haase 31 Jan '07

31 Jan '07

Hi, On the numpy page: http://numpy.scipy.org/ the is this text: """ Questions? Ask them at the numpy-discussion(a)lists.sourceforge.net mailing list """ Is sourceforge.net still used for numpy ? Also I think it should point to the list sign-up-page since only members are allowed to posts to the list - is this correct ? - Sebastian.

1 0