# [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
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 = []
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

> 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)
>