[Numpy-discussion] A Cython apply_along_axis function
Keith Goodman
kwgoodman@gmail....
Fri Dec 10 11:25:25 CST 2010
On Wed, Dec 1, 2010 at 6:07 PM, Keith Goodman <kwgoodman@gmail.com> wrote:
> On Wed, Dec 1, 2010 at 5:53 PM, David <david@silveregg.co.jp> wrote:
>
>> On 12/02/2010 04:47 AM, Keith Goodman wrote:
>>> It's hard to write Cython code that can handle all dtypes and
>>> arbitrary number of dimensions. The former is typically dealt with
>>> using templates, but what do people do about the latter?
>>
>> The only way that I know to do that systematically is iterator. There is
>> a relatively simple example in scipy/signal (lfilter.c.src).
>>
>> I wonder if it would be possible to add better support for numpy
>> iterators in cython...
>
> Thanks for the tip. I'm starting to think that for now I should just
> template both dtype and ndim.
I ended up templating both dtype and axis. For the axis templating I
used two functions: looper and loop_cdef.
LOOPER
Make a 3d loop template:
>>> loop = '''
.... for iINDEX0 in range(nINDEX0):
.... for iINDEX1 in range(nINDEX1):
.... amin = MAXDTYPE
.... for iINDEX2 in range(nINDEX2):
.... ai = a[INDEXALL]
.... if ai <= amin:
.... amin = ai
.... y[INDEXPOP] = amin
.... '''
Import the looper function:
>>> from bottleneck.src.template.template import looper
Make a loop over axis=0:
>>> print looper(loop, ndim=3, axis=0)
for i1 in range(n1):
for i2 in range(n2):
amin = MAXDTYPE
for i0 in range(n0):
ai = a[i0, i1, i2]
if ai <= amin:
amin = ai
y[i1, i2] = amin
Make a loop over axis=1:
>>> print looper(loop, ndim=3, axis=1)
for i0 in range(n0):
for i2 in range(n2):
amin = MAXDTYPE
for i1 in range(n1):
ai = a[i0, i1, i2]
if ai <= amin:
amin = ai
y[i0, i2] = amin
LOOP_CDEF
Define parameters:
>>> ndim = 3
>>> dtype = 'float64'
>>> axis = 1
>>> is_reducing_function = True
Import loop_cdef:
>>> from bottleneck.src.template.template import loop_cdef
Make loop initialization code:
>>> print loop_cdef(ndim, dtype, axis, is_reducing_function)
cdef Py_ssize_t i0, i1, i2
cdef int n0 = a.shape[0]
cdef int n1 = a.shape[1]
cdef int n2 = a.shape[2]
cdef np.npy_intp *dims = [n0, n2]
cdef np.ndarray[np.float64_t, ndim=2] y = PyArray_EMPTY(2, dims,
NPY_float64, 0)
More information about the NumPy-Discussion
mailing list