Scipy：2次元複素数積分の計算を高速化する

Question

scipy.integrateからdblquadを使用して2次元複素積分を繰り返し計算したいとします。評価の数が非常に多いので、私のコードの評価速度を上げたいと思います。

Dblquadは複雑な被積分関数を扱うことができないようです。したがって、Iは実部と虚部に複雑な積分を分割した：

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、およびLAM予めdefind変数です。私は、次のコマンドを使用半径raの円の面積の積分を評価するために：

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] 

プロファイラによる2つの被積分を、約35000回評価されます。合計計算には約1秒かかりますが、これは私が気にしているアプリケーションにとっては長すぎます。

私はPythonとScipyを使った科学計算の初心者です。評価スピードを向上させる方法を指摘するコメントに満足しています。 integrand_real関数とintegrand_complex関数のコマンドを書き直す方法はありますか？

Cythonのようなツールを使ってこれらの関数をコンパイルするのは意味がありますか？はいの場合：このアプリケーションに最も適したツールはどれですか？

Answer 1

あなたは以下を参照してください、Cythonを使用することにより、約10速度の要因を得ることができます。

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 

これはおそらく十分ではありません、と悪いニュースは、これはおそらくかなり であるということです適応インテグレーションに同じアルゴリズムを使用している場合は、C/Fortranのすべての速度に近い値（最大で2倍）を返します。 （scipy.integrate.quad はすでにFortranに入っています）

また、 の統合を行うには、さまざまな方法を検討する必要があります。これはいくつかの考えを必要とします---私の頭の上の から多くを提供することはできません。

または、 が評価されるまでの許容差を減らすことができます。

 
# 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

マルチプロセッシング（マルチスレッド）を考えましたか？単純な並列処理が答えかもしれないので、あなたは最終的な統合を行う必要はないようです。統合する必要があったとしても、最終的な統合を行う前に、スレッドの実行が完了するのを待つことができます。つまり、すべてのワーカーが完了するまでメインスレッドをブロックすることができます。

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