6

scipy.special次数 n と引数 x のベッセル関数が xjv(n,x)でベクトル化されていることに気付きました。

In [14]: import scipy.special as sp In [16]: sp.jv(1, range(3)) # n=1, [x=0,1,2] Out[16]: array([ 0., 0.44005059, 0.57672481])

しかし、対応する球面ベッセル関数のベクトル化された形式はありませんsp.sph_jn。

In [19]: sp.sph_jn(1,range(3)) 

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-19-ea59d2f45497> in <module>()
----> 1 sp.sph_jn(1,range(3)) #n=1, 3 value array

/home/glue/anaconda/envs/fibersim/lib/python2.7/site-packages/scipy/special/basic.pyc in sph_jn(n, z)
    262     """
    263     if not (isscalar(n) and isscalar(z)):
--> 264         raise ValueError("arguments must be scalars.")
    265     if (n != floor(n)) or (n < 0):
    266         raise ValueError("n must be a non-negative integer.")

ValueError: arguments must be scalars.

さらに、球ベッセル関数は、1 回のパスで N のすべての次数を計算します。したがって、n=5引数 にベッセル関数が必要な場合x=10、n=1,2,3,4,5 が返されます。実際には、jn とその導関数を 1 回のパスで返します。

In [21]: sp.sph_jn(5,10)
Out[21]: 
(array([-0.05440211,  0.07846694,  0.07794219, -0.03949584, -0.10558929,
        -0.05553451]),
 array([-0.07846694, -0.0700955 ,  0.05508428,  0.09374053,  0.0132988 ,
        -0.07226858]))

API にこの非対称性が存在するのはなぜですか。ベクトル化された、または少なくともより迅速に (つまり cython で) 球状ベッセル関数を返すライブラリを知っている人はいますか?

4

3 に答える 3

4

cython 関数を記述して計算を高速化できます。最初に行う必要があるのは、fortran 関数のアドレスを取得することです。Python でSPHJこれを行う方法は次のとおりです。

from scipy import special as sp
sphj = sp.specfun.sphj

import ctypes
addr = ctypes.pythonapi.PyCObject_AsVoidPtr(ctypes.py_object(sphj._cpointer))

次に、Cython で fortran 関数を直接呼び出すことができますprange()。計算を高速化するためにマルチコアを使用していることに注意してください。

%%cython -c-Ofast -c-fopenmp --link-args=-fopenmp
from cpython.mem cimport PyMem_Malloc, PyMem_Free
from cython.parallel import prange
import numpy as np
import cython
from cpython cimport PyCObject_AsVoidPtr
from scipy import special

ctypedef void (*sphj_ptr) (const int *n, const double *x, 
                            const int *nm, const double *sj, const double *dj) nogil

cdef sphj_ptr _sphj=<sphj_ptr>PyCObject_AsVoidPtr(special.specfun.sphj._cpointer)


@cython.wraparound(False)
@cython.boundscheck(False)
def cython_sphj2(int n, double[::1] x):
    cdef int count = x.shape[0]
    cdef double * sj = <double *>PyMem_Malloc(count * sizeof(double) * (n + 1))
    cdef double * dj = <double *>PyMem_Malloc(count * sizeof(double) * (n + 1))
    cdef int * mn    = <int *>PyMem_Malloc(count * sizeof(int))
    cdef double[::1] res = np.empty(count)
    cdef int i
    if count < 100:
        for i in range(x.shape[0]):
            _sphj(&n, &x[i], mn + i, sj + i*(n+1), dj + i*(n+1))
            res[i] = sj[i*(n+1) + n]    #choose the element you want here        
    else:
        for i in prange(count,  nogil=True):
            _sphj(&n, &x[i], mn + i, sj + i*(n+1), dj + i*(n+1))
            res[i] = sj[i*(n+1) + n]    #choose the element you want here
    PyMem_Free(sj)
    PyMem_Free(dj)
    PyMem_Free(mn)
    return res.base

sphj()比較するために、forloopを呼び出す Python 関数を次に示します。

import numpy as np
def python_sphj(n, x):
    sphj = special.specfun.sphj
    res = np.array([sphj(n, v)[1][n] for v in x])
    return res

以下は、10 個の要素に対する %timit の結果です。

x = np.linspace(1, 2, 10)
r1 = cython_sphj2(4, x)
r2 = python_sphj(4, x)
assert np.allclose(r1, r2)
%timeit cython_sphj2(4, x)
%timeit python_sphj(4, x)

結果：

10000 loops, best of 3: 21.5 µs per loop
10000 loops, best of 3: 28.1 µs per loop

100000 要素の結果は次のとおりです。

x = np.linspace(1, 2, 100000)
r1 = cython_sphj2(4, x)
r2 = python_sphj(4, x)
assert np.allclose(r1, r2)
%timeit cython_sphj2(4, x)
%timeit python_sphj(4, x)

結果：

10 loops, best of 3: 44.7 ms per loop
1 loops, best of 3: 231 ms per loop
于 2015-02-11T01:15:13.480 に答える