<html><head></head><body><div class="gmail_quote"><br>
"Jaime Fernández del Río" <jaime.frio@gmail.com> a écrit :<blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); 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><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.
</div></div>
<p style="margin-top: 2.5em; margin-bottom: 1em; border-bottom: 1px solid #000"></p><pre class="k9mail"><hr /><br />NumPy-Discussion mailing list<br />NumPy-Discussion@scipy.org<br /><a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br /></pre></blockquote></div><br clear="all">Hi,<br>
<br>
If you're recursive computation is always linear, like Fibonacci, your usecase falls in the scope of scipy.signal.lfilter which is compiled code based.<br>
<br>
For your last example, I think parameters would be : numerator b = [1] and denominator a=[1 -x]<br>
(can't check the code on my cell phone sorry)<br>
<br>
best,<br>
Pierre<br>
-- <br>
Pierre Haessig</body></html>