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

Reply via email to