[Numpy-discussion] Giving numpy the ability to multi-iterate excluding an axis

John Salvatier jsalvati at u.washington.edu
Wed Dec 22 03:44:28 EST 2010

This now makes sense to me, and I think it should work :D. This is all very
cool. This is going to do big things for cython and numpy.

Some hopefully constructive criticism:

When first reading through the API description, the way oa_ndim and oa_axes
work is not clear. I think your description would be clearer if you explain
what oa_ndim means (I gather something like "the number of axes over which
you wish to iterate"), currently it just says "These parameters let you
control in detail how the axes of the operand arrays get matched together
and iterated."

It's also not totally clear to me how offsetting works. What are the offsets
measured from? It seems like they are measured from another iterator, but
I'm not sure and I don't see how it gets that information.


On Tue, Dec 21, 2010 at 5:12 PM, Mark Wiebe <mwwiebe at gmail.com> wrote:

> On Mon, Dec 20, 2010 at 1:42 PM, John Salvatier <jsalvati at u.washington.edu
> > wrote:
>> A while ago, I asked a whether it was possible to multi-iterate over
>> several ndarrays but exclude a certain axis(
>> http://www.mail-archive.com/numpy-discussion@scipy.org/msg29204.html),
>> sort of a combination of PyArray_IterAllButAxis and PyArray_MultiIterNew. My
>> goal was to allow creation of relatively complex ufuncs that can allow
>> reduction or directionally dependent computation and still use broadcasting
>> (for example a moving averaging ufunc that can have changing averaging
>> parameters). I didn't get any solutions, which I take to mean that no one
>> knew how to do this.
>> I am thinking about trying to make a numpy patch with this functionality,
>> and I have some questions: 1) How difficult would this kind of task be for
>> someone with non-expert C knowledge and good numpy knowledge? 2) Does anyone
>> have advice on how to do this kind of thing?
> You may be able to do what you would like with the new iterator I've
> written.  In particular, it supports nesting multiple iterators by providing
> either pointers or offsets, and allowing you to specify any subset of the
> axes to iterate.  Here's how the code to do this in a simple 3D case might
> look, for making axis 1 the inner loop:
> PyArrayObject *op[2] = {a,b};
> npy_intp axes_outer[2] = {0,2}};
> npy_intp *op_axes[2];
> npy_intp axis_inner = 1;
> npy_int32 flags[2] = {NPY_ITER_READONLY, NPY_ITER_READONLY};
> NpyIter *outer, *inner;
> NpyIter_IterNext_Fn oiternext, iiternext;
> npy_intp *ooffsets;
> char **idataptrs;
> op_axes[0] = op_axes[1] = axes_outer;
> outer = NpyIter_MultiNew(2, op, NPY_ITER_OFFSETS,
>                            NPY_KEEPORDER, NPY_NO_CASTING, flags, NULL, 2,
> op_axes, 0);
> op_axes[0] = op_axes[1] = &axis_inner;
> inner = NpyIter_MultiNew(2, op, 0, NPY_KEEPORDER, NPY_NO_CASTING, flags,
> NULL, 1, op_axes, 0);
> oiternext = NpyIter_GetIterNext(outer);
> iiternext = NpyIter_GetIterNext(inner);
> ooffsets = (npy_intp *)NpyIter_GetDataPtrArray(outer);
> idataptrs = NpyIter_GetDataPtrArray(inner);
> do {
>    do {
>       char *a_data = idataptrs[0] + ooffsets[0], *b_data = idataptrs[0] +
> ooffsets[0];
>       /* Do stuff with the data */
>    } while(iiternext());
>    NpyIter_Reset(inner);
> } while(oiternext());
> NpyIter_Deallocate(outer);
> NpyIter_Deallocate(inner);
> Extending to more dimensions, or making both the inner and outer loops have
> multiple dimensions, isn't too crazy.  Is this along the lines of what you
> need?
> If you check out my code, note that it currently isn't exposed as NumPy API
> yet, but you can try a lot of things with the Python exposure.
> Cheers,
> Mark
> _______________________________________________
> 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/20101222/3e401b2b/attachment.html>

More information about the NumPy-Discussion mailing list