[Tutor] improving speed using and recalling C functions
Peter Otten
__peter__ at web.de
Sat Apr 12 12:21:22 CEST 2014
Gabriele Brambilla wrote:
> Ok guys, when I wrote that email I was excited for the apparent speed
> increasing (it was jumping the bottleneck for loop for the reason peter
> otten outlined).
> Now, instead the changes, the speed is not improved (the code still
> running from this morning and it's at one forth of the dataset).
>
> What can I do to speed it up?
Not as easy as I had hoped and certainly not as pretty, here's my
modification of the code you sent me. What makes it messy is that
I had to inline your kappa() function; my first attempt with
numpy.vectorize() didn't help much. There is still stuff in the
'for gammar...' loop that doesn't belong there, but I decided it
was time for me to stop ;)
Note that it may still be worthwhile to consult a numpy expert
(which I'm not!).
from scipy import stats
import matplotlib.pyplot as plt
from scipy import optimize
from matplotlib import colors, ticker, cm
import numpy as np
phamin = 0
phamax = 2*pi
obamin = 0
obamax = pi
npha = 100
nobs = 181
stepPHA = (phamax-phamin)/npha
stepOB = (obamax-obamin)/nobs
freq = 10
c = 2.9979*(10**(10))
e = 4.8032*(10**(-10))
hcut = 1.0546*(10**(-27))
eVtoErg = 1.6022*(10**(-12))
from math import *
import numpy as np
from scipy.interpolate import interp1d
kaparg = [
-3.0, -2.0, -1.52287875, -1.22184875, -1.0, -0.69897,
-0.52287875, -0.39794001, -0.30103, -0.22184875,
-0.15490196, 0.0, 0.30103, 0.60205999, 0.69897,
0.77815125, 0.90308999, 1.0]
kapval = [
-0.6716204 , -0.35163999, -0.21183163, -0.13489603,
-0.0872467 , -0.04431225, -0.03432803, -0.04335142,
-0.05998184, -0.08039898, -0.10347378, -0.18641901,
-0.52287875, -1.27572413, -1.66958623, -2.07314329,
-2.88941029, -3.7212464 ]
my_inter = interp1d(kaparg, kapval)
def LEstep(n):
Emin = 10**6
Emax = 5*(10**10)
Lemin = log10(Emin)
Lemax = log10(Emax)
stepE = (Lemax-Lemin)/n
return stepE, n, Lemin, Lemax
def mymain(stepENE, nex, Lemin, Lemax, freq):
eel = np.array(list(range(nex)))
eels = np.logspace(Lemin, Lemax, num=nex, endpoint=False)
rlc = c/(2*pi*freq)
sigmas = [1, 3, 5, 10, 30]
MYMAPS = [
np.zeros([npha, nobs, nex], dtype=float) for _ in sigmas]
alpha = '60_'
ALPHA = (1.732050808/c)*(e**2)
for count, my_line in enumerate(open('datasm0_60_5s.dat')):
myinternet = []
print('reading the line', count, '/599378')
my_parts = np.array(my_line.split(), dtype=float)
phase = my_parts[4]
zobs = my_parts[5]
rho = my_parts[6]
gmils = my_parts[7:12]
i = int((phase-phamin)/stepPHA)
j = int((zobs-obamin)/stepOB)
for gammar, MYMAP in zip(gmils, MYMAPS):
omC = (1.5)*(gammar**3)*c/(rho*rlc)
gig = omC*hcut/eVtoErg
omega = (10**(eel*stepENE+Lemin))*eVtoErg/hcut
x = omega/omC
kap = np.empty(x.shape)
sel = x >= 10.0
zsel = x[sel]
kap[sel] = 1.2533 * np.sqrt(zsel)*np.exp(-zsel)
sel = x < 0.001
zsel = x[sel]
kap[sel] = (2.1495 * np.exp(0.333333333 * np.log(zsel))
- 1.8138 * zsel)
sel = ~ ((x >= 10.0) | (x < 0.001))
zsel = x[sel]
result = my_inter(np.log10(zsel))
kap[sel] = 10**result
Iom = ALPHA*gammar*kap
P = Iom*(c/(rho*rlc))/(2*pi)
phps = P/(hcut*omega)
www = phps/(stepPHA*sin(zobs)*stepOB)
MYMAP[i,j] += www
for sigma, MYMAP in zip(sigmas, MYMAPS):
print(sigma)
filename = "_".join(str(p) for p in
["skymap", alpha, sigma, npha, phamin, phamax, nobs,
obamin, obamax, nex, Lemin, Lemax, '.dat']
)
x, y, z = MYMAP.shape
with open(filename, 'ab') as MYfile:
np.savetxt(
MYfile,
MYMAP.reshape(x*y, z, order="F").T,
delimiter=",", fmt="%s", newline=",\n")
if __name__ == "__main__":
if len(sys.argv)<=1:
stepENE, nex, Lemin, Lemax = LEstep(200)
elif len(sys.argv)<=2:
stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))
else:
stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))
freq=float(sys.argv[2])
mymain(stepENE, nex, Lemin, Lemax, freq)
For reference here is the original (with the loop over gmlis
instead of gmils):
> import sys
>
> from math import *
> from scipy import ndimage
> from scipy import stats
> import matplotlib.pyplot as plt
> from scipy import optimize
> from matplotlib import colors, ticker, cm
> import numpy as np
> import cProfile
> import pstats
>
> phamin=0
> phamax=2*pi
> obamin=0
> obamax=pi
> npha=100
> nobs=181
> stepPHA=(phamax-phamin)/npha
> stepOB=(obamax-obamin)/nobs
> freq=10
> c=2.9979*(10**(10))
> e=4.8032*(10**(-10))
> hcut=1.0546*(10**(-27))
> eVtoErg=1.6022*(10**(-12))
>
>
> from math import *
> import numpy as np
> from scipy.interpolate import interp1d
>
>
> def kappa(z):
> N=18
> kaparg = [-3.0, -2.0, -1.52287875, -1.22184875, -1.0, -0.69897, -0.52287875, -0.39794001, -0.30103, -0.22184875, -0.15490196, 0.0, 0.30103, 0.60205999, 0.69897, 0.77815125, 0.90308999, 1.0]
> kapval = [-0.6716204 , -0.35163999, -0.21183163, -0.13489603, -0.0872467 , -0.04431225, -0.03432803, -0.04335142, -0.05998184, -0.08039898, -0.10347378, -0.18641901, -0.52287875, -1.27572413,
-1.66958623, -2.07314329, -2.88941029, -3.7212464 ]
> zlog=log10(z)
> if z < 0.001:
> k = 2.1495 * exp (0.333333333 * log (z)) - 1.8138 * z
> return (k)
> elif z >= 10.0:
> k = 1.2533 * sqrt (z) * exp (-z)
> return (k)
> else:
> my_inter = interp1d(kaparg, kapval)
> my_z = np.array([zlog])
> result = my_inter(my_z)
> valuelog = result[0]
> k=10**valuelog
> return(k)
>
>
>
>
> def LEstep(n):
> Emin=10**6
> Emax=5*(10**10)
> Lemin=log10(Emin)
> Lemax=log10(Emax)
> stepE=(Lemax-Lemin)/n
> return (stepE, n, Lemin, Lemax)
>
>
> def mymain(stepENE, nex, Lemin, Lemax, freq):
>
>
> eel = list(range(nex))
> eels = np.logspace(Lemin, Lemax, num=nex, endpoint=False)
>
> indpha = list(range(npha))
> indobs = list(range(nobs))
> rlc = c/(2*pi*freq)
>
> #creating an empty 3D vector
> MYMAPS = [np.zeros([npha, nobs, nex], dtype=float), np.zeros([npha, nobs, nex], dtype=float), np.zeros([npha, nobs, nex], dtype=float), np.zeros([npha, nobs, nex], dtype=float), np.zeros([npha, nobs,
nex], dtype=float)]
>
>
> count=0
>
>
> alpha = '60_'
>
> for my_line in open('datasm0_60_5s.dat'):
> myinternet = []
> gmlis = []
> print('reading the line', count, '/599378')
> my_parts = [float(i) for i in my_line.split()]
> phase = my_parts[4]
> zobs = my_parts[5]
> rho = my_parts[6]
>
> gmils=[my_parts[7], my_parts[8], my_parts[9], my_parts[10], my_parts[11]]
>
> i = int((phase-phamin)/stepPHA)
> j = int((zobs-obamin)/stepOB)
>
> for gammar, MYMAP in zip(gmils, MYMAPS):
>
> omC = (1.5)*(gammar**3)*c/(rho*rlc)
> gig = omC*hcut/eVtoErg
>
> for w in eel:
> omega = (10**(w*stepENE+Lemin))*eVtoErg/hcut
> x = omega/omC
> kap = kappa(x)
> Iom = (1.732050808/c)*(e**2)*gammar*kap
> P = Iom*(c/(rho*rlc))/(2*pi)
> phps = P/(hcut*omega)
> www = phps/(stepPHA*sin(zobs)*stepOB)
> MYMAP[i,j,w] += www
>
> count = count + 1
>
>
>
> sigmas = [1, 3, 5, 10, 30]
>
> multis = zip(sigmas, MYMAPS)
>
> for sigma, MYMAP in multis:
>
> print(sigma)
> filename='skymap_'+alpha+'_'+str(sigma)+'_'+str(npha)+'_'+str(phamin)+'_'+str(phamax)+'_'+str(nobs)+'_'+str(obamin)+'_'+str(obamax)+'_'+str(nex)+'_'+str(Lemin)+'_'+str(Lemax)+'_.dat'
>
> MYfile = open(filename, 'a')
> for k in eel:
> for j in indobs:
> for i in indpha:
> A=MYMAP[i, j, k]
> stringa = str(A) + ','
> MYfile.write(stringa)
> accapo = '\n'
> MYfile.write(accapo)
>
> MYfile.close()
>
>
> if __name__ == "__main__":
> if len(sys.argv)<=1:
> stepENE, nex, Lemin, Lemax = LEstep(200)
> elif len(sys.argv)<=2:
> stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))
> else:
> stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))
> freq=float(sys.argv[2])
>
>
> #mymain(stepENE, nex, Lemin, Lemax, freq)
>
> #print('profile')
> cProfile.run('mymain(stepENE, nex, Lemin, Lemax, freq)', 'restats', 'time')
>
> p = pstats.Stats('restats')
> p.strip_dirs().sort_stats('name')
> p.sort_stats('time').print_stats(20)
>
More information about the Tutor
mailing list