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

[Mark Wiebe]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20101221/d4f4cf11/attachment.html>