Mailman 3 March 2006 - 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

Numpy x Matlab: some synthetic benchmarks
by Paulo Jose da Silva e Silva 11 Feb '07

11 Feb '07

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

recarray field names
by Erin Sheldon 08 Jul '06

08 Jul '06

Hi all- I'm new to numpy and I'm slowly weening myself off of IDL. So far the experience has been very positive. I use recarrays a lot. I often read these recarrays from fits files using pyfits (modified to work with numpy). I find myself doing the following more than I would like: if 'TAG' in rec._coldefs.names: .... It seems messy to be accessing this "hidden" attribute in this way. Is there a plan to add methods to more transparently do such things? e.g. def fieldnames(): return _coldefs.names def field_exists(fieldname): return fieldname.upper() in _coldefs.names def field_index(fieldname): if field_exists(fieldname): return _coldefs.names.index(fieldname.upper()) else: return -1 # or None maybe Thanks, Erin

5 13

numexpr thoughts
by Tim Hochberg 04 May '06

04 May '06

1. David already mentioned opcodes to copy from an integer array into a float register. Another idea would be to have an opcode to copy from a strided array into a register. This would save a copy for noncontiguous matrices and opens up the door to doing broadcasting. I think we'd have to tweak the main interpreter loop a bit, but I think it's doable. Depending how crunched we feel for opcode space, the casting array and this array could overlap (OP_CAST_AND_COPY for instance). 2. Something I noticed while writing tests is that the rearangement of operations for 'a/const' to '(1/const)*a' changes the results slightly (try const=3). I don't think I care -- I just thought I'd point it out. 3. It may simplify things to use copy_c to eliminate a bunch of the extra bytecodes for operating on functions of more than one argument. I need to check the timings on this -- I don't know how much of a speed hit using copy_c will cause. 4. At some point I plan to add complex operations. My original thought was to just duplicate the current, floating-point interpreter for complex numbers. However, if we have opcodes for casting, it might instead make sense to have one interpreter that works with both. This would be more efficient for mixed expressions since we wouldnt' have to case everything complex right away. It looks like there's enough opcode space, although I worry a bit whether making that main switch loop huge is going to slow things down. Probably not, but it's something to keep an eye on. 5. Same thoughts as 4., but for integer operations. Potentially the interpreter could work with longs, doubles and cdoubles, downcasting as appropriate on the way out. It's looking like numexpr could come pretty close to evaluating most numpy expressions if we put enough work into it. None of the required changes look like they should slow down the original, simple, fast case. However we should keep an eye on that as we go along. It's all very cool; I'd like to thank David for providing such a fun toy. -tim

4 12

Dumping record arrays
by Pierre GM 22 Apr '06

22 Apr '06

Folks, I'd like to dump/pickle some record arrays. The pickling works, the unpickling raises a ValueError (on my version of numpy 0.9.6). (cf below). Is this already corrected in the svn version ? Thx ########################################################################### # x1 = array([21,32,14]) x2 = array(['my','first','name']) x3 = array([3.1, 4.5, 6.2]) r = rec.fromarrays([x1,x2,x3], names='id, word, number') r.dump('dumper') rb=load('dumper') --------------------------------------------------------------------------- exceptions.ValueError Traceback (most recent call last) /home/backtopop/Work/workspace-python/pyflows/src/<ipython console> /usr/lib64/python2.4/site-packages/numpy/core/numeric.py in load(file) 331 if isinstance(file, type("")): 332 file = _file(file,"rb") --> 333 return _cload(file) 334 335 # These are all essentially abbreviations /usr/lib64/python2.4/site-packages/numpy/core/_internal.py in _reconstruct(subtype, shape, dtype) 251 252 def _reconstruct(subtype, shape, dtype): --> 253 return ndarray.__new__(subtype, shape, dtype) 254 255 ValueError: ('data-type with unspecified variable length', <function _reconstruct at 0x2aaaafcf1578>, (<class 'numpy.core.records.recarray'>, (0,), 'V'))

3 2

Multiple inheritance from ndarray
by Robert Lupton 18 Apr '06

18 Apr '06

I have a swig extension that defines a class that inherits from both a personal C-coded image struct (actImage), and also from Numeric's UserArray. This works very nicely, but I thought that it was about time to upgrade to numpy. The code looks like: from UserArray import * class Image(UserArray, actImage): def __init__(self, *args): actImage.__init__(self, *args) UserArray.__init__(self, self.getArray(), 'd', copy=False, savespace=False) I can't figure out how to convert this to use ndarray, as ndarray doesn't seem to have an __init__ method, merely a __new__. So what's the approved numpy way to handle multiple inheritance? I've a nasty idea that this is a python question that I should know the answer to, but I'm afraid that I don't... R

4 8

ndarray.fill and ma.array.filled
by Sasha 11 Apr '06

11 Apr '06

In an ideal world, any function that accepts ndarray would accept ma.array and vice versa. Moreover, if the ma.array has no masked elements and the same data as ndarray, the result should be the same. Obviously current implementation falls short of this goal, but there is one feature that seems to make this goal unachievable. This feature is the "filled" method of ma.array. Pydoc for this method reports the following: | filled(self, fill_value=None) | A numeric array with masked values filled. If fill_value is None, | use self.fill_value(). | | If mask is nomask, copy data only if not contiguous. | Result is always a contiguous, numeric array. | # Is contiguous really necessary now? That is not the best possible description ("filled" is "filled"), but the essence is that the result of a.filled(value) is a contiguous ndarray obtained from the masked array by copying non-masked elements and using value for masked values. I would like to propose to add a "filled" method to ndarray. I see several possibilities and would like to hear your opinion: 1. Make filled simply return self. 2. Make filled return a contiguous copy. 3. Make filled replace nans with the fill_value if array is of floating point type. Unfortunately, adding "filled" will result is a rather confusing situation where "fill" and "filled" both exist and have very different meanings. I would like to note that "fill" is a somewhat odd ndarray method. AFAICT, it is the only non-special method that mutates the array. It appears to be just a performance trick: the same result can be achived with "a[...] = ".

13 68

New patch for MA
by Pierre GM 06 Apr '06

06 Apr '06

Folks, You can find a new patch for MA on the wiki http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/ma-200603… along with a test suite. The 'clip' method should now work with array arguments. Were also added cumsum, cumprod, std, var and squeeze. I'll deal with flags, setflags, setfield, dump and others when I'll have a better idea of how it works -- which probably won't happen anytime soon, as I don't really have time to dig in the code for these functions. AAMOF, I'm more interested in checking/patching some other aspects of numpy for MA (eg, mlab...) Once again, please send me your comments and suggestions. Thx for everything P.

3 4

Re: [Numpy-discussion] first impressions with numpy
by Tim Hochberg 04 Apr '06

04 Apr '06

Sebastian Haase wrote: > Thanks Tim, > that's OK - I got the idea... > BTW, is there a (policy) reason that you sent the first email just to > me and not the mailing list !? No. Just clumsy fingers. Probably the same reason the functions got all garbled! > > I would really be more interested in comments to my first point ;-) > I think it's important that numpy will not be to cryptic and only for > "hackers", but nice to look at ... (hope you get what I mean ;-) Well, I think it's probably a good idea and it sounds like Travis like the idea " for some of the builtin types". I suspect that's code for "not types for which it doesn't make sense, like recarrays". > I'm have developed an "image analysis algorithm development" platform > (based on wxPython + PyShell) that more and more people in our lab are > using (instead of Matlab !) and I changed the default sys.displayhook > to print str(...) instead of repr(...) > mainly to get .3 instead of .2999999999998 > but seing int32 instead of <i4 would also be much "nicer". I agree that str should display something nicer. Repr should probably stay the same though. > Thanks for your great work .. Oh, I'm not doing much in the way of great work. I'm mostly just causing Travis headaches. -tim > Sebastian > > > Tim Hochberg wrote: > >> Tim Hochberg wrote: >> >>> Sebastian Haase wrote: >>> >>>> Hi, >>>> I'm a long time user of numarray. Now I downloaded numpy for the >>>> first time - and am quite excited to maybe soon being able to use >>>> things like weave ! >>>> Thanks for all the good work ! >>>> >>> [SNIP] >>> >>>> >>>> >>>> 2) >>>> This is probably more numarray related: >>>> Why does numarray.asarray( numpy.array([1]) ) return a numpy array, >>>> not a numarray ??? This is even true for numarray.array( >>>> numpy.array([1]) ) !! >>>> >>> I expect that numarray is grabbing the output of __array__() and >>> using that for asarray. What's happening in the second case is hard >>> to know, but I bet it's some snafu with assuming that anything that >>> isn't a numarray, but implements __array__ must be returning a new >>> numarray. That's just a guess. FWIW, I'm using the following two >>> functions to go back and forth between numarray and numpy. So far >>> they seem to work and they're faster (or actually work in the case >>> of asnumarray). You may want to rename them to suit your personal >>> preferences. If you use them, please let me know if you find any >>> problems. >>> >>> Regards, >>> >> I see this got thoroughly mangled somehow. If you want a copy, let me >> know and I'll just send you it as a file. >> >> -tim >> > > >

9 14

numpy error handling
by Tim Hochberg 02 Apr '06

02 Apr '06

I've just been looking at how numpy handles changing the behaviour that is triggered when there are numeric error conditions (overflow, underflow, etc.). If I understand it correctly, and that's a big if, I don't think I like it nearly as much as the what numarray has in place. It appears that numpy uses the two functions, seterr and geterr, to set and query the error handling. These set/read a secret variable stored in the local scope. I assume that the various ufuncs then examine that value to determine how to handle errors. The secret variable approach is a little clunky, but that's not what concerns me. What concerns me is that this approach is *only* useful for built in numpy functions and falls down if we call any user defined functions. Suppose we want to be warned on underflow. Setting this is as simple as: def func(*args): numpy.seterr(under='warn') # do stuff with args return result Since seterr is local to the function, we don't have to reset the error handling at the end, which is convenient. And, this works fine if all we are doing is calling numpy functions and methods. However, if we are calling a function of our own devising we're out of luck since the called function will not inherit the error settings that we have set. Thus we have no way to influence the error settings of functions downstream from us. Compare this to the numarray approach, where you write code like the following: def func(*args) numarray.Error.pushMode(underflow='warn') try: # do stuff with args return result finally: numarray.Error.popMode() At first glance this looks a lot more clunky. However, it's also a lot more useful. In this case, pushMode pushes a new error mode value onto a global[1] stack of mode values. Thus, any functions that func calls, will inherit it's error settings, unless they explicitly override them. As for the clunkiness factor, that's partly just the choice of camelCase versus smashwords, but also the try-finally. The second, however, will be addressed in Python 2.5 with the 'with' statement. At that point, we could have a numpy equivalent to the above numarray code that looked something like: def func(*args) with numpy.errors(underflow='warn'): # do stuff with args return result I believe we should rework the error handling to follow the numarray model before things get set in stone at 1.0. The 'with' statement support can be worked in a bit later if need be, although I suspect that will be easy. I also would prefer more verbose keys ala numarray (underflow, overflow, dicidebyzero and invalid) than those currently used by numpy (under, over, divide and invalid). And (will he never stop) I like numarrays defaults better here too: overflow='warn', underflow='ignore', dividebyzero='warn', invalid='warn'. Currently, numpy defaults to ignore for all cases. These last points are relatively minor though. Regards, -tim [1] I believe it's actually thread local, and if it's not it should be, but that detail is not important for our discussion here.

4 7