Question about scientific calculations in Python
Jason Orendorff
jason at jorendorff.com
Wed Mar 13 00:36:05 EST 2002
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()
More information about the Python-list
mailing list