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.gabrielebrambi...@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 <__pete...@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@python.org >> To unsubscribe or change subscription options: >> https://mail.python.org/mailman/listinfo/tutor >> > >
_______________________________________________ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: https://mail.python.org/mailman/listinfo/tutor