[Numpy-discussion] 3D array problem challenging in Python
Happyman
bahtiyor_zohidov@mail...
Wed Jan 2 19:27:26 CST 2013
Hi Chris,
Thanks a lot.
I did as you advised..but, unfortunately, I could not "negotiate" with "quad" function at all.
what do you think if quad function can get arrays or not?
Вторник, 1 января 2013, 19:53 -08:00 от Chris Barker - NOAA Federal <chris.barker@noaa.gov>:
>On Sun, Dec 30, 2012 at 6:40 PM, Happyman < bahtiyor_zohidov@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@noaa.gov
>_______________________________________________
>NumPy-Discussion mailing list
>NumPy-Discussion@scipy.org
>http://mail.scipy.org/mailman/listinfo/numpy-discussion
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20130103/b473c855/attachment.html
More information about the NumPy-Discussion
mailing list