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?
I'm trying to take baby steps towards writing an apply_along_axis function that takes as input a cython function, numpy array, and axis. I'm using the following numpy ticket as a guide but I'm really just copying and pasting without understanding:
http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx
Can anyone spot why I get a segfault on the call to nanmean_1d in apply_along_axis?
import numpy as np cimport numpy as np import cython
cdef double NAN = <double> np.nan ctypedef np.float64_t (*func_t)(void *buf, np.npy_intp size, np.npy_intp s)
def apply_along_axis(np.ndarray[np.float64_t, ndim=1] a, int axis):
cdef func_t nanmean_1d cdef np.npy_intp stride, itemsize cdef int ndim = a.ndim cdef np.float64_t out
itemsize = <np.npy_intp> a.itemsize
if ndim == 1: stride = a.strides[0] // itemsize # convert stride bytes --> items out = nanmean_1d(a.data, a.shape[0], stride) else: raise ValueError("Not yet coded")
return out
cdef np.float64_t nanmean_1d(void *buf, np.npy_intp n, np.npy_intp s): "nanmean of buffer." cdef np.float64_t *a = <np.float64_t *> buf # cdef np.npy_intp i, count = 0 cdef np.float64_t asum, ai if s == 1: for i in range(n): ai = a[i] if ai == ai: asum += ai count += 1 else: for i in range(n): ai = a[i*s] if ai == ai: asum += ai count += 1 if count > 0: return asum / count else: return NAN
On 12/01/2010 08:47 PM, 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?
What you typically do is to use the C-level iterator API. In fact there's a recent thread on cython-users that does exactly this ("How can I use PyArray_IterAllButAxis..."). Of course, make sure you take the comments of that thread into account (!).
I feel that is easier to work with than what you do below. Not saying it couldn't be easier, but it's not too bad once you get used to it.
Dag Sverre
I'm trying to take baby steps towards writing an apply_along_axis function that takes as input a cython function, numpy array, and axis. I'm using the following numpy ticket as a guide but I'm really just copying and pasting without understanding:
http://projects.scipy.org/numpy/attachment/ticket/1213/_selectmodule.pyx
Can anyone spot why I get a segfault on the call to nanmean_1d in apply_along_axis?
import numpy as np cimport numpy as np import cython
cdef double NAN =<double> np.nan ctypedef np.float64_t (*func_t)(void *buf, np.npy_intp size, np.npy_intp s)
def apply_along_axis(np.ndarray[np.float64_t, ndim=1] a, int axis):
cdef func_t nanmean_1d cdef np.npy_intp stride, itemsize cdef int ndim = a.ndim cdef np.float64_t out itemsize =<np.npy_intp> a.itemsize if ndim == 1: stride = a.strides[0] // itemsize # convert stride bytes --> items out = nanmean_1d(a.data, a.shape[0], stride) else: raise ValueError("Not yet coded") return outcdef np.float64_t nanmean_1d(void *buf, np.npy_intp n, np.npy_intp s): "nanmean of buffer." cdef np.float64_t *a =<np.float64_t *> buf # cdef np.npy_intp i, count = 0 cdef np.float64_t asum, ai if s == 1: for i in range(n): ai = a[i] if ai == ai: asum += ai count += 1 else: for i in range(n): ai = a[i*s] if ai == ai: asum += ai count += 1 if count> 0: return asum / count else: return NAN _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi Keith,
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...
cheers,
David
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.
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. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
I enthusiastically support better iterator support for cython
On Wed, Dec 1, 2010 at 6:09 PM, John Salvatier jsalvati@u.washington.edu wrote:
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. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
I enthusiastically support better iterator support for cython
I enthusiastically welcome contributions along this line.
- Robert
On 12/02/2010 08:17 AM, Robert Bradshaw wrote:
On Wed, Dec 1, 2010 at 6:09 PM, John Salvatier jsalvati@u.washington.edu wrote:
On Wed, Dec 1, 2010 at 6:07 PM, Keith Goodmankwgoodman@gmail.com wrote:
On Wed, Dec 1, 2010 at 5:53 PM, Daviddavid@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. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
I enthusiastically support better iterator support for cython
I enthusiastically welcome contributions along this line.
Me too :-)
I guess we're moving into more Cython-list territory, so let's move any follow-ups there (posting this one both places).
Just in case anybody is wondering what something like this could look like, here's a rough scetch complete with bugs. The idea would be to a) add some rudimentary support for using the yield keyword in Cython to make a generator function, b) inline the generator function if the generator is used directly in a for-loop. This should result in very efficient code, and would also be much easier to implement than a general purpose generator.
@cython.inline cdef array_iter_double(np.ndarray a, int axis=-1): cdef np.flatiter it ita = np.PyArray_IterAllButAxis(a, &axis) cdef Py_ssize_t stride = a.strides[axis], length = a.shape[axis], i while np.PyArray_ITER_NOTDONE(ita): for i in range(length): yield <double*>(np.PyArray_ITER_DATA(it) + )[i * stride])[0] # TODO: Probably yield indices as well np.PyArray_ITER_NEXT(it) # TODO: add faster special-cases for stride == sizeof(double)
# Use NumPy iterator API to sum all values of array with # arbitrary number of dimensions: cdef double s = 0, value for value in array_iter_double(myarray): s += value # at this point, the contents of the array_iter_double function is copied, # and "s += value" simply inserted everywhere "yield" occurs in the function
Dag Sverre
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)
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?
There are number of ways to do it. NumPy's C API has an iterator that returns an axis on demand. Mine just collects an array with pointers to the first element in each axis. The latter is more friendly to parallelization (OpenMP or Python threads with released GIL), which is why I wrote it, otherwise it has no advantage over NumPy's.
Sturla
participants (6)
-
Dag Sverre Seljebotn -
David -
John Salvatier -
Keith Goodman -
Robert Bradshaw -
Sturla Molden