how to solve this exercise with python numpy vectorization?

Question

how to solve this exercise 4.5 on page 2 with python Numpy vectorization?

Link to download:

https://dl.dropbox.com/u/92795325/Python%20Scripting%20for%20Computational%20Scien%20-%20H.P.%20%20Langtangen.pdf

I tried this with Python loop, but I need Vectorization version.

from numpy import *
import time

def fun1(x):
       return 2*x+1

def integ(a,b,n):
       t0 = time.time()
       h = (b-a)/n
       a1 = (h/2)*fun1(a)
       b1 = (h/2)*fun1(b)
       c1 = 0
       for i in range(1,n,1):
              c1 = fun1((a+i*h))+c1
       t1 = time.time()
       return a1+b1+h*c1, t1-t0

Hi Ardeshir, do you know what it means to "vectorize" the problem? — askewchan, Apr 05 '13 at 14:50
I understand, I was asking if you did :) If you do understand, have you tried to implement it? — askewchan, Apr 05 '13 at 15:43
I have already understand the "vectorize" in recursions, but I can't find out what I have to do with this exercise. I am beginner in Python and I want to see some more examples of vectorizing. Here is the Recursion I have understood: https://dl.dropbox.com/u/92795325/Recursion%20%26%20Vectorization.py The problem is that each argument of the new array must be the sum of the two arguments beside that! (very simple and clear) — Ardeshir DVD, Apr 06 '13 at 07:42

score 0 · Accepted Answer · answered Apr 06 '13 at 22:20

To "vectorize" using numpy, all this means is that instead of doing an explicit loop like,

for i in range(1, n):
    c = c + f(i)

Then instead you should make i into a numpy array, and simply take its sum:

i = np.arange(1,n)
c = i.sum()

And numpy automatically does the vectorization for you. The reason this is faster is because numpy loops are done in a better optimized way than a plain python loop, for a variety of reasons. Generally speaking, the longer the loop/array, the better the advantage. Here is your trapezoidal integration implemented:

import numpy as np

def f1(x):
    return 2*x + 1

# Here's your original function modified just a little bit:
def integ(f,a,b,n):
    h = (b-a)/n
    a1 = (h/2)*f(a)
    b1 = (h/2)*f(b)
    c1 = 0
    for i in range(1,n,1):
        c1 = f((a+i*h))+c1
    return a1 + b1 + h*c1

# Here's the 'vectorized' function:
def vinteg(f, a, b, n):
    h = (b-a) / n
    ab = 0.5 * h * (f(a)+f(b)) #only divide h/2 once

    # use numpy to make `i` a 1d array:
    i = np.arange(1, n) 
    # then, passing a numpy array to `f()` means that `f` returns an array
    c = f(a + h*i) # now c is a numpy array

    return ab + c.sum() # ab + np.sum(c) is equivalent

Here I will import what I named tmp.py into an ipython session for easier timing than using time.time:

import trap
f = trap.f1
a = 0
b = 100
n = 1000

timeit trap.integ(f, a, b, n)
#1000 loops, best of 3: 378 us per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 51.6 us per loop

Wow, seven times faster.

See if it helps much for smaller n

n = 10

timeit trap.integ(f, a, b, n)
#100000 loops, best of 3: 6 us per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 43.4 us per loop

Nope, much slower for small loops! What about very large n?

n = 10000

timeit trap.integ(f, a, b, n)
#100 loops, best of 3: 3.69 ms per loop

timeit trap.vinteg(f, a, b, n)
#10000 loops, best of 3: 111 us per loop

Thirty times faster!

Thanks for your help. I see... But I have some more question if you don't mind. First, is "vectorized version" always more efficient than loops--i.e. should I avoid using loops whenever possible? Second, what is the difference between loops and vectors that makes such big difference? Why there is not such thing in (for example) FORTRAN and C programming? — Ardeshir DVD, Apr 08 '13 at 16:47
@ArdeshirDVD First: Usually faster, but not always. Also it's also more readable than loops, which can more important. Second: The difference isn't that they're loops, but the way the data is stored and iterated over. Numpy builtin functions are mainly written and compiled in C which is faster. They don't exist in Fortran or C because F. and C don't need them: those are compiled languages, and it's the fact that python is an "interpreted" language that slows down its loops. — askewchan, Apr 08 '13 at 19:46
@ArdeshirDVD Since python is interpreted and not compiled, and because its variables are all dynamically typed, each iteration of the loop must check the types of the variable, which is wasteful for a loop over 10000 floats of identical type. See questions like http://stackoverflow.com/a/8097669/1730674 and http://stackoverflow.com/q/3033329/1730674 for a little discussion on why python is so slow running for some things — askewchan, Apr 08 '13 at 19:55

how to solve this exercise with python numpy vectorization?

1 Answers1