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

numpy distutils patch
by Jarrod Millman 21 Nov '07

21 Nov '07

Hey, I was looking at numpy/distutils/ccompiler.py and noticed that it has a fix for distutils.util.split_quoted. Here is the relevant code from split_quoted in numpy.distutils.ccompiler: ----------------------------------- def split_quoted(s): <snip> if _has_white_re.search(s[beg+1:end-1]): s = s[:beg] + s[beg+1:end-1] + s[end:] pos = m.end() - 2 else: # Keeping quotes when a quoted word does not contain # white-space. XXX: send a patch to distutils pos = m.end() <snip> ----------------------------------- Here is the relevant code from split_quoted in distutils.util: ----------------------------------- def split_quoted(s): <snip> s = s[:beg] + s[beg+1:end-1] + s[end:] pos = m.end() - 2 <snip> ----------------------------------- Does anyone know if a patch was ever submitted upstream? If not, is there any reason that a patch shouldn't be submitted? Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/

2 3

beginner question: rank-1 arrays
by Robin 16 Nov '07

16 Nov '07

Hi, I am trying to implement a project in scipy. I think I am getting somewhere finally. However in my code (I am converting from MATLAB) it is important to maintain 2d arrays, and keep the difference between row and column vectors. After working through some initial problems I think I am getting more of a picture of how things work in numpy. However I am finding my code littered with things like: np.array(np.r_[0:nterms],ndmin=2) (for a row vector) or np.array(np.r_[0:nterms],ndmin=2).T (for a column vector) Coming from matlab and being use to 0:10 for row or (0:10)' for column this seems a bit messy. Is there a better way of constructing row/column 2d arrays from a slice type range? Thanks, Robin

6 9

vectorizing loops
by Mathew Yeates 02 Nov '07

02 Nov '07

Anybody know of any tricks for handling something like z[0]=1.0 for i in range(100): out[i]=func1(z[i]) z[i+1]=func2(out[i]) ??

10 19

matrix problem: float to matrix power
by Alan G Isaac 01 Nov '07

01 Nov '07

>>> 1.0**numpy.array([1,2,3]) array([ 1., 1., 1.]) >>> 1.0**numpy.mat([1,2,3]) Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: unsupported operand type(s) for ** or pow(): 'float' and 'matrix' Why the restriction for matrices? Same question for matrices conformable for broadcasting: why not produce the equivalent of the array result? Cheers, Alan Isaac

4 7

numpy FFT memory accumulation
by Ray S 31 Oct '07

31 Oct '07

I am using fftRes = abs(fft.rfft(data_array[end-2**15:end])) to do running analysis on streaming data. The N never changes. It sucks memory up at ~1MB/sec with 70kHz data rate and 290 ffts/sec. (Interestingly, Numeric FFT accumulates much slower..) (Commenting out that one line stops memory growth.) What can one do to alleviate this? Can I del() some numpy object or such? It's a bit of an issue for a program that needs to run for weeks. It's purpose is simply to argmax() the largest bin, which always falls within a small range - do I have another, better option than fft? Cheers, Ray Schumacher Blue Cove Interactive

3 2