[Tutor] improving speed using and recalling C functions
Gabriele Brambilla
gb.gabrielebrambilla at gmail.com
Sat Apr 12 14:22:34 CEST 2014
Ok guys,
I'm not expert about profile but help me to look at it.
this one is for 715853 elements (to multiply by 5, and for each of this N*5
there is a loop of 200 times)
Sat Apr 12 04:58:50 2014 restats
9636507991 function calls in 66809.764 seconds
Ordered by: internal time
List reduced from 47 to 20 due to restriction <20>
ncalls tottime percall cumtime percall
filename:lineno(function)
1 13548.507 13548.507 66809.692 66809.692
skymapsI.py:44(mymain)
125800544 13539.337 0.000 15998.925 0.000
interpolate.py:394(_call_linear)
880603808 5353.382 0.000 5353.382 0.000
{numpy.core.multiarray.array}
715853000 4998.740 0.000 52861.634 0.000
instruments.py:10(kappa)
251601088 4550.940 0.000 4550.940 0.000 {method 'reduce'
of 'numpy.ufunc' objects}
125800544 4312.078 0.000 10163.614 0.000
interpolate.py:454(_check_bounds)
125800544 2944.126 0.000 14182.917 0.000
interpolate.py:330(__init__)
125800544 2846.577 0.000 29484.248 0.000
interpolate.py:443(_evaluate)
125800544 1665.852 0.000 6000.603 0.000 polyint.py:82(_set_yi)
125800544 1039.455 0.000 1039.455 0.000 {method 'clip' of
'numpy.ndarray' objects}
251601088 944.848 0.000 944.848 0.000 {method 'reshape' of
'numpy.ndarray' objects}
251601088 922.928 0.000 1651.218 0.000
numerictypes.py:735(issubdtype)
503202176 897.044 0.000 3434.768 0.000
numeric.py:392(asarray)
125800544 816.401 0.000 32242.481 0.000
polyint.py:37(__call__)
251601088 787.593 0.000 5338.533 0.000 _methods.py:31(_any)
125800544 689.779 0.000 1989.101 0.000
polyint.py:74(_reshape_yi)
125800544 638.946 0.000 638.946 0.000 {method
'searchsorted' of 'numpy.ndarray' objects}
125800544 606.778 0.000 2257.996 0.000
polyint.py:102(_set_dtype)
125800544 598.000 0.000 6598.602 0.000
polyint.py:30(__init__)
629002720 549.358 0.000 549.358 0.000 {issubclass}
looking at tottime it seems that skymaps mymain() and interpolate take the
same big amount of time...right?
So it's true that I have to slow down mymain() but interpolate is a problem
too!
do you agree with me?
Now I will read Peter Otten's code and run the new simulation with it
thanks
Gabriele
2014-04-12 6:21 GMT-04:00 Peter Otten <__peter__ at web.de>:
> 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)
> >
>
>
> _______________________________________________
> Tutor maillist - Tutor at python.org
> To unsubscribe or change subscription options:
> https://mail.python.org/mailman/listinfo/tutor
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/tutor/attachments/20140412/e292037b/attachment-0001.html>
More information about the Tutor
mailing list