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.
Bago On Fri, Oct 25, 2013 at 12:32 PM, Jaime Fernández del Río < jaime.f...@gmail.com> wrote: > 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 > >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion