배열 입력으로 파이썬 scipy quad를 호출하는 방법

Question

중첩 된 scipy.integrate.quad 호출을 사용하여 2 차원 적분을 통합했습니다. integrand는 numpy 함수로 이루어져 있습니다. 따라서 입력 배열을 전달하고 각 입력에 대해 한 번 호출하는 것보다 입력 배열을 전달하는 것이 훨씬 효율적입니다. numpy의 배열로 인해 ~ 2 배 빠른 속도입니다.

그러나 .... 필자가 구현 한 통합 패키지를 하나의 차원에만 통합하려는 경우 - 다른 차원에 대한 입력 배열을 사용하면 "scipy"쿼드 팩 패키지가 수행 할 수없는 것처럼 보입니다. 배열 된 입력을 처리하는 것은 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 함수를 통해 가능합니다. 나는이 문제를 오랫동안 가지고 있었고,이 vectorize 함수를 생각해 냈습니다.

이처럼 사용할 수 있습니다

vectorized_function = numpy.vectorize (your_function)

출력 = vectorized_function (your_array_input)

Answer 2

두려워 여기 내 음성 질문에 대답하고 있습니다. 나는 그것이 가능하다고 생각하지 않는다. 쿼드처럼 보이는 것은 뭔가 다른 것으로 쓰여진 라이브러리의 일종의 항구입니다 - 내부에서 어떻게 진행되는지 정의하는 라이브러리이기 때문에 라이브러리 자체를 다시 디자인하지 않으면 내가 원하는 것을 할 수 없을 것입니다.

여러 D 통합에 타이밍 문제가있는 사용자를 위해 최선의 방법은 전용 통합 라이브러리를 사용하고있는 것으로 나타났습니다. 나는 '쿠바 (cuba)'가 매우 효율적인 멀티 D 통합 루틴을 갖고있는 것으로 보였다.

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

이 루틴은 내가 그들에게 얘기를 꿀꺽 꿀꺽를 사용하여 종료 있도록 C로 작성 - 결국 또한 효율성을 위해 다 내 적분 다시 썼다 - 일 최대 부하를 질주 ....

을