如何使用数组输入调用 Python Sipy Quad

Question

我正在使用嵌套的scipy.integrate.quad调用来集成二维积分。 积分由 numpy 函数组成 - 因此向它传递输入数组比循环访问输入并为每个输入调用一次要高效得多 - 由于 numpy 的数组，它快 ~2 个数量级。

然而。。。。如果我想只在一个维度上集成我的积分 - 但在另一个维度上有一个输入数组，事情就会下降 - 似乎"Scipy"quadpack 包无法做 numpy 处理数组输入所做的任何事情。 有没有人看到过这个 - 或者找到了修复它的方法 - 或者我误解了它。 我从四边形得到的错误是：

Traceback (most recent call last):
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:Python27libsite-packagesscipyintegratequadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:Python27libsite-packagesscipyintegratequadpack.py", line 312, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
quadpack.error: Supplied function does not return a valid float.

我在下面放了一个我正在尝试做的事情的卡通版本 - 我实际上正在做的事情有一个更复杂的积分，但这是 gyst。

肉在顶部 - 底部正在做基准测试以显示我的观点。

import numpy as np
import time
from scipy.integrate import quad

def Integrand(x, y):
    '''
        Integrand
    '''
    return np.sin(x)*np.sin( y )
def Integrate_x(y):
    '''
        Integrate over x given (y)
    '''
    return quad(Integrand, 0, np.pi/2, args=(y))[0]

def fnIntegrate_x(ystart, yend, nsteps, ArrayInput = False):
    '''
    '''
    yarray = np.arange(ystart,yend, (yend - ystart)/float(nsteps))
    I = np.zeros(nsteps)
    if ArrayInput :
        I = Integrate_x(yarray)
    else :
        for i,y in enumerate(yarray) :
            I[i] = Integrate_x(y)
    return y, I


NCALLS_SET = 1000
NSETS = 10
SETS_t = np.zeros(NSETS)
for i in np.arange(NSETS) :
    XInputs = np.random.rand(NCALLS_SET, 2)
    t0 = time.time()
    for x in XInputs :
        Integrand(x[0], x[1])
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET
print "Benchmarking Integrand - Single Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
Benchmarking Integrand - Single Values:
NCALLS_SET:  1000
NSETS:  10
TimePerCall(s):  1.23999834061e-05 4.06987868647e-06
'''




NCALLS_SET = 1000
NSETS = 10
SETS_t = np.zeros(NSETS)
for i in np.arange(NSETS) :
    XInputs = np.random.rand(NCALLS_SET, 2)
    t0 = time.time()
    Integrand(XInputs[:,0], XInputs[:,1])
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET
print "Benchmarking Integrand - Array Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
Benchmarking Integrand - Array Values:
NCALLS_SET:  1000
NSETS:  10
TimePerCall(s):  2.00009346008e-07 1.26497018465e-07
'''






NCALLS_SET = 1000
NSETS = 100
SETS_t = np.zeros(NSETS)
for i in np.arange(NSETS) :

    t0 = time.time()
    fnIntegrate_x(0, 1, NCALLS_SET, False)
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET
print "Benchmarking fnIntegrate_x - Single Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
NCALLS_SET:  1000
NSETS:  100
TimePerCall(s):  0.000165750000477 8.61204306241e-07
TotalTime:  16.5750000477
'''




NCALLS_SET = 1000
NSETS = 100
SETS_t = np.zeros(NSETS)
for i in np.arange(NSETS) :

    t0 = time.time()
    fnIntegrate_x(0, 1, NCALLS_SET, True)
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET
print "Benchmarking fnIntegrate_x - Array Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
'''
****  Doesn't  work!!!! *****
Traceback (most recent call last):
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:UsersJPDocumentsPythonTestingQuadTestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:Python27libsite-packagesscipyintegratequadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:Python27libsite-packagesscipyintegratequadpack.py", line 312, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
quadpack.error: Supplied function does not return a valid float.
'''

Answer 1

可以通过numpy.vectorize函数实现。我有这个问题很长时间了，然后想到了这个矢量化功能。

你可以像这样使用它：

vectorized_function = numpy.vectorize(your_function)
output = vectorized_function(your_array_input)

Answer 2

怕我在这里用否定来回答我自己的问题。 我认为这是不可能的。 似乎 quad 是用其他东西编写的库的某种端口 - 因此它是内部的库定义事情是如何完成的 - 所以如果不重新设计库本身，就不可能做我想做的事情。

对于在多D集成上遇到时序问题的其他人，我发现最好的方法是使用专用的集成库。 我发现"cuba"似乎有一些非常有效的多D集成例程。

http://www.feynarts.de/cuba/

这些例程是用 c 编写的，所以我最终使用 SWIG 与它们交谈 - 最终为了提高效率，我用 c 重写了我的积分 - 这加快了加载速度......

Answer 3

使用 quadpy（我的一个项目）。它是完全矢量化的，因此可以处理任何形状的数组值函数，并且速度非常快。

Answer 4

我在将所有维度上从 -np.inf 到 np.inf 的概率密度函数集成时遇到了这个问题。

我通过创建一个接受 *args 的包装器函数、将 args 转换为 numpy 数组并集成包装器函数来修复它。

我认为使用 numpy 的矢量化仅集成所有值相等的子空间。

下面是一个示例：

from scipy.integrate import nquad
from scipy.stats import multivariate_normal
mean = [0., 0.]
cov = np.array([[1., 0.],
                [0., 1.]])
bivariate_normal = multivariate_normal(mean=mean, cov=cov)
def pdf(*args):
    x = np.array(args)
    return bivariate_normal.pdf(x)
integration_range = [[-18, 18], [-18, 18]]
nquad(pdf, integration_range)
Output: (1.000000000000001, 1.3429066352690133e-08)