Giving numpy the ability to multi-iterate excluding an axis
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? Best Regards, John
On Mon, Dec 20, 2010 at 1:42 PM, John Salvatier <jsalvati@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
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. John On Tue, Dec 21, 2010 at 5:12 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Mon, Dec 20, 2010 at 1:42 PM, John Salvatier <jsalvati@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Wed, Dec 22, 2010 at 12:44 AM, John Salvatier <jsalvati@u.washington.edu>wrote:
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."
Thanks, I've tried to clean up the description a bit.
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.
I added an example to the NEP to try to make it more clear, here's what I wrote: To help understand how the offsets work, here is a simple nested iteration example. Let's say our array a has shape (2, 3, 4), and strides (48, 16, 4). The data pointer for element (i, j, k) is at address PyArray_BYTES(a) + 48*i + 16*j + 4*k. Now consider two iterators with custom op_axes (0,1) and (2,). The first one will produce addresses like PyArray_BYTES(a) + 48*i + 16*j, and the second one will produce addresses likePyArray_BYTES(a) + 4*k. Simply adding together these values would produce invalid pointers. Instead, we can make the outer iterator produce offsets, in which case it will produce the values 48*i + 16*j, and its sum with the other iterator's pointer gives the correct data address. It's important to note that this will not work if any of the iterators share an axis. The iterator cannot check this, so your code must handle it. Additionally, taking a look at the ndarray strides documentation might help: http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.strides.ht... Cheers, Mark
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API? On Tue, Dec 21, 2010 at 5:12 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Mon, Dec 20, 2010 at 1:42 PM, John Salvatier <jsalvati@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@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Sat, Jan 1, 2011 at 11:23 AM, John Salvatier <jsalvati@u.washington.edu>wrote:
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API?
I've committed Python exposure for nested iteration to the new_iterator branch. In doing so, I also changed the mechanism in C. I found that it was simpler to expose to Python if I added a Reset function which gives new base data pointers, and this also simplifies C code using nested iterators. The Python code a = arange(2).reshape(2,1) b = arange(3).reshape(1,3) i, j = np.nested_iters([a,b], [[0],[1]]) for x in i: print "inner:" for y in j: print y[0], y[1] gives inner: 0 0 0 1 0 2 inner: 1 0 1 1 1 2 and C code for nested iteration looks something like this: NpyIter *iter1, *iter1; NpyIter_IterNext_Fn iternext1, iternext2; char **dataptrs1; /* * With the exact same operands, no copies allowed, and * no axis in op_axes used both in iter1 and iter2. * Buffering may be enabled for iter2, but not for iter1. */ iter1 = ...; iter2 = ...; iternext1 = NpyIter_GetIterNext(iter1); iternext2 = NpyIter_GetIterNext(iter2); dataptrs1 = NpyIter_GetDataPtrArray(iter1); do { NpyIter_ResetBasePointers(iter2, dataptrs1); do { /* Use the iter2 values */ } while (iternext2(iter2)); } while (iternext1(iter1)); Cheers, Mark
Wow, great! I'm excited to try this. I think your patch significantly increases the extendability of numpy. Is the C-API exposed currently? Can you use the API from Cython (meaning is the numpy.pxd file updated)? On Tue, Jan 4, 2011 at 12:04 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Sat, Jan 1, 2011 at 11:23 AM, John Salvatier <jsalvati@u.washington.edu
wrote:
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API?
I've committed Python exposure for nested iteration to the new_iterator branch. In doing so, I also changed the mechanism in C. I found that it was simpler to expose to Python if I added a Reset function which gives new base data pointers, and this also simplifies C code using nested iterators.
The Python code
a = arange(2).reshape(2,1) b = arange(3).reshape(1,3)
i, j = np.nested_iters([a,b], [[0],[1]]) for x in i: print "inner:" for y in j: print y[0], y[1]
gives
inner: 0 0 0 1 0 2 inner: 1 0 1 1 1 2
and C code for nested iteration looks something like this:
NpyIter *iter1, *iter1; NpyIter_IterNext_Fn iternext1, iternext2; char **dataptrs1;
/* * With the exact same operands, no copies allowed, and * no axis in op_axes used both in iter1 and iter2. * Buffering may be enabled for iter2, but not for iter1. */ iter1 = ...; iter2 = ...;
iternext1 = NpyIter_GetIterNext(iter1); iternext2 = NpyIter_GetIterNext(iter2); dataptrs1 = NpyIter_GetDataPtrArray(iter1);
do { NpyIter_ResetBasePointers(iter2, dataptrs1); do { /* Use the iter2 values */ } while (iternext2(iter2)); } while (iternext1(iter1));
Cheers, Mark
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, Jan 4, 2011 at 12:15 PM, John Salvatier <jsalvati@u.washington.edu>wrote:
Wow, great! I'm excited to try this. I think your patch significantly increases the extendability of numpy.
Is the C-API exposed currently? Can you use the API from Cython (meaning is the numpy.pxd file updated)?
The C-API isn't exposed yet, but that won't be too difficult since it's mostly a matter of adding all the functions to the arrays in the python setup files. I thought I might do that and look at plugging it into numexpr at the same time, since to be able to use the iterator's buffering and numexpr's multithreading together will require some small additions to the iterator. Cheers, Mark On Tue, Jan 4, 2011 at 12:04 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Sat, Jan 1, 2011 at 11:23 AM, John Salvatier < jsalvati@u.washington.edu> wrote:
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API?
I've committed Python exposure for nested iteration to the new_iterator branch. In doing so, I also changed the mechanism in C. I found that it was simpler to expose to Python if I added a Reset function which gives new base data pointers, and this also simplifies C code using nested iterators.
The Python code
a = arange(2).reshape(2,1) b = arange(3).reshape(1,3)
i, j = np.nested_iters([a,b], [[0],[1]]) for x in i: print "inner:" for y in j: print y[0], y[1]
gives
inner: 0 0 0 1 0 2 inner: 1 0 1 1 1 2
and C code for nested iteration looks something like this:
NpyIter *iter1, *iter1; NpyIter_IterNext_Fn iternext1, iternext2; char **dataptrs1;
/* * With the exact same operands, no copies allowed, and * no axis in op_axes used both in iter1 and iter2. * Buffering may be enabled for iter2, but not for iter1. */ iter1 = ...; iter2 = ...;
iternext1 = NpyIter_GetIterNext(iter1); iternext2 = NpyIter_GetIterNext(iter2); dataptrs1 = NpyIter_GetDataPtrArray(iter1);
do { NpyIter_ResetBasePointers(iter2, dataptrs1); do { /* Use the iter2 values */ } while (iternext2(iter2)); } while (iternext1(iter1));
Cheers, Mark
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Oh, and I'm not sure about Cython, since I've never looked into its details. I imagine Cython will want to short circuit some of the Python exposure code, since accessing the iterator values creates new array objects. -Mark On Tue, Jan 4, 2011 at 12:59 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Tue, Jan 4, 2011 at 12:15 PM, John Salvatier <jsalvati@u.washington.edu
wrote:
Wow, great! I'm excited to try this. I think your patch significantly increases the extendability of numpy.
Is the C-API exposed currently? Can you use the API from Cython (meaning is the numpy.pxd file updated)?
The C-API isn't exposed yet, but that won't be too difficult since it's mostly a matter of adding all the functions to the arrays in the python setup files. I thought I might do that and look at plugging it into numexpr at the same time, since to be able to use the iterator's buffering and numexpr's multithreading together will require some small additions to the iterator.
Cheers, Mark
On Tue, Jan 4, 2011 at 12:04 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Sat, Jan 1, 2011 at 11:23 AM, John Salvatier < jsalvati@u.washington.edu> wrote:
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API?
I've committed Python exposure for nested iteration to the new_iterator branch. In doing so, I also changed the mechanism in C. I found that it was simpler to expose to Python if I added a Reset function which gives new base data pointers, and this also simplifies C code using nested iterators.
The Python code
a = arange(2).reshape(2,1) b = arange(3).reshape(1,3)
i, j = np.nested_iters([a,b], [[0],[1]]) for x in i: print "inner:" for y in j: print y[0], y[1]
gives
inner: 0 0 0 1 0 2 inner: 1 0 1 1 1 2
and C code for nested iteration looks something like this:
NpyIter *iter1, *iter1; NpyIter_IterNext_Fn iternext1, iternext2; char **dataptrs1;
/* * With the exact same operands, no copies allowed, and * no axis in op_axes used both in iter1 and iter2. * Buffering may be enabled for iter2, but not for iter1. */ iter1 = ...; iter2 = ...;
iternext1 = NpyIter_GetIterNext(iter1); iternext2 = NpyIter_GetIterNext(iter2); dataptrs1 = NpyIter_GetDataPtrArray(iter1);
do { NpyIter_ResetBasePointers(iter2, dataptrs1); do { /* Use the iter2 values */ } while (iternext2(iter2)); } while (iternext1(iter1));
Cheers, Mark
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Cython just has interfaces to the C-API, I think. On Tue, Jan 4, 2011 at 1:01 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
Oh, and I'm not sure about Cython, since I've never looked into its details. I imagine Cython will want to short circuit some of the Python exposure code, since accessing the iterator values creates new array objects.
-Mark
On Tue, Jan 4, 2011 at 12:59 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Tue, Jan 4, 2011 at 12:15 PM, John Salvatier < jsalvati@u.washington.edu> wrote:
Wow, great! I'm excited to try this. I think your patch significantly increases the extendability of numpy.
Is the C-API exposed currently? Can you use the API from Cython (meaning is the numpy.pxd file updated)?
The C-API isn't exposed yet, but that won't be too difficult since it's mostly a matter of adding all the functions to the arrays in the python setup files. I thought I might do that and look at plugging it into numexpr at the same time, since to be able to use the iterator's buffering and numexpr's multithreading together will require some small additions to the iterator.
Cheers, Mark
On Tue, Jan 4, 2011 at 12:04 PM, Mark Wiebe <mwwiebe@gmail.com> wrote:
On Sat, Jan 1, 2011 at 11:23 AM, John Salvatier < jsalvati@u.washington.edu> wrote:
This thread is a bit old, but since it's not possible to use the C-API is possible to accomplish this same thing with the Python API?
I've committed Python exposure for nested iteration to the new_iterator branch. In doing so, I also changed the mechanism in C. I found that it was simpler to expose to Python if I added a Reset function which gives new base data pointers, and this also simplifies C code using nested iterators.
The Python code
a = arange(2).reshape(2,1) b = arange(3).reshape(1,3)
i, j = np.nested_iters([a,b], [[0],[1]]) for x in i: print "inner:" for y in j: print y[0], y[1]
gives
inner: 0 0 0 1 0 2 inner: 1 0 1 1 1 2
and C code for nested iteration looks something like this:
NpyIter *iter1, *iter1; NpyIter_IterNext_Fn iternext1, iternext2; char **dataptrs1;
/* * With the exact same operands, no copies allowed, and * no axis in op_axes used both in iter1 and iter2. * Buffering may be enabled for iter2, but not for iter1. */ iter1 = ...; iter2 = ...;
iternext1 = NpyIter_GetIterNext(iter1); iternext2 = NpyIter_GetIterNext(iter2); dataptrs1 = NpyIter_GetDataPtrArray(iter1);
do { NpyIter_ResetBasePointers(iter2, dataptrs1); do { /* Use the iter2 values */ } while (iternext2(iter2)); } while (iternext1(iter1));
Cheers, Mark
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
participants (2)
-
John Salvatier -
Mark Wiebe