Hi Chuck, Good questions. Responses inline below... Jason On Thu, Jun 10, 2010 at 8:26 AM, Charles R Harris <charlesr.harris@gmail.com

wrote:

On Wed, Jun 9, 2010 at 5:27 PM, Jason McCampbell < jmccampbell@enthought.com> wrote:

...

Hi everyone,

This is a follow-up to Travis's message on the re-factoring project from May 25th and the subsequent discussion. For background, I am a developer at Enthought working on the NumPy re-factoring project with Travis and Scott. The immediate goal from our perspective is to re-factor the core of NumPy into two architectural layers: a true core that is CPython-independent and an interface layer that binds the core to CPython.

A design proposal is now posted on the NumPy developer wiki: http://projects.scipy.org/numpy/wiki/NumPyRefactoring

The write-up is a fairly high-level description of what we think the split will look like and how to deal with issues such as memory management. There are also placeholders listed as 'TBD' where more investigation is still needed and will be filled in over time. At the end of the page there is a section on the C API with a link to a function-by-function breakdown of the C API and whether the function belongs in the interface layer, the core, or need to be split between the two. All functions listed as 'core' will continue to have an interface-level wrapper with the same name to ensure source-compatibility.

All of this, particularly the interface/core function designations, is a first analysis and in flux. The goal is to get the information out and elicit discussion and feedback from the community.

A few thoughts came to mind while reading the initial writeup.

1) How is the GIL handled in the callbacks.

How to handle the GIL still requires some thought. The cleanest way, IMHO, would is for the interface layer to release the lock prior to calling into the core and then each callback function in the interface is responsible for re-acquiring it. That's straightforward to define as a rule and should work well in general, but I'm worried about potential performance issues if/when a callback is called in a loop. A few optimization points is ok, but too many and it will just be a source of heisenbugs. One other option is to just use the existing release/acquire macros in NumPy and redirect them to the interface layer. Any app that isn't CPython would just leave those callback pointers NULL. It's less disruptive but leaves some very CPython-specific behavior in the core.

2) What about error handling? That is tricky to get right, especially in C and with reference counting.

The error reporting functions in the core will likely look a lot like the CPython functions - they seem general enough. The biggest change is the CPython ones take a PyObject as the error type. 99% of the errors reported in NumPy use one of a half-dozen pre-defined types that are easy to translate. There is at least one case where an object type (complex number) is dynamically and used as the type, but so far I believe it's only one case. The reference counting does get a little more complex because a core routine will need to decref the core object on error and the interface layer will need to similarly detect the error and potentially do it's own decref. Each layer is still responsible for it's own clean up, but there are now two opportunities to introduce leaks.

3) Is there a general policy as to how the reference counting should be handled in specific functions? That is, who does the reference incrementing/decrementing?

Both layers should implement the existing policy for the objects that it manages. Essentially a function can use it's caller's reference but needs to increment the count if it's going to store it. A new instance is returned with a refcnt of 1 and the caller needs to clean it up when it's no longer needed. But that means that if the core returns a new NpyArray instance to the interface layer, the receiving function in the interface must allocate a PyObject wrapper around it and set the wrapper's refcnt to 1 before returning it. Is that what you were asking? 4) Boost has some reference counted pointers, have you looked at them? C++

is admittedly a very different animal for this sort of application.

There is also need to replace the usage of PyDict and other uses of CPython for basic data structures that aren't present in C. Having access to C++ for this and reference counting would be nice, but has the potential to break builds for everyone who use the C API. I think it's worth discussing for the future but a bigger (and possibly more contentious) change than we are able to take on for this project.

Chuck

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

On Fri, Jun 11, 2010 at 1:18 AM, Sebastien Binet <binet@cern.ch> wrote:

On Thu, 10 Jun 2010 10:47:10 -0500, Jason McCampbell <jmccampbell@enthought.com> wrote:

...

...

4) Boost has some reference counted pointers, have you looked at them? C++ is admittedly a very different animal for this sort of application.

There is also need to replace the usage of PyDict and other uses of CPython for basic data structures that aren't present in C.  Having access to C++ for this and reference counting would be nice, but has the potential to break builds for everyone who use the C API.  I think it's worth discussing for the future but a bigger (and possibly more contentious) change than we are able to take on for this project.

for the dict part, a probably good enough replacement could be the hash table from:

http://c-algorithms.sourceforge.net/

My long quest for a usable set of basic structures in C lead me to basekit. It is BSD and the code is very simple, I like it very much: http://github.com/stevedekorte/basekit It is used in the IO programming language, which is in the same ballpark as python (very dynamic ala smalltalk) and uses a GC, which means the allocation strategy used in basekit is viable for VM which do not depend on reference counting. David

Den 10.06.2010 18:48, skrev Sturla Molden:

ctypes will also make porting to other Python implementations easier (or even other languages: Ruby, JacaScript) easier. Not to mention that it will make NumPy impervious to changes in the Python C API.

Linking is also easier with ctypes. I started getting a lot of linker errors when switching to 64-bit Python on Windows. Just throwing out all of Cython et al. in favour of plain C keeps the problem away. Sturla

On Jun 10, 2010, at 11:48 AM, Sturla Molden wrote:

I have a few radical suggestions:

There are some good ideas there. I suspect we can't address all of them in the course of this re-factoring effort, but I really appreciate you putting them out there, because they are useful things to consider.

1. Use ctypes as glue to the core DLL, so we can completely forget about refcounts and similar mess. Why put manual reference counting and error handling in the core? It's stupid.

If we could get to a single core DLL, then perhaps this would work. It may be very difficult to actually get to that point though from where we are because there is a lot to the current CPython interface that would have to be re-thought. Right now, it looks like there is still a need for an interface layer.

2. The core should be a plain DLL, loadable with ctypes. (I know David Cournapeau and Robert Kern is going to hate this.) But if Python can have a custom loader for .pyd files, so can NumPy for it's core DLL. For ctypes we just need to specify a fully qualified path to the DLL, which can be read from a config file or whatever.

This approach does not build a new Python type in compiled code. There are speed disadvantages to this --- especially for the numpy scalars.

5. Allow OpenMP pragmas in the core. If arrays are above a certain size, it should switch to multi-threading.

This is an interesting idea. -Travis

On Fri, Jun 11, 2010 at 7:25 AM, Sturla Molden <sturla@molden.no> wrote:

Den 10.06.2010 22:07, skrev Travis Oliphant:

...

...

2. The core should be a plain DLL, loadable with ctypes. (I know David Cournapeau and Robert Kern is going to hate this.) But if Python can have a custom loader for .pyd files, so can NumPy for it's core DLL. For ctypes we just need to specify a fully qualified path to the DLL, which can be read from a config file or whatever.

I am not sure why you say I would hate it - that's exactly what I would like numpy core to be (a library). Ctypes is an implementation detail. For reference counting, I am not sure we can get away with it unless we use something like LUA stack technique for the core. The big issue with reference counting is how to deal with non C VM - there has to be known good practices to make libraries reusable from say the CLI or the JVM, but I don't know them.

By the way: Is Cython mature enough to toss all C out the door?

Maybe not all, but I am convinced most of it could. Today, I would not have rewritten lfilter in C, for example, I would have just used cython. One issue with cython compared to C is that it makes compilation much slower altogether, at least with gcc. I have never been able to pin-point the exact issue (it does not seem like gcc should be that slow compiling the generated C). David

Den 11.06.2010 09:14, skrev Sebastien Binet:

it of course depends on the granularity at which you wrap and use numpy-core but tight loops calling ctypes ain't gonna be pretty performance-wise.

Tight loops in Python are never pretty. The purpose of vectorization with NumPy is to avoid tight Python loops.

I really think using cython would be a better option. one'd get the python2<-> python3 transition for "free".

Perhaps. Cython has nice syntax for PEP 3118 buffers. But there are downsides to this: - Cython is not C. We want the core to be a C library independent of Python. - Cython depends on CPython. - Compiling Cython generated C takes time. - Linkage can be a PITA (I use 64-bit Python now, and import libraries for gcc are missing.) Sturla

Thu, 10 Jun 2010 18:48:04 +0200, Sturla Molden wrote: [clip]

5. Allow OpenMP pragmas in the core. If arrays are above a certain size, it should switch to multi-threading.

Some places where Openmp could probably help are in the inner ufunc loops. However, improving the memory efficiency of the data access pattern is another low-hanging fruit for multidimensional arrays. I'm not sure how Openmp plays together with thread-safety; especially if used in outer constructs. Whether this is possible to do easily is not so clear, as Numpy uses e.g. iterator and loop objects, which cannot probably be shared between different openmp threads. So one would have to do some extra work to parallelize these parts manually. But probably if multithreading is desired, openmp sounds like the way to go... -- Pauli Virtanen

On Fri, Jun 11, 2010 at 6:56 AM, Sturla Molden <sturla@molden.no> wrote:

Also about array iterators in NumPy's C base (i.e. for doing something along an axis): we don't need those. There is a different way of coding which leads to faster code.

1. Collect an array of pointers to each subarray (e.g. using std::vector<dtype*> or dtype**) 2. Dispatch on the pointer array...

Do you have the code for this ? That's something I wanted to do, but never took the time to do. Faster generic iterator would be nice, but very hard to do in general.

Another thing I did when reimplementing lfilter was "copy-in copy-out" for strided arrays.

What is copy-in copy out ? I am not familiar with this term ? David

On Thu, Jun 10, 2010 at 6:27 PM, Sturla Molden <sturla@molden.no> wrote:

Den 11.06.2010 00:57, skrev David Cournapeau:

Do you have the code for this ? That's something I wanted to do, but never took the time to do. Faster generic iterator would be nice, but very hard to do in general.

/* this computes the start adress for every vector along a dimension (axis) of an ndarray */

typedef PyArrayObject ndarray;

inline char *datapointer(ndarray *a, npy_intp *indices) { char *data = a->data; int i; npy_intp j; for (i=0; i < a->nd; i++) { j = indices[i]; data += j * a->strides[i]; } return data; }

static double ** _get_axis_pointer(npy_intp *indices, int axis, ndarray *a, double **axis_ptr_array, int dim) { /* recursive axis pointer search for 4 dimensions or more */ npy_intp i; double *axis_ptr; if (dim == a->nd) { /* done recursing dimensions, compute address of this axis */ axis_ptr = (double *) datapointer(a, indices); *axis_ptr_array = axis_ptr; return (axis_ptr_array + 1); } else if (dim == axis) { /* use first element only */ indices[dim] = 0; axis_ptr_array = _get_axis_pointer(indices, axis, a, axis_ptr_array, dim+1); return axis_ptr_array; } else { /* iterate range of indices */ npy_intp length = a->dimensions[dim]; for (i=0; i < length; i++) { indices[dim] = i; axis_ptr_array = _get_axis_pointer(indices, axis, a, axis_ptr_array, dim+1); } return axis_ptr_array; } }

static double **get_axis_pointers(ndarray *a, int axis, npy_intp *count) { npy_intp indices[NPY_MAXDIMS]; double **out, **tmp; npy_intp i, j; j = 1; for (i=0; i < a->nd; i++) { if (i != axis) j *= a->dimensions[i]; } *count = j;

out = (double **) malloc(*count * sizeof(double*)); if (out == NULL) { *count = 0; return NULL; } tmp = out; switch (a->nd) { /* for one dimension we just return the data pointer */ case 1: *tmp = (double *)a->data; break; /* specialized for two dimensions */ case 2: switch (axis) { case 0: for (i=0; i<a->dimensions[1]; i++) *tmp++ = (double *)(a->data + i*a->strides[1]); break;

case 1: for (i=0; i<a->dimensions[0]; i++) *tmp++ = (double *)(a->data + i*a->strides[0]); break; } break; /* specialized for three dimensions */ case 3: switch (axis) { case 0: for (i=0; i<a->dimensions[1]; i++) for (j=0; j<a->dimensions[2]; j++) *tmp++ = (double *)(a->data + i*a->strides[1] + j*a->strides[2]); break; case 1: for (i=0; i<a->dimensions[0]; i++) for (j=0; j<a->dimensions[2]; j++) *tmp++ = (double *)(a->data + i*a->strides[0] + j*a->strides[2]); break;

case 2: for (i=0; i<a->dimensions[0]; i++) for (j=0; j<a->dimensions[1]; j++) *tmp++ = (double *)(a->data + i*a->strides[0] + j*a->strides[1]);

} break; /* four or more dimensions: use recursion */ default: for (i=0; i<a->nd; i++) indices[i] = 0; _get_axis_pointer(indices, axis, a, out, 0);

} done: return out;

}

Another thing I did when reimplementing lfilter was "copy-in copy-out" for strided arrays.

What is copy-in copy out ? I am not familiar with this term ?

Strided memory access is slow. So it often helps to make a temporary copy that are contiguous.

But for an initial refactoring it probably falls in the category of premature optimization. Another thing to avoid on the first go around is micro-optimization, as it tends to complicate the code and often doesn't do much for performance. <snip> Chuck

Den 11.06.2010 03:02, skrev Charles R Harris:

But for an initial refactoring it probably falls in the category of premature optimization. Another thing to avoid on the first go around is micro-optimization, as it tends to complicate the code and often doesn't do much for performance.

First, to refractor NumPy we must get rid of array iterators as Python objects. Retaining them as Python objects defies the whole purpose. Second, creating an array of pointers has the advantage of being more friendly to OpenMP and multi-core CPUs. We don't need to serialize access to an iterator, and we allow OpenMP to do load balancing. For things like FFT on vectors along an axis in a multidimensional array (e.g. in multi-dimensional FFTs), this is not a micro-optimization. Not to mention that it is easier to iterate over an array than use a complicated iterator object in C. Sturla

On 06/11/2010 09:27 AM, Sturla Molden wrote:

Strided memory access is slow. So it often helps to make a temporary copy that are contiguous.

Ah, ok, I did not know this was called copy-in/copy-out, thanks for the explanation. I agree this would be a good direction to pursue, but maybe out of scope for the first refactoring, cheers, David

A Friday 11 June 2010 02:27:18 Sturla Molden escrigué:

...

...

Another thing I did when reimplementing lfilter was "copy-in copy-out" for strided arrays.

What is copy-in copy out ? I am not familiar with this term ?

Strided memory access is slow. So it often helps to make a temporary copy that are contiguous.

In my experience, this technique will only work well with strided arrays if you are going to re-use the data of these temporaries in cache, or your data is unaligned. But if you are going to use the data only once (and this is very common in NumPy element-wise operations), this is rather counter- productive for strided arrays. For example, in numexpr, we made a lot of different tests comparing "copy-in copy-out" and direct access techniques for strided arrays. The result was that operations with direct access showed significantly better performance with strided arrays. On the contrary, for unaligned arrays the copy-in copy- out technique gave better results. Look at these times, where the arrays where unidimensional with a length of 1 million element each, but the results can be extrapolated to larger, multidimensional arrays (the original benchmark file is bench/vml_timing.py): -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Numexpr version: 1.3.2.dev169 NumPy version: 1.4.1rc2 Python version: 2.6.1 (r261:67515, Feb 3 2009, 17:34:37) [GCC 4.3.2 [gcc-4_3-branch revision 141291]] Platform: linux2-x86_64 AMD/Intel CPU? True VML available? True VML/MKL version: Intel(R) Math Kernel Library Version 10.1.0 Product Build 081809.14 for Intel(R) 64 architecture applications -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= To start with, times between numpy and numexpr are very similar for very simple expressions (except for unaligned arrays, where "copy-in copy-out" works pretty well for numexpr): ******************* Expression: i2 > 0 numpy: 0.0016 numpy strided: 0.0037 numpy unaligned: 0.0086 numexpr: 0.0016 Speed-up of numexpr over numpy: 0.9512 numexpr strided: 0.0039 Speed-up of numexpr over numpy: 0.964 numexpr unaligned: 0.0042 Speed-up of numexpr over numpy: 2.0598 When doing some basic operations (mind that there are no temporaries here, so numpy should be not in great disadvantage), direct access to strided data goes between 2x and 3x faster than numpy: ******************* Expression: f3+f4 numpy: 0.0060 numpy strided: 0.0176 numpy unaligned: 0.0166 numexpr: 0.0052 Speed-up of numexpr over numpy: 1.1609 numexpr strided: 0.0086 Speed-up of numexpr over numpy: 2.0584 numexpr unaligned: 0.0099 Speed-up of numexpr over numpy: 1.6785 ******************* Expression: f3+i2 numpy: 0.0060 numpy strided: 0.0176 numpy unaligned: 0.0176 numexpr: 0.0031 Speed-up of numexpr over numpy: 1.9137 numexpr strided: 0.0061 Speed-up of numexpr over numpy: 2.8789 numexpr unaligned: 0.0078 Speed-up of numexpr over numpy: 2.2411 Notice how, until now, absolute times in numexpr and strided arrays (using the direct technique) are faster than the unaligned case (copy-in copy-out). Also, when evaluating transcendental expressions (numexpr uses Intel's Vector Math Library, VML, here), direct access is again faster than NumPy: ******************* Expression: exp(f3) numpy: 0.0150 numpy strided: 0.0155 numpy unaligned: 0.0222 numexpr: 0.0030 Speed-up of numexpr over numpy: 5.0268 numexpr strided: 0.0081 Speed-up of numexpr over numpy: 1.9086 numexpr unaligned: 0.0066 Speed-up of numexpr over numpy: 3.3454 ******************* Expression: log(exp(f3)+1)/f4 numpy: 0.0486 numpy strided: 0.0563 numpy unaligned: 0.0639 numexpr: 0.0121 Speed-up of numexpr over numpy: 4.0332 numexpr strided: 0.0170 Speed-up of numexpr over numpy: 3.3067 numexpr unaligned: 0.0164 Speed-up of numexpr over numpy: 3.8833 However, now that I see the latter figures, I don't remember that we have checked whether a copy-in copy-out technique would work faster in combination with VML. By looking at the better absolute times in unaligned arrays, I'd say chances are that performance for the strided scenario *might* benefit from using copy-in/copy-out. Mmh, that's worth a try... -- Francesc Alted

Den 11.06.2010 10:17, skrev Pauli Virtanen:

...

1. Collect an array of pointers to each subarray (e.g. using std::vector<dtype*> or dtype**) 2. Dispatch on the pointer array...

This is actually what the current ufunc code does.

The innermost dimension is handled via the ufunc loop, which is a simple for loop with constant-size step and is given a number of iterations. The array iterator objects are used only for stepping through the outer dimensions. That is, it essentially steps through your dtype** array, without explicitly constructing it.

Yes, exactly my point. And because the iterator does not explicitely construct the array, it sucks for parallel programming (e.g. with OpenMP): - The iterator becomes a bottleneck to which access must be serialized with a mutex. - We cannot do proper work scheduling (load balancing)

We could in principle use OpenMP to parallelize the outer loop, as long as the inner ufunc part is implemented in C, and does not go back into Python.

Yes. And for that we need to explicitely construct a pointer array.

But if the problem is memory-bound, it is not clear that parallelization helps.

That is often the case, yes. Memory access patterns is one of the most important issues speed wise, particularly when working with strided array slices. It would also make sence to evaluate expressions like "y = b*x + a" without a temporary array for b*x. I know roughly how to do it, but don't have time to look at it before next year. (Yes I know about numexpr, I am talking about plain Python code.) Sturla

On 11 June 2010 11:12, Benjamin Root <ben.root@ou.edu> wrote:

On Fri, Jun 11, 2010 at 8:31 AM, Sturla Molden <sturla@molden.no> wrote:

...

It would also make sence to evaluate expressions like "y = b*x + a" without a temporary array for b*x. I know roughly how to do it, but don't have time to look at it before next year. (Yes I know about numexpr, I am talking about plain Python code.)

Sturla

If I may chime in here with my own experience with NumPy code...

I typically use older, "weaker" computers for my work. I am not doing real-time modeling or some other really advanced, difficult computations. For me, NumPy works "fast enough", even on an EeePC. My main issue is the one given above by Sturla. I find that NumPy's memory usage can go out-of-control very easily in long mathematical expressions. With a mix of constants and large matricies, each step in the order of operations seems to take up more memory. Often, I would run into a major slow-down from trashing the swap. This is fairly trivial to get around by operating over slices of the matrices at a time, but -- to me -- all of this talk about optimizing the speed of the operations without addressing the temporaries issue is like trying to tune-up the engine of a car without bothering to take the lead weights out of the trunk.

I should say, though, that I've gone through the process of removing all temporary allocation using ufunc output arguments (np.add(a,b,c)) only to discover that it didn't actually save any memory and it was slower. The nice thing about temporaries is that they're, well, temporary; they go away. On the other hand, since memory reads are very slow, optimizations that do more calculation per load/store could make a very big difference, eliminating temporaries as a side effect. Anne

Just my 2 cents.

Ben Root

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Fri, 11 Jun 2010 15:31:45 +0200, Sturla Molden wrote: [clip]

...

The innermost dimension is handled via the ufunc loop, which is a simple for loop with constant-size step and is given a number of iterations. The array iterator objects are used only for stepping through the outer dimensions. That is, it essentially steps through your dtype** array, without explicitly constructing it.

Yes, exactly my point. And because the iterator does not explicitely construct the array, it sucks for parallel programming (e.g. with OpenMP):

- The iterator becomes a bottleneck to which access must be serialized with a mutex. - We cannot do proper work scheduling (load balancing)

I don't necessarily agree: you can do for parallelized outer loop { critical section { p = get iterator pointer ++iterator } inner loop in region `p` } This does allow load balancing etc., as a free processor can immediately grab the next available slice. Also, it would be easier to implement with OpenMP pragmas in the current code base. Of course, the assumption here is that the outer iterator overhead is small compared to the duration of the inner loop. This must then be compared to the memory access overhead involved in the dtype** array. -- Pauli Virtanen

Fri, 11 Jun 2010 10:29:28 +0200, Hans Meine wrote: [clip]

Ideally, algorithms would get wrapped in between two additional pre-/postprocessing steps:

1) Preprocessing: After broadcasting, transpose the input arrays such that they become C order. More specifically, sort the strides of one (e.g. the first/the largest) input array decreasingly, and apply the same transposition to the other arrays (in + out if given).

2) Perform the actual algorithm, producing a C-order output array.

3) Apply the reverse transposition to the output array. [clip]

The post/preprocessing steps are fairly easy to implement in the ufunc machinery. The problem here was that this would result to a subtle semantic change, as output arrays from ufuncs would no longer necessarily be in C order. We need to explicitly to *decide* to make this change, before proceeding. I'd say that we should simply do this, as nobody should depend on this assumption. I think I there was some code ready to implement this shuffling. I'll try to dig it out and implement the shuffling. -- Pauli Virtanen

pe, 2010-06-11 kello 10:52 +0200, Hans Meine kirjoitti:

On Friday 11 June 2010 10:38:28 Pauli Virtanen wrote: [clip]

...

I think I there was some code ready to implement this shuffling. I'll try to dig it out and implement the shuffling.

That would be great!

Ullrich Köthe has implemented this for our VIGRA/numpy bindings: http://tinyurl.com/fast-ufunc At the bottom you can see that he basically wraps all numpy.ufuncs he can find in the numpy top-level namespace automatically.

Ok, here's the branch: http://github.com/pv/numpy-work/compare/master...feature;ufunc-memory-access... Some samples: (the reported times in braces are times without the optimization) x = np.zeros([100,100]) %timeit x + x 10000 loops, best of 3: 106 us (99.1 us) per loop %timeit x.T + x.T 10000 loops, best of 3: 114 us (164 us) per loop %timeit x.T + x 10000 loops, best of 3: 176 us (171 us) per loop x = np.zeros([100,5,5]) %timeit x.T + x.T 10000 loops, best of 3: 47.7 us (61 us) per loop x = np.zeros([100,5,100]).transpose(2,0,1) %timeit np.cos(x) 100 loops, best of 3: 3.77 ms (9 ms) per loop As expected, some improvement can be seen. There's also appears to be an additional 5 us (~ 700 inner loop operations it seems) overhead coming from somewhere; perhaps this can still be reduced. -- Pauli Virtanen

ti, 2010-06-15 kello 10:10 -0400, Anne Archibald kirjoitti:

Correct me if I'm wrong, but this code still doesn't seem to make the optimization of flattening arrays as much as possible. The array you get out of np.zeros((100,100)) can be iterated over as an array of shape (10000,), which should yield very substantial speedups. Since most arrays one operates on are like this, there's potentially a large speedup here. (On the other hand, if this optimization is being done, then these tests are somewhat deceptive.)

It does perform this optimization, and unravels the loop as much as possible. If all arrays are wholly contiguous, iterators are not even used in the ufunc loop. Check the part after /* Determine how many of the trailing dimensions are contiguous */ However, in practice it seems that this typically is not a significant win -- I don't get speedups over the unoptimized numpy code even for shapes (2,)*20 where you'd think that the iterator overhead could be important: x = np.zeros((2,)*20) %timeit np.cos(x) 10 loops, best of 3: 89.9 ms (90.2 ms) per loop This is actually consistent with the result x = np.zeros((2,)*20) y = x.flatten() %timeit x+x 10 loops, best of 3: 20.9 ms per loop %timeit y+y 10 loops, best of 3: 21 ms per loop that you can probably verify with any Numpy installation. This result may actually be because the Numpy ufunc inner loops are not especially optimized -- they don't use fixed stride sizes etc., and are so not much more efficient than using array iterators.

On the other hand, it seems to me there's still some question about how to optimize execution order when the ufunc is dealing with two or more arrays with different memory layouts. In such a case, which one do you reorder in favour of?

The axis permutation is chosen so that the sum (over different arrays) of strides (for each dimension) is decreasing towards the inner loops. I'm not sure, however, that this is the optimal choice. Things may depend quite a bit on the cache architecture here.

Is it acceptable to return freshly-allocated arrays that are not C-contiguous?

Yes, it returns non-C-contiguous arrays. The contiguity is determined so that the ufunc loop itself can access the memory with a single stride. This may cause some speed loss in some expressions. Perhaps there should be a subtle bias towards C-order? But I'm not sure this is worth the bother. -- Pauli Virtanen

ti, 2010-06-15 kello 11:35 -0400, Anne Archibald kirjoitti: [clip: Numpy ufunc inner loops are slow]

This is a bit strange. I think I'd still vote for including this optimization, since one hopes the inner loop will get faster at some point. (If nothing else, the code generator can probably be made to generate specialized loops for common cases).

Ok, apparently I couldn't read. Numpy already unrolls the ufunc loop for C-contiguous arrays, that's the reason there was no difference. This should be fairer: x = np.zeros((2,)*20).transpose([1,0]+range(2,20)) y = x.flatten() %timeit x+x 10 loops, best of 3: 122 ms per loop %timeit y+y 10 loops, best of 3: 19.9 ms per loop And with the loop order optimization: %timeit x+x 10 loops, best of 3: 20 ms per loop To get an idea how things look without collapsing the loops, one can engineer nastier arrays: x = np.zeros([50*5,50*3,50*2])[:-23:5, :-29:3, 31::2] %timeit x+x 100 loops, best of 3: 2.83 ms (3.59 ms) per loop [clip: minimize sum of strides]

I suspect that it isn't optimal. As I understand it, the key restriction imposed by the cache architecture is that an entire cache line - 64 bytes or so - must be loaded at once. For strides that are larger than 64 bytes I suspect that it simply doesn't matter how big they are. (There may be some subtle issues with cache associativity but this would be extremely architecture-dependent.) So I would say, first do any dimensions in which some or all strides are less than 64 bytes, starting from the smallest. After that, any order you like is fine.

This is all up to designing a suitable objective function. Suppose `s_{i,j}` is the stride of array `i` in dimension `j`. The order-by-sum-of-strides rule reads f_j = \sum_i |s_{i,j}| and an axis permutation which makes the sequence `f_j` decreasing is chosen. Other heuristics would use a different f_j. One metric that I used previously when optimizing the reductions, was the number of elements that fit in a cache line, as you suggested. The sum of this for all arrays should be maximized in the inner loops. So one could take the objective function f_j = - C_1 C_2 \sum_i CACHELINE / (1 + C_2 |s_{i,j}|) with some constants C_1 and C_2, which, perhaps, would give better results. Note that the stride can be zero for broadcasted arrays, so we need to add a cutoff for that. With a bias toward C order, one could then have f_j = - C_3 j - C_1 C_2 \sum_i CACHELINE / (1 + C_2 |s_{i,j}|) Now the question would be just choosing C_1, C_2, and C_3 suitably. Another optimization could be flipping negative strides. This sounds a bit dangerous, though, so one would need to think if it could e.g. break a += a[::-1] etc. [clip]

I'm more worried this may violate some users' assumptions. If a user knows they need an array to be in C order, really they should use ascontiguousarray. But as it stands it's enough to make sure that it's newly-allocated as the result of an arithmetic expression. Incautious users could suddenly start seeing copies happening, or even transposed arrays being passed to, C and Fortran code.

Yes, it's a semantic change, and we have to make it consciously (and conscientiously, to boot :). I believe f2py checked the contiguity of intent(inout) arguments? The same should go for most C code that accepts Numpy arrays, at least PyArray_IsContiguous should be checked before assuming it is true. But this would still mean that code that previously worked, could start throwing up exceptions. Personally, I believe this change is worth making, with suitable mention in the release notes.

It's also worth exploring common usage patterns to make sure that numpy still gravitates towards using C contiguous arrays for everything. I'm imagining a user who at some point adds a transposed array to a normal array, then does a long series of computations on the result. We want as many of those operations as possible to operate on contiguous arrays, but it's possible that an out-of-order array could propagate indefinitely, forcing all loops to be done with one array having large strides, and resulting in output that is stil out-of-order.

I think, at present, non-C-contiguous arrays will propagate indefinitely.

Some preference for C contiguous output is worth adding.

Suggestions for better heuristics are accepted, just state it in the form of an objective function :)

It would also be valuable to build some kind of test harness to track the memory layout of the arrays generated in some "typical" calculations as well as the ordering of the loop used.

More generally, this problem is exactly what ATLAS is for - finding cache-optimal orders for linear algebra. So we shouldn't expect it to be simple.

Yes, I do not think any of us is expecting it to be simple. I don't think we can aim for the optimal solution, since it is ill-defined, but only for one that is "good enough in practice". -- Pauli Virtanen

On Tue, Jun 15, 2010 at 9:37 AM, David Cournapeau <cournape@gmail.com>wrote:

On Wed, Jun 16, 2010 at 12:16 AM, Pauli Virtanen <pav@iki.fi> wrote:

...

ti, 2010-06-15 kello 10:10 -0400, Anne Archibald kirjoitti:

...

Correct me if I'm wrong, but this code still doesn't seem to make the optimization of flattening arrays as much as possible. The array you get out of np.zeros((100,100)) can be iterated over as an array of shape (10000,), which should yield very substantial speedups. Since most arrays one operates on are like this, there's potentially a large speedup here. (On the other hand, if this optimization is being done, then these tests are somewhat deceptive.)

It does perform this optimization, and unravels the loop as much as possible. If all arrays are wholly contiguous, iterators are not even used in the ufunc loop. Check the part after

/* Determine how many of the trailing dimensions are contiguous */

However, in practice it seems that this typically is not a significant win -- I don't get speedups over the unoptimized numpy code even for shapes

(2,)*20

where you'd think that the iterator overhead could be important:

I unfortunately don't have much time to look into the code ATM, but tests should be run with different CPU. When I implemented the neighborhood iterator, I observed significant (somtimes several tens of %) differences - the gcc version also matters,

That's a common problem with trying to optimize at that level, things become architecture and compiler dependent. Reminds me a bit of an experiment where the experimenter was using genetic optimization to design a circuit on a chip and the optimal design ended up taking advantage of some stray capacitance. Chuck

On Thu, Jun 10, 2010 at 6:48 PM, Sturla Molden <sturla@molden.no> wrote:

I have a few radical suggestions:

1. Use ctypes as glue to the core DLL, so we can completely forget about refcounts and similar mess. Why put manual reference counting and error handling in the core? It's stupid.

I totally agree, I thought that the refactoring was supposed to provide simple data structures and simple algorithms to perform the C equivalents of sin,cos,exp, dot, +,-,*,/, dot, inv, ... Let me explain at an example what I expected: In the core C numpy library there would be new "numpy_array" struct with attributes numpy_array->buffer numpy_array->dtype numpy_array->ndim numpy_array->shape numpy_array->strides numpy_array->owndata etc. that replaces the current PyArrayObject which contains Python C API stuff: typedef struct PyArrayObject { PyObject_HEAD char *data; /* pointer to raw data buffer */ int nd; /* number of dimensions, also called ndim */ npy_intp *dimensions; /* size in each dimension */ npy_intp *strides; /* bytes to jump to get to the next element in each dimension */ PyObject *base; /* This object should be decref'd upon deletion of array */ /* For views it points to the original array */ /* For creation from buffer object it points to an object that shold be decref'd on deletion */ /* For UPDATEIFCOPY flag this is an array to-be-updated upon deletion of this one */ PyArray_Descr *descr; /* Pointer to type structure */ int flags; /* Flags describing array -- see below*/ PyObject *weakreflist; /* For weakreferences */ void *buffer_info; /* Data used by the buffer interface */ } PyArrayObject; Example: -------------- If one calls the following Python code x = numpy.zeros((N,M,K), dtype=float) the memory allocation would be done on the Python side. Calling a ufunc like y = numpy.sin(x) would first allocate the memory for y on the Python side and then call a C function a la numpy_unary_ufunc( double (*fcn_ptr)(double), numpy_array *x, numpy_array *y) If y is already allocated, one would call y = numpy.sin(x, out = y) Similarly z = x*y would first allocate the memory for z and then call a C function a la numpy_binary_ufunc( double (*fcn_ptr)(double, double), numpy_array *x, numpy_array *y, numpy_array *z) similarly other functions like dot: z = dot(x,y, out = z) would simply call a C function a la numpy_dot( numpy_array *x, numpy_array *y, numpy_array *z) If one wants to use numpy functions on the C side only, one would use the numpy_array struct manually. I.e. one has to do the memory management oneself in C. Which is perfectly ok since one is just interested in using the algorithms.

2. The core should be a plain DLL, loadable with ctypes. (I know David Cournapeau and Robert Kern is going to hate this.) But if Python can have a custom loader for .pyd files, so can NumPy for it's core DLL. For ctypes we just need to specify a fully qualified path to the DLL, which can be read from a config file or whatever.

That would probably the easiest way.

3. ctypes will take care of all the mess regarding the GIL. Again there is no need to re-invent the wheel. ctypes.CDLL releases the GIL when calling into C, and re-acquires before making callbacks to Python. ctypes.PyDLL keeps the GIL for legacy library code that are not threadsafe.

I have not much experience with parallelizing codes. If one implements the algorithms side effect free it should be quite easy to parallelize them on the C level, not?

ctypes will also make porting to other Python implementations easier (or even other languages: Ruby, JacaScript) easier. Not to mention that it will make NumPy impervious to changes in the Python C API.

4. Write the core in C++ or Fortran (95/2003), not C. ANSI C (aka C89). Use std::vector<> instead of malloc/free for C++, and possibly templates for generics on various dtypes.

The only reason I see for C++ is the possibility to use meta programming which is very ill-designed. I'd rather like to see some simple code preprocessing on C code than C++ template meta programming. And it should be possible to avoid mallocs in the C code, not?

5. Allow OpenMP pragmas in the core. If arrays are above a certain size, it should switch to multi-threading.

Sturla

Just my 2 cents, Sebastian

Den 10.06.2010 15:26, skrev Charles R Harris:

A few thoughts came to mind while reading the initial writeup.

1) How is the GIL handled in the callbacks. 2) What about error handling? That is tricky to get right, especially in C and with reference counting. 3) Is there a general policy as to how the reference counting should be handled in specific functions? That is, who does the reference incrementing/decrementing? 4) Boost has some reference counted pointers, have you looked at them? C++ is admittedly a very different animal for this sort of application.

Chuck

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Den 12.06.2010 15:57, skrev David Cournapeau:

Anything non trivial will require memory allocation and object ownership conventions. If the goal is interoperation with other languages and vm, you may want to use something else than plain malloc, to interact better with the allocation strategies of the host platform (reference counting, garbage collection).

We can allocate memory on the Python side e.g. using Python's bytearray, ctypes array, or a Python string. We don't need the convenstion of "ownership" if we use a Python object to store the data buffer (e.g. bytearray). It owns itself, and commit suicide when garbage collected. Sturla

Den 13.06.2010 02:39, skrev David Cournapeau:

But the point is to get rid of the python dependency, and if you don't allow any api call to allocate memory, there is not much left to implement in the core.

Memory allocation is platform dependent. A CPython version could use bytearray, other platforms could e.g. use arrays from a gc managed heap (IronPython, Jython). Because memory management is platform dependent, it does not naturally belong in the core. Having ndarrays allocate buffers with malloc() and "own" buffers they've allocated complicate things terribly. The purpose of "ownership" is to know which ndarray is to call free() on the buffer. Why go through all that pain? It's just a duplication of Python's garbage collection. Re-inventing the wheel is stupid. Let buffers be Python objects that clean themselves up. If NumPy does not allocate memory on it's own, there will be no leaks due to errors in NumPy. There is still work to do in the core, i.e. the computational loops in array operators, broadcasting, ufuncs, copying data between buffers, etc. The C functions in the core would then be called with the output array already allocated. Sturla

Den 13.06.2010 05:47, skrev David Cournapeau:

This only works in simple cases. What do you do when you don't know the output size ?

First: If you don't know, you don't know. Then you're screwed and C is not going to help. Second: If we cannot figure out how much to allocate before starting to use C (very unlikely), we can always use a callback to Python. Or we can have two C functions, the first determining the amount of allocation, the second doing the computation.

How do you deal with reallocation ?

That is platform dependent. A buffer object could e.g. have a realloc method.

How do you deal with ufunc and broadcasting buffering without allocating in the library ?

We allocate whatever we need on the Python side, possibly with a callback to Python if that's what we need. Sturla

David Cournapeau wrote:

...

In the core C numpy library there would be new "numpy_array" struct with attributes

numpy_array->buffer

Anything non trivial will require memory allocation and object ownership conventions.

I totally agree -- I've been thinking for a while about a core array data structure that you could use from C/C++, and would interact well with numpy and Python -- it would be even better if it WAS numpy. I was thinking that at the root of it would be a "data_block" object (the buffer in the above), that would have a reference counting system. It would be its own system, but hopefully be able to link to Python's easily when used with Python. we're talking C here, so trying to have a full fledged memory management would be practically inventing another language, but if the just the data blocks of memory could be managed, that would allow things like sub-arrays, and slices, and views to work well. I as partly inspired by how useless C++ std::valarrays appeared to be -- a slice of a valarray is not the same as a valarray, because the original valrray is managing the memory (and I may have that wrong) -- what this said to me is that there needs to be a central system managing the memory blocks, while the user code can manage the higher level objects. If the C core were left with absolutely no memory management, then I think we'd have a next to useless system in raw C. -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker@noaa.gov

On Sat, Jun 12, 2010 at 11:38 AM, Dag Sverre Seljebotn < dagss@student.matnat.uio.no> wrote:

Christopher Barker wrote:

...

David Cournapeau wrote:

...

...

In the core C numpy library there would be new "numpy_array" struct with attributes

numpy_array->buffer

...

Anything non trivial will require memory allocation and object ownership conventions.

I totally agree -- I've been thinking for a while about a core array data structure that you could use from C/C++, and would interact well with numpy and Python -- it would be even better if it WAS numpy.

I was thinking that at the root of it would be a "data_block" object (the buffer in the above), that would have a reference counting system. It would be its own system, but hopefully be able to link to Python's easily when used with Python.

I think taking PEP 3118, strip out the Python-specific things, and then add memory management conventions, would be a good starting point.

Simply a simple header file/struct definition and specification, which could in time become a de facto way of exporting multidimensional array data between C libraries, between Fortran and C and so on (Kurt Smith's fwrap could easily be adapted to support it). The C-NumPy would then be a library on top of this spec (mainly ufuncs operating on such structs).

The memory management conventions needs some thought, as you say, because of slices -- but a central memory allocator is not good enough because one would often be accessing memory that's allocated with other purposes in mind (and should not be deallocated, or deallocated in a special way). So refcounts + deallocator callback seems reasonable.

(Not that I'm involved in this, just my 2 cents.)

This is more the way I see things, except I would divide the bottom layer into two parts, views and memory. The memory can come from many places -- memmaps, user supplied buffers, etc. -- but we should provide a simple reference counted allocator for the default. The views correspond more to PEP 3118 and simply provide data types, dimensions, and strides, much as arrays do now. However, I would confine the data types to those available in C with a bit extra information as to precision, because. Object arrays would be a special case of pointer arrays (void pointer arrays?) and structured arrays/Unicode might be a special case of char arrays. The more complicated dtypes would then be built on top of those. Some things just won't be portable, pointers in particular, but such is life. As to languages, I think we should stay with C. C++ has much to offer for this sort of thing but would be quite a big jump and maybe not as universal as C. FORTRAN is harder to come by than C and older versions didn't have such things as unsigned integers. I really haven't used FORTRAN since the 77 version, so haven't much idea what the modern version looks like, but I do suspect we have more C programmers than FORTRAN programmers, and adding a language translation on top of a design refactoring is just going to complicate things. Chuck

If I could, I would like to throw out another possible feature that might need to be taken into consideration for designing the implementation of numpy arrays. One thing I found somewhat lacking -- if that is the right term -- is a way to convolve a numpy array with an arbitrary windowing element. I managed to do some fancy indexing approaches, but they never did seem efficient. So -- if it is at all possible -- I would like to see a design that not only vectorizes well, but could also be used be used efficiently in arbitrary (or "fancy") indexing schemes as well. Maybe what I am asking is technically infeasible... I don't know. Just something to consider. I'll get back to doing documentation... Ben Root

Charles Harris wrote:

On Sat, Jun 12, 2010 at 11:38 AM, Dag Sverre Seljebotn < dagss@student.matnat.uio.no> wrote:

...

Christopher Barker wrote:

...

David Cournapeau wrote:

...

...

In the core C numpy library there would be new "numpy_array" struct with attributes

numpy_array->buffer

...

Anything non trivial will require memory allocation and object ownership conventions.

I totally agree -- I've been thinking for a while about a core array data structure that you could use from C/C++, and would interact well with numpy and Python -- it would be even better if it WAS numpy.

I was thinking that at the root of it would be a "data_block" object (the buffer in the above), that would have a reference counting system. It would be its own system, but hopefully be able to link to Python's easily when used with Python.

I think taking PEP 3118, strip out the Python-specific things, and then add memory management conventions, would be a good starting point.

Simply a simple header file/struct definition and specification, which could in time become a de facto way of exporting multidimensional array data between C libraries, between Fortran and C and so on (Kurt Smith's fwrap could easily be adapted to support it). The C-NumPy would then be a library on top of this spec (mainly ufuncs operating on such structs).

The memory management conventions needs some thought, as you say, because of slices -- but a central memory allocator is not good enough because one would often be accessing memory that's allocated with other purposes in mind (and should not be deallocated, or deallocated in a special way). So refcounts + deallocator callback seems reasonable.

(Not that I'm involved in this, just my 2 cents.)

This is more the way I see things, except I would divide the bottom layer into two parts, views and memory. The memory can come from many places -- memmaps, user supplied buffers, etc. -- but we should provide a simple reference counted allocator for the default. The views correspond more to PEP 3118 and simply provide data types, dimensions, and strides, much as arrays do now. However, I would confine the data types to those available in C with a bit extra information as to precision, because. Object arrays would be a special case of pointer arrays (void pointer arrays?) and structured arrays/Unicode might be a special case of char arrays. The more complicated dtypes would then be built on top of those. Some things just won't be portable, pointers in particular, but such is life.

As to languages, I think we should stay with C. C++ has much to offer for this sort of thing but would be quite a big jump and maybe not as universal as C. FORTRAN is harder to come by than C and older versions didn't have such things as unsigned integers. I really haven't used FORTRAN since the 77 version, so haven't much idea what the modern version looks like, but I do suspect we have more C programmers than FORTRAN programmers, and adding a language translation on top of a design refactoring is just going to complicate things.

I'm not sure how Fortran got into this, but if it was from what I wrote: I certainly didn't suggest any part of NumPy is written in Fortran. Sorry for causing confusion. What I meant: If there's an "ndarray.h" which basically just contains the minimum C version of PEP 3118 for passing around and accessing arrays. Without any runtime dependencies -- of course, writing code using it will be much easier when using the C-NumPy runtime library, but for simply exporting data from an existing library one could do without the runtime dependency. (To be more precise about fwrap: In fwrap there's a need to pass multidimensional, flexible-size array data back and forth between (autogenerated) Fortran and (autogenerated) C. It currently defines its own structs (or extra arguments, but boils down to the same thing...) for this purpose, but if an "ndarray.h" could create a standard for passing array data then that would be a natural choice instead.) Dag Sverre

On Sat, Jun 12, 2010 at 3:57 PM, David Cournapeau <cournape@gmail.com> wrote:

On Sat, Jun 12, 2010 at 10:27 PM, Sebastian Walter <sebastian.walter@gmail.com> wrote:

...

On Thu, Jun 10, 2010 at 6:48 PM, Sturla Molden <sturla@molden.no> wrote:

...

I have a few radical suggestions:

1. Use ctypes as glue to the core DLL, so we can completely forget about refcounts and similar mess. Why put manual reference counting and error handling in the core? It's stupid.

I totally agree, I  thought that the refactoring was supposed to provide simple data structures and simple algorithms to perform the C equivalents of sin,cos,exp, dot, +,-,*,/, dot, inv, ...

Let me explain at an example what I expected:

In the core C numpy library there would be new  "numpy_array" struct with attributes

 numpy_array->buffer  numpy_array->dtype  numpy_array->ndim  numpy_array->shape  numpy_array->strides  numpy_array->owndata etc.

that replaces the current  PyArrayObject which contains Python C API stuff:

typedef struct PyArrayObject {        PyObject_HEAD        char *data;             /* pointer to raw data buffer */        int nd;                 /* number of dimensions, also called ndim */        npy_intp *dimensions;       /* size in each dimension */        npy_intp *strides;          /* bytes to jump to get to the                                   next element in each dimension */        PyObject *base;         /* This object should be decref'd                                   upon deletion of array */                                /* For views it points to the original array */                                /* For creation from buffer object it points                                   to an object that shold be decref'd on                                   deletion */                                /* For UPDATEIFCOPY flag this is an array                                   to-be-updated upon deletion of this one */        PyArray_Descr *descr;   /* Pointer to type structure */        int flags;              /* Flags describing array -- see below*/        PyObject *weakreflist;  /* For weakreferences */        void *buffer_info;      /* Data used by the buffer interface */ } PyArrayObject;

Example: --------------

If one calls the following Python code x = numpy.zeros((N,M,K), dtype=float) the memory allocation would be done on the Python side.

Calling a ufunc like y = numpy.sin(x) would first allocate the memory for y on the Python side and then call a C function a la numpy_unary_ufunc( double (*fcn_ptr)(double), numpy_array *x, numpy_array *y)

If y is already allocated, one would call y = numpy.sin(x, out = y)

Similarly z = x*y would first allocate the memory for z and then call a C function a la numpy_binary_ufunc( double (*fcn_ptr)(double, double), numpy_array *x, numpy_array *y, numpy_array *z)

similarly other functions like dot: z = dot(x,y, out = z)

would simply call a C function a la numpy_dot( numpy_array *x, numpy_array *y, numpy_array *z)

If one wants to use numpy functions on the C side only, one would use the numpy_array struct manually. I.e. one has to do the memory management oneself in C. Which is perfectly ok since one is just interested in using the algorithms.

Anything non trivial will require memory allocation and object ownership conventions. If the goal is interoperation with other languages and vm, you may want to use something else than plain malloc, to interact better with the allocation strategies of the host platform (reference counting, garbage collection).

I'm just saying that the "host platform" could do the memory management and not libnumpy. I.e. libnumpy could be just a collection of algorithms. Reimplementing half of the Python C API somehow doesn't feel right to me. Those users who like to use C++ could write a class with methods that internally call the libnumpy functions: -------------- example code ----------------- class Array{ numpy_array *_array; public: const Array operator+(Array &rhs) const { Array retval( ... arguments for the right type and dimensions of the output...); numpy_add((*this)->_array, rhs->_array, retval->_array); return retval; } }; -------------- end code ----------------- I.e. let C++ do all the memory management and type inference but the numpy core C API does the number crunching. In other languages (Python, Ruby, R, whatever) one would implement a similar class. I cannot speak for others, but something about these lines is what I'd love to see since it would make it relatively easy to use numpy functionality even in existing C/C++/R/Ruby codes. Sebastian

...

The only reason I see for C++ is the possibility to use meta programming which is very ill-designed. I'd rather like to see some simple code preprocessing on C code than C++ template meta programming.

I don't think anyone is seriously considering changing languages. Especially if interoperation is desired, C is miles ahead of C++ anyway.

David _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Den 15.06.2010 18:30, skrev Sturla Molden:

A very radical solution would be to get rid of all C, and go for a "pure Python" solution. NumPy could build up a text string with OpenCL code on the fly, and use the OpenCL driver as a "JIT compiler" for fast array expressions. Most GPUs and CPUs will support OpenCL, and thus there will be no need for a compiled language like C or Fortran for fast computation in the near future.

Paradoxically, with OpenCL, Python is going to be better for fast numerical computation than C or Fortran, because Python is better at generating an manipulating dynamic text. Sturla