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 <[email protected]>:
> 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 - [email protected]
> To unsubscribe or change subscription options:
> https://mail.python.org/mailman/listinfo/tutor
>
_______________________________________________
Tutor maillist - [email protected]
To unsubscribe or change subscription options:
https://mail.python.org/mailman/listinfo/tutor
