python - simps - sympy integrate
Scipy: acelerar el cálculo de una integral compleja en 2D (2)
Quiero calcular repetidamente una integral compleja de dos dimensiones usando dblquad de scipy.integrate. Como el número de evaluaciones será bastante alto, me gustaría aumentar la velocidad de evaluación de mi código.
Dblquad no parece ser capaz de manejar integrandos complejos. Por lo tanto, he dividido el integrando complejo en una parte real y una parte imaginaria:
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, y lam son variables definidas de antemano. Para evaluar la integral sobre el área de un círculo con radio ra, uso los siguientes comandos:
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]
Según el perfilador, los dos integrandos se evalúan unas 35000 veces. El cálculo total tarda aproximadamente un segundo, que es demasiado largo para la aplicación que tengo en mente.
Soy un principiante de la informática científica con Python y Scipy y estaría contento con los comentarios que señalan formas de mejorar la velocidad de evaluación. ¿Hay formas de volver a escribir los comandos en las funciones integrand_real e integrand_complex que podrían conducir a mejoras de velocidad significativas?
¿Tendría sentido compilar esas funciones usando herramientas como Cython? En caso afirmativo, ¿qué herramienta se ajusta mejor a esta aplicación?
¿Has considerado el multiprocesamiento (multihilo)? Parece que no es necesario realizar una integración final (en todo el conjunto), por lo que la respuesta podría ser un procesamiento en paralelo simple. Incluso si tuviera que integrar, puede esperar a que se ejecuten los hilos para finalizar el cálculo antes de realizar la integración final. Es decir, puede bloquear el hilo principal hasta que todos los trabajadores hayan terminado.
Puedes ganar un factor de aproximadamente 10 en velocidad usando Cython, mira abajo:
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
Esto probablemente no sea suficiente, y la mala noticia es que probablemente sea bastante cercano (tal vez un factor de 2 como máximo) a todo-en-C / velocidad de Fortran --- si se usa el mismo algoritmo para la integración adaptativa. (scipy.integrate.quad ya está en Fortran).
Para llegar más lejos, necesitaría considerar diferentes formas de realizar la integración. Esto requiere pensar un poco, no puedo ofrecer mucho desde lo más alto de mi cabeza ahora.
Alternativamente, puede reducir la tolerancia hasta la cual se evalúa la integral.
# 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