<div dir="ltr"><br><div class="gmail_extra"><br><div class="gmail_quote">On Wed, Nov 25, 2015 at 12:21 AM, Sebastian Berg <span dir="ltr"><<a href="mailto:sebastian@sipsolutions.net" target="_blank">sebastian@sipsolutions.net</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="HOEnZb"><div class="h5">On Di, 2015-11-24 at 16:49 -0800, Eli Bendersky wrote:<br>
><br>
><br>
> On Mon, Nov 23, 2015 at 2:09 PM, Sebastian Berg<br>
> <<a href="mailto:sebastian@sipsolutions.net">sebastian@sipsolutions.net</a>> wrote:<br>
>         On Mo, 2015-11-23 at 13:31 -0800, Eli Bendersky wrote:<br>
>         > Hello,<br>
>         ><br>
>         ><br>
>         > I'm trying to understand the buffering done by the Numpy<br>
>         iterator<br>
>         > interface (the new post 1.6-one) when running ufuncs on<br>
>         arrays that<br>
>         > require broadcasting. Consider this simple case:<br>
>         ><br>
>         > In [35]: m = np.arange(16).reshape(4,4)<br>
>         > In [37]: n = np.arange(4)<br>
>         > In [39]: m + n<br>
>         > Out[39]:<br>
>         > array([[ 0,  2,  4,  6],<br>
>         >        [ 4,  6,  8, 10],<br>
>         >        [ 8, 10, 12, 14],<br>
>         >        [12, 14, 16, 18]])<br>
>         ><br>
>         ><br>
><br>
><br>
>         On first sight this seems true. However, there is one other<br>
>         point to<br>
>         consider here. The inner ufunc loop can only handle a single<br>
>         stride. The<br>
>         contiguous array `n` has to be iterated as if it had the<br>
>         strides<br>
>         `(0, 8)`, which is not the strides of the contiguous array `m`<br>
>         which can<br>
>         be "unrolled" to 1-D. Those effective strides are thus not<br>
>         contiguous<br>
>         for the inner ufunc loop and cannot be unrolled into a single<br>
>         ufunc<br>
>         call!<br>
><br>
>         The optimization (which might kick in a bit more broadly<br>
>         maybe), is thus<br>
>         that the number of inner loop calls is minimized, whether that<br>
>         is worth<br>
>         it, I am not sure, it may well be that there is some worthy<br>
>         optimization<br>
>         possible here.<br>
>         Note however, that this does not occur for large inner loop<br>
>         sizes<br>
>         (though I think you can find some "bad" sizes):<br>
><br>
>         ```<br>
>         In [18]: n = np.arange(40000)<br>
><br>
>         In [19]: m = np.arange(160000).reshape(4,40000)<br>
><br>
>         In [20]: o = m + n<br>
>         Iterator: Checking casting for operand 0<br>
>         op: dtype('int64'), iter: dtype('int64')<br>
>         Iterator: Checking casting for operand 1<br>
>         op: dtype('int64'), iter: dtype('int64')<br>
>         Iterator: Checking casting for operand 2<br>
>         op: <null>, iter: dtype('int64')<br>
>         Iterator: Setting allocated stride 1 for iterator dimension 0<br>
>         to 8<br>
>         Iterator: Setting allocated stride 0 for iterator dimension 1<br>
>         to 320000<br>
>         Iterator: Copying inputs to buffers<br>
>         Iterator: Expanding inner loop size from 8192 to 40000 since<br>
>         buffering<br>
>         wasn't needed<br>
>         Any buffering needed: 0<br>
>         Iterator: Finished copying inputs to buffers (buffered size is<br>
>         40000)<br>
>         ```<br>
><br>
><br>
> The heuristic in the code says that if we can use a single stride and<br>
> that's larger than the buffer size (which I assume is the default<br>
> buffer size, and can change) then it's "is_onestride" and no buffering<br>
> is done.<br>
><br>
><br>
> So this led me to explore around this threshold (8192 items by default<br>
> on my machine), and indeed we can notice funny behavior:<br>
><br>
> In [51]: %%timeit n = arange(8192); m =<br>
> np.arange(8192*40).reshape(40,8192)<br>
> o = m + n<br>
>    ....:<br>
> 1000 loops, best of 3: 274 µs per loop<br>
><br>
> In [52]: %%timeit n = arange(8292); m =<br>
> np.arange(8292*40).reshape(40,8292)<br>
> o = m + n<br>
>    ....:<br>
> 1000 loops, best of 3: 229 µs per loop<br>
><br>
><br>
> So, given this, it's not very clear why the "optimization" kicks in.<br>
> Buffering for small sizes seems like a mistake.<br>
><br>
<br>
</div></div>I am pretty sure it is not generally a mistake. Consider the case of an<br>
10000x3 array (note that shrinking the buffer can have great advantage<br>
though, I am not sure if this is done usually).<br>
If you have (10000, 3) + (3,) arrays, then the ufunc outer loop would<br>
have 10000x overhead. Doing the buffering (which I believe has some<br>
extra code to be faster), will lower this to a few ufunc inner loop<br>
calls.<br>
I have not timed it, but I would be a bit surprised if it was not faster<br>
in this case at least. Even calling a C function (and looping) for an<br>
inner loop of 3 elements, should be quite a bit of overhead, and my<br>
guess is more overhead than the buffering.<span class="HOEnZb"><font color="#888888"><br></font></span></blockquote><div><br></div><div>Yes, that's a good point for arrays shaped like this. I guess all this leaves us with is a realization that the heuristic *could* be tuned somewhat for arrays where the inner dimension is large - as in the case I demonstrated above it's nonsensical to have a computation be 20% faster when the array size increases over an arbitrary threshold.<br><br></div><div>I'll see if I can find some time to dig more into this and figure out where the knobs to tweak the heuristic are.<br><br></div><div>Thanks for the enlightening discussion, Sebastian<br></div><div><br></div><div>Eli<br></div><div><br><br><br> </div><div><br><br><br><br> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="HOEnZb"><font color="#888888">
- Sebastian<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
><br>
> Eli<br>
><br>
><br>
><br>
><br>
>         Anyway, feel free to have a look ;). The code is not the most<br>
>         read one<br>
>         in NumPy, and it would not surprise me a lot if you can find<br>
>         something<br>
>         to tweak.<br>
><br>
>         - Sebastian<br>
><br>
><br>
>         ><br>
>         > If I instrument Numpy (setting NPY_IT_DBG_TRACING and such),<br>
>         I see<br>
>         > that when the add() ufunc is called, 'n' is copied into a<br>
>         temporary<br>
>         > buffer by the iterator. The ufunc then gets the buffer as<br>
>         its data.<br>
>         ><br>
>         ><br>
>         > My question is: why is this buffering needed? It seems<br>
>         wasteful, since<br>
>         > no casting is required here, no special alignment problems<br>
>         and also<br>
>         > 'n' is contiguously laid out in memory. It seems that it<br>
>         would be more<br>
>         > efficient to just use 'n' in the ufunc instead of passing in<br>
>         the<br>
>         > buffer. What am I missing?<br>
>         ><br>
>         ><br>
>         > Thanks in advance,<br>
>         ><br>
>         > Eli<br>
>         ><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>
><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>
><br>
><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>
</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></div>