Mailman 3 February 2006 - 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

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

Multiple inheritance from ndarray
by Robert Lupton 18 Apr '06

18 Apr '06

I have a swig extension that defines a class that inherits from both a personal C-coded image struct (actImage), and also from Numeric's UserArray. This works very nicely, but I thought that it was about time to upgrade to numpy. The code looks like: from UserArray import * class Image(UserArray, actImage): def __init__(self, *args): actImage.__init__(self, *args) UserArray.__init__(self, self.getArray(), 'd', copy=False, savespace=False) I can't figure out how to convert this to use ndarray, as ndarray doesn't seem to have an __init__ method, merely a __new__. So what's the approved numpy way to handle multiple inheritance? I've a nasty idea that this is a python question that I should know the answer to, but I'm afraid that I don't... R

4 8

Numpy and PEP 343
by Tim Hochberg 10 Mar '06

10 Mar '06

<pie-in-the-sky> An idea that has popped up from time to time is delaying evalution of a complicated expressions so that the result can be computed more efficiently. For instance, the matrix expression: a = b*c + d*e results in the creation of two, potentially large, temporary matrices and also does a couple of extra loops at the C level than the equivalent expression implemented in C would. The general idea has been to construct some sort of psuedo-object, when the numerical operations are indicated, then do the actual numerical operations at some later time. This would be very problematic if implemented for all arrays since it would quickly become impossible to figure out what was going on, particularly with view semantics. However, it could result in large performance improvements without becoming incomprehensible if implemented in small enough chunks. A "straightforward" approach would look something like: numpy.begin_defer() # Now all numpy operations (in this thread) are deferred a = b*c + d*e # 'a' is a special object that holds pointers to # 'b', 'c', 'd' and 'e' and knows what ops to perform. numpy.end_defer() # 'a' performs the operations and now looks like an array Since 'a' knows the whole series of operations in advance it can perform them more efficiently than would be possible using the basic numpy machinery. Ideally, the space for 'a' could be allocated up front, and all of the operations could be done in a single loop. In practice the optimization might be somewhat less ambitious, depending on how much energy people put into this. However, this approach has some problems. One is the syntax, which clunky and a bit unsafe (a missing end_defer in a function could cause stuff to break very far away). The other is that I suspect that this sort of deferred evaluation makes multiple views of an array even more likely to bite the unwary. The syntax issue can be cleanly addressed now that PEP 343 (the 'with' statement) is going into Python 2.5. Thus the above would look like: with numpy.deferral(): a = b*c + d*e Just removing the extra allocation of temporary variables can result in 30% speedup for this case[1], so the payoff would likely be large. On the down side, it could be quite a can of worms, and would likely require a lot of work to implement. Food for thought anyway. </pie-in-the-sky> -tim [1] from timeit import Timer print Timer('a = b*c + d*e', 'from numpy import arange;b=c=d=e=arange(100000.)').timeit(10000) print Timer('a = b*c; multiply(d,e,temp); a+=temp', 'from numpy import arange, zeros, multiply;' 'b=c=d=e=arange(100000.);temp=zeros([100000], dtype=float)').timeit(10000) => 94.8665989672 62.6143562939

13 33

Matrix times scalar is wacky
by Bill Baxter 08 Mar '06

08 Mar '06

Multiplying a matrix times a scalar seems to return junk for some reason: >>> A = numpy.asmatrix(numpy.rand(1,2)) >>> A matrix([[ 0.30604211, 0.98475225]]) >>> A * 0.2 matrix([[ 6.12084210e-002, 7.18482614e-290]]) >>> 0.2 * A matrix([[ 6.12084210e-002, 7.18482614e-290]]) >>> numpy.__version__ '0.9.5' --billyb

6 9

Thoughts on an ndarray super-class
by Travis Oliphant 05 Mar '06

05 Mar '06

The bigndarray class is going to disappear (probably in the next release of NumPy). It was a stop-gap measure as the future of 64-bit fixes in Python was unclear. Python 2.5 will have removed the 64-bit limitations that led to the bigndarray and so it will be removed. I have been thinking, however, of replacing it with a super-class that does not define the dimensions or strides. In other words, the default array would be just a block of memory. The standard array would inherit from the default and add dimension and strides pointers. I was thinking that this might make it easier for sub-classes using fixed-sized dimensions and strides. I'm not sure if that would actually be useful, but since I was thinking about the disappearance of the bigndarray, I thought I would ask for comments. -Travis

4 4

Matlab page on scipy wiki
by Bill Baxter 03 Mar '06

03 Mar '06

I added some content to the "NumPy/SciPy for Matlab users" page on the scipy wiki. But my knowledge of NumPy/SciPy isn't sufficient to fill in the whole chart of equivalents that I laid out. If folks who know both could browse by and maybe fill in a blank or two, that would be great. I think this will be a helpful "getting started" page for newbies to NumPy coming from matlab, like me. One of the most frustrating things is when you sit down and can't figure out how to do the most basic things that do in your sleep in another environment (like making a column vector). So hopefully this page will help. The URL is : http://www.scipy.org/Wiki/NumPy_for_Matlab_Addicts Thanks, Bill Baxter

17 39

setting package path
by Stefan van der Walt 02 Mar '06

02 Mar '06

In numpytest.py, set_package_path is provided for handling path changes while doing unit tests. It reads def set_package_path(level=1): """ Prepend package directory to sys.path. set_package_path should be called from a test_file.py that satisfies the following tree structure: <somepath>/<somedir>/test_file.py Then the first existing path name from the following list <somepath>/build/lib.<platform>-<version> <somepath>/.. is prepended to sys.path. ... However, the line that supposedly calculates "somepath/.." is d = os.path.dirname(os.path.dirname(os.path.abspath(testfile))) which calculates "somepath". Which is wrong: the docstring, the code or my interpretation? Stéfan

2 2

Re: [SciPy-user] Table like array
by Travis Oliphant 01 Mar '06

01 Mar '06

Michael Sorich wrote: > Hi, > > I am looking for a table like array. Something like a 'data frame' > object to those familiar with the statistical languages R and Splus. > This is mainly to hold and manipulate 2D spreadsheet like data, which > tends to be of relatively small size (compared to what many people > seem to use numpy for), heterogenous, have column and row names, and > often contains missing data. You could subclass the ndarray to produce one of these fairly easily, I think. The missing data item could be handled by a mask stored along with the array (or even in the array itself). Or you could use a masked array as your core object (though I'm not sure how it handles the arbitrary (i.e. record-like) data-types yet). Alternatively, and probably the easiest way to get started, you could just create your own table-like class and use simple 1-d arrays or 1-d masked arrays for each of the columns --- This has always been a way to store record-like tables. It really depends what you want the data-frames to be able to do and what you want them to "look-like." > A RecArray seems potentially useful, as it allows different fields to > have different data types and holds the name of the field. However it > doesn't seem easy to manipulate the data. Or perhaps I am simply > having difficulty finding documentation on there features. Adding a new column/field means basically creating a new array with a new data-type and copying data over into the already-defined fields. Data-types always have a fixed number of bytes per item. What those bytes represent can be quite arbitrary but it's always fixed. So, it is always "more work" to insert a new column. You could make that seamless in your table class so the user doesn't see it though. You'll want to thoroughly understand the dtype object including it's attributes and methods. Particularly the fields attribute of the dtype object. > eg > adding a new column/field (and to a lesser extent a new row/record) to > the recarray Adding a new row or record is actually similar because once an array is created it is usually resized by creating another array and copying the old array into it in the right places. > Changing the field/column names > make a new table by selecting a subset of fields/columns. (you can > select a single field/column, but not multiple). Right. So far you can't select multiple columns. It would be possible to add this feature with a little-bit of effort if there were a strong demand for it, but it would be much easier to do it in your subclass and/or container class. 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? > It would also be nice for the table to be able to deal easily with > masked data (I have not tried this with recarray yet) and perhaps also > to be able to give the rows/records unique ids that could be used to > select the rows/records (in addition to the row/record index), in the > same way that the fieldnames can select the fields. Adding fieldnames to the "rows" is definitely something that a subclass would be needed for. I'm not sure how you would even propose to select using row names. Would you also use getitem semantics? > Can anyone comment on this issue? Particularly whether code exists for > this purpose, and if not ideas about how best to go about developing > such a Table like array (this would need to be limited to python > programing as my ability to program in c is very limited). I don't know of code that already exists for this, but I don't think it would be too hard to construct your own data-frame object. I would probably start with an implementation that just used standard arrays of a particular type to represent the internal columns and then handle the indexing using your own over-riding of the __getitem__ and __setitem__ special methods. This would be the easiest to get working, I think. -Travis

3 2

subclassing ndaray
by Colin J. Williams 28 Feb '06

28 Feb '06

I have a subclass Bar, a 1-dim array which has some methods and some attributes. One of the attributes is a view of the Bar to permit different shaping. Suppose that 'a' is an instance of 'Bar', which has a method 'show' and a view attribute 'v'. a ^ 15 returns a Bar instance, with its methods but without the attributes. I am attempt to change this, Bar has a method __xor__, see below: def __xor__(self, other): ''' Exclusive or: __xor__(x, y) => x ^ y . ''' z= 1 << this loops to the recursion limit result= ArrayType.__xor__(self, other) n= self.n result.n= n result.rowSize= self.rowSize result.show= self.show result.v= _n.reshape(result.view(), (n*n, n*n)) return result Could anyone suggest a workaround please? Colin W.

5 10