Scipy: Beschleunigung der Berechnung eines komplexen 2D-Integrals

Question

Ich möchte wiederholt ein zweidimensionales komplexes Integral mit dblquad von scipy.integrate berechnen. Da die Anzahl der Auswertungen sehr hoch sein wird, möchte ich die Auswertungsgeschwindigkeit meines Codes erhöhen.

Dblquad scheint nicht in der Lage zu sein, komplexe Integrale zu behandeln. Somit habe ich den Komplex Integrand in eine reale aufzuspalten und Imaginärteil:

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 und lam sind Variablen defind im Voraus. Um das Integral über die Fläche eines Kreises mit dem Radius ra bewerten ich die folgenden Befehle verwenden:

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] 

Nach dem Profiler die beiden Integra werden etwa 35000 mal bewertet. Die Gesamtberechnung dauert etwa eine Sekunde, was für die Anwendung, die mir in den Sinn kommt, zu lang ist.

Ich bin ein Anfänger zum wissenschaftlichen Rechnen mit Python und Scipy und würde mich über Kommentare freuen, die Wege zur Verbesserung der Auswertungsgeschwindigkeit aufzeigen. Gibt es Möglichkeiten, die Befehle in den Funktionen integrand_real und integrand_complex neu zu schreiben, was zu signifikanten Geschwindigkeitsverbesserungen führen könnte?

Wäre es sinnvoll, diese Funktionen mit Tools wie Cython zu kompilieren? Wenn ja: Welches Tool passt am besten zu dieser Anwendung?

Answer 1

Sie können einen Faktor von etwa 10 in der Geschwindigkeit gewinnen, indem Cython verwenden, siehe unten:

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 

Dies ist wahrscheinlich nicht genug ist, und die schlechte Nachricht ist, dass dies wahrscheinlich ist ganz schließen Sie (vielleicht höchstens Faktor 2) an alles-in-C/Fortran-Geschwindigkeit --- wenn Sie denselben Algorithmus für adaptive Integration verwenden. (scipy.integrate.quad selbst ist bereits in Fortran.)

Um weiter zu kommen, müssten Sie verschiedene Möglichkeiten, die Integration zu tun. Dies erfordert einiges Denken --- kann nicht viel von die Spitze meines Kopfes jetzt bieten.

Alternativ können Sie die Toleranz reduzieren, bis zu der das Integral ausgewertet wird.

 
# 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

Haben Sie Multiprozessing (Multithreading) in Betracht gezogen? Es scheint, dass Sie keine endgültige Integration (über den gesamten Satz) benötigen, so dass die einfache parallele Verarbeitung die Antwort sein könnte. Selbst wenn Sie integrieren mussten, können Sie darauf warten, dass Threads ausgeführt werden, um die Berechnung abzuschließen, bevor Sie die endgültige Integration durchführen. Das heißt, Sie können den Hauptthread blockieren, bis alle Worker fertig sind.

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