Hi, I like how jn works because I can array as arg (I use jn on 3D array). This is not the case for sph_jn & lpmn :-( I tried to bypass it, with something like def foo(m,n): x = arange(0,Lx,dx) y = arange(0,Ly,dy) z = arange(0,Lz,dz) tab = zeros((len(x), len(y), len(z)), dtype='f') for i in range(0,len(x)): for j in range(0,len(y)): for k in range(0,len(z)): tab[i,j,k] = cos(theta(x[i],y[j],z[k])) # tab[i,j,k] = Ymn_theta_p(x[i],y[j],z[k],m,n)*sph_jn(n,r(x[i],y[j],z[k]))[n][0]*rho(x[i],y[j],z[k]) return (tab) but knowing my arrays have more than hundred of thousand cells, it is _very_ slow. I guess sph_jn & lpmn are written like this because they return arrays (0 to n order, derivative, etc). I wish to have Legendre & spherical Bessel functions which could accept array as arg. I only need the last order (n) and the derivative (which could be called by another func) for the same order. How could I write this ? Or is there another way to do the same thing ? Ex. : In [1]: from scipy.special import * In [2]: x=arange(0,1,0.1) In [3]: print jn(1,x) [ 0. 0.04993753 0.09950083 0.14831882 0.19602658 0.24226846 0.28670099 0.32899574 0.36884205 0.40594955] In [4]: print sph_jn(1,x) --------------------------------------------------------------------------- exceptions.ValueError Traceback (most recent call last) /home/fred/<ipython console> /usr/local/lib/python2.4/site-packages/scipy/special/basic.py in sph_jn(n, z) 212 """ 213 if not (isscalar(n) and isscalar(z)): --> 214 raise ValueError, "arguments must be scalars." 215 if (n!= floor(n)) or (n<0): 216 raise ValueError, "n must be a non-negative integer." ValueError: arguments must be scalars. Cheers, -- Fred.
I hope someone else will find a better answer to your question, but here is something you could try if you are really hungry for speed. make a copy of the function sphj in scipy/lib/special/specfun/specfun.f and rename it tweak it so it takes a vector argument. Loop over all inputs in vector return vector of answer (I'd like to help, but my knowledge of fortran is pretty limited.) make a copy of the sphj function interface in scipy/lib/special/specfun.pyf and rename it tweak it so it accepts arrays as input and returns an array, something like subroutine sphj_1(n,x,nx, nm, sj,dj) ! in :specfun:specfun.f integer intent(in), check(n>=1) :: n double precision intent(in), dimension(nx), depend(nx) :: x integer intent(in) :: nx integer intent(out) :: nm double precision intent(out),dimension(nx),depend(nx) :: sj double precision intent(out),dimension(nx),depend(nx) :: dj end subroutine sph_1j Recompile scipy. You could also modify sphj in basic.py to do the loop there, but I doubt you'd gain much in terms of speed up. I had a similar problem once, and after coding the loop inside the f code, I got a speedup of about 6x compared to the python loop. Good luck. David 2006/7/27, fred <fredantispam@free.fr>:
Hi,
I like how jn works because I can array as arg (I use jn on 3D array).
This is not the case for sph_jn & lpmn :-(
I tried to bypass it, with something like
def foo(m,n): x = arange(0,Lx,dx) y = arange(0,Ly,dy) z = arange(0,Lz,dz) tab = zeros((len(x), len(y), len(z)), dtype='f') for i in range(0,len(x)): for j in range(0,len(y)): for k in range(0,len(z)): tab[i,j,k] = cos(theta(x[i],y[j],z[k])) # tab[i,j,k] =
Ymn_theta_p(x[i],y[j],z[k],m,n)*sph_jn(n,r(x[i],y[j],z[k]))[n][0]*rho(x[i],y[j],z[k]) return (tab)
but knowing my arrays have more than hundred of thousand cells, it is _very_ slow.
I guess sph_jn & lpmn are written like this because they return arrays (0 to n order, derivative, etc).
I wish to have Legendre & spherical Bessel functions which could accept array as arg. I only need the last order (n) and the derivative (which could be called by another func) for the same order.
How could I write this ?
Or is there another way to do the same thing ?
Ex. :
In [1]: from scipy.special import *
In [2]: x=arange(0,1,0.1)
In [3]: print jn(1,x) [ 0. 0.04993753 0.09950083 0.14831882 0.19602658 0.24226846 0.28670099 0.32899574 0.36884205 0.40594955]
In [4]: print sph_jn(1,x)
--------------------------------------------------------------------------- exceptions.ValueError Traceback (most recent call last)
/home/fred/<ipython console>
/usr/local/lib/python2.4/site-packages/scipy/special/basic.py in sph_jn(n, z) 212 """ 213 if not (isscalar(n) and isscalar(z)): --> 214 raise ValueError, "arguments must be scalars." 215 if (n!= floor(n)) or (n<0): 216 raise ValueError, "n must be a non-negative integer."
ValueError: arguments must be scalars.
Cheers,
-- Fred. _______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
David Huard a écrit :
I hope someone else will find a better answer to your question, but here is something you could try if you are really hungry for speed.
[snip] Ok, thanks for the hint, I'm working on it. However, I would like to see how jn(n,x) is written. I did not find nothing about it (no "def jn:"), except the fact that jn(n,x) is written in C, in cephes/. But I did not find out how the interface C/python is written, as in fortran (source specfun.f -> specfun.pyf -> basic.py). I'm going on the bad way ? Any suggestion are welcome. Cheers, -- Fred.
fred wrote:
David Huard a écrit :
I hope someone else will find a better answer to your question, but here is something you could try if you are really hungry for speed.
[snip]
Ok, thanks for the hint, I'm working on it.
However, I would like to see how jn(n,x) is written. I did not find nothing about it (no "def jn:"), except the fact that jn(n,x) is written in C, in cephes/. But I did not find out how the interface C/python is written, as in fortran (source specfun.f -> specfun.pyf -> basic.py). I'm going on the bad way ?
It is a ufunc, so it is created by calling PyUFunc_FromFuncAndData in Lib/special/_cephesmodule.c . The actual implementation of the math-bits is in Lib/special/cephes/jn.c . -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Robert Kern a écrit :
It is a ufunc, so it is created by calling PyUFunc_FromFuncAndData in Lib/special/_cephesmodule.c . Hmm, it's a little too obscure for me :-(
Does it mean that the solution proposed by David is the only one ? I don't like it because you have to pass the dimension of x as arg and because of the loop (thanks anyway, David, you showed me an interesting way ;-) Cheers, -- Fred.
fred wrote:
Robert Kern a écrit :
It is a ufunc, so it is created by calling PyUFunc_FromFuncAndData in Lib/special/_cephesmodule.c . Hmm, it's a little too obscure for me :-(
Does it mean that the solution proposed by David is the only one ? I don't like it because you have to pass the dimension of x as arg and because of the loop (thanks anyway, David, you showed me an interesting way ;-)
It's probably the easiest one. There is a reason that the current forms of the spherical harmonic functions are not ufuncs like jn is: they can't be shoved into the ufunc framework since they return arrays for scalar inputs. Note that in the f2py interface that David wrote, you won't have to pass in the dimension of x as an argument from Python; you just have to write the underlying FORTRAN subroutine to do so. The f2py interface takes care of inferring the dimensions of your inputs for you. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco
Robert Kern a écrit : > fred wrote: > >>Robert Kern a écrit : >> >> >>>It is a ufunc, so it is created by calling PyUFunc_FromFuncAndData in >>>Lib/special/_cephesmodule.c . >> >>Hmm, it's a little too obscure for me :-( >> >>Does it mean that the solution proposed by David is the only one ? >>I don't like it because you have to pass the dimension of x as arg and >>because of the loop (thanks anyway, David, you showed me an interesting >>way ;-) > > > It's probably the easiest one. There is a reason that the current forms of the > spherical harmonic functions are not ufuncs like jn is: they can't be shoved > into the ufunc framework since they return arrays for scalar inputs. Yes, this is the reason I wrote mylpmn for instance, which returns - only the Legendre function and not its derivative ; - only for order m and degree n. > Note that in the f2py interface that David wrote, you won't have to pass in the > dimension of x as an argument from Python; you just have to write the underlying > FORTRAN subroutine to do so. The f2py interface takes care of inferring the > dimensions of your inputs for you. Ok, I did not really read it completely. Thanks to have pointed me on that. Cheers, -- Fred.
Robert Kern a écrit : [snip]
It's probably the easiest one. There is a reason that the current forms of the spherical harmonic functions are not ufuncs like jn is: they can't be shoved into the ufunc framework since they return arrays for scalar inputs.
marsu[pts/0]:~/{5}/> ipython Python 2.4.3 (#1, Jul 25 2006, 18:55:45) Type "copyright", "credits" or "license" for more information. IPython 0.7.2 -- An enhanced Interactive Python. ? -> Introduction to IPython's features. %magic -> Information about IPython's 'magic' % functions. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [1]: from scipy.special import * In [2]: help(lpmn) Help on function lpmn in module scipy.special.basic: lpmn(m, n, z) Associated Legendre functions of the first kind, Pmn(z) and its derivative, Pmn'(z) of order m and degree n. Returns two arrays of size (m+1,n+1) containing Pmn(z) and Pmn'(z) for all orders from 0..m and degrees from 0..n. z can be complex. In [3]: help(jn) Help on ufunc object: jn = class ufunc(__builtin__.object) | Optimized functions make it possible to implement arithmetic with arrays efficiently | | Methods defined here: | | __call__(...) | x.__call__(...) <==> x(...) | | __repr__(...) | x.__repr__() <==> repr(x) | | __str__(...) | x.__str__() <==> str(x) | | accumulate(...) | | outer(...) | | reduce(...) | | reduceat(...) .../... How can I get info about jn functions ?? Cheers, -- Fred.
fred wrote:
Robert Kern a écrit :
[snip]
It's probably the easiest one. There is a reason that the current forms of the spherical harmonic functions are not ufuncs like jn is: they can't be shoved into the ufunc framework since they return arrays for scalar inputs.
marsu[pts/0]:~/{5}/> ipython Python 2.4.3 (#1, Jul 25 2006, 18:55:45) Type "copyright", "credits" or "license" for more information.
IPython 0.7.2 -- An enhanced Interactive Python. ? -> Introduction to IPython's features. %magic -> Information about IPython's 'magic' % functions. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more.
In [1]: from scipy.special import *
In [2]: help(lpmn) Help on function lpmn in module scipy.special.basic:
lpmn(m, n, z) Associated Legendre functions of the first kind, Pmn(z) and its derivative, Pmn'(z) of order m and degree n. Returns two arrays of size (m+1,n+1) containing Pmn(z) and Pmn'(z) for all orders from 0..m and degrees from 0..n.
z can be complex.
In [3]: help(jn) Help on ufunc object:
jn = class ufunc(__builtin__.object) | Optimized functions make it possible to implement arithmetic with arrays efficiently | | Methods defined here: | | __call__(...) | x.__call__(...) <==> x(...) | | __repr__(...) | x.__repr__() <==> repr(x) | | __str__(...) | x.__str__() <==> str(x) | | accumulate(...) | | outer(...) | | reduce(...) | | reduceat(...) .../...
How can I get info about jn functions ??
Cheers,
You can try ipython Python 2.4.1 (#1, Sep 12 2005, 23:33:18) Type "copyright", "credits" or "license" for more information. IPython 0.7.3.svn -- An enhanced Interactive Python. ? -> Introduction to IPython's features. %magic -> Information about IPython's 'magic' % functions. help -> Python's own help system. object? -> Details about 'object'. ?object also works, ?? prints more. In [1]: from scipy import * In [2]: special.jn? Type: ufunc Base Class: <type 'numpy.ufunc'> String Form: <ufunc 'jn'> Namespace: Interactive Docstring: y = jn(x1,x2) y=jn(n,x) returns the Bessel function of integer order n at x. In [3]: info(special.jn) y = jn(x1,x2) y=jn(n,x) returns the Bessel function of integer order n at x. Nils
fred wrote:
Nils Wagner a écrit :
You can try
[snip]
Ok, but why help() does not work on jn() ? Because it's a ufunc ?
For some reason Python's help does not print the docstring of an instance but only of it's type (i.e. the general ufunc). I think this is broken behavior. IPython does the right thing. So does scipy's info function. -Travis
Robert Kern a écrit : [snip]
Note that in the f2py interface that David wrote, you won't have to pass in the dimension of x as an argument from Python; you just have to write the underlying FORTRAN subroutine to do so. The f2py interface takes care of inferring the dimensions of your inputs for you.
That works quite fine, thanks to all you have helped/answered to me : I wrote lpmn3D() and sphjn3D() functions that fit my needs. Thanks again. Cheers, -- Fred.
participants (5)
-
David Huard -
fred -
Nils Wagner -
Robert Kern -
Travis Oliphant