<div dir="ltr">Jaime brought up the same issue recently, along with some other issues for ufunc.reduceat:<div><a href="https://mail.scipy.org/pipermail/numpy-discussion/2016-March/075199.html">https://mail.scipy.org/pipermail/numpy-discussion/2016-March/075199.html</a><br></div><div><br></div><div>I completely agree with both of you that the current behavior for empty slices is very strange and should be changed to remove the special case. Nathaniel Smith voiced the same opinion on the GitHub issue [1].</div><div><br></div><div>I think the path forward here (as Nathaniel writes) is pretty clear:</div><div>1. Start issuing a FutureWarning about a future change.</div><div>2. Fix this in a release or two.</div><div><br></div><div>[1] <a href="https://github.com/numpy/numpy/issues/834">https://github.com/numpy/numpy/issues/834</a></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Jul 29, 2016 at 11:42 AM, Erik Brinkman <span dir="ltr"><<a href="mailto:erik.brinkman@umich.edu" target="_blank">erik.brinkman@umich.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi,<div><br></div><div>The behavior of a ufuncs reduceat on empty slices seems a little strange, and I wonder if there's a reason behind it / if there's a route to potentially changing it. First, I'll go into a little background.</div><div><br></div><div>I've been making a lot of use of ufuncs <font face="monospace">reduceat</font> functionality on staggered arrays. In general, I'll have "<font face="monospace">n</font>" arrays, each with size "<font face="monospace">s[n]</font>" and I'll store them in one array "<font face="monospace">x</font>", such that "<font face="monospace">s.sum() == x.size</font>". <font face="monospace">reduceat</font> is great because I use</div><div><font face="monospace"><br></font></div><div><font face="monospace">ufunc.reduceat(x, np.insert(s[:-1].cumsum(), 0, 0))</font></div><div><br></div><div>to get some summary information about each array. However, reduceat seems to behave strangely for empty slices. To make things concrete, let's assume:</div><div><br></div><div><font face="monospace">import numpy as np</font></div><div><font face="monospace">s = np.array([3, 0, 2])</font></div><div><span style="font-family:monospace;line-height:18px">x = np.arange(s.sum())</span><font face="monospace"><br></font></div><div><font face="monospace">inds = np.insert(s[:-1].cumsum(), 0, 0)</font></div><div><font face="monospace"># [0, 3, 3]</font></div><div><font face="monospace">np.add.reduceat(x, inds)</font></div><div><font face="monospace"># [3, 3, 7] not [3, 0, 7]</font></div><div><font face="monospace"># This is distinct from</font></div><div><font face="monospace">np.fromiter(map(np.add.reduce, np.array_split(x, inds[1:])), x.dtype, s.size - 1)<br></font></div><div><font face="monospace"># [3, 0, 7] what I wanted</font></div><div><br></div><div><span style="line-height:1.5">The </span><a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.reduceat.html" style="line-height:1.5" target="_blank">current documentation</a><span style="line-height:1.5"> on reduceat first states:</span></div><div><br></div><div>For i in <font face="monospace">range(len(indices))</font>, reduceat computes <font face="monospace">ufunc.reduce(a[indices[i]:indices[i+1]])</font><br></div><div><br></div><div>That would suggest the outcome, that I expected. However, in the examples section it goes into a bunch of examples which contradict that statement and instead suggest that the actual algorithm is more akin to:</div><div><br></div><div><font face="monospace">ufunc.reduce(a[indices[i]:indices[i+1]]) if indices[i+1] > indices[i] else <span style="line-height:1.5">indices[i]</span><br></font></div><div><br></div><div>Looking at the <a href="https://github.com/numpy/numpy/blob/2af06c804931aae4b30bb3349bc60271b0b65381/numpy/core/src/umath/ufunc_object.c#L3697" target="_blank">source</a>, it seems like it's copying indices[i], and then while there are more elements to process it keeps reducing, resulting in this unexpected behavior. It seems like the proper thing to do would be start with <font face="monospace">ufunc.identity</font>, and then reduce. This is slightly less performant than what's implemented, but more "correct." There could, of course, just be a switch to copy the identity only when the slice is empty.</div><div><br></div><div>Is there a reason it's implemented like this? Is it just for performance, or is this strange behavior <i>useful</i> somewhere? It seems like "fixing" this would be bad because you'll be changing somewhat documented functionality in a backwards incompatible way. What would the best approach to "fixing" this be? add another function "reduceat_"? add a flag to reduceat to do the proper thing for empty slices?</div><div><br></div><div>Finally, is there a good way to work around this? I think for now I'm just going to mask out the empty slices and use insert to add them back in, but if I'm missing an obvious solution, I'll look at that too. I need to mask them out because, <font face="monospace">np.add.reduceat(x, 5)</font> would ideally return 0, but instead it throws an error since 5 is out of range...</div><div><br></div><div>Thanks for indulging my curiosity,</div><div>Erik</div></div>
<br>_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="https://mail.scipy.org/mailman/listinfo/numpy-discussion" rel="noreferrer" target="_blank">https://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
<br></blockquote></div><br></div>