Scipy: Tăng tốc tính toán tích phân phức hợp 2D

Question

Tôi muốn lặp lại phép tính tích phân phức tạp hai chiều bằng cách sử dụng dblquad từ scipy.integrate. Vì số lượng đánh giá sẽ khá cao, tôi muốn tăng tốc độ đánh giá của mã của mình.

Dblquad dường như không thể xử lý các tích phân phức tạp. Do đó, tôi đã chia phần tích phân phức tạp thành một phần thực và một phần tưởng tượng:

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, và lam là các biến defind. Để đánh giá không thể thiếu trên diện tích hình tròn với bán kính ra tôi sử dụng các lệnh sau:

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] 

Theo profiler hai integrands được đánh giá khoảng 35.000 lần. Tổng số tính toán mất khoảng một giây, quá dài đối với ứng dụng tôi có trong đầu.

Tôi là người mới bắt đầu sử dụng máy tính khoa học với Python và Scipy và sẽ hài lòng với các nhận xét chỉ ra cách cải thiện tốc độ đánh giá. Có cách viết lại các lệnh trong các hàm integrand_real và integrand_complex có thể dẫn đến các cải tiến tốc độ đáng ngờ không?

Bạn có nên biên dịch các chức năng đó bằng các công cụ như Cython không? Nếu có: Công cụ nào phù hợp nhất với ứng dụng này?

Answer 1

Bạn có thể đạt được một yếu tố khoảng 10 về tốc độ bằng cách sử dụng Cython, xem dưới đây:

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 

Đây có lẽ là không đủ, và tin xấu là điều này có lẽ là khá đóng (có thể là hệ số tối đa là 2) với mọi tốc độ C/Fortran --- nếu sử dụng cùng một thuật toán để tích hợp thích ứng. (scipy.integrate.quad chính nó đã có ở Fortran.)

Để tiếp tục, bạn cần cân nhắc các cách khác nhau để thực hiện tích hợp . Điều này đòi hỏi một số suy nghĩ --- không thể cung cấp nhiều từ đỉnh đầu của tôi bây giờ.

Ngoài ra, bạn có thể giảm dung sai cho đến khi số tích phân được đánh giá.

 
# 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

Bạn đã xem xét đa xử lý (đa luồng) chưa? Dường như bạn không cần phải làm một sự tích hợp cuối cùng (trên toàn bộ tập hợp) để xử lý song song đơn giản có thể là câu trả lời. Ngay cả khi bạn đã có để tích hợp, bạn có thể chờ đợi cho chạy chủ đề để kết thúc tính toán trước khi thực hiện hội nhập cuối cùng. Đó là, bạn có thể chặn các chủ đề chính cho đến khi tất cả các công nhân đã hoàn thành.

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