[Numpy-discussion] Using NumPy iterators in C extension
Mon Jun 6 10:10:32 CDT 2011
On Mon, Jun 6, 2011 at 4:06 AM, Irwin Zaid <email@example.com>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
>>> import numpy as np
>>> a = np.zeros((2,3,4))
>>> it = np.nditer([a,None],
... op_axes=[range(a.ndim) + [-1], None],
>>> b = it.operands
>>> print a.shape, b.shape, it.shape
(2, 3, 4) (2, 3, 4, 5) (24,)
> Thanks in advance for your help! I can provide more information if
> something is unclear.
> NumPy-Discussion mailing list
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the NumPy-Discussion