Mailman 3 January 2009 - NumPy-Discussion

example reading binary Fortran file
by David Froger 13 Feb '13

13 Feb '13

Hy, My question is about reading Fortran binary file (oh no this question again...) Until now, I was using the unpack module like that : def lread(f,fourBeginning,fourEnd,*tuple): from struct import unpack """Reading a Fortran binary file in litte-endian""" if fourBeginning: f.seek(4,1) for array in tuple: for elt in xrange(array.size): transpose(array).flat[elt] = unpack(array.dtype.char,f.read(array.itemsize))[0] if fourEnd: f.seek(4,1) After googling, I read that fopen and npfille was deprecated, and we should use numpy.fromfile and ndarray.tofile, but despite of the documentaion, the cookbook, the mailling list and google, I don't succed in making a simple example. Considering the simple Fortran code below what is the Python script to read the four arrrays? What about if my pc is little endian and the file big endian? I think it will be a good idea to put the Fortran writting-arrays code and the Python reading-array script in the cookbook and maybe a page to help people comming from Fortran to start with Python ? Best, David Froger program makeArray implicit none integer,parameter:: nx=10,ny=20 real(4),dimension(nx,ny):: ux,uy,p integer :: i,j open(11,file='uxuyp.bin',form='unformatted') do i = 1,nx do j = 1,ny ux(i,j) = real(i*j) uy(i,j) = real(i)/real(j) p (i,j) = real(i) + real(j) enddo enddo write(11) ux,uy write(11) p close(11) end program makeArray

8 22

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

Test failures on 2.6
by T J 19 Nov '10

19 Nov '10

Hi, I'm getting a couple of test failures with Python 2.6, Numpy 1.2.0, Nose 0.10.4: nose version 0.10.4 ..........................................................................................................................................................................................................................................................................................................................................................................................................................................F................K........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................./share/home/me/usr/lib/python2.6/site-packages/numpy/lib/tests/test_io.py:68: SyntaxWarning: assertion is always true, perhaps remove parentheses? assert(c.readlines(), ........................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................................./share/home/me/usr/lib/python2.6/site-packages/numpy/ma/tests/test_core.py:1315: SyntaxWarning: assertion is always true, perhaps remove parentheses? assert(store._mask, True) /home/me/usr/lib/python2.6/site-packages/numpy/ma/tests/test_core.py:1322: SyntaxWarning: assertion is always true, perhaps remove parentheses? assert(store._mask, True) /home/me/usr/lib/python2.6/site-packages/numpy/ma/tests/test_core.py:1989: SyntaxWarning: assertion is always true, perhaps remove parentheses? assert(test.mask, [0,1,0,0,0,0,0,0,0,0]) ...............................................E................................................................................................................................................................................ ====================================================================== ERROR: Tests the min/max functions with explicit outputs ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/me/usr/lib/python2.6/site-packages/numpy/ma/tests/test_core.py", line 653, in test_minmax_funcs_with_output result = npfunc(xm,axis=0,out=nout) File "/home/me/usr/lib/python2.6/site-packages/numpy/core/fromnumeric.py", line 1525, in amin return amin(axis, out) File "/home/me/usr/lib/python2.6/site-packages/numpy/ma/core.py", line 2978, in min np.putmask(out, newmask, np.nan) ValueError: cannot convert float NaN to integer ====================================================================== FAIL: test_umath.TestComplexFunctions.test_against_cmath ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/me/usr/lib/python2.6/site-packages/nose-0.10.4-py2.6.egg/nose/case.py", line 182, in runTest self.test(*self.arg) File "/home/me/usr/lib/python2.6/site-packages/numpy/core/tests/test_umath.py", line 268, in test_against_cmath assert abs(a - b) < atol, "%s %s: %s; cmath: %s"%(fname,p,a,b) AssertionError: arcsin 2: (1.57079632679-1.31695789692j); cmath: (1.57079632679+1.31695789692j) ---------------------------------------------------------------------- Ran 1726 tests in 8.856s FAILED (KNOWNFAIL=1, errors=1, failures=1) <nose.result.TextTestResult run=1726 errors=1 failures=1>

2 1

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

Is anyone knowledgeable about dll deployment on windows ?
by David Cournapeau 01 Dec '09

01 Dec '09

Hi, My head hurt trying to understand dll management with windows. I wanted to find a sane way to use dll for numscons, but I can't see how to do this, so I was wondering if anyone on this ML had any deep knowledge on how to install dll, and reuse them with python extensions ? The problem is the following: - you have some dll (say MKL) in a path (say C:\program files\MKL). - how do you tell another dll (more exactly a python extension, a .pyd) to look for mkl dll into the MKL path ? It seems that there is no rpath-like capability on windows, and the only way to load dll is to put them into some windows directory, or to put the dll path into PATH. Is it really the only way ? cheers, David

8 25

inplace matrix multiplication
by Frédéric Bastien 29 Apr '09

29 Apr '09

Hi, I would like to know how I can make a call to the blas function gemm in numpy. I need a multiply and accumulate for matrix and I don't want to allocate a new matrix each time I do it. thanks for your time Frederic Bastien

4 6

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

Automatic differentiation (was Re: second-order gradient)
by Rob Clewley 11 Mar '09

11 Mar '09

> Maybe we should focus on writing a decent 'deriv' function then. I > know Konrad Hinsen's Scientific had a derivatives package > (Scientific.Functions.Derivatives) that implemented automatic > differentiation: > > http://en.wikipedia.org/wiki/Automatic_differentiation That would be great, but wouldn't that be best suited as a utility requiring Sympy? You'll want to take advantage of all sorts of symbolic classes, especially for any source code transformation approach. IMO Hinsen's implementation isn't a very efficient or attractive solution to AD given the great existing C/C++ codes out there. Maybe we should be looking to provide a python interface to an existing open source package such as ADOL-C, but I'm all in favour of a new pure python approach too. What would be perfect is to have a single interface to a python AD package that would support a faster implementation if the user wished to install a C/C++ package, otherwise would default to a pure python equivalent. -Rob

4 4

N-D array interface page is out of date
by Andrew Straw 09 Mar '09

09 Mar '09

Hi, I just noticed that the N-D array interface page is outdated and doesn't mention the buffer interface that is standard with Python 2.6 and Python 3.0: http://numpy.scipy.org/array_interface.shtml This page is linked to from http://numpy.scipy.org/ I suggest, at the minimum, modifying the page with really annoying blinking red letters at the top (or other suitable warning) that this is deprecated at that people should use http://www.python.org/dev/peps/pep-3118/ instead.

4 8

NumPy-Discussion January 2009