Mailman 3 July 2005 - 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

Suggested Addition to the new array PEP
by Chris Barker 12 Aug '05

12 Aug '05

Hi all, but particularly Travis, I just took the time to read the new array PEP: http://numeric.scipy.org/PEP.txt It looks great. I do have one suggestion for an addition method: change_type(NewType) # or some other name: Change the type that the bytes in the array represent, without changing the data at. Dimensions will change if the new type has a different itemsize than the old. This would be the equivalent of: A = Numeric.fromstring(A.tostring(), NewType) but without the extra copies of the data. When would anyone want this? I've written a function that does exactly this with Numeric arrays in order to efficiently slice, dice and re-arange a bunch of binary data files. The files consist of a bunch of C structs written sequentially to a file. In this case the structs are: (long, long, char) I have a file which is N of these structs written sequentially to a file, I can do: NumBytes = 9 # number of bytes in a record # Read the file data = fromstring(file.read(NumBytes*NumTimesteps*NumLEs),UnsignedInt8) file.close() data.shape = (NumTimesteps,NumLEs,NumBytes) # extract LE data: LEs = data[:,:,:8].copy() # extract the flags flags = data[:,:,8].copy() del data # here I can create the Int32 array without making a copy: # changetype returns a rank-1 array changetype(LEs,Int32) # It now gets reshaped to fit what it represents. LEs.shape = (NumTimesteps,NumLEs,2) I think it would only work in the general case with a contiguous array. All this really requires is to re-set the attributes of the array. with Numeric, the code is (for contiguous arrays): NewDescription = PyArray_DescrFromType(typecode); NewElementSize = NewDescription->elsize; if (TotalBytes % NewElementSize > 0){ PyErr_SetString (PyExc_ValueError, "The input array is the wrong size for the requested type"); return NULL; } Array->descr = NewDescription; // does the old descr need to be freed?? how?? Array->nd = 1; Array->strides[0] = NewElementSize; Array->dimensions[0] = TotalBytes / NewElementSize; The core By the way, I'm also not seeing tostring() or tobuffer() listed in the PEP -Chris -- Christopher Barker, Ph.D. Oceanographer NOAA/OR&R/HAZMAT (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker(a)noaa.gov

2 1

numarray: PyArrayObject supports number protocols?
by Edward C. Jones 04 Aug '05

04 Aug '05

I just saw in the docs, Section 14.2.4: It should be noted that unlike earlier versions of numarray, the present PyArrayObject structure is a first class python object, with full support for the number protocols in C. Does this mean that I can add two numarrays in C using "PyNumber_Add"?

3 2

Re: [Numpy-discussion] LAPACK/BLAS for Fedora core 4 i386 ?
by Chris Barker 29 Jul '05

29 Jul '05

Jochen Küpper wrote: > "Chris Barker" <Chris.Barker(a)noaa.gov> writes: >>Do you mean this? : > yes OK, so I've done what I can. > To find out you might have to look at the LAPACK functions you are > using and compare them against the ones provided by ATLAS. Well, I'm testing using LinearAlgebra.solve_linear_equations() from Numeric (and numarray). The reason I'm wondering is that on my last machine, a 1ghz PIII laptop running Gentoo, using the Gentoo provided atlas/lapack, I got a 6-7 times speed up over lapack_lite. On my new machine, a 2GHz Pentium M laptop running FC4, I got about a 2 times speed-up, which is nice, but not nearly as impressive. Another possible issue is that I used pre-complied binaries form the atlas site, which are a bit old, maybe I should compile myself. Any thoughts? -Chris -- Christopher Barker, Ph.D. Oceanographer NOAA/OR&R/HAZMAT (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker(a)noaa.gov

4 4

New e-mail address
by Pieter van Beek 28 Jul '05

28 Jul '05

*English follows below* Geachte relatie, Ooit (misschien recent, misschien al weer jaren geleden) heeft u e-mail verstuurd naar mijn voormalige e-mail adres: _vanbeek(a)q-factor.nl_. Bij deze wil ik u mededelen dat mijn e-mail adres is gewijzigd in pieter(a)djinnit.com. Gelieve deze wijziging door te voeren in uw adresbestand. Met vriendelijke groet, Pieter van Beek ------------------------------------------------------------------------ Dear relation, Once (maybe recently, maybe years ago) you have been sending mail to my former e-mail address: _vanbeek(a)q-factor.nl_. Hereby I'd like to inform you that my e-mail address has changed into pieter(a)djinnit.com. Please change my e-mail address in your address book. Kind regards, Pieter van beek

1 0

Interfacing to other array classes
by Joachim Dahl 28 Jul '05

28 Jul '05

I am one the authors of 'cvxopt', a Python package that includes simple matrix classes aimed at numerical linear algebra for convex solvers, and I saw the discussion on the list regarding interfacing to such classes. I think the idea of an advanced interface or buffer protocol to exchange information between simple array types without using (or knowing) the API of the different classes is great. What is the status of the new array interface? When it is ready for testing, I will definitely try it out. - Joachim

1 0

LAPACK/BLAS for Fedora core 4 i386 ?
by Chris Barker 28 Jul '05

28 Jul '05

Hi all, Can anyone recommend the best way to get a native LAPACK installed on Fedora core 4? I'm quite surprised that I can't find an atlas rpm with yum. While I'm at it, when I was using Gentoo, it came with a nice atlas blas/lapack. Unfortunately, the atlas lapack does not include all of lapack. It has everything Numeric/numarray needs, but not some stuff i need for another project (Banded Matrix solvers). Does anyone have a suggestion for how to add the addional lapack stuff I need, while still using atlas stuff where possible? Frankly, I don't understand why atlas doesn't just include all of lapack, using generic versions for anything they haven't optimized, it would be a lot easier to get one stop shopping. -Chris

3 4

Re[2]: [Numpy-discussion] convolve
by Alan G Isaac 27 Jul '05

27 Jul '05

On Wed, 27 Jul 2005, René Bastian apparently wrote: > Outcommenting > if len(data) < len(kernel): > data, kernel = kernel, data > in convolve.py gives the correct result. The author is not listed for this numarray module, so we cannot directly ask why this behavior was chosen. You might use Numeric.convolve instead of numarray.convolve.Convolve.convolve: Numeric gives the expected result here. fwiw, Alan Isaac

3 7

design issues - octave 'incompatibilities'
by Soeren Sonnenburg 27 Jul '05

27 Jul '05

Dear all, I am new to numarray and as I am trying to use it day-by-day I am now wondering about certain numeric/numarray design issues. Please forgive me / correct me as I probably misunderstood certain issues. -- Why did you choose row-major instead of column major format as practiced by R/octave/matlab... Doesn't that mean performance problems as fortran/blas is optimized if you work with the transposed version ? -- why do vectors have no 'orientation', i.e. there are only row but no column vectors (or why do you treat matrices/vectors differently, i.e. have vectors at all as a separate type) numarray: a=array([[1,2,3],[4,5,6],[7,8,9]]) >>> a[0,:] array([1, 2, 3]) >>> a[:,0] array([1, 4, 7]) vs. octave: a=[1,2,3;4,5,6;7,8,9]; >> a(1,:) ans = 1 2 3 >> a(:,1) ans = 1 4 7 -- How can one obtain submatrices from full matrices: numarray gives only elements (which is very, very rarely needed and should IMHO rather be some extra function): >>> i=[0,1] >>> j=[0,2] >>> a[i,j] array([1, 6]) vs octave: >> i=[1,2]; >> j=[1,3]; >> a(i,j) ans = 1 3 4 6 all one can do in numarray currently is as awkward as this: >>> a.take(i,0).take(j,1) array([[1, 3], [4, 6]]) also mixing slice and variable does not work a[:,i] - hmmhh seems that is a major issue for me... -- why don't you allow iterating over the whole matrix via a single index ? i.e. a[3] Traceback (most recent call last): File "<stdin>", line 1, in ? IndexError: Index out of range vs octave >> a(4) ans = 2 Are there more elegant ways to do this ? Which issues are likely to be fixed in the future ? Best regards, Soeren

6 19

2-D FFT input/ouput storage
by Daniel Sheltraw 27 Jul '05

27 Jul '05

Hello All Now that I have Numeric loading properly at the Python command line I went in search of some docs on the 2-D FFT. How are the input and output to the forward and inverse 2-D FFTs stored? Is it simply a 2-D version of the 1-D storage? When I do "print FFT.inverse_fft2d.__doc__" I get no helpful info about storage. Thanks, Daniel

2 1