Ctrl+K

特殊函數 (scipy.special)#

scipy.special 套件的主要功能是定義數學物理中眾多的特殊函數。可用的函數包括 airy、elliptic、bessel、gamma、beta、hypergeometric、parabolic cylinder、mathieu、spheroidal wave、struve 和 kelvin 函數。還有一些低階統計函數，這些函數不適合一般使用，因為 stats 模組提供了更簡單的介面來使用這些函數。這些函數大多數可以接受陣列引數並傳回陣列結果，遵循與 Numerical Python 中其他數學函數相同的廣播規則。這些函數中許多也接受複數作為輸入。如需包含單行描述的可用函數完整列表，請輸入 >>> help(special). 每個函數也都有自己的說明文件，可以使用 help 存取。如果您沒有看到您需要的函數，請考慮編寫它並貢獻到程式庫中。您可以使用 C、Fortran 或 Python 編寫函數。請在程式庫的原始碼中尋找每種類型函數的範例。

實數階貝索函數 (jv, jn_zeros)#

貝索函數是貝索微分方程式的一系列解，具有實數或複數階 alpha

\[x^2 \frac{d^2 y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0\]

除了其他用途外，這些函數還出現在波動傳播問題中，例如薄鼓皮的振動模式。以下是以邊緣固定的圓形鼓皮為例

>>> from scipy import special
>>> import numpy as np
>>> def drumhead_height(n, k, distance, angle, t):
...    kth_zero = special.jn_zeros(n, k)[-1]
...    return np.cos(t) * np.cos(n*angle) * special.jn(n, distance*kth_zero)
>>> theta = np.r_[0:2*np.pi:50j]
>>> radius = np.r_[0:1:50j]
>>> x = np.array([r * np.cos(theta) for r in radius])
>>> y = np.array([r * np.sin(theta) for r in radius])
>>> z = np.array([drumhead_height(1, 1, r, theta, 0.5) for r in radius])

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_axes(rect=(0, 0.05, 0.95, 0.95), projection='3d')
>>> ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap='RdBu_r', vmin=-0.5, vmax=0.5)
>>> ax.set_xlabel('X')
>>> ax.set_ylabel('Y')
>>> ax.set_xticks(np.arange(-1, 1.1, 0.5))
>>> ax.set_yticks(np.arange(-1, 1.1, 0.5))
>>> ax.set_zlabel('Z')
>>> plt.show()
"This code generates a 3-D representation of the vibrational modes on a drum head viewed at a three-quarter angle. A circular region on the X-Y plane is defined with a Z value of 0 around the edge. Within the circle a single smooth valley exists on the -X side and a smooth peak exists on the +X side. The image resembles a yin-yang at this angle."

特殊函數的 Cython 綁定 (scipy.special.cython_special)#

SciPy 也為 special 中許多函數的純量、類型化版本提供了 Cython 綁定。以下 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 的說明，請參閱 Cython 文件。）在範例中，函數 csc.gamma 的運作方式基本上與其 ufunc 對應物 gamma 類似，儘管它採用 C 類型作為引數，而不是 NumPy 陣列。請特別注意，該函數經過重載以支援實數和複數引數；編譯時會選取正確的變體。函數 csc.sici 的運作方式與 sici 略有不同；對於 ufunc，我們可以寫成 ai, bi = sici(x)，而在 Cython 版本中，多個傳回值會作為指標傳遞。將其視為類似於使用輸出陣列呼叫 ufunc 可能會有幫助：sici(x, out=(si, ci))。

使用 Cython 綁定有兩個潛在優勢

  • 它們避免了 Python 函數的額外負擔

  • 它們不需要 Python 全域直譯器鎖 (GIL)

以下章節討論如何利用這些優勢來潛在地加速您的程式碼，當然，始終應該先分析程式碼，以確保額外的努力是值得的。

避免 Python 函數的額外負擔#

對於 special 中的 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 微秒才能執行。顯然，這個範例是人為設計的：可以直接呼叫 special.jv(np.arange(100), 1) 並獲得與 cython_tight_loop 一樣快的結果。重點是，如果 Python 函數的額外負擔在您的程式碼中變得顯著，那麼 Cython 綁定可能會很有用。

釋放 GIL#

通常需要在許多點評估特殊函數，並且通常評估可以很容易地並行化。由於 Cython 綁定不需要 GIL，因此可以使用 Cython 的 prange 函數輕鬆地並行執行它們。例如，假設我們想要計算亥姆霍茲方程式的基本解

\[\Delta_x G(x, y) + k^2G(x, y) = \delta(x - y),\]

其中 \(k\) 是波數，\(\delta\) 是狄拉克 delta 函數。已知在二維中，唯一（輻射）解是

\[G(x, y) = \frac{i}{4}H_0^{(1)}(k|x - y|),\]

其中 \(H_0^{(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

（有關在 Cython 中編譯並行程式碼的說明，請參閱此處。）如果上述 Cython 程式碼在檔案 test.pyx 中，那麼我們可以編寫一個非正式基準測試，比較該函數的並行版本和序列版本

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 秒。

不在 scipy.special 中的函數#

有些函數未包含在 special 中，因為使用 NumPy 和 SciPy 中的現有函數可以輕鬆實作。為了防止重新發明輪子，本節提供了幾個此類函數的實作，希望這些實作說明如何處理類似的函數。在所有範例中，NumPy 都以 np 匯入，而 special 則以 sc 匯入。

二元熵函數

def binary_entropy(x):
    return -(sc.xlogy(x, x) + sc.xlog1py(1 - x, -x))/np.log(2)

在 [0, 1] 上的矩形步階函數

def step(x):
    return 0.5*(np.sign(x) + np.sign(1 - x))

平移和縮放可用於取得任意步階函數。

斜坡函數

def ramp(x):
    return np.maximum(0, x)
在本頁中