[SciPy] Special functions

Special functions

 

패키지의 주요 기능 scipy.special은 수학 물리학의 수많은 특수 기능의 정의이다. 사용 가능한 함수에는 airy, elliptic, bessel, gamma, beta, hypergeometric, parabolic cylinder, mathieu, spheroidal wave, struve 및 kelvin이 포함된다. 또한, 이러한 기능에 대한 더 쉬운 인터페이스가 stats 모듈에서 제공되기 때문에 일반적인 사용을 위한 것이 아닌 일부 저수준 통계 기능이 있다. 이러한 함수의 대부분은 Numerical Python의 다른 수학 함수와 동일한 브로드캐스팅 규칙에 따라 배열 인수를 취하고 배열 결과를 반환할 수 있다. 이러한 함수 중 다수는 복소수를 입력으로 받아들인다. 각 기능에는 도움말을 사용하여 액세스할 수 있는 자체 문서도 있다. 필요한 함수가 표시되지 않으면 함수를 작성하여 라이브러리에 기여하는 것을 고려해야 한다. C, Fortran 또는 Python으로 함수를 작성할 수 있다.

 

Bessel functions of real order (jv, jn_zeros)

 

Bessel 함수는 실수 또는 복소수 차수 알파를 사용하는 Bessel의 미분 방정식에 대한 솔루션 제품군이다.

 

from scipy import special
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt

x = np.arange(-5,5,0.1)
y = np.arange(-5,5,0.1)
X,Y = np.meshgrid(x,y)
Z = X*np.exp(-X**2 - Y**2)

fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111, projection='3d')

# Plot a 3D surface
ax.plot_surface(X, Y, Z)

plt.show()

 

Cython Bindings for Special Functions

 

SciPy는 또한 스칼라에 대한 Cython 바인딩을 제공하며 특수한 많은 기능의 유형이 지정된 버전이다.

 

cimport scipy.special.cython_special as csc

cdef:
    double x = 1
    double complex z = 1 + 1j
    double si, ci, rgam
    double complex cgam

rgam = csc.gamma(x)
print(rgam)
cgam = csc.gamma(z)
print(cgam)
csc.sici(x, &si, &ci)
print(si, ci)

 

Cython 바인딩을 사용하면 두 가지 잠재적 이점이 있다.

 

Python 함수 오버헤드를 방지한다. Python Global Interpreter Lock (GIL)이 필요하지 않다

 

Python 함수 오버헤드 방지

 

특수한 ufunc의 경우 벡터화, 즉 함수에 배열을 전달하여 Python 함수 오버헤드를 방지한다. 일반적으로 이 접근 방식은 꽤 잘 작동하지만 때로는 예를 들어, 자체 ufunc를 구현할 때 루프 내부의 스칼라 입력에 대한 특수 함수를 호출하는 것이 더 편리하다. 이 경우 Python 함수 오버헤드가 중요해질 수 있다.

 

import scipy.special as sc
cimport scipy.special.cython_special as csc

def python_tight_loop():
    cdef:
        int n
        double x = 1

    for n in range(100):
        sc.jv(n, x)

def cython_tight_loop():
    cdef:
        int n
        double x = 1

    for n in range(100):
        csc.jv(n, x)

 

한 컴퓨터 python_tight_loop에서 실행하는 데 약 131마이크로초가 걸렸고 cython_tight_loop은 실행하는 데 약 18.2마이크로초가 걸렸다. 이처럼 빠르게 호출하고 결과를 얻을 수 있다. 요점은 Python 함수 오버헤드가 코드에서 중요해지면 Cython 바인딩이 유용할 수 있다는 것이다.

 

GIL

 

종종 여러 지점에서 특수 기능을 평가해야 하며 일반적으로 평가는 사소하게 병렬화할 수 있다. Cython 바인딩에는 GIL이 필요하지 않으므로 Cython의 prange 기능을 사용하여 병렬로 쉽게 실행할 수 있다. 예를 들어, Helmholtz 방정식에 대한 기본 솔루션을 계산하려고 한다고 가정한다.

 

 

k는 파수이고 δ는 Dirac 델타 함수이다. 2차원에서 고유 (발산) 솔루션은 다음과 같다고 알려져 있다.

 

 

H0(1)는 제1종 한켈 함수, 즉 함수 hankel1이다. 다음 예제는 이 함수를 병렬로 계산하는 방법을 보여준다.

 

from libc.math cimport fabs
cimport cython
from cython.parallel cimport prange

import numpy as np
import scipy.special as sc
cimport scipy.special.cython_special as csc

def serial_G(k, x, y):
    return 0.25j*sc.hankel1(0, k*np.abs(x - y))

@cython.boundscheck(False)
@cython.wraparound(False)
cdef void _parallel_G(double k, double[:,:] x, double[:,:] y,
                      double complex[:,:] out) nogil:
    cdef int i, j

    for i in prange(x.shape[0]):
        for j in range(y.shape[0]):
            out[i,j] = 0.25j*csc.hankel1(0, k*fabs(x[i,j] - y[i,j]))

def parallel_G(k, x, y):
    out = np.empty_like(x, dtype='complex128')
    _parallel_G(k, x, y, out)
    return out

import timeit
import numpy as np
from test import serial_G, parallel_G

def main():
    k = 1
    x, y = np.linspace(-100, 100, 1000), np.linspace(-100, 100, 1000)
    x, y = np.meshgrid(x, y)

    def serial():
        serial_G(k, x, y)

    def parallel():
        parallel_G(k, x, y)

    time_serial = timeit.timeit(serial, number=3)
    time_parallel = timeit.timeit(parallel, number=3)
    print("Serial method took {:.3} seconds".format(time_serial))
    print("Parallel method took {:.3} seconds".format(time_parallel))

if __name__ == "__main__":
    main()

 

한 쿼드 코어 컴퓨터에서 직렬 방식은 1.29초, 병렬 방식은 0.29초가 걸렸다.

 

https://docs.scipy.org/doc/scipy/tutorial/special.html

Special functions (scipy.special) — SciPy v1.9.3 Manual