Ctrl+K

特殊函数 (scipy.special)#

scipy.special 包的主要特点是定义了大量的数学物理特殊函数。可用的函数包括 airy、elliptic、bessel、gamma、beta、hypergeometric、抛物柱面、mathieu、球状波、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\) 是狄拉克δ函数。已知在二维空间中，唯一的（辐射）解是

\[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)
在本页