<div dir="ltr"><div><div><br></div>This behavior seems to depend on the order in which elements of the arrays are processes. That seems like a dangerous thing to rely on, the main reason I can thing of that someone would want to change the loop order is to implement parallel ufuncs.<br>
<br></div>Bago<br><br></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Fri, Oct 25, 2013 at 12:32 PM, Jaime Fernández del Río <span dir="ltr"><<a href="mailto:jaime.frio@gmail.com" target="_blank">jaime.frio@gmail.com</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">I recently came up with a way of vectorizing some recursive sequence calculations. While it works, I am afraid it is relying on implementation details potentially subject to change. The basic idea is illustrated by this function, calculating the first n items of the Fibonacci sequence:<div>
<br></div><div><p style="margin:0px">def fibonacci(n):</p>
<p style="margin:0px"> fib = np.empty(n, dtype=np.intp)</p>
<p style="margin:0px"> fib[:2] = 1</p>
<p style="margin:0px"> np.add(fib[:-2], fib[1:-1], out=fib[2:])</p>
<p style="margin:0px"> return fib</p><p style="margin:0px"><br></p><p style="margin:0px">>>> fibonacci(10)</p><p style="margin:0px">array([ 1, 1, 2, 3, 5, 8, 13, 21, 34, 55], dtype=int64)</p><p style="margin:0px">
<br></p><p style="margin:0px">I believe that the biggest issue that could break this is if the ufunc decided to buffer the arrays, as this is relying on the inputs and outputs of np.add sharing the same memory.</p><p style="margin:0px">
<br></p><p style="margin:0px">You can use the same idea to do more complicated things, for instance to calculate the items of the sequence:</p><p style="margin:0px"><br></p><p style="margin:0px">f[0] = a[0]</p><p style="margin:0px">
f[n] = a[n] + x * f[n-1]</p><p style="margin:0px"><br></p><p style="margin:0px">from numpy.lib.stride_tricks import as_strided</p><p style="margin:0px">from numpy.core.umath_tests import inner1d</p><p style="margin:0px">
<br>
</p><p style="margin:0px">def f(a, x):</p><p style="margin:0px"> out = np.array(a, copy=True, dtype=np.double)</p><p style="margin:0px"> n = len(a)</p><p style="margin:0px"> out_view = as_strided(out, shape=(n-1, 2), strides=out.strides*2)</p>
<p style="margin:0px"> inner1d(out_view, [x, 1], out=out[1:])</p><p style="margin:0px">
</p><p style="margin:0px"> return out</p><p style="margin:0px"><br></p><p style="margin:0px">>>> f(np.arange(10), 0.1)</p><p style="margin:0px">array([ 0. , 1. , 2.1 , 3.21 , 4.321 ,</p>
<p style="margin:0px"> 5.4321 , 6.54321 , 7.654321 , 8.7654321 , 9.87654321])</p><div><br></div><div>Again, I think buffering is the clearest danger of doing something like this, as the first input and output must share the same memory for this to work. That this is a real concern is easy to see: since `inner1d` only has loops registered for long ints and double floats:</div>
<div><br></div><div><div>>>> inner1d.types</div><div>['ll->l', 'dd->d']</div></div><div><br></div><div>the above function `f` doesn't work if the `out` array is created, e.g. as np.float32, since there will be buffering happening because of the type casting.</div>
<div><br></div><div>So I have two questions:</div><div><br></div><div>1. Is there some other reason, aside from buffering, that could go wrong, or change in a future release?</div><div>2. As far as buffering is concerned, I thought of calling `np.setbuffer(1)` before any of these functions, but it complains and requests that the value be a multiple of 16. Is there any other way to ensure that the data is fetched from an updated array in every internal iteration?</div>
<div><br></div><div>Thanks!</div><span class="HOEnZb"><font color="#888888"><div><br></div><div>Jaime</div><div><br></div>-- <br>(\__/)<br>( O.o)<br>( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.
</font></span></div></div>
<br>_______________________________________________<br>
NumPy-Discussion mailing list<br>
<a href="mailto:NumPy-Discussion@scipy.org">NumPy-Discussion@scipy.org</a><br>
<a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br>
<br></blockquote></div><br></div>