11

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

Dblquad は複雑な被積分関数を処理できないようです。したがって、複素被積分関数を実部と虚部に分割しました。

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 はあらかじめ定義された変数です。半径 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 などのツールを使用してこれらの関数をコンパイルすることは理にかなっていますか? はいの場合: このアプリケーションに最適なツールはどれですか?

4

3 に答える 3

14

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 に含まれています。)

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

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

# Python で行う
#
# >>> import pyximport; pyximport.install(reload_support=True)
# >>> cythonmodule をインポート

cimport numpy を np として
cimport cython

"complex.h" からの cdef extern:
    double complex csqrt(double complex z) nogil
    double complex cexp(double complex z) nogil
    double creal(ダブル コンプレックス z) ノーギル
    double cimag(double complex z) ノーギル

libc.math cimport sqrt から

scipy.integrateインポートdblquadから

cdef クラス パラメータ:
    cdef public double lam、y0、k、zxp、z、ra

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

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

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

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

def doit(lam, y0, k, zxp, z, ra):
    p = Params(ラム=ラム, y0=y0, k=k, zxp=zxp, z=z, ra=ra)
    rr, err = dblquad(integrand_real, -ra, ra, ラムダ x: -ymax(x, p), ラムダ x: ymax(x, p), args=(p,))
    ri, err = dblquad(integrand_imag, -ra, ra, ラムダ x: -ymax(x, p), ラムダ x: ymax(x, p), args=(p,))
    rr + 1j*ri を返す
于 2013-04-03T17:37:48.530 に答える
4

マルチプロセッシング (マルチスレッド) を検討したことがありますか? （セット全体にわたって）最終的な統合を行う必要がないように思われるため、単純な並列処理が答えになる可能性があります。統合する必要があったとしても、最終的な統合を行う前に、実行中のスレッドが計算を終了するのを待つことができます。つまり、すべてのワーカーが完了するまでメイン スレッドをブロックできます。

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

于 2013-04-03T16:47:40.783 に答える