"Jaime Fernández del Río" <jaime.f...@gmail.com> a écrit :
>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:
>
>def fibonacci(n):
>
> fib = np.empty(n, dtype=np.intp)
>
> fib[:2] = 1
>
> np.add(fib[:-2], fib[1:-1], out=fib[2:])
>
> return fib
>
>
>>>> fibonacci(10)
>
>array([ 1, 1, 2, 3, 5, 8, 13, 21, 34, 55], dtype=int64)
>
>
>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.
>
>
>You can use the same idea to do more complicated things, for instance
>to
>calculate the items of the sequence:
>
>
>f[0] = a[0]
>
>f[n] = a[n] + x * f[n-1]
>
>
>from numpy.lib.stride_tricks import as_strided
>
>from numpy.core.umath_tests import inner1d
>
>
>def f(a, x):
>
> out = np.array(a, copy=True, dtype=np.double)
>
> n = len(a)
>
> out_view = as_strided(out, shape=(n-1, 2), strides=out.strides*2)
>
> inner1d(out_view, [x, 1], out=out[1:])
>
> return out
>
>
>>>> f(np.arange(10), 0.1)
>
>array([ 0. , 1. , 2.1 , 3.21 , 4.321 ,
>
> 5.4321 , 6.54321 , 7.654321 , 8.7654321 , 9.87654321])
>
>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:
>
>>>> inner1d.types
>['ll->l', 'dd->d']
>
>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.
>
>So I have two questions:
>
>1. Is there some other reason, aside from buffering, that could go
>wrong,
>or change in a future release?
>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?
>
>Thanks!
>
>Jaime
>
>--
>(\__/)
>( O.o)
>( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus
>planes
>de dominación mundial.
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>NumPy-Discussion mailing list
>NumPy-Discussion@scipy.org
>http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi,
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.
For your last example, I think parameters would be : numerator b = [1] and
denominator a=[1 -x]
(can't check the code on my cell phone sorry)
best,
Pierre
--
Pierre Haessig
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
