카테고리 보관물: scicomp

scicomp

Matlab의 통합이 Scipy에서 통합 성능을 능가하는 이유는 무엇입니까? 상당히 유사합니다). 파이썬 import numpy as np from

matlab이 숫자 통합과 Scipy를 처리하는 방식에 약간의 좌절이 있습니다. 아래 테스트 코드에서 다음과 같은 차이점이 있습니다.

  1. Matlab의 버전은 내 파이썬 버전 보다 평균 24 배 빠릅니다 !
  2. Matlab의 버전은 경고없이 정수를 계산할 수 있으며 파이썬은 반환합니다 nan+nanj

언급 된 두 가지 점과 관련하여 파이썬에서 동일한 성능을 얻으려면 어떻게해야합니까? 문서에 따르면 두 방법 모두 적분을 근사하기 위해 “전역 적응 구적법”을 사용해야합니다.

아래는 두 버전의 코드입니다 (파이썬은 복잡한 정수를 처리 할 수 ​​있도록 정수 함수를 작성해야하지만 상당히 유사합니다).

파이썬

import numpy as np
from scipy import integrate
import time

def integral(integrand, a, b,  arg):
    def real_func(x,arg):
        return np.real(integrand(x,arg))
    def imag_func(x,arg):
        return np.imag(integrand(x,arg))
    real_integral = integrate.quad(real_func, a, b, args=(arg))
    imag_integral = integrate.quad(imag_func, a, b, args=(arg))
    return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)


def f_integrand(s, omega):
    sigma = np.pi/(np.pi+2)
    xs = np.exp(-np.pi*s/(2*sigma))
    x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+np.sqrt(1-xs**2))
    x2 = 1-2*sigma/np.pi*(1-xs)
    zeta = x2+x1*1j
    Vc = 1/(2*sigma)
    theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
    t1 = 1/np.sqrt(1+np.tan(theta)**2)
    t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
    return np.real((t1-1j*t2)/np.sqrt(zeta**2-1))*np.exp(1j*omega*s/Vc);

t0 = time.time()
omega = 10
result = integral(f_integrand, 0, np.inf, omega)
print time.time()-t0
print result

MATLAB

function [ out ] = f_integrand( s, omega )
    sigma = pi/(pi+2);
    xs = exp(-pi.*s./(2*sigma));
    x1 = -2*sigma./pi.*(log(xs./(1+sqrt(1-xs.^2)))+sqrt(1-xs.^2));
    x2 = 1-2*sigma./pi.*(1-xs);
    zeta = x2+x1*1j;
    Vc = 1/(2*sigma);
    theta =  -1*asin(exp(-pi./(2.0.*sigma).*s));
    t1 = 1./sqrt(1+tan(theta).^2);
    t2 = -1./sqrt(1+1./tan(theta).^2);
    out = real((t1-1j.*t2)./sqrt(zeta.^2-1)).*exp(1j.*omega.*s./Vc);
end

t=cputime;
omega = 10;
result = integral(@(s) f_integrand(s,omega),0,Inf)
time_taken = cputime-t

답변

이 질문에는 두 가지 매우 다른 하위 질문이 있습니다. 첫 번째 문제 만 다룰 것입니다.

Matlab의 버전은 내 파이썬 버전 보다 평균 24 배 빠릅니다 !

두 번째는 주관적입니다. 나는 사용자에게 적분에 문제가 있음을 알리는 것이 좋은 것이며이 SciPy 동작은 Matlab보다 성능이 뛰어나며 침묵을 유지하고 어떻게 든 Matlab 엔지니어들만 알고있는 방식으로 내부적으로 처리하려고 시도한다고 말합니다. 최고라고 결정했습니다.

NaN waring을 피하기 위해 통합 범위를 0 에서 30으로 ( 0 에서 np.inf 대신 ) 변경하고 JIT 컴파일을 추가했습니다. 솔루션을 벤치마킹하기 위해 통합을 300 회 반복했습니다. 결과는 랩톱에서 나왔습니다.

JIT 컴파일이없는 경우 :

$ ./test_integrate.py
34.20992112159729
(0.2618828053067563+0.24474506983644717j)

JIT 컴파일

$ ./test_integrate.py
0.8560323715209961
(0.261882805306756+0.24474506983644712j)

이 방법으로 두 줄의 코드를 추가 하면 비 JIT 버전에 비해 약 40 배 의 Python 코드 속도가 향상 됩니다. 더 나은 비교를 제공 할 랩톱에는 Matlab이 없지만 24/40 = 0.6보다 PC에 잘 맞으면 JIT가있는 Python은 Matlab보다 두 배 빠릅니다 .이 특정 사용자 알고리즘. 전체 코드 :

#!/usr/bin/env python3
import numpy as np
from scipy import integrate
from numba import complex128,float64,jit
import time

def integral(integrand, a, b,  arg):
    def real_func(x,arg):
        return np.real(integrand(x,arg))
    def imag_func(x,arg):
        return np.imag(integrand(x,arg))
    real_integral = integrate.quad(real_func, a, b, args=(arg))
    imag_integral = integrate.quad(imag_func, a, b, args=(arg))
    return real_integral[0] + 1j*imag_integral[0]

vintegral = np.vectorize(integral)


@jit(complex128(float64, float64), nopython=True, cache=True)
def f_integrand(s, omega):
    sigma = np.pi/(np.pi+2)
    xs = np.exp(-np.pi*s/(2*sigma))
    x1 = -2*sigma/np.pi*(np.log(xs/(1+np.sqrt(1-xs**2)))+np.sqrt(1-xs**2))
    x2 = 1-2*sigma/np.pi*(1-xs)
    zeta = x2+x1*1j
    Vc = 1/(2*sigma)
    theta =  -1*np.arcsin(np.exp(-np.pi/(2.0*sigma)*s))
    t1 = 1/np.sqrt(1+np.tan(theta)**2)
    t2 = -1/np.sqrt(1+1/np.tan(theta)**2)
    return np.real((t1-1j*t2)/np.sqrt(zeta**2-1))*np.exp(1j*omega*s/Vc);

t0 = time.time()
omega = 10
for i in range(300):
    #result = integral(f_integrand, 0, np.inf, omega)
    result = integral(f_integrand, 0, 30, omega)
print (time.time()-t0)
print (result)

PC의 차이점을 보려면 @jit 줄을 주석으로 처리하십시오.

답변

때때로 통합 기능을 JIT 할 수 없습니다. 이 경우 다른 통합 방법을 사용하는 것이 해결책입니다.

나는 추천한다 scipy.integrate.romberg (ref) .
romberg복잡한 함수를 통합하고 함수를 배열로 평가할 수 있습니다.

답변