(How)Scipy能否有效地集成具有数组值参数的函数(无循环)

Question

我想集成一个使用高效(矢量化/并行化(方法获取数组参数的函数。

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

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]))

然而，我的实际函数需要处理更大的数组输入，这导致计算持续时间过长，令人无法接受。

我尝试了以下方面的变体，并了解到(正如对这个问题的评论中所指出的(，scipy.integrate.quadrature的vec_func=True选项实际上并不能使向量值自变量作为参数传递到正在集成的函数中。[除此之外：这与MATLABintegral函数有很大不同，选项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在不使用循环的情况下集成数组参数函数？

Answer 1

矢量化quad_vec将在scipy 1.4中发布。

Answer 2

quadry(我的一个项目(已经将计算矢量化。