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] whereas the code: for i in range(0,n-1): y[i+1] = y[i] + u[i+1] 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