ok Peter Otten code works (very fast),
and this is the profile
Sat Apr 12 11:15:39 2014 restats
92834776 function calls in 6218.782 seconds
Ordered by: internal time
List reduced from 41 to 20 due to restriction <20>
ncalls tottime percall cumtime percall filename:lineno(function)
1 5301.641 5301.641 6218.763 6218.763 skymapsIO.py:50(mymain)
3489985 380.469 0.000 452.478 0.000
interpolate.py:394(_call_linear)
3489985 98.186 0.000 227.229 0.000
interpolate.py:454(_check_bounds)
6979970 96.567 0.000 96.567 0.000 {method 'reduce' of
'numpy.ufunc'objects}
3489985 44.853 0.000 738.135 0.000 interpolate.py:443(_evaluate)
7677978 41.010 0.000 41.010 0.000 {numpy.core.multiarray.array}
5 40.430 8.086 40.621 8.124 npyio.py:882(savetxt)
3489985 26.952 0.000 26.952 0.000 {method 'clip' of
'numpy.ndarray'objects}
3489985 24.749 0.000 24.749 0.000 {method 'searchsorted' of
'numpy.ndarray' objects}
3489985 22.457 0.000 828.238 0.000 polyint.py:37(__call__)
6979970 19.720 0.000 116.287 0.000 _methods.py:31(_any)
6979980 15.330 0.000 35.092 0.000 numeric.py:392(asarray)
3489985 14.847 0.000 45.039 0.000 polyint.py:57(_prepare_x)
3489990 12.904 0.000 12.904 0.000 {method 'reshape' of
'numpy.ndarray' objects}
6979970 12.757 0.000 129.044 0.000 {method 'any' of
'numpy.ndarray' objects}
3489985 11.624 0.000 11.624 0.000 {method 'astype' of
'numpy.ndarray' objects}
3489985 10.077 0.000 10.077 0.000 {numpy.core.multiarray.empty}
3489985 9.945 0.000 22.607 0.000 polyint.py:63(_finish_y)
3489985 7.051 0.000 7.051 0.000 {method 'ravel' of
'numpy.ndarray' objects}
697998 6.746 0.000 6.746 0.000 {zip}
So I think that in this way it's ok.
Thank you all very much,
Gabriele
p.s: I didn't know this way to write: is there a tutorial for this kind of
operations?
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
2014-04-12 9:34 GMT-04:00 Gabriele Brambilla <gb.gabrielebrambi...@gmail.com
>:
> 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
