如何使用数组输入调用Python Scipy Quad函数

Question

我正在使用嵌套的scipy.integrate.quad调用来积分一个二维被积函数。该被积函数由numpy函数组成 - 因此将其传递给一个输入数组比循环遍历每个输入并调用一次更有效率 - 由于numpy的数组，速度提高了大约两个数量级。

然而……如果我想在另一个维度上积分我的被积函数 - 但是使用另一个维度上的输入数组时，事情就会出问题 - 似乎'scipy' quadpack包无法处理数组输入所做的任何操作。是否有其他人看到过这种情况 - 或者找到了解决方法 - 还是我误解了它。我从quad获得的错误是：

Traceback (most recent call last):
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.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.

我在下面放了一个卡通图，展示了我想做的事情 - 实际上我所做的有更复杂的被积函数，但这是要点。

关键在于顶部，底部正在进行基准测试以显示我的观点。

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:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.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.

'''

- JPH

Answer 1

可以通过numpy.vectorize函数实现。我长时间以来一直遇到这个问题，后来想到了这个vectorize函数。

你可以像这样使用它：

vectorized_function = numpy.vectorize(your_function)

output = vectorized_function(your_array_input)

- Me Man

---

很酷，我会在某个时候尝试这个 - 已经有一段时间没有看过这段代码了 - 但将来可能会有用。 - JPH

Answer 2

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

- Nico Schlömer

---

谢谢！很好知道。虽然我已经有一段时间没看这个了，但希望对其他人有用。 - JPH

Answer 3

恐怕我要用否定的方式回答自己的问题。我认为这是不可能的。似乎quad是某种其他语言库的移植，因此内部的库定义了如何完成任务，所以想要做到我想要的事情可能需要重新设计库本身。

对于其他在多次D积分中遇到时间问题的人，我发现使用专门的积分库是最好的方法。我发现“古巴”似乎有一些非常高效的多次D积分例程。

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

这些例程是用C语言编写的，所以我最终使用了SWIG与它们进行通信——最后我还将我的被积函数重写成了C语言，从而提高了效率，速度得到了大幅度的提升...

- JPH

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)

- John

---

你被禁止评论是有原因的，这并不意味着你可以在答案部分发表评论。请删除此评论。 - Rob

我找到了他的问题的解决方案，所以我不会删除它。 - John