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

setting the attributes of an array of object
by jelle 16 Oct '07

16 Oct '07

Hi, I'm wondering whether i can re-write the following idiom with numpy arrays: for i in some_list: i.some_attr = some_value it would be wonderful if one was able to write this idiom as arr[just_these].some_attr = some_value or setattr(arr[just_these], 'some_attr', some_value) since often expensive loops through lists of object could be avoided. any thoughts on this are much appreciated, thanks, -jelle

4 3

adding field to rec array
by John Hunter 05 Oct '07

05 Oct '07

I have a record array r and I want to add a new field to it. I have been looking at setfield but I am not sure how to use it for this purpose. Eg # r is some npy record array N = len(r) x = npy.zeros(N) # add array of floats x to r with dtype name 'jdh' and type '<f8' Any suggestions? JDH

3 4

indexing bug?
by Nadia Dencheva 03 Oct '07

03 Oct '07

Is indexing with floats really allowed in numpy? >>> a=numpy.array([1,2,3,4]) >>> a[2.99] 3 Nadia Dencheva

9 12

Re: [Numpy-discussion] indexing bug?
by Nadezhda Dencheva 29 Sep '07

29 Sep '07

>> This should return an error and not silently truncate to int. > >Why do you say that? I see now this was a Friday afternoon lack of imagination :) The current behaviour is consistent and well >defined: > >On a somewhat related note, you may also be interested in the PEP at > >http://docs.python.org/whatsnew/pep-357.html > Interesting. Thanks for explaining this. Nadia

1 0

ANN: maskedarray
by Pierre GM 27 Sep '07

27 Sep '07

All, The latest version of maskedarray has just been released on the scipy SVN sandbox. This version fixes the inconsistencies in filling (see below) and introduces some minor modifications for optimization purposes (see below as well). Many thanks to Eric Firing and Matt Knox for the fruitful discussions at the origin of this release! In addition, a bench.py file has been introduced, to compare the speed of numpy.ma and maskedarray. Once again, thanks to Eric for his first draft. Please feel free to try it and send me some feedback. Modifications: * Consistent filling ! In numpy.ma, the division of array A by array B works in several steps: - A is filled w/ 0 - B is filled w/ 1 - A/B is computed - the output mask is updated as the combination of A.mask, B.mask and the domain mask (B==0) The problems with this approach are that (i) it's not useful to fill A and B beforehand if the values will be masked anyway; (ii) nothing prevents infs to show up, as the domain is taken into account at the end only. In this latest version of maskedarray, the same division is decomposed as: - a copy of B._data is filled with 1 with the domain (B==0) - the division of A._data by this copy is computed - the output mask is updated as the combination of A.mask, B.mask and the domain mask (B==0). Prefilling on the domain avoids the presence of nans/infs. However, this comes with the price of making some functions and methods slower than their numpy.ma counterparts, as you'll be able to observe for sqrt and log with the bench.py file. An alternative would be to avoid filling at all, at the risk of leaving nans and infs. * masked_invalid / fix_invalid Two new functions are introduced. masked_invalid(x) masks x where x is nan or inf. fix_invalid(x) returns (a copy of) x, where invalid values (nans & infs) are replaced by fill_value. * No mask shrinking Following Paul Dubois and Sasha's example, I eventually had to get rid of the semi-automatic shrinking of the mask in __getitem__, which appeared to be a major bottleneck. In other words, one can end up with an array full of False instead of nomask, which may slow things down a bit. You can force a mask back to nomask with the new shrink_mask method. *_sharedmask Here again, I followed Paul and Sasha's ideas and reintroduce the _sharedmask flag to prevent inadequate propagation of the mask. When creating a new array with x=masked_array(data, mask=m), x._mask is initially a reference to m and x._sharedmask is True. When x is modified, x._mask is copied to prevent a propagation back to m.

1 0

matplotlib - representation of nan values in 2D
by Dirk Zickermann 26 Sep '07

26 Sep '07

Dear pylab group, for the represenation of 2d measurement data I use the contourplot function from matplotlib. Some points in this map are not measurabel, therefore I get a non numerical value (nan) output. >From this data I want to generate a map and a histogram plot. This works fine, as long as I use regular matrix/array data structure without any voids. My questions are: (1) How can I make use of plotting data with NAN values as contour and asigns such values eg as black points? (2) How can I use the matplotlib hist() function with this data, that also include such missing data points? Maybe there is an easy workaround for this. Thanks a lot for your support, Dirk (python2.51/matplotlib-0.90.1, win32) my code: ... import matplotlib ... my2dData=[[1,2,3,4,5.0 ,NaN,7,8,9,10],[10,9,8,7,6,5,4,3,2,1]] ... figure(1) imshow(my2dData) pylab.show() ..

2 2