如何調用python scipy quad與數組輸入

Question

我正在使用嵌套的scipy.integrate.quad調用來集成2維積分。被積函數由numpy函數組成 - 所以傳遞給它一個輸入數組要比效率更高，而不是循環輸入併爲每個輸入調用一次 - 由於numpy的數組，它的速度要快2個數量級。

但是，如果我想將我的integrand只整合到一個維度上 - 但是在其他維度上輸入的數據會下降 - 看起來像'scipy'quadpack軟件包無法做到無論numpy如何處理陣列輸入。有沒有其他人看到過這種情況 - 或者找到了解決方法 - 或者我誤解了它。我從四得到的錯誤是：

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. 

我已經把的我嘗試下面做一個卡通版 - 什麼我實際上做了一個比較複雜的積但這是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:\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. 

''' 

Answer 1

可以通過numpy.vectorize函數。我有這個問題很長一段時間，然後來到這個向量化函數。

，你可以使用它像這樣：

vectorized_function = numpy.vectorize（your_function）

輸出= vectorized_function（your_array_input）

Answer 2

害怕我在這裏回答我的問題。我不認爲這是可能的。看起來像quad是某種庫的某種端口，它是用別的方式編寫的 - 因此它是內部的庫，它定義瞭如何完成工作 - 所以如果不重新設計庫本身，可能無法做我想做的事。

對於有多個D集成的時間問題的其他人，我發現最好的方法是使用專用的集成庫。我發現「古巴」似乎有一些非常高效的多維集成程序。

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

這些程序用C語言編寫，所以我結束了使用痛飲與他們交談 - 最終也爲了提高效率重新寫我在C積 - 這加快東西負載....