Mailman 3 February 2008 - 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

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

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

failure building numpy using icc
by Glen W. Mabey 15 May '08

15 May '08

Hello, I'm using svn numpy and get the following error upon executing /usr/local/bin/python2.5 setup.py config --noisy --cc=/opt/intel/cce/10.0.025/bin/icc --compiler=intel --fcompiler=intel build_clib build_ext I see: conv_template:> build/src.linux-x86_64-2.5/numpy/core/src/scalartypes.inc Traceback (most recent call last): File "setup.py", line 96, in <module> setup_package() File "setup.py", line 89, in setup_package configuration=configuration ) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/core.py", line 184, in setup return old_setup(**new_attr) File "/usr/local/lib/python2.5/distutils/core.py", line 151, in setup dist.run_commands() File "/usr/local/lib/python2.5/distutils/dist.py", line 974, in run_commands self.run_command(cmd) File "/usr/local/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/command/build_ext.py", line 56, in run self.run_command('build_src') File "/usr/local/lib/python2.5/distutils/cmd.py", line 333, in run_command self.distribution.run_command(command) File "/usr/local/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/command/build_src.py", line 130, in run self.build_sources() File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/command/build_src.py", line 147, in build_sources self.build_extension_sources(ext) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/command/build_src.py", line 252, in build_extension_sources sources = self.template_sources(sources, ext) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/command/build_src.py", line 359, in template_sources outstr = process_c_file(source) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/conv_template.py", line 185, in process_file % (sourcefile, process_str(''.join(lines)))) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/conv_template.py", line 150, in process_str newstr[sub[0]:sub[1]], sub[4]) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/conv_template.py", line 117, in expand_sub % (line, template_re.sub(namerepl, substr))) File "/home/gmabey/src/DiamondBack/Diamondback/src/numpy-20080228_svn/numpy/distutils/conv_template.py", line 113, in namerepl return names[name][thissub[0]] KeyError: 'PREFIX' And I do not see any errors when building the same svn version with gcc (on a different machine). I've unsuccessfully tried to follow that backtrace of functions to figure out exactly what is going on. Any hints/suggestions? Thanks, Glen Mabey

3 9

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

ndarray.__mul__ is too gready?
by Pearu Peterson 24 Apr '08

24 Apr '08

Hi, Say, one defines a class A that does not like to have numpy arrays in the left-hand-side of an operation, say *, but in the rhs it is expected. For an example, consider A as a huge customized matrix that defines only matrix * vector where vector can be numpy array and scalar * matrix where scalar is python float or int. Users may accidentally hit vector * matrix and this should raise an exception. But when vector is numpy array, then ndarray.__mul__ calls element * A for each array element and here A.__rmul__ gets only the types of array elements (that are floats, for instance) as rhs and hence is not able to decide that it should return NotImplemented. Instead, it computes scalar * matrix and the final result will be array([element1*matrix, element2*matrix, ..]) which is not desired already because of the memory usage and should have raised an exception instead. Here is a small example illustrating the problem: class A: def __rmul__(self, other): if isinstance(other, ndarray): return NotImplemented if isinstance(other, (int, float)): return self # pretend other==1 for an example return NotImplemented def __repr__(self): return 'A()' >>> array([1,2,3]) * A() array([A(), A(), A()], dtype=object) # expected TypeError Is there any workaround to this problem? In other situations/applications the above behavior of ndarray.__mul__ is a feature but is there any way to avoid this feature? Thanks, Pearu

2 1

PCA on set of face images
by devnew＠gmail.com 14 Apr '08

14 Apr '08

hi guys I have a set of face images with which i want to do face recognition using Petland's PCA method.I gathered these steps from their docs 1.represent matrix of face images data 2.find the adjusted matrix by substracting the mean face 3.calculate covariance matrix (cov=A* A_transpose) where A is from step2 4.find eigenvectors and select those with highest eigenvalues 5.calculate facespace=eigenvectors*A when it comes to implementation i have doubts as to how i should represent the matrix of face images? using PIL image.getdata() i can make an array of each greyscale image. Should the matrix be like each row contains an array representing an image? That will make a matrix with rows=numimages and columns=numpixels cavariancematrix =A *A_transpose will create a square matrix of shape(numimages,numimages) Using numpy.linalg.eigh(covariancematrix) will give eigenvectors of same shape as the covariance matrix. I would like to know if this is the correct way to do this..I have no big expertise in linear algebra so i would be grateful if someone can confirm the right way of doing this RoyG

4 8

numpyx.pyx (recent svn) works?
by Neal Becker 09 Apr '08

09 Apr '08

I tried numpyx.pyx with cython-0.9.6.12. Here's what I got: In [2]: import numpyx In [3]: numpyx.test() -=-=-=-=-=-=-=-=-=-= printing array info for ndarray at 0x0 print number of dimensions: 0 address of strides: 0x0 strides: memory dump: -1e-30 -=-=-=-=-=-=-=-=-=-= -=-=-=-=-=-=-=-=-=-= --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/nbecker/cython/<ipython console> in <module>() /home/nbecker/cython/numpyx.pyx in numpyx.test() 94 95 for arr in [arr1,arr2,arr3,arr4,arr5]: ---> 96 print arr 97 print_array_info(arr) 98 /home/nbecker/cython/numpyx.pyx in numpyx.print_array_info() 12 13 print '-='*10 ---> 14 print 'printing array info for ndarray at 0x%0lx'%(<c_python.Py_intptr_t>arr,) 15 print 'print number of dimensions:',arr.nd 16 print 'address of strides: 0x%0lx'%(<c_python.Py_intptr_t>arr.strides,) TypeError: only length-1 arrays can be converted to Python scalars

8 24

NumPy-Discussion February 2008