Mailman 3 January 2007 - NumPy-Discussion

Should non ufunc numpy functions behave like ufunc regarding casting to output argument ?
by David Cournapeau Jan. 16, 2007

Jan. 16, 2007

Hi, I am trying to add support for out argument to one C function using numpy API (still the clip function). I was wondering about the expected behaviour when out does not have the "expected" type. For example, using again the clip function (but the question is not specific to this function) In [1]: import numpy In [2]: a = numpy.linspace(0, 10, 101) In [3]: b = numpy.zeros(a.shape, dtype = numpy.int32) In [4]: print a.dtype float64 In [5]: a.clip(0.1, 0.5, b) Should this be equivalent to b = a.clip(0.1, 0.5); b = b.astype(numpy.int32) (ie, the casting is done at the end, similar to an ufunc) ? cheers, David

5 7

Best way to compare two implementations of the same numpy function ?
by David Cournapeau Jan. 16, 2007

Jan. 16, 2007

Hi, I am finalizing my patch to SVN numpy to improve the clip function. I managed to modify the code so that the generated C sources are generated through the genapi system in numpy/core, etc.... and it compiles fine. The only problem I have left is "where" to expose the function to numpy ? I would like to have a new function to be able to test against the old version (I have a few tens of tests to compare results between current numpy clip and my function), but I am not sure how to do it so it can be commited as a patch to numpy SVN (if the patch is found acceptable by numpy developers, of cours), cheers, David

1 0

Can PyArray_CastTo be used to "subcast" an array of float to int
by David Cournapeau Jan. 16, 2007

Jan. 16, 2007

Hi, I would like to convert a double array to an *existing* array of int32 using the C numpy API. I tried using PyArray_CastTo, but this gives me an array full of zeros instead of converting the items to int. I am doing something like the following to convert a PyArrayObject in_o of double: ndescr = PyArray_DescrNewFromType(NPY_INT32); /* should be equivalent to astype */ ar = (PyArrayObject*)PyArray_FromArray((PyArrayObject*)in_o, ndescr, NPY_FORCECAST); if (ar == NULL) { goto clean_ndescr; } /* Using C Api for casting */ dims = PyArray_DIMS(in_o); Py_INCREF(ndescr); casted = (PyArrayObject*)PyArray_Zeros(((PyArrayObject*)in_o)->nd, dims, ndescr, 0); if (casted == NULL) { goto clean_ar; } st = PyArray_CastTo((PyArrayObject*)in_o, casted); if (st) { goto clean_casted; } When using this code in python with a = numpy.linspace(0, 10, 101) as input, I got array([0, 0, 1, 1, 2, 2, 3, 3, 4, 5] for ar, and array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]) for array casted. Can't PyArray_CastTo be used to 'unsafely' cast an array, or is there something I am missing ? cheers, David

1 1

Come to PyCon if you can
by Travis Oliphant Jan. 15, 2007

Jan. 15, 2007

Greetings, This is a little plug for coming to PyCon Feb 22-25. I am going to be there from Thursday night (Feb. 22) until Saturday evening (Feb. 25), and would like to meet with anyone interested to discuss some of the issues surrounding NumPy/SciPy and the PEPs we are working on for Python. I would stay for the sprints but I can't this year. If others are interested in using the sprints for SciPy/NumPy development, then we can still coordinate efforts. I will be presenting a talk on Understanding NumPy which is targeted to the beginner/intermediate audience on Saturday morning. There are other talks as well that have a distinct scientific computing flavor that look very interesting. Hopefully we can organize a BOF session as well during the conference (on Friday evening preferrably). Here are more details: http://us.pycon.org/TX2007/ Early bird registration ends Jan 15 (today!) Best regards, -Travis O.

2 1

rounding all element of an array with more than 1d
by Giorgio Luciano Jan. 15, 2007

Jan. 15, 2007

I've posted also to scipy group but with no luck. At first I thought it was an easy task but... I've to round the numbers inside of a 2d array but nothing seems to work. the only message I've found in the archive is an old function custom made for doing that.. wouldn't be useful to have it also included something that do this task in numpy ? here's and example of the error using round Traceback (most recent call last): File "<pyshell#25>", line 1, in <module> round(vp) TypeError: only length-1 arrays can be converted to Python scalars Any feedback appreciated

2 2

Fast clip for native types, 2d version
by David Cournapeau Jan. 15, 2007

Jan. 15, 2007

Hi, First, I wanted to thank everybody who helped me to clarify many points concerning memory layout of numpy arrays; I think I have a much clearer idea of the way numpy arrays behave at the C level. I've used all those informations to correct my initial implementation of clip to improve the clip function for common cases: it speeds up things only for native endianness, and scalar min and max (both contiguous and non contiguous cases). I've attached the new version code (only for one type to avoid too big emails; you have to dl the archive to actually compile the implementation); the whole package with tests + profiling script is there: http://www.ar.media.kyoto-u.ac.jp/members/david/archives/fastclip.tgz If this looks Ok, I will prepare a patch against current numpy, with the C sources being generated by numpy.distutils instead of the tool I am using now (autogen) Now, to improve other cases (mainly implementing an in-place clip function + non scalar min/max), there are some clarifications needed, mainly related to broadcast rules the current clip implementation which seems to break numpy conventions: 1: the old implementation returns an array which has the same endianness than the input array. This is a bit odd, because when the input is byte swapped, the returned array is still byte swapped, which seems to be against numpy convention. Here is some code which seem odd to me (code assumes little endian machine) a = numpy.random.randn(3, 2) b = a.astype(a.dtype.newbyteorder('>')) c = b.copy() assert a.dtype.isnative assert not b.dtype.isnative assert not c.dtype.isnative # Endianness behaviour of basic operation with numpy arrays print (a + b).dtype.isnative #one arg is non native -> returns native print (b + c).dtype.isnative # both args not native -> returns native # Now, what's happening endian-wise with clip: print numpy.clip(a, -0.5, 0.5).dtype.isnative # everything native -> returns native print numpy.clip(b, -0.5, 0.5).dtype.isnative # input array non native -> returns non native print numpy.clip(b, a, 0.5).dtype.isnative # input array non native, native array min -> returns native The fact that the output's endianness depends on min/max arguments being arrays or not does not seem really coherent ? 2: the old implementation does not upcast the input array. If the input is int32, and min/max are float32, the function fails; if input is float32, and min/max float64, the output is still float32. Again, this seems against the expected numpy behaviour ? 3: the old implementation supports clipping with complex arrays. I don't see any obvious meaningful implementation of clipping in those cases (using the module to compare them ?) If breaking those oddities is allowed, this would make the improvements much simpler to code, cheers, David

2 5

About extending NumPy with c
by asaf david Jan. 13, 2007

Jan. 13, 2007

Hello i want to write a c function that receive and process a numpy array. to do so, i'm following this tutorial http://numpy.scipy.org/numpydoc/numpy-13.html as a start, i copy pasted the "A Simple Example" from the tutorial and compiled the module as usual with no errors/warnings. however, when trying to run the function on a legal numpy array, i'm getting an error message: "python.exe has encountered a problem and needs to close, send, dont send etc' i've narrowed it down and it happens when running this line if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &array)) any ideas how to solve it? thanks here's my complete c file: #include <python.h> #include <numeric/arrayobject.h> static PyObject *trace(PyObject *self, PyObject *args) { PyArrayObject *array; double sum; int i, n; if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &array)) return NULL; if (array->nd != 2 || array->descr->type_num != PyArray_DOUBLE) { PyErr_SetString(PyExc_ValueError, "array must be two-dimensional and of type float"); 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]); return PyFloat_FromDouble(sum); } static PyMethodDef foo_methods[] = { { "trace", trace, METH_VARARGS, NULL}, { NULL, NULL, 0, NULL } }; PyMODINIT_FUNC initfoo() { Py_InitModule3("foo", foo_methods, "bla bla"); }

2 1

Need help for implementing a fast clip in numpy (was slow clip)
by David Cournapeau Jan. 13, 2007

Jan. 13, 2007

Hi, I am finally implementing a C function to replace the current slow implementation of clip in python as promised a few weeks ago. The idea is to implement it as the following: def clip(input, min, max): a = input.copy() putmask(a, a <= min, min) putmask(a, a >= max, max) return a I don't have that much experience in writing general numpy functions, so I was wondering of other people could advise me on the following points. - should I care only about numeric types, eg something where PyArray_ISNUMBER returns true ? - If input is an array of type T1, min a scalar of type T2 and max a scalar of type T3, which type should be the reference, assuming they are all basic numeric types (float, int, etc...). Eg, if T1 is integer, T2 and T3 are float, what am I supposed to do ? Should I use the basic rule as told in numpy doc, and in that case, and do something like this to determine common datatype (pseudo code): int typenum; for each arg in arglist: typenum = PyArray_ObjectType(arg, typenum) - PyArray_PutMask does not assume anything on the datatype of the input array nor this alignement, and treat its binary content as a buffer to char*. Do I need to do the same, or can I just use "normal" pointer arithmetic on arrays which elements can be casted to basic C types ? cheers, David

9 29

Can numpy compute ln(640320**3 + 744)/163**.5 to 30 places?
by Dick Moores Jan. 13, 2007

Jan. 13, 2007

Can python compute ln(640320**3 + 744)/163**.5 to 30 places? I ran across this the other day. It's value is said to be the same as pi to about 30 places. I'd like to know if it is possible to do this in python with numpy. A friend pointed out that the Windows Calculator gives about 30 significant digits. I tried it and got 3.1415926535897932384626433832797. Pi to 32 places is 3.14159265358979323846264338327950, so the calculator is accurate to 31 significant digits (including the 3.). So, how about numpy? Thanks, Dick Moores Win XP Pro Python 2.5 yet to install numpy

1 0

Overload binary operators with custom array type
by Adam Jenkins Jan. 12, 2007

Jan. 12, 2007

I have written my own array class, let's call it "myarray", which behaves like a numpy.ndarray. I would like to support binary operators with mixed operand types, i.e.: lhs + rhs where one of the operands is a myarray object, and the other is a numpy.ndarray object. The problem is that I want the operation in this case to be handled by my add method, not by numpy's add method. If lhs is a myarray object, then this is no problem, since myarray.__add__(lhs, rhs) will be called. However if lhs is a numpy array and rhs is a myarray, then numpy.ndarray.__add__(lhs, rhs) will be called, whereas what I would like to happen would be for myarray.__radd__(rhs, lhs) to be called. Is there a way I can get numpy to do the latter? Python will only call the rhs.__radd__ method if lhs.__add__ returns NotImplemented, but as far as I can tell, the numpy.ndarray.__add__ method never returns NotImplemented. It seems that for "lhs + rhs", where lhs is a numpy array, numpy does one of two things: 1. if rhs has an __array__ method, than it calls rsh.__array__(), which is expected to return a numpy.ndarray object. Then it performs the binary operation with the result of __array__ for the rhs. 2. If rhs doesn't have an __array__ method, then it treats rhs as some kind of scalar, and tries to add each element of lhs to rhs. Under no circumstances does numpy just assume that it doesn't know what to do with rhs and return NotImplemented. My question is, is there a way I can get the effect of having myarray.__radd__(rhs, lhs) in this case, or in some way get numpy to invoke my code to handle this operation rather than just assuming rhs is a scalar? As for why I don't just implement myarray.__array__(), the reason is because it would be expensive. A myarray object doesn't actually contain the array data; rather it contains a handle referring to an array which resides in another process. The server process which actually contains the arrays can even be on another machine. So converting a myarray object to a numpy.ndarray involves transferring all the array data to the python process. For various reasons, for my application it's preferable that I handle binary operators with mixed numpy/myarray operand types myself, rather than converting the myarray object to a numpy.ndarray. Thanks in advance for any suggestions you can give. Adam

3 4