Question about scientific calculations in Python

[Jason Orendorff]

I benchmarked it; results:

  plain Python: 8.557 sec
  optimized Python: 2.325 sec
  Numeric Python: 0.753 sec

This is for 1000 numbers in both s_vector and r_vector.

## Jason Orendorff    http://www.jorendorff.com/



from Numeric import *
from random import random
from time import clock

def f1(s_vector, f_s_vector, r_vector, atoms):
    n = 0
    i_vector = []
    for s in s_vector:
        sum = 0
        for r in r_vector:
            x = 2 * pi * r * s
            sum = sum + (2 * (sin(x))/x)
        intensity = atoms * pow(f_s_vector[n], 2) * (1 + sum / atoms)
        i_vector.append(intensity)
        n = n + 1
    return i_vector

def f2(s_vector, f_s_vector, r_vector, atoms):
    # Precompute the 2*pi*r, which is a constant inside of the s loop
    precomputed_r1_vector = [2 * pi * r for r in r_vector]
    local_sin = math.sin  # local variables have fast lookup
    i_vector = []
    i_vector_append = i_vector.append

    # using the zip gets rid of the extra n+1 and [n], at the expense
    # of making a new list.  I suspect it's faster, but I'm not sure
    for s, f_s_vector in zip(s_vector, f_s_vector):
        sum = 0.0
        for scaled_r in precomputed_r1_vector:
            x = scaled_r * s
            sum = sum + 2 * (local_sin(x)/x)

        # Did some algebra here to merge the first and last terms,
        # which gets rid of the division.  Also replaced the function call
        # to 'pow' with the equivalent use of **2
        i_vector_append((f_s_vector**2) * (atoms + sum))

    return i_vector

def f3(S, FS, R, atoms):
    Rcol = R[:, NewAxis]
    x = 2 * pi * Rcol * S
    y = 2 * sin(x) / x
    sums = add.reduce(y)
    I = atoms * FS**2 * (1 + sums/atoms)
    return I

def timefn(msg, fn, *args):
    t0 = clock()
    result = fn(*args)
    t1 = clock()
    print msg + ": %.3f sec" % (t1-t0)

def measure():
    S = array([random()-.5 for i in range(1000)])
    R = array([random()-.5 for i in range(1000)])
    FS = array([random()*0.5 for i in range(1000)])
    atoms = 1000

    timefn("plain Python",     f1, S, FS, R, atoms)
    timefn("optimized Python", f2, S, FS, R, atoms)
    timefn("Numeric Python",   f3, S, FS, R, atoms)

if __name__ == '__main__':
    measure()