0

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

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

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

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

我尝试了以下变化，并获得了知识（如对此问题的评论中所述），vec_func=True选项scipy.integrate.quadrature实际上并不能将向量值参数作为参数传递到正在集成的函数中。[旁白：这使得它与 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 集成数组参数函数而不使用循环？

4

2 回答 2

1

Vectorizedquad_vec将在 scipy 1.4 发布时提供。

于 2019-11-05T13:33:22.903 回答
1

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

于 2019-11-05T21:59:11.803 回答