Mailman 3 February 2007 - NumPy-Discussion

fill_value in masked arrays
by Keith Goodman Feb. 11, 2007

Feb. 11, 2007

The doc strings for MaskedArray and ma.argsort do not seem consistent. The doc string for MaskedArray says "The fill_value is not used for computation within this module." But the doc string for ma.argsort says "Treating masked values as if they have the value fill_value, return sort indices for sorting along given axis." I'm trying to argsort matrices that have missing values.

2 1

Numpy x Matlab: some synthetic benchmarks
by Paulo Jose da Silva e Silva Feb. 11, 2007

Feb. 11, 2007

Hello, Travis asked me to benchmark numpy versus matlab in some basic linear algebra operations. Here are the resuts for matrices/vectors of dimensions 5, 50 and 500: Operation x'*y x*y' A*x A*B A'*x Half 2in2 Dimension 5 Array 0.94 0.7 0.22 0.28 1.12 0.98 1.1 Matrix 7.06 1.57 0.66 0.79 1.6 3.11 4.56 Matlab 1.88 0.44 0.41 0.35 0.37 1.2 0.98 Dimension 50 Array 9.74 3.09 0.56 18.12 13.93 4.2 4.33 Matrix 81.99 3.81 1.04 19.13 14.58 6.3 7.88 Matlab 16.98 1.94 1.07 17.86 0.73 1.57 1.77 Dimension 500 Array 1.2 8.97 2.03 166.59 20.34 3.99 4.31 Matrix 17.95 9.09 2.07 166.62 20.67 4.11 4.45 Matlab 2.09 6.07 2.17 169.45 2.1 2.56 3.06 Obs: The operation Half is actually A*x using only the lower half of the matrix and vector. The operation 2in2 is A*x using only the even indexes. Of course there are many repetitions of the same operation: 100000 for dim 5 and 50 and 1000 for dim 500. The inner product is number of repetitions is multiplied by dimension (it is very fast). The software is numpy svn version 1926 Matlab 6.5.0.180913a Release 13 (Jun 18 2002) Both softwares are using the *same* BLAS and LAPACK (ATLAS for sse). As you can see, numpy array looks very competitive. The matrix class in numpy has too much overhead for small dimension though. This overhead is very small for medium size arrays. Looking at the results above (specially the small dimensions ones, for higher dimensions the main computations are being performed by the same BLAS) I believe we can say: 1) Numpy array is faster on usual operations but outerproduct (I believe the reason is that the dot function uses the regular matrix multiplication to compute outer-products, instead of using a special function. This can "easily" changes). In particular numpy was faster in matrix times vector operations, which is the most usual in numerical linear algebra. 2) Any operation that involves transpose suffers a very big penalty in numpy. Compare A'*x and A*x, it is 10 times slower. In contrast Matlab deals with transpose quite well. Travis is already aware of this and it can be probably solved. 3) When using subarrays, numpy is a slower. The difference seems acceptable. Travis, can this be improved? Best, Paulo Obs: Latter on (in a couple of days) I may present less synthetic benchmarks (a QR factorization and a Modified Cholesky).

13 27

Function numpy.get_include() has a side effet
by Nicolas Grilly Feb. 10, 2007

Feb. 10, 2007

Hello, I have a package with a distutils script. This script calls numpy.get_include(), in order to get numpy's header files location. But the function seems to have a side effect by silently patching distutils. Some classes and functions of standard library distutils are replaced by numpy's ones. For example, before calling numpy.get_include(), the function new_compiler comes from standard library: >>> from distutils.ccompiler import new_compiler >>> new_compiler.__module__ 'distutils.ccompiler' But after calling numpy.get_include(), the function comes from numpy: >>> import numpy >>> numpy.get_include() >>> from distutils.ccompiler import new_compiler >>> new_compiler.__module__ 'numpy.distutils.ccompiler' This is because the code of numpy.get_include() imports numpy.distutils, and this triggers the patching mechanism. This is a problem because the behavior of distutils is different before and after the call to numpy.get_include(), and it breaks my distutils script. Is it possible to remove this side effect? Thanks, Nicolas Grilly

1 0

averaging array, with fill values
by Evan Mason Feb. 9, 2007

Feb. 9, 2007

Hi, I want to get the mean over an array with 3 dimensions. The array has, say, dimensions 2x3x4 and I want the result to be 3x4; so the averaging is to be done over the 1st dimension. Also, there are some fill values, and I want these to be excluded from the calculation. I know I can use the ma module to locate the fill values, but after that I am unsure how to exclude them, and how to do the averaging. In [150]: x = rand(2,3,4) In [151]: x[0,2,1]=999 # fill value In [152]: x[1,2,3]=999 In [153]: x[0,1,3]=999 In [154]: x Out[154]: array([[[ 3.01880915e-02, 5.77085271e-02, 7.59176038e-01, 4.15271486e-01], [ 5.48643693e-01, 3.84995126e-01 , 5.01683678e-01, 9.99000000e+02], [ 5.72779004e-01, 9.99000000e+02, 4.18143934e-01, 5.84781674e-01]], [[ 8.90443158e-01, 3.76986788e-01, 8.81270409e-01, 5.19094405e-01], [ 8.12944573e-01, 3.89858156e-01, 5.99219891e-01, 9.99000000e+02], [ 2.31215256e-01, 5.93222883e-01, 2.45004093e-01, 9.18647954e-01]]]) In [155]: Thanks in advance for any help. -Evan

1 0

Re: [Numpy-discussion] Request for porting pycdf to NumPy
by Daran Rife Feb. 9, 2007

Feb. 9, 2007

Yep, I was referring to: http://pysclint.sourceforge.net/pycdf/ Regarding the issue of deciding which netCDF interface to adopt as the "standard", although it is unlikely we'll ever get consensus on this, I have tried several of the netCDF interfaces, including the one in Scientific, and have found that pycdf is the "cleanest"; it is very easy to use and understand. Daran -- Travis Oliphant wrote: > Is pycdf the same as pynetcdf? Where is it located? > I presume it's this one: http://pysclint.sourceforge.net/pycdf/

5 4

numarray argmax problem
by Otto Tronarp Feb. 9, 2007

Feb. 9, 2007

Hi, I have a problem with numarray.argmax, the following code import numarray as N import sys c = N.zeros((10, ), N.Float64) while 1: print 'B: ', sys.getrefcount(None) l = N.argmax(c) print 'A: ', sys.getrefcount(None) print Dies with: Fatal Python error: deallocating None Aborted I'm using numarray 1.5.2 Any chance that there is an easy fix for this? I know I should consider switching to numpy, but I don't have the time right now to do that. Otto

2 1

bug or feature?
by humufr＠yahoo.fr Feb. 9, 2007

Feb. 9, 2007

I have a big problem with numpy, numarray and Numeric (all version) If I'm using the script at the bottom, I obtain these results: var1 before function is [3 4 5] var2 before function is 1 var1 after function must be [3 4 5] is [ 9 12 15] <------ problem var2 after function must be 1 is 1 var3 must be the [9 12 15] is [ 9 12 15] var4 must be the 'toto' is toto I'm very surprised by the line noted. I always thinking that the input variable didn't change the variable itself outside the function. It's the comportement for var2 but var1 is changed and it's a big problem (at least for me). The only object in python with this behavior are the numeric object (Numeric, numarray or numpy), with list or other kind of object I have the expected result (the var1 before to go inside the function) I can't keep the input variable if I'm not doing a copy before to call a function... Is it normal and so do I have to do a copy of the input data each time I'm calling a function? Thanks Nicolas #!/usr/bin/env python import numpy print "numpy version ", numpy.__version__ def test(var1,var2): #print "var1 input function is",var1 #print "var2 input function is",var2 var1 *=3 var2 = 'toto' return var1,var2 var1=numpy.array([3,4,5]) var2=1 print "var1 before function is ",var1 print "var2 before function is ",var2 var3,var4=test(var1,var2) print "var1 after function must be [3 4 5] is ",var1 print "var2 after function must be 1 is ", var2 print "var3 must be the [9 12 15] is ", var3 print "var4 must be the 'toto' is ", var4

5 5

Comparing x and x.view
by Pierre GM Feb. 8, 2007

Feb. 8, 2007

All, I want to compare whether two arrays point to the same data. I've been using 'is' so far, but I'm wondering whether it's the right approach. If x is a plain ndarray, `x is x`, and `x is not x.view()`. I understand the second one (I think so...), `x` and `x.view `are two different Python objects. However, `x.__array_interface__ == x.view().__array_interface__`, which means that the underlying memory hasn't been modified at all, right ? In other terms, the data hasn't been copied, it's just being accessed slightly differently. So, when I'm using `x is y` to test whether some data has been copied, I should in fact compare the __array_interface__s, shouldn't I ? Sorry for the poor phrasing, and thanks a lot for your forthcoming inputs. P.

4 6

Latest Array-Interface PEP
by Travis Oliphant Feb. 8, 2007

Feb. 8, 2007

I'm attaching my latest extended buffer-protocol PEP that is trying to get the array interface into Python. Basically, it is a translation of the numpy header files into something as simple as possible that can still be used to describe a complicated block of memory to another user. My purpose is to get feedback and criticisms from this community before display before the larger Python community. -Travis PEP: <unassigned> Title: Extending the buffer protocol to include the array interface Version: $Revision: $ Last-Modified: $Date: $ Author: Travis Oliphant <oliphant(a)ee.byu.edu> Status: Draft Type: Standards Track Created: 28-Aug-2006 Python-Version: 2.6 Abstract This PEP proposes extending the tp_as_buffer structure to include function pointers that incorporate information about the intended shape and data-format of the provided buffer. In essence this will place an array interface directly into Python. Rationale Several extensions to Python utilize the buffer protocol to share the location of a data-buffer that is really an N-dimensional array. However, there is no standard way to exchange the additional N-dimensional array information so that the data-buffer is interpreted correctly. The NumPy project introduced an array interface (http://numpy.scipy.org/array_interface.shtml) through a set of attributes on the object itself. While this approach works, it requires attribute lookups which can be expensive when sharing many small arrays. One of the key reasons that users often request to place something like NumPy into the standard library is so that it can be used as standard for other packages that deal with arrays. This PEP provides a mechanism for extending the buffer protocol (which already allows data sharing) to add the additional information needed to understand the data. This should be of benefit to all third-party modules that want to share memory through the buffer protocol such as GUI toolkits, PIL, PyGame, CVXOPT, PyVoxel, PyMedia, audio libraries, video libraries etc. Proposal Add bf_getarrview and bf_relarrview function pointers to the buffer protocol to allow objects to share a view on a memory pointer including information about accessing it as an N-dimensional array. Add the TP_HAS_ARRAY_BUFFER flag to types that define this extended buffer protocol. Also a few additionsl C-API calls should perhaps be added to Python to facilitate creating new PyArrViewObjects. Specification: static PyObject* bf_getarrayview (PyObject *obj) This function must return a new reference to a PyArrViewObject which contains the details of the array information exposed by the object. If failure occurs, then NULL is returned and an exception set. static int bf_relarrayview(PyObject *obj) If not NULL then this will be called when the object returned by bf_getarrview is destroyed so that the underlying object can be aware when acquired "views" are released. The object that defines bf_getarrview should not re-allocate memory (re-size itself) while views are extant. A 0 is returned on success and a -1 and an error condition set on failure. The PyArrayViewObject has the structure typedef struct { PyObject_HEAD void *data; /* pointer to the beginning of data */ int nd; /* the number of dimensions */ Py_ssize_t *shape; /* c-array of size nd giving shape */ Py_ssize_t *strides; /* SEE BELOW */ PyObject *base; /* the object this is a "view" of */ PyObject *format; /* SEE BELOW */ int flags; /* SEE BELOW */ } PyArrayViewObject; strides -- a c-array of size nd providing the striding information which is the number of bytes to skip to get to the next element in that dimension. format -- a Python data-format object (PyDataFormatObject) which contains information about how each item in the array should be interpreted. flags -- an integer of flags. PYARR_WRITEABLE is the only flag that must be set appropriately by types. Other flags: PYARR_ALIGNED, PYARR_C_CONTIGUOUS, PYARR_F_CONTIGUOUS, and PYARR_NOTSWAPPED can all be determined from the rest of the PyArrayViewObject using the UpdateFlags C-API. The PyDataFormatObject has the structure typedef struct { PyObject_HEAD PySimpleformat primitive; /* basic primitive type */ int flags; /* byte-order, isaligned */ int itemsize; /* SEE BELOW */ int alignment; /* SEE BELOW */ PyObject *extended; /* SEE BELOW */ } PyDataFormatObject; enum Pysimpleformat {PY_BIT='1', PY_BOOL='?', PY_BYTE='b', PY_SHORT='h', PY_INT='i', PY_LONG='l', PY_LONGLONG='q', PY_UBYTE='B', PY_USHORT='H', PY_UINT='I', PY_ULONG='L', PY_ULONGLONG='Q', PY_FLOAT='f', PY_DOUBLE='d', PY_LONGDOUBLE='g', PY_CFLOAT='F', PY_CDOUBLE='D', PY_CLONGDOUBLE='G', PY_OBJECT='O', PY_CHAR='c', PY_UCS2='u', PY_UCS4='w', PY_FUNCPTR='X', PY_VOIDPTR='V'}; Each of these simple formats has a special character code which can be used to identify this primitive in a nested python list. flags -- flags for the data-format object. Specified masks are PY_NATIVEORDER PY_BIGENDIAN PY_LITTLEENDIAN PY_IGNORE itemsize -- the total size represented by this data-format in bytes unless the primitive is PY_BIT in which case it is the size in bits. For data-formats that are simple 1-d arrays of the underlying primitive, this total size can represent more than one primitive (with extended still NULL). alignment -- For the primitive types this is offsetof(struct {char c; type v;},v) extended -- NULL if this is a primitive data-type or no additional information is available. If primitive is PY_FUNCPTR, then this can be a tuple with >=1 element: (args, {dim0, dim1, dim2, ...}). args -- A list (of at least length 2) of data-format objects specifying the input argument formats with the last argument specifying the output argument data-format (use None for void inputs and/or outputs). For other primitives, this can be a tuple with >=2 elements: (names, fields, {dim0, dim1, dim2, ...}) Use None for both names and fields if they should be ignored. names -- An ordered list of string or unicode objects giving the names of the fields for a structure data-format. fields -- a Python dictionary with ordered-keys given by the list in names. Each entry in the dictionary is a 3-tuple containing (data-format-object, offset, meta-information) where meta-information is Py_None if there is no meta-information. Offset is given in bytes from the start of the record or in bits if PY_BIT is the primitive. Any additional entries in the extended tuple (dim0, dim1, etc.) are interpreted as integers which specify that this data-format is an array of the given shape of the fundamental data-format specified by the remainder of the DataFormat Object. The dimensions are specified so that the last-index is always assumed to vary the fastest (C-order). The constructor of a PyArrViewObject allocates the memory for shape and strides and the destructor frees that memory. The constructor of a PyDataFormatObject allocates the objects it needs for fields, names, and shape. C-API void PyArrayView_UpdateFlags(PyObject *view, int flags) /* update the flags on the array view object provided */ PyDataFormatObject *Py_NewSimpleFormat(Pysimpleformat primitive) /* return a new primitive data-format object */ PyDataFormatObject *Py_DataFormatFromCType(PyObject *ctype) /* return a new data-format object from a ctype */ int Py_GetPrimitiveSize(Pysimpleformat primitive) /* return the size (in bytes) of the provided primitive */ PyDataFormatObject *Py_AlignDataFormat(PyObject *format) /* take a data-format object and construct an aligned data-format object where all fields are aligned on appropriate boundaries for the compiler */ Discussion The information provided in the array view object is patterned after the way a multi-dimensional array is defined in NumPy -- including the data-format object which allows a variety of descriptions of memory depending on the need. Reference Implementation Supplied when the PEP is accepted. Copyright This document is placed in the public domain.

15 50

getting indices for array positions
by Christian Meesters Feb. 8, 2007

Feb. 8, 2007

Hi This questions might seem stupid, but I didn't get a clever solution myself, or found one in the archives, the cookbook, etc. . If I overlooked something, please give a pointer. Well, if I have an 1D array like [ 0. , 0.1, 0.2, 0.3, 0.4, 0.5] ,a scalar like 0.122 and want to retrieve the index postion of the closest value of the scalar in the array: Is there any fast method to get this? Right now I've implemented the following method: def _get_value_index(value, a): mindiff = 1e20 index = 0 for intensity, temp_index in zip(a, xrange(a.shape[0])): diff = abs(intensity - value) #closer to given value? if diff <= mindiff: mindiff = diff index = temp_index return index It works, but is akward and takes too much time (I've no benchmark), if the array is long and the method called often within a different function. But it should help clarify the problem. TIA Christian

4 6