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

[Chris Barker - NOAA Federal]

On Sun, Dec 30, 2012 at 6:40 PM, Happyman <bahtiyor_zohidov at mail.ru> wrote:

> 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)

in this case, you can do things like:

m = np.where(m.imag < 0.0, np.conj(m), m)

to vectorize.




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

here I use:

xs = np.asarray(xs, dtype-the_dtype_you_want)

it is essentially a no-op if xs is already an array, and will convert
it if it isn't.

> q=[ ]
> for i,s in enumerate(xlist.flat):
>
>               if float(s)==0.0: # To avoid a singularity at x=0
>                          q.append(0.0)

again, look to use np.where, or "fancy indexing":

ind = xs == 0.0
q[xs==0.0] = 0.0

>                          q.append(((L*L)/(2*pi) * (c * (an.real + bn.real
> )).sum()))

even if you do need the loop -- pre-allocate the result array (with
np.zeros() ), and then put stuf in it -- it will should be faster than
using a list.

> 3) Problem: 1)  I used "try" to avoid whether 'D' is singular or not!!! IS
> there better way beside this?

The other option is an if test -- try is faster if it's a rare
occurrence, slower if it's common.

> def F2(a,s):
> for i,d in enumerate(Dslist.flat):    # IS there any wayy to avoid from the
> loop here in this case???

see above.

note that using the where() or fancy indexing does mean you need to go
through the loop multiple times, but still probably much faster then
looping. For full-on speed for this sort of thing, Cython is a nice
option.

-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