Scipy: Przyspieszanie obliczania całki zespolonej 2D

Question

Chcę wielokrotnie obliczyć dwuwymiarową całkę zespoloną za pomocą dblquad z scipy.integrate. Ponieważ liczba ocen będzie dość wysoka, chciałbym zwiększyć szybkość oceny mojego kodu.

Dblquad wydaje się nie być w stanie obsłużyć złożonych elementów integrujących. Tak więc, mam podzielić złożony całki w prawdziwym i wyimaginowanym części:

def integrand_real(x, y): 
    R1=sqrt(x**2 + (y-y0)**2 + z**2) 
    R2=sqrt(x**2 + y**2 + zxp**2) 
    return real(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1)) 

def integrand_imag(x,y): 
    R1=sqrt(x**2 + (y-y0)**2 + z**2) 
    R2=sqrt(x**2 + y**2 + zxp**2) 
    return imag(exp(1j*k*(R1-R2)) * (-1j*z/lam/R2/R1**2) * (1+1j/k/R1)) 

y0, z, ZXP, k, i lam są zmiennymi defind wyprzedzeniem. Aby ocenić całki na obszarze okręgu o promieniu ra używam następujące polecenia:

from __future__ import division 
from scipy.integrate import dblquad 
from pylab import * 

def ymax(x): 
    return sqrt(ra**2-x**2) 

lam = 0.000532 
zxp = 5. 
z = 4.94 
k = 2*pi/lam 
ra = 1.0 

res_real = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x)) 
res_imag = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x), lambda x: ymax(x)) 
res = res_real[0]+ 1j*res_imag[0] 

Według profilera dwa podcałkowych oceniane są około 35.000 razy. Całkowite obliczenie zajmuje około jednej sekundy, co jest zbyt długie dla aplikacji, którą mam na myśli.

Jestem początkującym informatykiem naukowym z Python i Scipy i byłbym szczęśliwy z komentarzy wskazujących sposoby poprawy szybkości oceny. Czy istnieją sposoby przepisywania poleceń w funkcjach integrand_real i integrand_complex, które mogą prowadzić do znacznych ulepszeń prędkości?

Czy ma sens kompilowanie tych funkcji przy użyciu narzędzi takich jak Cython? Jeśli tak: które narzędzie najlepiej pasuje do tej aplikacji?

Answer 1

można uzyskać współczynnik około 10 prędkości za pomocą Cython, patrz poniżej:

In [87]: %timeit cythonmodule.doit(lam=lam, y0=y0, zxp=zxp, z=z, k=k, ra=ra) 
1 loops, best of 3: 501 ms per loop 
In [85]: %timeit doit() 
1 loops, best of 3: 4.97 s per loop 

Jest to prawdopodobnie nie wystarczy, a najgorsze jest to, że jest to prawdopodobnie dość zamknij (być może współczynnik 2) do wszystkiego-w-C/prędkości Fortranu --- jeśli używasz tego samego algorytmu do integracji adaptacyjnej. (scipy.integrate.quad sam jest już w Fortranie.)

Aby uzyskać więcej informacji, należy rozważyć różne sposoby integracji z . To wymaga trochę myślenia - nie mogę teraz zaoferować większej wartości od .

Alternatywnie można zmniejszyć tolerancję, do której oszacowano całkę .

 
# Do in Python 
# 
# >>> import pyximport; pyximport.install(reload_support=True) 
# >>> import cythonmodule 

cimport numpy as np 
cimport cython 

cdef extern from "complex.h": 
    double complex csqrt(double complex z) nogil 
    double complex cexp(double complex z) nogil 
    double creal(double complex z) nogil 
    double cimag(double complex z) nogil 

from libc.math cimport sqrt 

from scipy.integrate import dblquad 

cdef class Params: 
    cdef public double lam, y0, k, zxp, z, ra 

    def __init__(self, lam, y0, k, zxp, z, ra): 
     self.lam = lam 
     self.y0 = y0 
     self.k = k 
     self.zxp = zxp 
     self.z = z 
     self.ra = ra 

@cython.cdivision(True) 
def integrand_real(double x, double y, Params p): 
    R1 = sqrt(x**2 + (y-p.y0)**2 + p.z**2) 
    R2 = sqrt(x**2 + y**2 + p.zxp**2) 
    return creal(cexp(1j*p.k*(R1-R2)) * (-1j*p.z/p.lam/R2/R1**2) * (1+1j/p.k/R1)) 

@cython.cdivision(True) 
def integrand_imag(double x, double y, Params p): 
    R1 = sqrt(x**2 + (y-p.y0)**2 + p.z**2) 
    R2 = sqrt(x**2 + y**2 + p.zxp**2) 
    return cimag(cexp(1j*p.k*(R1-R2)) * (-1j*p.z/p.lam/R2/R1**2) * (1+1j/p.k/R1)) 

def ymax(double x, Params p): 
    return sqrt(p.ra**2 + x**2) 

def doit(lam, y0, k, zxp, z, ra): 
    p = Params(lam=lam, y0=y0, k=k, zxp=zxp, z=z, ra=ra) 
    rr, err = dblquad(integrand_real, -ra, ra, lambda x: -ymax(x, p), lambda x: ymax(x, p), args=(p,)) 
    ri, err = dblquad(integrand_imag, -ra, ra, lambda x: -ymax(x, p), lambda x: ymax(x, p), args=(p,)) 
    return rr + 1j*ri 

Answer 2

Czy brałeś pod uwagę proces wieloprocesowy (wielowątkowość)? Wygląda na to, że nie musisz wykonywać końcowej integracji (w całym zestawie), więc prostym przetwarzaniem równoległym może być odpowiedź. Nawet jeśli musiałeś się zintegrować, możesz poczekać na uruchomienie wątków, aby zakończyć obliczenia przed ostateczną integracją. Oznacza to, że możesz zablokować główny wątek, aż wszyscy pracownicy zakończą pracę.

http://docs.python.org/2/library/multiprocessing.html