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

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 27 Apr '08

27 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

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

05 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

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

23 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

Re: [Numpy-discussion] [C++-sig] Overloading sqrt(5.5)*myvector
by Bruce Sherwood 08 Jan '08

08 Jan '08

Roman Yakovenko wrote: > On Dec 29, 2007 7:47 AM, Bruce Sherwood <Bruce_Sherwood(a)ncsu.edu> wrote: > >> I realized belatedly that I should upgrade from Boost 1.33 to 1.34. >> Alas, that didn't cure my problem. >> > Can you post small and complete example of what you are trying to achieve? > I don't have a "small and complete" example available, but I'll summarize from earlier posts. VPython (vpython.org) has its own vector class to mimic the properties of 3D vectors used in physics, in the service of easy creation of 3D animations. There is a beta version which imports numpy and uses it internally; the older production version uses Numeric. Boost python and thread libraries are used to connect the C++ VPython code to Python. There is operator overloading that includes scalar*vector and vector*scalar, both producing vector. With Numeric, sqrt produced a float, which was a scalar for the operator overloading. With numpy, sqrt produces a numpy.float64 which is caught by vector*scalar but not by scalar*vector, which means that scalar*vector produces an ndarray rather than a vector, which leads to a big performance hit in existing VPython programs. The overloading and Boost code is the same in the VPython/Numeric and VPython/numpy versions. I don't know whether the problem is with numpy or with Boost or with the combination of the two. Here is the relevant part of the vector class: inline vector operator*( const double s) const throw() { return vector( s*x, s*y, s*z); } and here is the free function for right multiplication: inline vector operator*( const double& s, const vector& v) { return vector( s*v.x, s*v.y, s*v.z); } Maybe the problem is in the Boost definitions: py::class_<vector>("vector", py::init< py::optional<double, double, double> >()) .def( self * double()) .def( double() * self) Left multiplication is fine, but right multiplication isn't. A colleague suggested the following Boost declarations but cautioned that he wasn't sure of the syntax for referring to operator, and indeed this doesn't compile: .def( "__mul__", &vector::operator*(double), "Multiply vector times scalar") .def( "__rmul__", &operator*(const double&, const vector&), "Multiply scalar times vector") I would really appreciate a Boost or numpy expert being able to tell me what's wrong (if anything) with these forms. However, I may have a useful workaround as I described in a post to the numpy discussion list. A colleague suggested that I do something like this for sqrt and other such mathematical functions: def sqrt(x): try: return mathsqrt(x) except TypeError: return numpysqrt(x) That is, first try the simple case of a scalar argument, handled by the math module sqrt, and only use the numpy sqrt routine in the case of an array argument. Even with the overhead of the try/except machinery, one gets must faster square roots for scalar arguments this way than with the numpy sqrt. Bruce Sherwood

5 6

Numpy or Boost problem
by Bruce Sherwood 07 Jan '08

07 Jan '08

I'm not sure whether this is a Numpy problem or a Boost problem, so I'm posting to both communities. In old Numeric, type(sqrt(5.5)) was float, but in numpy, type(sqrt(5.5)) is numpy.float64. This leads to a big performance hit in calculations in a beta version of VPython, using the VPython 3D "vector" class, compared with the old version that used Numeric (VPython is a 3D graphics module for Python; see vpython.org). Operator overloading of the VPython vector class works fine for vector*sqrt(5.5) but not for sqrt(5.5)*vector. The following free function catches 5.5*vector but fails to catch sqrt(5.5)*vector, whose type ends up as numpy.ndarray instead of the desired vector, with concomitant slow conversions in later vector calculations: inline vector operator*( const double& s, const vector& v) { return vector( s*v.x, s*v.y, s*v.z); } I've thrashed around on this, including trying to add this: inline vector operator*( const npy_float64& s, const vector& v) { return vector( s*v.x, s*v.y, s*v.z); } But the compiler correctly complains that this is in conflict with the version of double*vector, since in fact npy_float64 is actually double. It's interesting and presumably meaningful to the knowledgeable (not me) that vector*sqrt(5.5) yields a vector, even though the overloading speaks of double, not a specifically numpy name: inline vector operator*( const double s) const throw() { return vector( s*x, s*y, s*z); } VPython uses Boost, and the glue concerning vectors includes the following: py::class_<vector>("vector", py::init< py::optional<double, double, double> >()) .def( self * double()) .def( double() * self) As far as I can understand from the Boost Python documentation, this is the proper way to specify the left-hand and right-hand overloadings. But do I have to add something like .def( npy_float64() * self)? Help would be much appreciated. Bruce Sherwood

5 5

NumPy-Discussion December 2007