On Mon, Jun 6, 2011 at 4:06 AM, Irwin Zaid <[email protected]>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 > [email protected] > http://mail.scipy.org/mailman/listinfo/numpy-discussion >
_______________________________________________ NumPy-Discussion mailing list [email protected] http://mail.scipy.org/mailman/listinfo/numpy-discussion
