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