[Numpy-discussion] Using NumPy iterators in C extension

[Mark Wiebe]

On Mon, Jun 6, 2011 at 4:06 AM, Irwin Zaid <irwin.zaid at physics.ox.ac.uk>wrote:

> Hi all,
>
> I am writing a C extension for NumPy that implements a custom function.
> I have a minor question about iterators, and their order of iteration,
> that I would really appreciate some help with. Here goes...
>
> My function takes a sequence of N-dimensional input arrays and a single
> (N+1)-dimensional output array. For each input array, it iterates over
> the N dimensions of that array and the leading N dimensions of the
> output array. It stores a result in the trailing dimension of the output
> array that is cumulative for the entire sequence of input arrays.
>
> Crucially, I need to be sure that the input and output elements always
> match during an iteration, otherwise the output array becomes incorrect.
> Meaning, is the kth element visited in the output array during a single
> iteration the same for each input array?
>
> So far, I have implemented this function using two iterators created
> with NpyIter_New(), one for an N-dimensional input array and one for the
> (N+1)-dimensional output array. I call NpyIter_RemoveAxis() on the
> output array iterator. I am using NPY_KEEPORDER for the NPY_ORDER flag,
> which I suspect is not what I want.
>
> How should I guarantee that I preserve order? Should I be using a
> different order flag, like NPY_CORDER? Or should I be creating a
> multi-array iterator with NpyIter_MultiNew for a single input array and
> the output array, though I don't see how to do that in my case.
>

Using a different order flag would make things consistent, but the better
solution is to use NpyIter_AdvancedNew, and provide 'op_axes' to describe
how your two arrays relate. Here's the Python equivalent of what I believe
you want to do, in this case allowing the iterator to allocate the output
for you:

>>> import numpy as np
>>>
>>> a = np.zeros((2,3,4))
>>>
>>> it = np.nditer([a,None],
...                 flags=['multi_index'],
...                 op_flags=[['readonly'],['writeonly','allocate']],
...                 op_axes=[range(a.ndim) + [-1], None],
...                 itershape=(a.shape+(5,)))
>>> it.remove_axis(it.ndim-1)
>>> it.remove_multi_index()
>>>
>>> b = it.operands[1]
>>>
>>> print a.shape, b.shape, it.shape
(2, 3, 4) (2, 3, 4, 5) (24,)

Cheers,
Mark


> Thanks in advance for your help! I can provide more information if
> something is unclear.
>
> Best,
>
> Irwin
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20110606/ed364023/attachment.html>