Używam zagnieżdżonego scipy.integrate.Quad połączenia, aby zintegrować integransę 2 wymiarową. Integransa jest wykonana z funkcji NUMPY - więc jest o wiele bardziej wydajny, aby przekazać go tablicę wejść - niż do pętli przez wejścia i nazwać go raz na każdy - jest ~ 2 rzędami wielkości szybciej ze względu na macierzy numpy.

Jednakże .... Jeśli chcę zintegrować moją integransę nad jednym wymiarem - ale z tablicą wejść przez inne wymiar rzeczy upadają - wydaje się, że pakiet Quadpack "Scipy" nie jest w stanie robić tego, co to jest Numpy robi, aby poradzić sobie z tabelami wejściowymi. Czy ktoś inny to widział - i znalazł sposób naprawienia - czy jestem nieporozumieniem. Błąd, który dostaję z quadem to:

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.

Umieściłem wersję kreskówki tego, co próbuję zrobić poniżej - co faktycznie robię, ma bardziej skomplikowaną integransę, ale jest to GYST.

Mięso jest na górze - dno robi benchmarking, aby pokazać mój punkt.

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.

'''
3
JPH 10 sierpień 2012, 00:00

3 odpowiedzi

Najlepsza odpowiedź

Jest to możliwe dzięki funkcji numpy.Vectorize. Miałem ten problem przez długi czas, a potem podszedł do tej funkcji Vectorize.

Możesz go użyć w ten sposób:

vectorized_function = numpy.vectorize(your_function)

output = vectorized_function(your_array_input)
1
Gabriel 14 listopad 2018, 20:11

Boi się, że odpowiadam na moje własne pytanie z negatywnym tutaj. Nie sądzę, że to możliwe. Wygląda na to, że Quad jest jakiś rodzaj portu biblioteki napisanej w czymś innym - jako taka jest biblioteka w środku, która określa, jak się robi rzeczy - więc prawdopodobnie nie jest możliwe, aby zrobić to, czego chciałem bez przeprojektowania samego biblioteki.

Dla innych osób z kwestiami czasowymi w wielu integracji D, znalazłem najlepszy sposób korzystania z dedykowanej biblioteki integracji. Odkryłem, że "Kuba" wydawała się mieć pewne dość wydajne procedury integracji Multi D.

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

Procedury te są napisane w C, więc skończyło się za pomocą SWIG, aby porozmawiać z nimi - i ostatecznie również do wydajności ponownie napisał moją integransę w C - które przyspieszyły ładunki ....

2
JPH 3 marzec 2013, 18:02

Użyj Quadpy (projekt mojego). Jest on w pełni wektoryzowany, więc może obsługiwać funkcje wyceniane w tabliczce dowolnego kształtu i robi to bardzo szybko.

1
Nico Schlömer 12 listopad 2019, 15:46