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