Dag Sverre Seljebotn skrev:

Nitpick: This will fail on large arrays. I guess numpy.npy_intp is the right type to use in this case?

Yup. You are right. Thanks. Sturla

Sturla Molden skrev:

By the way, here is a more polished version, does it look ok?

No it doesn't... Got to keep the GIL for the general case (sorting object arrays). Fixing that. SM

Hello Sturla, I had a quick look at your code. Looks fine. A few notes... In "select" you should replace numpy with np. In "_median" how can you, if n==2, use s[] if s is not defined? What if n==1? Also, I think when returning an empty array, it should be of the same type you would get in the other cases. You could replace _median with the following. Best, Luca def _median(x, inplace): assert(x.ndim == 1) n = x.shape[0] if n > 2: k = n >> 1 s = select(x, k, inplace=inplace) if n & 1: return s[k] else: return 0.5*(s[k]+s[:k].max()) elif n == 0: return np.empty(0, dtype=x.dtype) elif n == 2: return 0.5*(x[0]+x[1]) else: # n == 1 return x[0]

Sturla Molden wrote:

Dag Sverre Seljebotn skrev:

...

Nitpick: This will fail on large arrays. I guess numpy.npy_intp is the right type to use in this case?

By the way, here is a more polished version, does it look ok?

http://projects.scipy.org/numpy/attachment/ticket/1213/generate_qselect.py http://projects.scipy.org/numpy/attachment/ticket/1213/quickselect.pyx

I didn't look at the algorithm, but the types look OK (except for the gil as you say). Comments: a) Is the cast to numpy.npy_intp really needed? I'm pretty sure shape is defined as numpy.npy_intp*. b) If you want higher performance with contiguous arrays (which occur a lot as inplace=False is default I guess) you can do np.ndarray[T, ndim=1, mode="c"] to tell the compiler the array is contiguous. That doubles the number of function instances though...

Cython needs something like Java's generics by the way :-)

Yes, we all long for that. It will come as soon as somebody volunteers I suppose -- it shouldn't be all that difficult, but I don't think any of the existing devs will be up for it any time soon. Dag Sverre

On Wed, 2 Sep 2009, Dag Sverre Seljebotn wrote:

Sturla Molden wrote:

...

Dag Sverre Seljebotn skrev:

...

Nitpick: This will fail on large arrays. I guess numpy.npy_intp is the right type to use in this case?

By the way, here is a more polished version, does it look ok?

http://projects.scipy.org/numpy/attachment/ticket/1213/generate_qselect.py http://projects.scipy.org/numpy/attachment/ticket/1213/quickselect.pyx

I didn't look at the algorithm, but the types look OK (except for the gil as you say). Comments:

a) Is the cast to numpy.npy_intp really needed? I'm pretty sure shape is defined as numpy.npy_intp*. b) If you want higher performance with contiguous arrays (which occur a lot as inplace=False is default I guess) you can do

np.ndarray[T, ndim=1, mode="c"]

to tell the compiler the array is contiguous. That doubles the number of function instances though...

...

Cython needs something like Java's generics by the way :-)

Yes, we all long for that. It will come as soon as somebody volunteers I suppose -- it shouldn't be all that difficult, but I don't think any of the existing devs will be up for it any time soon.

Danilo's C++ project has some baby steps in that direction, though it'll need to be expanded quite a bit to handle this. - Robert

I am very interested in this result. I have wanted to know how to do an apply_along_axis function for a while now. On Tue, Nov 30, 2010 at 11:21 AM, Keith Goodman wrote:

On Tue, Sep 1, 2009 at 2:37 PM, Sturla Molden wrote:

...

Dag Sverre Seljebotn skrev:

...

Nitpick: This will fail on large arrays. I guess numpy.npy_intp is the right type to use in this case?

By the way, here is a more polished version, does it look ok?

http://projects.scipy.org/numpy/attachment/ticket/1213/generate_qselect.py

...

http://projects.scipy.org/numpy/attachment/ticket/1213/quickselect.pyx

This is my favorite numpy/scipy ticket. So I am happy that I can contribute in a small way by pointing out a bug. The search for the k-th smallest element is only done over the first k elements (that's the bug) instead of over the entire array. Specifically "while l < k" should be "while l < r".

I added a median function to the Bottleneck package: https://github.com/kwgoodman/bottleneck

Timings:

...

...

import bottleneck as bn arr = np.random.rand(100, 100) timeit np.median(arr) 1000 loops, best of 3: 762 us per loop timeit bn.median(arr) 10000 loops, best of 3: 198 us per loop

What other functions could be built from a selection algorithm?

nanmedian scoreatpercentile quantile knn select others?

But before I add more functions to the package I need to figure out how to make a cython apply_along_axis function. For the first release I am hand coding the 1d, 2d, and 3d cases. Boring to write, hard to maintain, and doesn't solve the nd case.

Does anyone have a cython apply_along_axis that takes a cython reducing function as input? The ticket has an example but I couldn't get it to run. If no one has one (the horror!) I'll begin to work on one sometime after the first release. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

Hi, On Tue, Nov 30, 2010 at 11:35 AM, Keith Goodman wrote:

On Tue, Nov 30, 2010 at 11:25 AM, John Salvatier wrote:

...

I am very interested in this result. I have wanted to know how to do an

My first thought was to write the reducing function like this

cdef np.float64_t namean(np.ndarray[np.float64_t, ndim=1] a):

but cython doesn't allow np.ndarray in a cdef.

Sorry for the ill-considered hasty reply, but do you mean that this: import numpy as np cimport numpy as cnp cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a): return np.nanmean(a) # just a placeholder is not allowed? It works for me. Is it a cython version thing? (I've got 0.13), See you, Matthew

I think last time I looked into how to apply a function along an axis I thought that the PyArray_IterAllButAxis would not work for that task ( http://docs.scipy.org/doc/numpy/reference/c-api.array.html#PyArray_IterAllBu...), but I think perhaps I misunderstood it. I'm looking into how to use it. On Tue, Nov 30, 2010 at 12:06 PM, Keith Goodman wrote:

On Tue, Nov 30, 2010 at 11:58 AM, Matthew Brett wrote:

...

Hi,

On Tue, Nov 30, 2010 at 11:35 AM, Keith Goodman wrote:

...

On Tue, Nov 30, 2010 at 11:25 AM, John Salvatier wrote:

...

I am very interested in this result. I have wanted to know how to do an

My first thought was to write the reducing function like this

cdef np.float64_t namean(np.ndarray[np.float64_t, ndim=1] a):

but cython doesn't allow np.ndarray in a cdef.

Sorry for the ill-considered hasty reply, but do you mean that this:

import numpy as np cimport numpy as cnp

cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a): return np.nanmean(a) # just a placeholder

is not allowed? It works for me. Is it a cython version thing? (I've got 0.13),

Oh, that's nice! I'm using 0.11.2. OK, time to upgrade. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

@Keith Goodman I think I figured it out. I believe something like the following will do what you want, iterating across one axis specially, so it can apply a median function along an axis. This code in particular is for calculating a moving average and seems to work (though I haven't checked my math). Let me know if you find any problems. def ewma(a, d, int axis = -1): out = np.empty(a.shape, dtype) cdef np.flatiter ita, ito ita = np.PyArray_IterAllButAxis(a, &axis) ito = np.PyArray_IterAllButAxis(out, &axis) cdef int i cdef int axis_length = a.shape[axis] cdef int a_axis_stride = a.strides[axis]/a.itemsize cdef int o_axis_stride = out.strides[axis]/out.itemsize cdef double avg = 0.0 cdef double weight = 1.0 - np.exp(-d) while np.PyArray_ITER_NOTDONE(ita): avg = 0.0 for i in range(axis_length): avg += (np.PyArray_ITER_DATA (ita))[i * a_axis_stride ] * weight + avg * (1 - weight) (np.PyArray_ITER_DATA (ito))[i * o_axis_stride ] = avg np.PyArray_ITER_NEXT(ita) np.PyArray_ITER_NEXT(ito) return out On Tue, Nov 30, 2010 at 9:46 PM, Felix Schlesinger wrote:

...

...

import numpy as np cimport numpy as cnp

...

cdef cnp.float64_t namean(cnp.ndarray[cnp.float64_t, ndim=1] a): return np.nanmean(a) # just a placeholder

...

is not allowed? It works for me. Is it a cython version thing? (I've got 0.13),

Oh, that's nice! I'm using 0.11.2. OK, time to upgrade.

Oh wow, does that mean that http://trac.cython.org/cython_trac/ticket/177 is fixed? I couldn't find anything in the release notes about that, but it would be great news. Does the cdef function acquire and hold the buffer?

Felix

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

On Wed, Sep 2, 2009 at 1:25 PM, Chad Netzer wrote:

On Mon, Aug 31, 2009 at 9:06 PM, Sturla Molden wrote:

...

We recently has a discussion regarding an optimization of NumPy's median to average O(n) complexity. After some searching, I found out there is a selection algorithm competitive in speed with Hoare's quick select. It has the advantage of being a lot simpler to implement. In plain Python:

...

Chad, you can continue to write quick select using NumPy's C quick sort in numpy/core/src/_sortmodule.c.src. When you are done, it might be about 10% faster than this. :-)

I was sick for a bit last week, so got stalled on my version, but I'll be working on it this weekend. I'm going for a more general partition function, that could have slightly more general use cases than just a median. Nevertheless, its good to see there could be several options, hopefully at least one of which can be put into numpy.

By the way, as far as I can tell, the above algorithm is exactly the same idea as a non-recursive Hoare (ie. quicksort) selection: Do the partition, then only proceed to the sub-partition that must contain the nth element. My version is a bit more general, allowing partitioning on a range of elements rather than just one, but the concept is the same. The numpy quicksort already does non recursive sorting.

I'd also like to, if possible, have a specialized 2D version, since image media filtering is one of my interests, and the C version works on 1D (raveled) arrays only.

There are special hardwired medians for 2,3,5,9 elements, which covers a lot of image processing. They aren't in numpy, though ;) David has implemented a NeighborhoodIter that could help extract the elements if you want to deal with images. Chuck

On Thu, Sep 3, 2009 at 00:09, Sturla Molden wrote:

Chad Netzer skrev:

...

I'd also like to, if possible, have a specialized 2D version, since image media filtering is one of my interests, and the C version works on 1D (raveled) arrays only. I agree. NumPy (or SciPy) could have a select module similar to the sort module. If the select function takes an axis argument similar to the sort functions, only a small change to the current np.median would needed.

Take a look at this:

http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx

Here is a select function that takes an axis argument. There are specialized versions for 1D, 2D, and 3D. Input can be contiguous or not. For 4D and above, axes are found by recursion on the shape array. Thus it should be fast regardless of dimensions.

When he is talking about 2D, I believe he is referring to median filtering rather than computing the median along an axis. I.e., replacing each pixel with the median of a specified neighborhood around the pixel. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco

Chad Netzer wrote:

But Charles Harris's earlier suggestion of some hard coded medians for common filter template sizes (ie 3x3, 5x5, etc.) may be a nice addition to scipy, especially if it can be generalized somewhat to other filters.

For 2D images try looking into PIL : ImageFilter.MedianFilter Cheers, Jon

On Thu, Sep 3, 2009 at 12:14 AM, Sturla Molden wrote:

Chad Netzer skrev:

...

That's right, Robert. Basically, I meant doing a median on a square (or rectangle) "view" of an array, without first having to ravel(), thus generally saving a copy. But actually, since my selection based median overwrites the source array, it may not save a copy anyway.

Avoiding copies of timy buffers is futile optimization.

QuickSelect is overkill for tiny buffers like common filter kernels. Insertion sort is fine.

Shell sort is nice for sizes in the awkward range. Chuck