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

Unpleasant behavior with poly1d and numpy scalar multiplication
by Fernando Perez 15 Feb '10

15 Feb '10

Hi all, consider this little script: from numpy import poly1d, float, float32 p=poly1d([1.,2.]) three=float(3) three32=float32(3) print 'three*p:',three*p print 'three32*p:',three32*p print 'p*three32:',p*three32 which produces when run: In [3]: run pol1d.py three*p: 3 x + 6 three32*p: [ 3. 6.] p*three32: 3 x + 6 The fact that multiplication between poly1d objects and numbers is: - non-commutative when the numbers are numpy scalars - different for the same number if it is a python float vs a numpy scalar is rather unpleasant, and I can see this causing hard to find bugs, depending on whether your code gets a parameter that came as a python float or a numpy one. This was found today by a colleague on numpy 1.0.4.dev3937. It feels like a bug to me, do others agree? Or is it consistent with a part of the zen of numpy I've missed thus far? Thanks, f

6 17

Compact way of performing array math with specified result type?
by Russell E. Owen 08 Feb '10

08 Feb '10

I often find myself doing simple math on sequences of numbers (which might or might not be numpy arrays) where I want the result (and thus the inputs) coerced to a particular data type. I'd like to be able to say: numpy.divide(seq1, seq2, dtype=float) but ufuncs don't allow on to specify a result type. So I do this instead: numpy.array(seq1, dtype=float) / numpy.array(seq2, dtype=float) Is there a more compact solution (without having to create the result array first and supply it as an argument)? -- Russell

4 3

Slicing/selection in multiple dimensions simultaneously
by Mike Ressler 17 Mar '09

17 Mar '09

The following seems to be a wart: is it expected? Set up a 10x10 array and some indexing arrays: a=arange(100) a.shape=(10,10) q=array([0,2,4,6,8]) r=array([0,5]) Suppose I want to extract only the "even" numbered rows from a - then print a[q,:] <works - output deleted> Every fifth column: print a[:,r] <works - output deleted> Only the even rows of every fifth column: print a[q,r] --------------------------------------------------------------------------- <type 'exceptions.ValueError'> Traceback (most recent call last) /.../.../.../<ipython console> in <module>() <type 'exceptions.ValueError'>: shape mismatch: objects cannot be broadcast to a single shape But, this works: print a[q,:][:,r] [[ 0 5] [20 25] [40 45] [60 65] [80 85]] So why does the a[q,r] form have problems? Thanks for your insights. Mike -- mike.ressler(a)alum.mit.edu

8 16

Request for advice: project to get NumPy working in IronPython
by Giles Thomas 26 Apr '08

26 Apr '08

Hi, At Resolver Systems, we have a product that is written in IronPython - the .NET Python implementation - and allows users to use that language to script a spreadsheet-like interface. Because they're using IronPython, they can access their existing .NET objects and libraries, which has worked out really well for us and for them. But there's an increasing number of users who would like access to CPython C extensions - in particular, NumPy. We *don't* want to just hack something together that would just allow users of our own product to access NumPy; what we'd really like to do is work on an OSS solution that will allow all IronPython users to access NumPy transparently, just as if they were using CPython. My question is, is this something that the readers of this list would like to get involved in, either in terms of contibuting, or even in terms on advice? We're under no illusions; this is a significant undertaking. And we're very keen on hearing from those who know NumPy better about how this should be structured, even if people don't have the time to write code for it. In particular: * How do people feel about a source-code-compatible solution, where (perhaps) we would maintain a project that basically provided an alternative makefile for NumPy - or, even better, we could work with the NumPy developers to contribute a .NET/IronPython package. * Would it be better to try for some kind of "binary compatibility", where we'd write some kind of "glue" that sat between the existing C extension .pyd files and the IronPython engine? Our gut feeling is that this would be much more work, but we might be missing something. * What should the work's relationship to the NumPy project be? Any thoughts, comments, or advice would be very much appreciated. Regards, Giles -- Giles Thomas MD & CTO, Resolver Systems Ltd. giles.thomas(a)resolversystems.com +44 (0) 20 7253 6372 We're hiring! http://www.resolversystems.com/jobs/ 17a Clerkenwell Road, London EC1M 5RD, UK VAT No.: GB 893 5643 79 Registered in England and Wales as company number 5467329. Registered address: 843 Finchley Road, London NW11 8NA, UK

12 16

ndarray.__mul__ is too gready?
by Pearu Peterson 24 Apr '08

24 Apr '08

Hi, Say, one defines a class A that does not like to have numpy arrays in the left-hand-side of an operation, say *, but in the rhs it is expected. For an example, consider A as a huge customized matrix that defines only matrix * vector where vector can be numpy array and scalar * matrix where scalar is python float or int. Users may accidentally hit vector * matrix and this should raise an exception. But when vector is numpy array, then ndarray.__mul__ calls element * A for each array element and here A.__rmul__ gets only the types of array elements (that are floats, for instance) as rhs and hence is not able to decide that it should return NotImplemented. Instead, it computes scalar * matrix and the final result will be array([element1*matrix, element2*matrix, ..]) which is not desired already because of the memory usage and should have raised an exception instead. Here is a small example illustrating the problem: class A: def __rmul__(self, other): if isinstance(other, ndarray): return NotImplemented if isinstance(other, (int, float)): return self # pretend other==1 for an example return NotImplemented def __repr__(self): return 'A()' >>> array([1,2,3]) * A() array([A(), A(), A()], dtype=object) # expected TypeError Is there any workaround to this problem? In other situations/applications the above behavior of ndarray.__mul__ is a feature but is there any way to avoid this feature? Thanks, Pearu

2 1

Problem with numpy.linalg.eig?
by Peter Creasey 04 Apr '08

04 Apr '08

Hi all, The following code calling numpy v1.0.4 fails to terminate on my machine, which was not the case with v1.0.3.1 from numpy import arange, float64 from numpy.linalg import eig a = arange(13*13, dtype = float64) a.shape = (13,13) a = a%17 eig(a) Regards, Peter

6 19

planet.scipy.org
by Jarrod Millman 01 Apr '08

01 Apr '08

Hey, I just wanted to announce that we now have a NumPy/SciPy blog aggregator thanks to Gaël Varoquaux: http://planet.scipy.org/ Feel free to contact me if you have a blog that you would like included. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

12 24

Problems with long
by Tom Johnson 12 Mar '08

12 Mar '08

Hi, I'm having some troubles with long. >>> from numpy import log >>> log(8463186938969424928L) 43.5822574833 >>> log(10454852688145851272L) <type 'exceptions.AttributeError'>: 'long' object has no attribute 'log' Thoughts?

3 3

Faster array version of ndindex
by Jonathan Taylor 22 Feb '08

22 Feb '08

I was needing an array representation of ndindex since ndindex only gives an iterator but array(list(ndindex)) takes too long. There is prob some obvious way to do this I am missing but if not feel free to include this code which is much faster. In [252]: time a=np.array(list(np.ndindex(10,10,10,10,10,10))) CPU times: user 11.61 s, sys: 0.09 s, total: 11.70 s Wall time: 11.82 In [253]: time a=ndtuples(10,10,10,10,10,10) CPU times: user 0.32 s, sys: 0.21 s, total: 0.53 s Wall time: 0.60 def ndtuples(*dims): """Fast implementation of array(list(ndindex(*dims))).""" # Need a list because we will go through it in reverse popping # off the size of the last dimension. dims = list(dims) # N will keep track of the current length of the indices. N = dims.pop() # At the beginning the current list of indices just ranges over the # last dimension. cur = np.arange(N) cur = cur[:,np.newaxis] while dims != []: d = dims.pop() # This repeats the current set of indices d times. # e.g. [0,1,2] -> [0,1,2,0,1,2,...,0,1,2] cur = np.kron(np.ones((d,1)),cur) # This ranges over the new dimension and 'stretches' it by N. # e.g. [0,1,2] -> [0,0,...,0,1,1,...,1,2,2,...,2] front = np.arange(d).repeat(N)[:,np.newaxis] # This puts these two together. cur = np.column_stack((front,cur)) N *= d return cur

4 4