[Numpy-discussion] 3D array problem challenging in Python

[Happyman]

 Hi Chris

Actually, I think the request by me was not so properly explained, I'm afraid (sorry about that!!!)
Let me explain precisely,

Last time I calculated matrix which is from only one file. it took exactly:  237.713999987 seconds. Let's say approx. 4 minute!

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:
R(rows, cols) and obtained all processed data onto R matrix.

HERE, below I put some of my codes which are connected with each other. Also I showed F1() and F2()

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 !!!
==================================
The codes here:
=================================================================
import numpy as np
import scipy.special as ss
from scipy.special import sph_jnyn,sph_jn,jv,yv
from scipy import integrate
import time
import os

---------------------------
1) Problem: no problem in this F0() function
---------------------------
Inputs: m=5+0.4j  - complex number as an example!
            x= one value - float!
---------------------------
#This function returns an, bn coefficients I don't want it to be vectorized because it is already done. it is working well!
def F0(m, x):
nmax = np.round(2.0+x+4.0*x**(1.0/3.0))
mx = m * x
j_x,jd_x,y_x,yd_x = ss.sph_jnyn(nmax, x)        #    sph_jnyn   -    from scipy special functions    
j_x = j_x[1:]
jd_x = jd_x[1:]
y_x = y_x[1:]
yd_x = yd_x[1:]
h1_x = j_x + 1.0j*y_x
h1d_x = jd_x + 1.0j*yd_x
j_mx,jd_mx = ss.sph_jn(nmax, mx)        #    sph_jn    -    from scipy special functions    
j_mx = j_mx[1:]
jd_mx = jd_mx[1:]
j_xp = j_x + x*jd_x
j_mxp = j_mx + mx*jd_mx
h1_xp = h1_x + x*h1d_x
m2 = m * m
an = (m2 * j_mx * j_xp - j_x * j_mxp)/(m2 * j_mx * h1_xp - h1_x * j_mxp)
bn = (j_mx * j_xp - j_x * j_mxp)/(j_mx * h1_xp - h1_x * j_mxp)

return an, bn

--------------------------------------
2) Problem:    1) To avoid loop
                       2) To return values from the function (below) no matter whether 'a' array or scalar!
--------------------------------------
  m=5+0.4j  - for example
  L = 30 - for example 
  a - array(one dimensional)
-------------------------------------- 
def F1(m,L,a):
xs = pi * a / L
if(m.imag < 0.0):
             m = conj(m)

# Want to make sure we can accept single arguments or arrays
try:
     xs.size
     xlist = xs
except:
     xlist = array(xs)

q=[ ]
for i,s in enumerate(xlist.flat):
              if float(s)==0.0: # To avoid a singularity at x=0
                         q.append(0.0)
               else:
                         x = np.round(s,7)
                         an,bn = F0(m,x)
                         n = arange(1.0, an.size + 1.0)
                         c = 2.0 * n + 1.0
                         q.append(((L*L)/(2*pi) * (c * (an.real + bn.real )).sum()))
return array(q)

-----------------------------
3) Problem: 1)  I used "try" to avoid whether 'D' is singular or not!!! IS there better way beside this?
                    2)  IS there any wayy to avoid from the loop here in this case???
-----------------------------
 a - array(one dimensional the same as above)
 s - array(one dimensional)
----------------------------- 
def F2(a,s):

try:                  # I used "try" to avoid whether 'D' is singular or not!!! IS there better way beside this?
      a.size
      Dslist = a
except:
      Dslist = np.array(a)

K=zeros( size( Dslist ) )

for i,d in enumerate(Dslist.flat):    # IS there any wayy to avoid from the loop here in this case???
              if float(d)==0.0:
                          K[i]=0.0
              else:
                          K[i] = np.exp(s**(-0.45)) * d)
return K

----------------------
4) Problem: F3 is depends on F1, F2 
----------------------
def F3(m=20,L=30):
                     F_file=loadtxt ( ' filename')     #  Instead, here we can create file, size is 120x150 will be 18000 values as I explained above.

                    val = [integrate.quad(lambda x: F1(m,L,x)*F2(x,i), 0.0, 7.0) for i in F_file]
                    return array(val)

F3 - what I really tried to get more efficient result, unfortunately, it took about one month for it!!! nothing was got! 

I really need help, I am stuck....
Воскресенье, 30 декабря 2012, 16:35 -08:00 от Chris Barker - NOAA Federal <chris.barker at noaa.gov>:
>On Sun, Dec 30, 2012 at 3:41 AM, Happyman < bahtiyor_zohidov at mail.ru > wrote:
>> nums=32
>> rows=120
>> cols=150
>>
>> for k in range(0,nums):
>>           for i in range(0,rows):
>>                      for j in range(0,cols):
>>                                       if float ( R[ k ] [ i ] [ j ] ) ==
>> 0.0:
>
>why the float() -- what data type is R?
>
>>                                       else:
>>                                                  val11[ i ] [ j ], val22[ i
>> ][ j ] = integrate.quad( lambda x :  F1(x)*F2(x) , 0 , pi)
>
>this is odd -- Do F1 and F2 depend on i,j, or k somehow? or are you
>somehow integerting over the k-dimension? In which case, I'm guessing
>that integration is you time killer anyway -- do some profiling to
>know for sure.
>
>-Chris
>
>-- 
>
>Christopher Barker, Ph.D.
>Oceanographer
>
>Emergency Response Division
>NOAA/NOS/OR&R            (206) 526-6959   voice
>7600 Sand Point Way NE   (206) 526-6329   fax
>Seattle, WA  98115       (206) 526-6317   main reception
>
>Chris.Barker at noaa.gov
>_______________________________________________
>NumPy-Discussion mailing list
>NumPy-Discussion at scipy.org
>http://mail.scipy.org/mailman/listinfo/numpy-discussion

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20121231/5449ee8b/attachment.html>