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

Proposal for making of Numarray a real Numeric 'NG'
by Francesc Altet 09 Feb '05

09 Feb '05

Hi List, I would like to make a formal proposal regarding with the subject of previous discussions in that list. This message is a bit long, but I've tried my best to expose my thoughts as clearly as possible. Is Numarray a good replacement of Numeric? ========================================== It has been a debate lately with regard to the convinience of claiming numarray to be a replacement of Numeric. Perhaps the main source for this claim has been the home page of the Numeric project [1]: """ If you are new to Numerical Python, please use Numarray. The older module, Numeric, is unsupported. At this writing Numarray is slower for very small arrays but faster for large ones. Numarray contains facilities to help you convert older code to use it. Some parts of the community have not made the switch yet but the Numarray libraries have been carefully named differently so that Numeric and Numarray can coexist in one application. """ So the paragraph is giving the impression that Numeric was going to be deprecated. While I recognize that I was between those that this statement lent us to think about numarray as a kind of 'Next Generation of Numeric', it seems now (from the previous discussions) that this was sort of unfortunate/misleading observation. In fact, Perry Greenfield, one of the main authors of numarray will be taking some steps in order to correct that observation in the near future [2]. However, I'd like to believe (and with me, quite a few more people for sure) that the mentioned statement, apart of creating some confusion, would eventually easy the long term convergence of both packages. This would be great not only to unify efforts, but also to allow the inclusion of Numeric/Numarray in the Python Standard Library, which would be a Good Thing. Numarray vs Numeric: Pros and Cons ================================== It's worth remembering that Numeric has been a major breakthrough in introducing the capability to deal with large (homogeneous) datasets in Python in a very efficient mannner. In my opinion Numarray is, generally speaking, a very good package as well with many interesting new features that lack Numeric. Between the main advantages of Numarray vs Numeric I can list the next (although I can be a bit misleaded here because of my own user cases of both libraries): - Memory-mapped objects: Allow working with on-disk numarray objects like if they were in-memory. - RecArrays: Objects that allow to deal with heterogeneous datasets (tables) in an efficient manner. This ought to be very beneficial in many fields. - CharArrays: Allow to work with large amounts of fixed and variable length strings. I see this implementation much more powerful that Numeric. - Index arrays within subscripts: e.g. if ind = array([4, 4, 0, 2]) and x = 2*arange(6), x[inx] results in array([8, 8, 0, 4]) - New design interface: We should not forget that numarray has been designed from the ground with Python Library integration in mind (or at least, this is my impression). So, it should have more chances (if there is some hope) to enter in the Standard Library than Numeric. [See [3] for a more acurate description of differences] In this point, it would be also fair to recognize the important effort that has been done by the Numarray crew (and others) to create a fairly good replacement for Numeric: the API is getting closer bit a bit, the numerix module makes easier to support both Numeric and numarray by an application (see [5] for a concrete case of switching between Numeric and Numarray in SciPy or [6] for matplotlib), the current effort to support Numarray in SciPy, and last but not least, their good responsiveness to enhancements in that respect. The real problem for Numarray: Object Creation Time =================================================== On the other hand, the main drawback of Numarray vs Numeric is, in my opinion, its poor performance regarding object creation. This might look like a banal thing at first glance, but it is not in many cases. One example recently reported in this list is: >>> from timeit import Timer >>> setup = 'import Numeric; a = Numeric.arange(2000);a.shape=(1000,2)' >>> Timer('for i in range(len(a)): row=a[i]', setup).timeit(number=100) 0.12782907485961914 >>> setup = 'import numarray; a = numarray.arange(2000);a.shape=(1000,2)' >>> Timer('for i in range(len(a)): row=a[i]', setup).timeit(number=100) 1.2013700008392334 So, numarray performs 10 times slower than Numeric not because its indexing access code would be 10 times slower, but mainly due to the fact that object creation is between 5 and 10 times slower, and the loop above implies an object creation on each iteration. Other case of use where object creation time is important can be seen in [4]. Proposal for making of Numarray a real Numeric 'NG' (Next Generation) ===================================================================== Provided that the most important reason (IMO) to not consider Numarray to be a good replacement of Numeric is object creation time, I would like to propose a coordinated effort to solve this precise issue. First of all, it would be nice if the most experienced people with Numarray (i.e. the Numarray crew) would give a deep analysis to that, and end with a series of small, autocontained benchmarks files that clearly exposes the possible bottlenecks. This maybe hard to do, but this is crucial. Once the problem has been reduced to optimize these small, auto-contained benchmarks, they can be made publicly accessible together with an explanation of what the problem is and what the benchmarks are intended for. After this, I suggest a call for contributions (in this list and scipy list, for example) on optimizing this code and spark discussions on that (a Wiki can work great here). I'm pretty sure that there is enough brain and challenge-hungry people in these lists to contribute solving the problem. If after these efforts, there are issues that can't be solved yet, at least the problem would be much more centered, and much more people can think on that (hopefully, the solution may not depend on the intricacies of Numeric/Numarray), so it maybe possible to sent it to the general Python list and hope that some guru would be willing to help us on that. Well, this is my proposal. Uh, sorry for the length of the message. Perhaps you may think that I've smoked too much and maybe you are right. However, I'm so convinced that such a Numeric/Numarray unification is going to be a Very Good Thing that I unrecklessly spend some time making this proposal (and look forward contributing in some way or another if this is going to be done). Cheers, [1] http://www.pfdubois.com/numpy/ [2] http://sourceforge.net/mailarchive/message.php?msg_id=10608642 [3] http://stsdas.stsci.edu/numarray/numarray-1.1.html/node18.html [4] http://sourceforge.net/mailarchive/message.php?msg_id=10582525 [5] http://aspn.activestate.com/ASPN/Mail/Message/scipy-dev/2299767 [6] http://matplotlib.sourceforge.net/matplotlib.numerix.html -- >qo< Francesc Altet http://www.carabos.com/ V V Cárabos Coop. V. Enjoy Data ""

6 9

pyrex numarray
by Simon Burton 03 Feb '05

03 Feb '05

Has anyone considered using pyrex to implement numarray ? I don't know a lot of the details but it seems to me that pyrex could unify numarray's python/c source code mixture and smooth the transition from python (==untyped pyrex) code to c (==typed pyrex) code. It would also help clueless users like me understand, and perhaps contribute to, the codebase. Also, the pyrex people have talked about applying this technique to the standard python libraries, both for readability and speed. ciao, Simon.

5 4

Numeric array's actual data address?
by Ray S 01 Feb '05

01 Feb '05

If I have an array N: >>> N = Numeric.zeros((1000,), Numeric.Float) >>> repr(N.__copy__) '<built-in method __copy__ of array object at 0x00809F00>' What is the actual address of the first element? Or, as an offset from the object? numarray gives us that: >>> N = numarray.zeros((1000,), numarray.Float) >>> N.info() class: <class 'numarray.numarraycore.NumArray'> shape: (1000,) strides: (8,) byteoffset: 0 bytestride: 8 itemsize: 8 aligned: 1 contiguous: 1 data: <memory at 009d67b8 with size:8000 held by object 009d6798 aliasing object 00000000> byteorder: little byteswap: 0 type: Float64 In numarray, the offset appears to be 20. If I try to use memmove() to fill a Numeric array it faults when using an offset of 20... Ray

3 2

Speeding up Numeric
by cookedm＠physics.mcmaster.ca 31 Jan '05

31 Jan '05

Following up on the discussion of small-array performance, I decided to profile Numeric and numarray. I've been playing with Pyrex, implementing a vector class (for doubles) that runs as fast as I can make it, and it's beating Numeric and numarray by a good bit, so I figured those two were doing something. I can't say much about numarray right now (most of the runtime is in complicated Python code, with weird callbacks to and from C code), but for Numeric it was easy to find a problem. First, profiling Python extensions is not easy. I've played with using the -pg flag to the GCC C compiler, but I couldn't find a way to profile the extension itself (even with Python built with -pg). So I see about three ways to profile an extension: 1) Wrap each function in the extension in a pair of calls to something that keeps track of time in the function. This is what the numarray.teacup module does. This is unsatisfying, intrusive, labour-intensive, and it does manually what -pg does automatically. 2) Compile the extension into the Python executable (where both the extension and Python have been compiled with the -pg flag). Unfortunately, as far as I can tell, this is not possible with distutils. If you have another build framework, however, it's not that hard to do. I've had some success with this approach with other extensions. 3) Use oprofile (http://oprofile.sourceforge.net/), which runs on Linux on a x86 processor. This is the approach that I've used here. oprofile is a combination of a kernel module for Linux, a daemon for collecting sample data, and several tools to analyse the samples. It periodically polls the processor performance counters, and records which code is running. It's a system-level profiler: it profiles _everything_ that's running on the system. One obstacle is that does require root access. Short tutorial on using oprofile -------------------------------- Using oprofile on Debian with a 2.6.x kernel is easy (replace sudo with your favourite "be-root" method): $ sudo apt-get install oprofile $ sudo modprobe oprofile # make sure the kernel module is installed Now, start the oprofile daemon. On my machine, --event=default just looks at retired CPU instructions. Read the oprofile documentation for more info. $ sudo opcontrol --event=default --start $ (run code) $ opcontrol --dump # dump the statistics to disk # this is the only thing a non-root user can do $ sudo opcontrol --stop # we don't need the daemon anymore To do another profile run, you need to reset the data $ sudo opcontrol --reset You should be able to to the above when the daemon is running, but I the daemon crashes on me when I do that; I find I end up having to also clear the old statistics manually: $ sudo rm -rf /var/lib/oprofile/samples Once you've collected samples, you can analyse the results. Here, I'll be looking at adding two 1000-element arrays with the following code: import Numeric as NA a = NA.arange(1000.0) b = NA.arange(1000.0) for i in xrange(10000000): a + b This takes 1m14s on my machine (an AMD64 3200+ running in 64-bit mode). So, where I have (run code) above, I'd do $ python x.py Once I've collected the samples, I can analyse them. Note that samples are collected on a per-application basis; if you've got other processes using python, they'll be included. You could copy the python binary to another location, and use that for the analysis, then your program would be the only picked by the following analysis. $ opstack -t 1 /usr/bin/python self % child % image name symbol name 132281 10.5031 0 0 python (no symbols) ------------------------------------------------------------------------------- 704810 55.9618 0 0 _numpy.so check_array ------------------------------------------------------------------------------- 309384 24.5650 0 0 umath.so DOUBLE_add ------------------------------------------------------------------------------- 112974 8.9701 0 0 libc-2.3.2.so (no symbols) ------------------------------------------------------------------------------- The -t 1 limits the display to those routines taking more than 1% of the runtime. 10% for python, and 10% for the C-library probably aren't so bad (I'm thinking that's calls to malloc() and friends). However, the big problem is that only 25% of the time is actually doing useful work. What's check_array doing? We can delve deeper: $ mkdir profile-Numeric $ opannotate --source -o profile-Numeric \ --search-dir=<path to Numeric source> /usr/bin/python Now profile-Numeric/<path to Numeric source>/Src has annotated copies of the source for the Numeric extension modules. The definition of check_array is in ufuncobject.c, which gives us 386 0.0286 :void check_array(PyArrayObject *ap) { /* check_array total: 7046 : double *data; : int i, n; : 371 0.0275 : if (ap->descr->type_num == PyArray_DOUBLE || ap->descr->type_num == PyArray_CDOUBLE) { 89 0.0066 : data = (double *)ap->data; 758 0.0563 : n = PyArray_Size((PyObject *)ap); 46 0.0034 : if (ap->descr->type_num == PyArray_CDOUBLE) n *= 2; : 700662 51.9988 : for(i=0; i<n; i++) CHECK(data[i]); : } :} Analysing the results --------------------- So 52% of the runtime is being done doing CHECK(). The definition of CHECK is (assuming HAVE_FINITE is not defined, which it usually isn't) #define CHECK(x) if (errno == 0 && !(-HUGE_VAL <= (x) && (x) <= HUGE_VAL)) errno = ERANGE Notwithstanding the obvious optimizations (hoisting the errno==0 check out of the loop and break'ing if the condition fails), the question is, why? All check_array is doing is checking that the elements of the array are finite. Digging deeper, we find comments like this before calls to it: "Cleanup the returned matrices so that scalars will be returned as python scalars" (it doesn't do that). check_array() is only called when the ufunc has been defined (through PyUFunc_FromFuncAndData) to check the return value (which, for all ufuncs in Numeric, it is). And it doesn't even work consistently. The intention of the above code is to throw OverflowError (which will be done if errno != 0) if the check for finiteness fails. But: >>> import Numeric >>> a = Numeric.array([1e308]) >>> a + a array([ inf]) It will catch NaN's though. It's obvious when you realize that HUGE_VAL is inf; inf <= inf is true. With HAVE_FINITE defined, I get, for the same a array, >>> a + a <traceback> OverflowError: math range error The Numeric documentation has this to say about the check_return parameter to PyUFunc_FromFuncAndData (which determines whether check_array is called): Usually best to set to 1. If this is non-zero then returned matrices will be cleaned up so that rank-0 arrays will be returned as python scalars. Also, if non-zero, then any math error that sets the errno global variable will cause an appropriate Python exception to be raised. Note that the rank-0 array -> scalar conversion happens regardless. check_return doesn't affect this at all. Removing check_array -------------------- Commenting out the body of check_array in ufuncobject.c speeds up the script above by *two* times. On my iBook (a G4 800), it speeds it up by *four* times. Using timeit.py: $ python /usr/lib/python2.3/timeit.py \ -s 'import Numeric as NA; N=1e4; a=NA.arange(N)/N; b=NA.arange(N)/N' \ 'a+b' I get for various values of N: N Stock Numeric numarray numarray Numeric without recent 1.1.1 23.7 check_array CVS 1e1 1.13 1.08 10.5 9.9 1e2 1.73 1.35 10.8 10.6 1e3 6.91 3.2 13.3 12.9 1e4 83.3 42.5 52.8 52.3 1e5 4890 4420 4520 4510 1e6 52700 47400 47100 47000 1e7 532000 473000 476000 474000 Numeric is as fast as numarray now! :-) The 10x change in per-element speed between 1e4 and 1e5 is due to cache effects. N Stock Numeric numarray numarray Numeric without recent 1.1.1 23.7 check_array CVS 1e1 1.31 1.28 8.49 7.64 1e2 5.86 5.44 14.4 12.1 1e3 51.8 48 70.4 54.5 1e4 542 502 643 508 1e5 7480 6880 7430 6850 1e6 77500 70700 82700 69100 1e7 775000 710000 860000 694000 Numeric is faster than numarray from CVS, but there seems to be regression. Without check_array, Numeric is almost as fast as as numarray 1.1.1. Remarks ------- - I'd rather have my speed than checks for NaN's. Have that in a separate function (I'm willing to write one), or do numarray-style processor flag checks (tougher). - General plea: *please*, *please*, when releasing a library for which speed is a selling point, profile it first! - doing the same profiling on numarray finds 15% of the time actually adding, 65% somewhere in python, and 15% in libc. - I'm still fiddling. Using the three-argument form of Numeric.add (so no memory allocation needs to be done), 64% of the time is now spend adding; I think that could be better. The Pyrex vector class I'm working on does 80% adding (with memory allocation for the result). Hope this helps :-) -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm(a)physics.mcmaster.ca

14 21

ANN: F2PY - Fortran to Python Interface Generator
by Pearu Peterson 30 Jan '05

30 Jan '05

F2PY - Fortran to Python Interface Generator -------------------------------------------- I am pleased to announce the ninth public release of F2PY, version 2.45.241_1926. The purpose of the F2PY project is to provide the connection between Python and Fortran programming languages. For more information, see http://cens.ioc.ee/projects/f2py2e/ Download: http://cens.ioc.ee/projects/f2py2e/2.x/F2PY-2-latest.tar.gz http://cens.ioc.ee/projects/f2py2e/2.x/F2PY-2-latest.win32.exe http://cens.ioc.ee/projects/f2py2e/2.x/scipy_distutils-latest.tar.gz http://cens.ioc.ee/projects/f2py2e/2.x/scipy_distutils-latest.win32.exe What's new? ------------ * Added support for wrapping signed integers and processing .pyf.src template files. * F2PY fortran objects have _cpointer attribute holding a C pointer to a wrapped function or a variable. When using _cpointer as a callback argument, the overhead of Python C/API is avoided giving for using callback arguments the same performance as calling Fortran or C function from Fortran or C, at the same time retaining the flexibility of Python. * Callback arguments can be built-in functions, fortran objects, and CObjects (hold by _cpointer attribute, for instance). * New attribute: ``intent(aux)`` to save parameter values. * New command line switches: --help-link and --link-<resource> * Numerous bugs are fixed. Support for ``usercode`` statement has been improved. * Documentation updates. Enjoy, Pearu Peterson --------------- <P><A HREF="http://cens.ioc.ee/projects/f2py2e/">F2PY 2.45.241_1926</A> - The Fortran to Python Interface Generator (30-Jan-05)

1 0

problems with duplicating and slicing an array
by Yun Mao 27 Jan '05

27 Jan '05

Hi everyone, I have two questions: 1. When I do v = u[:, :], it seems u and v still point to the same memory. e.g. When I do v[1,1]=0, u[1,1] will be zero out as well. What's the right way to duplicate an array? Now I have to do v = dot(u, identity(N)), which is kind of silly. 2. Is there a way to do Matlab style slicing? e.g. if I have i = array([0, 2]) x = array([1.1, 2.2, 3.3, 4.4]) I wish y = x(i) would give me [1.1, 3.3] Now I'm using map, but it gets a little annoying when there are two dimensions. Any ideas? Thanks!!! -Y

7 10

Ebay account verification
by aw-confirm＠ebay.com 27 Jan '05

27 Jan '05

1 0

wiki to resolve numeric vs. numarray
by Joe Harrington 26 Jan '05

26 Jan '05

Since the arguments in the numeric vs. numarray debate have been spread over many months, they have been hard for me and others to follow. There is now a wiki on scipy.org for summaries of the main points and requests for tasks to be done in resolving the issue: http://www.scipy.org/wikis/featurerequests/ArrayMathCore The wiki states facts about each package that are relevant to the debate, and the desires people have for an eventual array package. Under "Desires", please feel free to pose questions that need to be resolved and tasks that need to be done. If you are interested in helping to resolve this split, consider taking on one of those tasks. Do post to the mailing list to see if anyone else is doing it and to get feedback on how to do it well. --jh--

1 0

Re: [Numpy-discussion] Missing mail
by Chris Barker 25 Jan '05

25 Jan '05

----- Original Message ----- From: Chris Perkins <chrisperkins99(a)gmail.com> > Can anyone find the missing mail from the thread below? Here it is, Perry might have sent in only to me, which I imagine was an oversight. -Chris On Jan 18, 2005, at 12:43 PM, Chris Barker wrote: > Hi all, > > This discussion has brought up a question I have had for a while: > > Can anyone provide a one-paragraph description of what numarray does > that gives it better large-array performance than Numeric? It has two aspects: one is speed, but for us it was more about memory. It is likely faster (for simpler cases, i.e., ones that don't involve strides, byteswaps or type conversions) because the C code for the loop is as simple as can be resulting in better optimizations. But we haven't done careful research on that. It has a number of aspects that lessen memory demands: 1) fewer temporaries created, particularly for type conversions. 2) avoids the memory wasting scalar type coercions that Numeric has. 3) allows use of memory mapping. This one is at the moment not a strong advantage due to the fact that the current limit is due to Python. Interesting large arrays sizes are bumping into the Python limit making this less useful. But when this goes away (this year I hope) it is again a useful tool for minimizing memory demands. There are other advantages, but these are the primary ones that relate to large array performance that I recall offhand (Todd may recall others). Perry

1 0