（如何）Scipy 可以有效地将函数与数组值参数集成（没有循环）？答案

【问题标题】：(How) Can Scipy integrate functions with array valued arguments efficiently (without loops)?（如何）Scipy 可以有效地将函数与数组值参数集成（没有循环）？
【发布时间】：2020-03-01 07:41:03
【问题描述】：

我想使用高效（矢量化/并行化）方法集成一个接受数组参数的函数。

我可以使用代码中不需要的循环来获得所需的输出，如下面的（降低复杂性）示例：

import numpy as np
from scipy import integrate

c = np.array([1, 2])
r = np.array([2, 1])

def fun(p, c):
    p1 = np.array(np.zeros(c.shape))
    mask = np.array(np.sqrt(np.pi/(2*c)) < 1)
    p1[mask] = np.arccos(np.sqrt(np.pi/(2*c[mask])))

    d = np.zeros(c.shape)
    mask = np.abs(p) <= p1
    d[mask] = 1/(np.pi/(2*c[mask]**2) + np.cos(p))
    mask = np.logical_and(np.abs(p) > p1, np.abs(p) <= np.pi/2)
    d[mask] = 1/(np.pi/(2*c[mask]**2) +
                     ((np.cos(p1[mask]) - np.cos(p))/2))
    return(d)

def intgd(p, r, c):
    A = np.ones((np.size(r), np.size(c)))

    s = np.sin(r) - np.sin(p)
    A[s != 0] = np.sin(c[s != 0]*s[s != 0])/(c[s != 0]*s[s != 0])

    return 1/fun(p, c)**2*(A**2)

res = np.zeros((np.size(r), np.size(c)))
for ii in range(0, np.size(r)):
    for jj in range(0, np.size(c)):
        res[ii, jj], err = integrate.quad(intgd, -np.pi/2, np.pi/2,
                                           epsabs=1e-10, limit=100,
                                           args=(r[ii], c[jj]))

但是，我的实际函数需要处理更大的数组输入，这会导致计算持续时间过长。

我尝试了以下变体，并获得了以下知识（如 in comments on this question 所述），scipy.integrate.quadrature 的 vec_func=True 选项实际上并不能将向量值参数作为参数传递给函数被整合。 [除此之外：这使它与 MATLAB integral 函数完全不同，ArrayValued, true 选项确实启用了该功能，这导致积分评估更快，显然是并行化的。]

import numpy as np
from scipy import integrate

c = np.array([1, 2], ndmin=2)
r = np.array([2, 1])
r = r[:, np.newaxis]

def fun(p, c):
    p1 = np.zeros(c.shape)
    mask = np.array(np.sqrt(np.pi/(2*c)) < 1, ndmin=2)
    p1[mask] = np.arccos(np.sqrt(np.pi/(2*c[mask])))

    d = np.zeros(c.shape)
    mask = np.abs(p) <= p1
    d[mask] = 1/(np.pi/(2*c[mask]**2) + np.cos(p))
    mask = np.logical_and(np.abs(p) > p1, np.abs(p) <= np.pi/2)
    d[mask] = 1/(np.pi/(2*c[mask]**2) +
                     ((np.cos(p1[mask]) - np.cos(p))/2))
    return(d)

def intgd(p, r, c):
    A = np.ones((np.size(r), np.size(c)))

    c_bcr = np.broadcast_to(c, (np.size(r), np.size(c)))
    r_bcc = np.broadcast_to(r, (np.size(r), np.size(c)))

    s = np.sin(r_bcc) - np.sin(p)
    A[s != 0] = np.sin(c_bcr[s != 0]*s[s != 0])/(c_bcr[s != 0]*s[s != 0])

    return 1/fun(p, c)**2*(A**2)

res, err = integrate.quadrature(intgd, -np.pi/2, np.pi/2,
                                         args=(r, c), tol=1e-10, vec_func=True)

如何使用 Scipy 集成数组参数函数而不使用循环？

【问题讨论】：

您能简单描述一下您要避免的“循环”吗？什么变量？如果我没看错quadrature，它集成了一个标量值函数。 vec_func=True 表示它可以一次处理集成变量的多个值，并为每个值返回一个值。但它仍然是标量整合，而不是多维整合。 nquad 用于多维整合。
谢谢。它循环了我想要避免的输入数组参数的值。我认为 nquad 适用于要评估多个积分变量的多个积分；此处并非如此 - 只有一个积分变量：在示例中为 p。
即r和c的值的不同积分（广播与否）？
在示例中，p 是积分的变量，而 r 和 c 作为附加参数传入。由于 r 和 c 是值数组，因此它们在嵌套的 for 循环中循环，以生成形状为 (np.size(r), np.size(c)) 的输出，其中包含对 p 的积分结果对于 r 和 c 的每个组合。我不想遍历 r 和 c，我想要一个集成函数或方法，它将这些变量识别为数组值参数，并以有效的并行方式评估 r 和 c 的所有组合在 p 上的集成。
我认为您找到了关于该主题的最佳链接。我在回答关于solve_ivp、stackoverflow.com/questions/54991056/… 的类似问题时也发现了它。 scipy 没有 MATLAB 提供的所有 wiz-bang 功能。

标签： python python-3.x scipy numerical-integration quad

【解决方案1】：

矢量化 quad_vec 将在 scipy 1.4 发布时提供。

【讨论】：

感谢您的提醒，我们将非常欢迎！

【解决方案2】：

quadpy（我的一个项目）具有矢量化计算。

【讨论】：

感谢@Nico 的建议，我查看了quadpy，它看起来是一个非常好的库。不过，对我来说，我认为它的学习曲线可能有点太陡了，也许是因为我不是真正的数学家。我发现很难理解如何应用它 - 我可以看到需要了解解决方案域的几何形状，但实际上这超出了我的能力范围。您能否建议如何将其应用于上述问题？
用scipy.integrate.quad代替quapy.quad。