I want to repeatedly calculate a two-dimensional complex integral using dblquad from scipy.integrate. As the number of evaluations will be quite high I would like to increase the evaluation speed of my code.
Dblquad does not seem to be able to handle complex integrands. Thus, I have split the complex integrand into a real and an imaginary part:
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, and lam are variables defind in advance. To evaluate the integral over the area of a circle with radius ra I use the following commands:
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]
According to the profiler the two integrands are evaluated about 35000 times. The total calculation takes about one second, which is too long for the application I have in mind.
I am a beginner to scientific computing with Python and Scipy and would be happy about comments that point out ways of improving the evaluation speed. Are there ways of rewriting the commands in the integrand_real and integrand_complex functions that could lead to siginficant speed improvements?
Would it make sense to compile those functions using tools like Cython? If yes: Which tool would best fit this application?
You can gain a factor of about 10 in speed by using Cython, see below:
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
This is probably not enough, and the bad news is that this is probably quite close (maybe factor of 2 at most) to everything-in-C/Fortran speed --- if using the same algorithm for adaptive integration. (scipy.integrate.quad itself is already in Fortran.)
To get further, you'd need to consider different ways to do the integration. This requires some thinking --- can't offer much from the top of my head now.
Alternatively, you can reduce the tolerance up to which the integral is evaluated.
# 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