<div dir="ltr">ok Peter Otten code works (very fast),<div>and this is the profile</div><div><div>Sat Apr 12 11:15:39 2014    restats</div><div><br></div><div>         92834776 function calls in 6218.782 seconds</div><div><br>
</div><div>   Ordered by: internal time</div><div>   List reduced from 41 to 20 due to restriction <20></div><div><br></div><div>   ncalls  tottime  percall  cumtime  percall filename:lineno(function)</div><div>        1 5301.641 5301.641 6218.763 6218.763 skymapsIO.py:50(mymain)</div>
<div>  3489985  380.469    0.000  452.478    0.000 interpolate.py:394(_call_linear)</div><div>  3489985   98.186    0.000  227.229    0.000 interpolate.py:454(_check_bounds)</div><div>  6979970   96.567    0.000   96.567    0.000 {method 'reduce' of 'numpy.ufunc'objects}</div>
<div>  3489985   44.853    0.000  738.135    0.000 interpolate.py:443(_evaluate)</div><div>  7677978   41.010    0.000   41.010    0.000 {numpy.core.multiarray.array}</div><div>        5   40.430    8.086   40.621    8.124 npyio.py:882(savetxt)</div>
<div>  3489985   26.952    0.000   26.952    0.000 {method 'clip' of 'numpy.ndarray'objects}</div><div>  3489985   24.749    0.000   24.749    0.000 {method 'searchsorted' of 'numpy.ndarray' objects}</div>
<div>  3489985   22.457    0.000  828.238    0.000 polyint.py:37(__call__)</div><div>  6979970   19.720    0.000  116.287    0.000 _methods.py:31(_any)</div><div>  6979980   15.330    0.000   35.092    0.000 numeric.py:392(asarray)</div>
<div>  3489985   14.847    0.000   45.039    0.000 polyint.py:57(_prepare_x)</div><div>  3489990   12.904    0.000   12.904    0.000 {method 'reshape' of 'numpy.ndarray' objects}</div><div>  6979970   12.757    0.000  129.044    0.000 {method 'any' of 'numpy.ndarray' objects}</div>
<div>  3489985   11.624    0.000   11.624    0.000 {method 'astype' of 'numpy.ndarray' objects}</div><div>  3489985   10.077    0.000   10.077    0.000 {numpy.core.multiarray.empty}</div><div>  3489985    9.945    0.000   22.607    0.000 polyint.py:63(_finish_y)</div>
<div>  3489985    7.051    0.000    7.051    0.000 {method 'ravel' of 'numpy.ndarray' objects}</div><div>   697998    6.746    0.000    6.746    0.000 {zip}</div></div><div><br></div><div>So I think that in this way it's ok.</div>
<div><br></div><div>Thank you all very much,</div><div><br></div><div>Gabriele</div><div><br></div><div>p.s: I didn't know this way to write: is there a tutorial for this kind of operations? </div><div><span style="font-family:arial,sans-serif;font-size:13px">            kap = np.empty(x.shape)</span><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            sel = x >= 10.0</span><br style="font-family:arial,sans-serif;font-size:13px"><span style="font-family:arial,sans-serif;font-size:13px">            zsel = x[sel]</span><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            kap[sel] = 1.2533 * np.sqrt(zsel)*np.exp(-zsel)</span><br style="font-family:arial,sans-serif;font-size:13px"><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            sel = x < 0.001</span><br style="font-family:arial,sans-serif;font-size:13px"><span style="font-family:arial,sans-serif;font-size:13px">            zsel = x[sel]</span><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            kap[sel] = (2.1495 * np.exp(0.333333333 * np.log(zsel))</span><br style="font-family:arial,sans-serif;font-size:13px"><span style="font-family:arial,sans-serif;font-size:13px">                        - 1.8138 * zsel)</span><br style="font-family:arial,sans-serif;font-size:13px">
<br style="font-family:arial,sans-serif;font-size:13px"><span style="font-family:arial,sans-serif;font-size:13px">            sel = ~ ((x >= 10.0) | (x < 0.001))</span><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            zsel = x[sel]</span><br style="font-family:arial,sans-serif;font-size:13px"><span style="font-family:arial,sans-serif;font-size:13px">            result = my_inter(np.log10(zsel))</span><br style="font-family:arial,sans-serif;font-size:13px">
<span style="font-family:arial,sans-serif;font-size:13px">            kap[sel] = 10**result</span><br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">2014-04-12 9:34 GMT-04:00 Gabriele Brambilla <span dir="ltr"><<a href="mailto:gb.gabrielebrambilla@gmail.com" target="_blank">gb.gabrielebrambilla@gmail.com</a>></span>:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><p>Ok, i just run Peter's code and it seems really faster...I hope to don't mistake this time!</p><div class="">

<p>Thanks</p>
<p>Gabriele</p>
<p>sent from Samsung Mobile</p>
</div><div class="gmail_quote">Il giorno 12/apr/2014 08:22, "Gabriele Brambilla" <<a href="mailto:gb.gabrielebrambilla@gmail.com" target="_blank">gb.gabrielebrambilla@gmail.com</a>> ha scritto:<div><div class="h5">
<br type="attribution"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr">Ok guys, <div>I'm not expert about profile but help me to look at it.</div><div>this one is for 715853 elements (to multiply by 5, and for each of this N*5 there is a loop of 200 times)</div><div><br></div>


<div><div><div>Sat Apr 12 04:58:50 2014    restats</div><div><br></div><div>         9636507991 function calls in 66809.764 seconds</div><div><br></div><div>   Ordered by: internal time</div><div>   List reduced from 47 to 20 due to restriction <20></div>


<div><br></div><div>   ncalls        tottime    percall        cumtime   percall      filename:lineno(function)</div><div>        1       13548.507 13548.507   66809.692 66809.692   skymapsI.py:44(mymain)</div><div>125800544 13539.337   0.000        15998.925    0.000       interpolate.py:394(_call_linear)</div>


<div>880603808 5353.382    0.000         5353.382    0.000       {numpy.core.multiarray.array}<br></div><div>715853000 4998.740    0.000         52861.634  0.000        instruments.py:10(kappa)</div><div>251601088 4550.940    0.000         4550.940    0.000    {method 'reduce' of 'numpy.ufunc' objects}</div>


<div>125800544 4312.078    0.000        10163.614    0.000      interpolate.py:454(_check_bounds)</div><div>125800544 2944.126    0.000        14182.917    0.000 interpolate.py:330(__init__)<br></div><div>125800544 2846.577    0.000         29484.248    0.000 interpolate.py:443(_evaluate)</div>


<div>125800544 1665.852    0.000        6000.603    0.000 polyint.py:82(_set_yi)</div><div>125800544 1039.455    0.000         1039.455    0.000 {method 'clip' of 'numpy.ndarray' objects}</div><div>251601088  944.848    0.000           944.848    0.000 {method 'reshape' of 'numpy.ndarray' objects}</div>


<div>251601088  922.928    0.000          <a href="tel:1651.218%20%C2%A0%20%C2%A00.000" value="+16512180000" target="_blank">1651.218    0.000</a> numerictypes.py:735(issubdtype)</div><div>503202176  897.044    0.000           3434.768    0.000 numeric.py:392(asarray)</div>

<div>125800544  816.401    0.000         32242.481    0.000 polyint.py:37(__call__)</div>
<div>251601088  787.593    0.000          5338.533    0.000 _methods.py:31(_any)</div><div>125800544  689.779    0.000          1989.101    0.000 polyint.py:74(_reshape_yi)</div><div>125800544  638.946    0.000          638.946    0.000 {method 'searchsorted' of 'numpy.ndarray' objects}</div>


<div>125800544  606.778    0.000         2257.996    0.000 polyint.py:102(_set_dtype)</div><div>125800544  598.000    0.000            6598.602    0.000 polyint.py:30(__init__)</div><div>629002720  549.358    0.000           549.358    0.000 {issubclass}</div>


</div><div><br></div></div><div><br></div><div>looking at tottime it seems that skymaps mymain() and interpolate take the same big amount of time...right?</div><div><br></div><div>So it's true that I have to slow down mymain() but interpolate is a problem too!</div>


<div><br></div><div>do you agree with me? </div><div><br></div><div>Now I will read Peter Otten's code and run the new simulation with it</div><div><br></div><div>thanks</div><div><br></div><div>Gabriele</div></div><div class="gmail_extra">


<br><br><div class="gmail_quote">2014-04-12 6:21 GMT-04:00 Peter Otten <span dir="ltr"><<a href="mailto:__peter__@web.de" target="_blank">__peter__@web.de</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">


<div>Gabriele Brambilla wrote:<br>
<br>
> Ok guys, when I wrote that email I was excited for the apparent speed<br>
> increasing (it was jumping the bottleneck for loop for the reason peter<br>
> otten outlined).<br>
> Now, instead the changes, the speed is not improved (the code still<br>
> running from this morning and it's at one forth of the dataset).<br>
><br>
> What can I do to speed it up?<br>
<br>
</div>Not as easy as I had hoped and certainly not as pretty, here's my<br>
modification of the code you sent me. What makes it messy is that<br>
I had to inline your kappa() function; my first attempt with<br>
numpy.vectorize() didn't help much. There is still stuff in the<br>
'for gammar...' loop that doesn't belong there, but I decided it<br>
was time for me to stop ;)<br>
<br>
Note that it may still be worthwhile to consult a numpy expert<br>
(which I'm not!).<br>
<br>
from scipy import stats<br>
import matplotlib.pyplot as plt<br>
from scipy import optimize<br>
from matplotlib import colors, ticker, cm<br>
import numpy as np<br>
<br>
phamin = 0<br>
phamax = 2*pi<br>
obamin = 0<br>
obamax = pi<br>
npha = 100<br>
nobs = 181<br>
stepPHA = (phamax-phamin)/npha<br>
stepOB = (obamax-obamin)/nobs<br>
freq = 10<br>
c = 2.9979*(10**(10))<br>
e = 4.8032*(10**(-10))<br>
hcut = 1.0546*(10**(-27))<br>
eVtoErg = 1.6022*(10**(-12))<br>
<br>
from math import *<br>
import numpy as np<br>
from scipy.interpolate import interp1d<br>
<br>
kaparg = [<br>
    -3.0, -2.0, -1.52287875, -1.22184875, -1.0, -0.69897,<br>
     -0.52287875, -0.39794001, -0.30103, -0.22184875,<br>
     -0.15490196,  0.0, 0.30103, 0.60205999,  0.69897,<br>
     0.77815125,  0.90308999,  1.0]<br>
<br>
kapval = [<br>
    -0.6716204 , -0.35163999, -0.21183163, -0.13489603,<br>
     -0.0872467 , -0.04431225, -0.03432803, -0.04335142,<br>
     -0.05998184, -0.08039898, -0.10347378, -0.18641901,<br>
     -0.52287875, -1.27572413, -1.66958623, -2.07314329,<br>
     -2.88941029, -3.7212464 ]<br>
<br>
my_inter = interp1d(kaparg, kapval)<br>
<br>
def LEstep(n):<br>
    Emin = 10**6<br>
    Emax = 5*(10**10)<br>
    Lemin = log10(Emin)<br>
    Lemax = log10(Emax)<br>
    stepE = (Lemax-Lemin)/n<br>
    return stepE, n, Lemin, Lemax<br>
<div><br>
def mymain(stepENE, nex, Lemin, Lemax, freq):<br>
</div>    eel = np.array(list(range(nex)))<br>
<div>    eels = np.logspace(Lemin, Lemax, num=nex, endpoint=False)<br>
<br>
</div>    rlc = c/(2*pi*freq)<br>
<div><br>
    sigmas = [1, 3, 5, 10, 30]<br>
</div>    MYMAPS = [<br>
        np.zeros([npha, nobs, nex], dtype=float) for _ in sigmas]<br>
<br>
    alpha = '60_'<br>
    ALPHA = (1.732050808/c)*(e**2)<br>
    for count, my_line in enumerate(open('datasm0_60_5s.dat')):<br>
        myinternet = []<br>
<div>        print('reading the line', count, '/599378')<br>
</div>        my_parts = np.array(my_line.split(), dtype=float)<br>
<div>        phase = my_parts[4]<br>
        zobs = my_parts[5]<br>
        rho = my_parts[6]<br>
<br>
</div>        gmils = my_parts[7:12]<br>
<div><br>
        i = int((phase-phamin)/stepPHA)<br>
        j = int((zobs-obamin)/stepOB)<br>
<br>
</div><div>        for gammar, MYMAP in zip(gmils, MYMAPS):<br>
<br>
</div><div>            omC = (1.5)*(gammar**3)*c/(rho*rlc)<br>
            gig = omC*hcut/eVtoErg<br>
<br>
</div>            omega = (10**(eel*stepENE+Lemin))*eVtoErg/hcut<br>
            x = omega/omC<br>
<br>
            kap = np.empty(x.shape)<br>
            sel = x >= 10.0<br>
            zsel = x[sel]<br>
            kap[sel] = 1.2533 * np.sqrt(zsel)*np.exp(-zsel)<br>
<br>
            sel = x < 0.001<br>
            zsel = x[sel]<br>
            kap[sel] = (2.1495 * np.exp(0.333333333 * np.log(zsel))<br>
                        - 1.8138 * zsel)<br>
<br>
            sel = ~ ((x >= 10.0) | (x < 0.001))<br>
            zsel = x[sel]<br>
            result = my_inter(np.log10(zsel))<br>
            kap[sel] = 10**result<br>
<br>
            Iom = ALPHA*gammar*kap<br>
<div>            P = Iom*(c/(rho*rlc))/(2*pi)<br>
            phps = P/(hcut*omega)<br>
            www =  phps/(stepPHA*sin(zobs)*stepOB)<br>
</div>            MYMAP[i,j] += www<br>
<br>
    for sigma, MYMAP in zip(sigmas, MYMAPS):<br>
        print(sigma)<br>
        filename = "_".join(str(p) for p in<br>
            ["skymap", alpha, sigma, npha, phamin, phamax, nobs,<br>
            obamin, obamax, nex, Lemin, Lemax, '.dat']<br>
            )<br>
<br>
        x, y, z = MYMAP.shape<br>
        with open(filename, 'ab') as MYfile:<br>
            np.savetxt(<br>
                MYfile,<br>
                MYMAP.reshape(x*y, z, order="F").T,<br>
                delimiter=",", fmt="%s", newline=",\n")<br>
<div><br>
if __name__ == "__main__":<br>
    if len(sys.argv)<=1:<br>
        stepENE, nex, Lemin, Lemax = LEstep(200)<br>
    elif len(sys.argv)<=2:<br>
        stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))<br>
    else:<br>
        stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))<br>
        freq=float(sys.argv[2])<br>
<br>
    mymain(stepENE, nex, Lemin, Lemax, freq)<br>
<br>
<br>
</div>For reference here is the original (with the loop over gmlis<br>
instead of gmils):<br>
<br>
> import sys<br>
><br>
> from math import *<br>
> from scipy import ndimage<br>
> from scipy import stats<br>
> import matplotlib.pyplot as plt<br>
> from scipy import optimize<br>
> from matplotlib import colors, ticker, cm<br>
> import numpy as np<br>
> import cProfile<br>
> import pstats<br>
><br>
> phamin=0<br>
> phamax=2*pi<br>
> obamin=0<br>
> obamax=pi<br>
> npha=100<br>
> nobs=181<br>
> stepPHA=(phamax-phamin)/npha<br>
> stepOB=(obamax-obamin)/nobs<br>
> freq=10<br>
> c=2.9979*(10**(10))<br>
> e=4.8032*(10**(-10))<br>
> hcut=1.0546*(10**(-27))<br>
> eVtoErg=1.6022*(10**(-12))<br>
><br>
><br>
> from math import *<br>
> import numpy as np<br>
> from scipy.interpolate import interp1d<br>
><br>
><br>
> def kappa(z):<br>
>     N=18<br>
>     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]<br>
>     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,<br>
-1.66958623, -2.07314329, -2.88941029, -3.7212464 ]<br>
>     zlog=log10(z)<br>
>     if z < 0.001:<br>
>         k = 2.1495 * exp (0.333333333 * log (z)) - 1.8138 * z<br>
>         return (k)<br>
>     elif z >= 10.0:<br>
>         k = 1.2533 * sqrt (z) * exp (-z)<br>
>         return (k)<br>
>     else:<br>
>         my_inter = interp1d(kaparg, kapval)<br>
>         my_z = np.array([zlog])<br>
>         result = my_inter(my_z)<br>
>         valuelog = result[0]<br>
>         k=10**valuelog<br>
>         return(k)<br>
<div>><br>
><br>
><br>
><br>
> def LEstep(n):<br>
>     Emin=10**6<br>
>     Emax=5*(10**10)<br>
>     Lemin=log10(Emin)<br>
>     Lemax=log10(Emax)<br>
>     stepE=(Lemax-Lemin)/n<br>
>     return (stepE, n, Lemin, Lemax)<br>
><br>
><br>
</div><div>> def mymain(stepENE, nex, Lemin, Lemax, freq):<br>
><br>
><br>
</div><div>>     eel = list(range(nex))<br>
>     eels = np.logspace(Lemin, Lemax, num=nex, endpoint=False)<br>
><br>
>     indpha = list(range(npha))<br>
>     indobs = list(range(nobs))<br>
>     rlc = c/(2*pi*freq)<br>
><br>
</div>>     #creating an empty 3D vector<br>
<div>>     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,<br>

nex], dtype=float)]<br>
><br>
><br>
</div>>     count=0<br>
><br>
><br>
>     alpha = '60_'<br>
><br>
>     for my_line in open('datasm0_60_5s.dat'):<br>
<div>>         myinternet = []<br>
>         gmlis = []<br>
>         print('reading the line', count, '/599378')<br>
>         my_parts = [float(i) for i in my_line.split()]<br>
>         phase = my_parts[4]<br>
>         zobs = my_parts[5]<br>
>         rho = my_parts[6]<br>
><br>
</div><div>>         gmils=[my_parts[7], my_parts[8], my_parts[9], my_parts[10], my_parts[11]]<br>
><br>
</div><div>>         i = int((phase-phamin)/stepPHA)<br>
>         j = int((zobs-obamin)/stepOB)<br>
><br>
</div><div>>         for gammar, MYMAP in zip(gmils, MYMAPS):<br>
><br>
</div><div>>             omC = (1.5)*(gammar**3)*c/(rho*rlc)<br>
>             gig = omC*hcut/eVtoErg<br>
><br>
</div><div>>             for w in eel:<br>
>                 omega = (10**(w*stepENE+Lemin))*eVtoErg/hcut<br>
>                 x = omega/omC<br>
</div>>                 kap = kappa(x)<br>
<div>>                 Iom = (1.732050808/c)*(e**2)*gammar*kap<br>
>                 P = Iom*(c/(rho*rlc))/(2*pi)<br>
>                 phps = P/(hcut*omega)<br>
</div><div>>                 www =  phps/(stepPHA*sin(zobs)*stepOB)<br>
>                 MYMAP[i,j,w] += www<br>
><br>
>         count = count + 1<br>
><br>
><br>
><br>
</div><div>>     sigmas = [1, 3, 5, 10, 30]<br>
><br>
</div>>     multis = zip(sigmas, MYMAPS)<br>
><br>
>     for sigma, MYMAP in multis:<br>
><br>
>         print(sigma)<br>
>         filename='skymap_'+alpha+'_'+str(sigma)+'_'+str(npha)+'_'+str(phamin)+'_'+str(phamax)+'_'+str(nobs)+'_'+str(obamin)+'_'+str(obamax)+'_'+str(nex)+'_'+str(Lemin)+'_'+str(Lemax)+'_.dat'<br>



<div>><br>
>         MYfile = open(filename, 'a')<br>
>         for k in eel:<br>
>             for j in indobs:<br>
>                 for i in indpha:<br>
>                     A=MYMAP[i, j, k]<br>
>                     stringa = str(A) + ','<br>
>                     MYfile.write(stringa)<br>
>             accapo = '\n'<br>
>             MYfile.write(accapo)<br>
><br>
>         MYfile.close()<br>
><br>
><br>
</div><div>> if __name__ == "__main__":<br>
>     if len(sys.argv)<=1:<br>
>         stepENE, nex, Lemin, Lemax = LEstep(200)<br>
>     elif len(sys.argv)<=2:<br>
>         stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))<br>
>     else:<br>
>         stepENE, nex, Lemin, Lemax = LEstep(int(sys.argv[1]))<br>
>         freq=float(sys.argv[2])<br>
><br>
><br>
</div>> #mymain(stepENE, nex, Lemin, Lemax, freq)<br>
><br>
> #print('profile')<br>
> cProfile.run('mymain(stepENE, nex, Lemin, Lemax, freq)', 'restats', 'time')<br>
<div>><br>
> p = pstats.Stats('restats')<br>
> p.strip_dirs().sort_stats('name')<br>
</div>> p.sort_stats('time').print_stats(20)<br>
<div><div>><br>
<br>
<br>
_______________________________________________<br>
Tutor maillist  -  <a href="mailto:Tutor@python.org" target="_blank">Tutor@python.org</a><br>
To unsubscribe or change subscription options:<br>
<a href="https://mail.python.org/mailman/listinfo/tutor" target="_blank">https://mail.python.org/mailman/listinfo/tutor</a><br>
</div></div></blockquote></div><br></div>
</blockquote></div></div></div>
</blockquote></div><br></div>