Hi, On Sat, Apr 7, 2012 at 1:27 AM, Sameer Grover <sameer.grove...@gmail.com>wrote:
> On Saturday 07 April 2012 02:51 AM, Francesco Barale wrote: > > Hello Sameer, > > > > Thank you very much for your reply. My goal was to try to speed up the > loop > > describing the accumulator. In the (excellent) book I was mentioning in > my > > initial post I could find one example that seemed to match what I was > trying > > to do. Basically, it is said that a loop of the following kind: > > > > n = size(u)-1 > > for i in xrange(1,n,1): > > u_new[i] = u[i-1] + u[i] + u[i+1] > > > > should be equivalent to: > > > > u[1:n] = u[0:n-1] + u[1:n] + u[i+1] > This example is correct. > > What I was trying to point out was that the single expression y[1:n] = > y[0:n-1] + u[1:n] will iterate over the array but will not accumulate. > It will add y[0:n-1] to u[1:n] and assign the result to y[1:n]. > > For example, > y = [1,2,3,4] > u = [5,6,7,8] > Then y[0:n-1] = [1,2,3] and u[1:n]=[6,7,8] > > The statement y[1:n] = y[0:n-1] + u[1:n] implies > y[1:n] = [6+1,7+2,8+3] = [7,9,11] > yielding y = [1,7,9,11] > FWIIFO, if assignment in loop like this ever makes any sense (which I doubt), > > whereas the code: > > for i in range(0,n-1): > y[i+1] = y[i] + u[i+1] > then it can be captured in a function, like In []: def s0(y, u): ..: yn= y.copy() ..: for k in xrange(y.size- 1): ..: yn[k+ 1]= yn[k]+ u[k+ 1] ..: return yn ..: and now this function can be easily transformed to utilize cumsum, like In []: def s1(y, u): ..: un= u.copy() ..: un[0]= y[0] ..: return cumsum(un) ..: thus In []: y, u= rand(1e5), rand(1e5) In []: allclose(s0(y, u), s1(y, u)) Out[]: True and definitely this transformation will outperform a plain python loop In []: timeit s0(y, u) 10 loops, best of 3: 122 ms per loop In []: timeit s1(y, u) 100 loops, best of 3: 2.16 ms per loop In []: 122/ 2.16 Out[]: 56.48148148148148 My 2 cents, -eat > > will accumulate and give y = [1,7,14,22] > > Sameer > > Am I missing something? > > > > Regards, > > Francesco > > > > > > Sameer Grover wrote: > >> On Saturday 07 April 2012 12:14 AM, francesco82 wrote: > >>> Hello everyone, > >>> > >>> After reading the very good post > >>> > http://technicaldiscovery.blogspot.com/2011/06/speeding-up-python-numpy-cython-and.html > >>> and the book by H. P. Langtangen 'Python scripting for computational > >>> science' I was trying to speed up the execution of a loop on numpy > arrays > >>> being used to describe a simple difference equation. > >>> > >>> The actual code I am working on involves some more complicated > equations, > >>> but I am having the same exact behavior as described below. To test the > >>> improvement in speed I wrote the following in vect_try.py: > >>> > >>> #!/usr/bin/python > >>> import numpy as np > >>> import matplotlib.pyplot as plt > >>> > >>> dt = 0.02 #time step > >>> time = np.arange(0,2,dt) #time array > >>> u = np.sin(2*np.pi*time) #input signal array > >>> > >>> def vect_int(u,y): #vectorized function > >>> n = u.size > >>> y[1:n] = y[0:n-1] + u[1:n] > >>> return y > >>> > >>> def sc_int(u,y): #scalar function > >>> y = y + u > >>> return y > >>> > >>> def calc_vect(u, func=vect_int): > >>> out = np.zeros(u.size) > >>> for i in xrange(u.size): > >>> out = func(u,out) > >>> return out > >>> > >>> def calc_sc(u, func=sc_int): > >>> out = np.zeros(u.size) > >>> for i in xrange(u.size-1): > >>> out[i+1] = sc_int(u[i+1],out[i]) > >>> return out > >>> > >>> To verify the execution time I've used the timeit function in Ipython: > >>> > >>> import vect_try as vt > >>> timeit vt.calc_vect(vt.u) --> 1000 loops, best of 3: 494 us per loop > >>> timeit vt.calc_sc(vt.u) -->10000 loops, best of 3: 92.8 us per loop > >>> > >>> As you can see the scalar implementation looping one item at the time > >>> (calc_sc) is 494/92.8~5.3 times faster than the vectorized one > >>> (calc_vect). > >>> > >>> My problem consists in the fact that I need to iterate the execution of > >>> calc_vect in order for it to operate on all the elements of the input > >>> array. > >>> If I execute calc_vect only once, it will only operate on the first > slice > >>> of > >>> the vectors leaving the remaining untouched. My understanding was that > >>> the > >>> vector expression y[1:n] = y[0:n-1] + u[1:n] was supposed to iterate > over > >>> all the array, but that's not happening for me. Can anyone tell me > what I > >>> am > >>> doing wrong? > >>> > >>> Thanks! > >>> Francesco > >>> > >> 1. By vectorizing, we mean replacing a loop with a single expression. In > >> your program, both the scalar and vector implementations (calc_vect and > >> calc_sc) have a loop each. This isn't going to make anything faster. The > >> vectorized implementation is just a convoluted way of achieving the same > >> result and is slower. > >> > >> 2. The expression y[1:n] = y[0:n-1] + u[1:n] is /not/ equivalent to the > >> following loop > >> > >> for i in range(0,n-1): > >> y[i+1] = y[i] + u[i+1] > >> > >> It is equivalent to something like > >> > >> z = np.zeros(n-1) > >> for i in range(0,n-1): > >> z[i] = y[i] + u[i+1] > >> y[1:n] = z > >> > >> i.e., the RHS is computed in totality and then assigned to the LHS. This > >> is how array operations work even in other languages such as Matlab. > >> > >> 3. I personally don't think there is a simple/obvious way to vectorize > >> what you're trying to achieve. > >> > >> Sameer > >> > >> _______________________________________________ > >> 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 >
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
