<HTML><BODY>Hi Chris<br><br>Actually, I think the request by me was not so properly explained, I'm afraid (sorry about that!!!)<br>Let me explain precisely,<br><br>Last time I calculated matrix which is from only one file. it took exactly: 237.713999987 seconds. Let's say approx. 4 minute!<br><br>BUT, IF I use it for more than 30 files(if for one file it takes 4 minute, then 30 files will be 120 minute????) each of which has 18 000 values(this can be as an example manually any kind of array). I have to process turn by turn. In other words, I created ,after taking some advice from python users, of course, two dimensional matrix:<br>R(rows, cols) and obtained all processed data onto R matrix.<br><br>HERE, below I put some of my codes which are connected with each other. Also I showed F1() and F2()<br><br>Again the same problem here I want to optimize my codes in order to avoid "Loop" as well as to get quick response as much as possible. BUT, it seems really confusing, would be great to get help from Python programmers !!!<br>==================================<br>The codes here:<br>=================================================================<br><p>import numpy as np<br>import scipy.special as ss</p><p>from scipy.special import sph_jnyn,sph_jn,jv,yv<br>from scipy import integrate</p><p>import time<br>import os<br><br>---------------------------<br>1) Problem: no problem in this F0() function<br>---------------------------<br>Inputs: m=5+0.4j - complex number as an example!<br> x= one value - float!<br>---------------------------<br>#This function returns an, bn coefficients I don't want it to be vectorized because it is already done. it is working well!<br></p><p>def F0(m, x):</p><p>nmax = np.round(2.0+x+4.0*x**(1.0/3.0))<br> mx = m * x</p><p> j_x,jd_x,y_x,yd_x = ss.sph_jnyn(nmax, x) # sph_jnyn - from scipy special functions </p><p> j_x = j_x[1:]<br> jd_x = jd_x[1:]<br> y_x = y_x[1:]<br> yd_x = yd_x[1:]</p><p> h1_x = j_x + 1.0j*y_x<br> h1d_x = jd_x + 1.0j*yd_x</p><p> j_mx,jd_mx = ss.sph_jn(nmax, mx) # sph_jn - from scipy special functions <br> j_mx = j_mx[1:]<br> jd_mx = jd_mx[1:]</p><p> j_xp = j_x + x*jd_x<br> j_mxp = j_mx + mx*jd_mx<br> h1_xp = h1_x + x*h1d_x</p><p>m2 = m * m<br> an = (m2 * j_mx * j_xp - j_x * j_mxp)/(m2 * j_mx * h1_xp - h1_x * j_mxp)<br> bn = (j_mx * j_xp - j_x * j_mxp)/(j_mx * h1_xp - h1_x * j_mxp)<br><br>return an, bn<br><br>--------------------------------------<br>2) Problem: 1) To avoid loop<br> 2) To return values from the function (below) no matter whether 'a' array or scalar!<br>--------------------------------------<br> m=5+0.4j - for example<br> L = 30 - for example <br> a - array(one dimensional)<br>-------------------------------------- </p><p>def F1(m,L,a):</p><p>xs = pi * a / L<br> if(m.imag < 0.0):<br> m = conj(m)<br><br># Want to make sure we can accept single arguments or arrays<br>try:<br> xs.size<br> xlist = xs<br>except:<br> xlist = array(xs)<br><br>q=[ ]<br>for i,s in enumerate(xlist.flat):</p><p> if float(s)==0.0: # To avoid a singularity at x=0<br> q.append(0.0)<br> else:<br> x = np.round(s,7)<br> an,bn = F0(m,x)<br> n = arange(1.0, an.size + 1.0)<br> c = 2.0 * n + 1.0<br> q.append(((L*L)/(2*pi) * (c * (an.real + bn.real )).sum()))</p><p>return array(q)<br><br>-----------------------------<br>3) Problem: 1) I used "try" to avoid whether 'D' is singular or not!!! IS there better way beside this?<br> 2) IS there any wayy to avoid from the loop here in this case???<br>-----------------------------<br> a - array(one dimensional the same as above)<br> s - array(one dimensional)<br>----------------------------- </p><p>def F2(a,s):<br><br>try: # I used "try" to avoid whether 'D' is singular or not!!! IS there better way beside this?<br> a.size<br> Dslist = a<br>except:<br> Dslist = np.array(a)<br><br>K=zeros( size( Dslist ) )<br><br>for i,d in enumerate(Dslist.flat): # IS there any wayy to avoid from the loop here in this case???<br> if float(d)==0.0:<br> K[i]=0.0<br> else:<br> K[i] = np.exp(s**(-0.45)) * d)</p><p>return K<br><br>----------------------<br>4) Problem: F3 is depends on F1, F2 <br>----------------------<br></p><p>def F3(m=20,L=30):</p><p> F_file=loadtxt ( ' filename') # Instead, here we can create file, size is 120x150 will be 18000 values as I explained above.<br><br> val = [integrate.quad(lambda x: F1(m,L,x)*F2(x,i), 0.0, 7.0) for i in F_file]</p><p> return array(val)<br><br>F3 - what I really tried to get more efficient result, unfortunately, it took about one month for it!!! nothing was got! <br><br>I really need help, I am stuck....</p><br>Воскресенье, 30 декабря 2012, 16:35 -08:00 от Chris Barker - NOAA Federal <chris.barker@noaa.gov>:<br>
<blockquote style="border-left:1px solid #0857A6; margin:10px; padding:0 0 0 10px;" class="mailru-blockquote">
<div id=""><div class="js-helper js-readmsg-msg">
<style type="text/css"></style>
<div id="style_13569141930000000023" class="mr_read__body">
<base target="_self" href="https://e.mail.ru/cgi-bin/">
<div id="style_13569141930000000023_BODY">On Sun, Dec 30, 2012 at 3:41 AM, Happyman <<a href="sentmsg?compose&To=bahtiyor_zohidov@mail.ru">bahtiyor_zohidov@mail.ru</a>> wrote:<br>> nums=32<br>> rows=120<br>> cols=150<br>><br>> for k in range(0,nums):<br>> for i in range(0,rows):<br>> for j in range(0,cols):<br>> if float ( R[ k ] [ i ] [ j ] ) ==<br>> 0.0:<br><br>why the float() -- what data type is R?<br><br>> else:<br>> val11[ i ] [ j ], val22[ i<br>> ][ j ] = integrate.quad( lambda x : F1(x)*F2(x) , 0 , pi)<br><br>this is odd -- Do F1 and F2 depend on i,j, or k somehow? or are you<br>somehow integerting over the k-dimension? In which case, I'm guessing<br>that integration is you time killer anyway -- do some profiling to<br>know for sure.<br><br>-Chris<br><br>-- <br><br>Christopher Barker, Ph.D.<br>Oceanographer<br><br>Emergency Response Division<br>NOAA/NOS/OR&R (206) 526-6959 voice<br>7600 Sand Point Way NE (206) 526-6329 fax<br>Seattle, WA 98115 (206) 526-6317 main reception<br><br><a href="sentmsg?compose&To=Chris.Barker@noaa.gov">Chris.Barker@noaa.gov</a><br>_______________________________________________<br>NumPy-Discussion mailing list<br><a href="sentmsg?compose&To=NumPy%2dDiscussion@scipy.org">NumPy-Discussion@scipy.org</a><br><a href="http://mail.scipy.org/mailman/listinfo/numpy-discussion" target="_blank">http://mail.scipy.org/mailman/listinfo/numpy-discussion</a><br></div>
<base target="_self" href="https://e.mail.ru/cgi-bin/">
</div>
</div>
</div>
</blockquote>
<br></BODY></HTML>