[Numpy-discussion] A Cython apply_along_axis function

[Dag Sverre Seljebotn]

On 12/02/2010 08:17 AM, Robert Bradshaw wrote:
> On Wed, Dec 1, 2010 at 6:09 PM, John Salvatier
> <jsalvati at u.washington.edu>  wrote:
>    
>> On Wed, Dec 1, 2010 at 6:07 PM, Keith Goodman<kwgoodman at gmail.com>  wrote:
>>      
>>> On Wed, Dec 1, 2010 at 5:53 PM, David<david at 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 at 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