[Tutor] improving speed using and recalling C functions
Gabriele Brambilla
gb.gabrielebrambilla at gmail.com
Sat Apr 12 15:34:08 CEST 2014
Ok, i just run Peter's code and it seems really faster...I hope to don't
mistake this time!
Thanks
Gabriele
sent from Samsung Mobile
Il giorno 12/apr/2014 08:22, "Gabriele Brambilla" <
gb.gabrielebrambilla at gmail.com> ha scritto:
> 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.000numerictypes.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/691f0cd6/attachment-0001.html>
More information about the Tutor
mailing list