코드의 Cython 최적화

저는 Cython으로 파이썬 입자 추적 코드의 성능을 높이기 위해 고심하고 있습니다.코드의 Cython 최적화

분명히

from scipy.integrate import odeint 
import numpy as np 
from numpy import sqrt, pi, sin, cos 
from time import time as Time 
import multiprocessing as mp 
from functools import partial 

cLight = 299792458. 
Dim = 6 

class Integrator: 
    def __init__(self, ring): 
     self.ring = ring 

    def equations(self, X, s): 
     dXds = np.zeros(Dim) 

     E, B = self.ring.getEMField([X[0], X[2], s], X[4]) 

     h = 1 + X[0]/self.ring.ringRadius 
     p_s = np.sqrt(X[5]**2 - self.ring.particle.mass**2 - X[1]**2 - X[3]**2) 
     dtds = h*X[5]/p_s 
     gamma = X[5]/self.ring.particle.mass 
     beta = np.array([X[1], X[3], p_s])/X[5] 

     dXds[0] = dtds*beta[0] 
     dXds[2] = dtds*beta[1] 
     dXds[1] = p_s/self.ring.ringRadius + self.ring.particle.charge*(dtds*E[0] + dXds[2]*B[2] - h*B[1]) 
     dXds[3] = self.ring.particle.charge*(dtds*E[1] + h*B[0] - dXds[0]*B[2]) 
     dXds[4] = dtds 
     dXds[5] = self.ring.particle.charge*(dXds[0]*E[0] + dXds[2]*E[1] + h*E[2]) 
     return dXds 

    def odeSolve(self, X0, sRange): 
     sol = odeint(self.equations, X0, sRange) 
     return sol 

class Ring: 
    def __init__(self, particle): 
     self.particle = particle 
     self.ringRadius = 7.112 
     self.magicB0 = self.particle.magicMomentum/self.ringRadius 

    def getEMField(self, pos, time): 
     x, y, s = pos 
     theta = (s/self.ringRadius*180/pi) % 360 
     r = sqrt(x**2 + y**2) 
     arg = 0 if r == 0 else np.angle(complex(x/r, y/r)) 
     rn = r/0.045 

     k2 = 37*24e3 
     k10 = -4*24e3 

     E = np.zeros(3) 
     B = np.array([ 0, self.magicB0, 0 ]) 

     for i in range(4): 
      if ((21.9+90*i < theta < 34.9+90*i or 38.9+90*i < theta < 64.9+90*i) and (-0.05 < x < 0.05 and -0.05 < y < 0.05)): 
       E = np.array([ k2*x/0.045 + k10*rn**9*cos(9*arg), -k2*y/0.045 -k10*rn**9*sin(9*arg), 0]) 
       break 
     return E, B 

class Particle: 
    def __init__(self): 
     self.mass = 105.65837e6 
     self.charge = 1. 
     self.gm2 = 0.001165921 

     self.magicMomentum = self.mass/sqrt(self.gm2) 
     self.magicEnergy = sqrt(self.magicMomentum**2 + self.mass**2) 
     self.magicGamma = self.magicEnergy/self.mass 
     self.magicBeta = self.magicMomentum/(self.magicGamma*self.mass) 


def runSimulation(nParticles, tEnd): 
    particle = Particle() 
    ring = Ring(particle) 
    integrator = Integrator(ring) 

    Xs = np.array([ np.array([45e-3*(np.random.rand()-0.5)*2, 0, 0, 0, 0, particle.magicEnergy]) for i in range(nParticles) ]) 
    sRange = np.arange(0, tEnd, 1e-9)*particle.magicBeta*cLight 

    ode = partial(integrator.odeSolve, sRange=sRange) 

    t1 = Time() 

    pool = mp.Pool() 
    sol = np.array(pool.map(ode, Xs)) 

    t2 = Time() 
    print ("%.3f sec" %(t2-t1)) 

    return t2-t1

, 가장 시간이 많이 걸리는 과정이 odeSolve() 및 클래스 통합의 방정식()로 정의 ODE를 통합한다 :

은 여기 내 순수 파이썬 코드입니다. 또한 Ring 클래스의 getEMField() 메소드는 해결 과정에서 equations() 메소드만큼 호출됩니다. 나는 사이 썬를 사용하여 최대 속도 (~ 20 배 이상 10 배)의 상당한 금액을 얻기 위해 노력했지만 난 단지 다음과 같은 사이 썬 스크립트에 의해 속도 업 ~ 1.5 배 수준 가지고 :

import cython 
import numpy as np 
cimport numpy as np 
from libc.math cimport sqrt, pi, sin, cos 

from scipy.integrate import odeint 
from time import time as Time 
import multiprocessing as mp 
from functools import partial 

cdef double cLight = 299792458. 
cdef int Dim = 6 

@cython.boundscheck(False) 
cdef class Integrator: 
    cdef Ring ring 

    def __init__(self, ring): 
     self.ring = ring 

    cpdef np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] equations(self, 
        np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] X, 
        double s): 
     cdef np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] dXds = np.zeros(Dim) 
     cdef double h, p_s, dtds, gamma 
     cdef np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] beta, E, B 

     E, B = self.ring.getEMField([X[0], X[2], s], X[4]) 

     h = 1 + X[0]/self.ring.ringRadius 
     p_s = np.sqrt(X[5]*X[5] - self.ring.particle.mass*self.ring.particle.mass - X[1]*X[1] - X[3]*X[3]) 
     dtds = h*X[5]/p_s 
     gamma = X[5]/self.ring.particle.mass 
     beta = np.array([X[1], X[3], p_s])/X[5] 

     dXds[0] = dtds*beta[0] 
     dXds[2] = dtds*beta[1] 
     dXds[1] = p_s/self.ring.ringRadius + self.ring.particle.charge*(dtds*E[0] + dXds[2]*B[2] - h*B[1]) 
     dXds[3] = self.ring.particle.charge*(dtds*E[1] + h*B[0] - dXds[0]*B[2]) 
     dXds[4] = dtds 
     dXds[5] = self.ring.particle.charge*(dXds[0]*E[0] + dXds[2]*E[1] + h*E[2]) 
     return dXds 

    cpdef np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] odeSolve(self, 
       np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] X0, 
       np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] sRange): 
     sol = odeint(self.equations, X0, sRange) 
     return sol 

@cython.boundscheck(False) 
cdef class Ring: 
    cdef Particle particle 
    cdef double ringRadius 
    cdef double magicB0 

    def __init__(self, particle): 
     self.particle = particle 
     self.ringRadius = 7.112 
     self.magicB0 = self.particle.magicMomentum/self.ringRadius 

    cpdef tuple getEMField(self, 
        list pos, 
        double time): 
     cdef double x, y, s 
     cdef double theta, r, rn, arg, k2, k10 
     cdef np.ndarray[np.double_t, ndim=1, negative_indices=False, mode="c"] E, B 

     x, y, s = pos 
     theta = (s/self.ringRadius*180/pi) % 360 
     r = sqrt(x*x + y*y) 
     arg = 0 if r == 0 else np.angle(complex(x/r, y/r)) 
     rn = r/0.045 

     k2 = 37*24e3 
     k10 = -4*24e3 

     E = np.zeros(3) 
     B = np.array([ 0, self.magicB0, 0 ]) 

     for i in range(4): 
      if ((21.9+90*i < theta < 34.9+90*i or 38.9+90*i < theta < 64.9+90*i) and (-0.05 < x < 0.05 and -0.05 < y < 0.05)): 
       E = np.array([ k2*x/0.045 + k10*rn**9*cos(9*arg), -k2*y/0.045 -k10*rn**9*sin(9*arg), 0]) 
       #E = np.array([ k2*x/0.045, -k2*y/0.045, 0]) 
       break 
     return E, B 

cdef class Particle: 
    cdef double mass 
    cdef double charge 
    cdef double gm2 

    cdef double magicMomentum 
    cdef double magicEnergy 
    cdef double magicGamma 
    cdef double magicBeta 

    def __init__(self): 
     self.mass = 105.65837e6 
     self.charge = 1. 
     self.gm2 = 0.001165921 

     self.magicMomentum = self.mass/sqrt(self.gm2) 
     self.magicEnergy = sqrt(self.magicMomentum**2 + self.mass**2) 
     self.magicGamma = self.magicEnergy/self.mass 
     self.magicBeta = self.magicMomentum/(self.magicGamma*self.mass) 

def runSimulation(nParticles, tEnd): 
    particle = Particle() 
    ring = Ring(particle) 
    integrator = Integrator(ring) 

    #nParticles = 5 
    Xs = np.array([ np.array([45e-3*(np.random.rand()-0.5)*2, 0, 0, 0, 0, particle.magicEnergy]) for i in range(nParticles) ]) 
    sRange = np.arange(0, tEnd, 1e-9)*particle.magicBeta*cLight 

    ode = partial(integrator.odeSolve, sRange=sRange) 

    t1 = Time() 

    pool = mp.Pool() 
    sol = np.array(pool.map(ode, Xs)) 

    t2 = Time() 
    print ("%.3f sec" %(t2-t1)) 

    return t2-t1

내가 얻기 위해 무엇을해야을 Cython의 최대 효과? (나는 Cython 대신 Numba를 시도했는데 실제로 Numba의 성능 향상은 엄청났습니다 (~ 20x 속도 향상).하지만 Numba를 Python 클래스 인스턴스와 함께 사용하는 것은 매우 어려웠고 Numba 대신 Cython을 사용하기로 결정했습니다.

출처

2017-12-15 hopeso

병목 현상을 찾기 위해 코드를 벤치마킹 해 보셨습니까? 빠른 리드 쓰루를 통해, Cython이나 Numba가 많은 속도 향상을 제공 할 수 있다는 것을 즉시 알 수 없습니다. 대부분의 작업은 이미 벡터화 된 방식으로 수행되고 있습니다. 먼저 [line profiler] (https://jakevdp.github.io/PythonDataScienceHandbook/01.07-timing-and-profiling.html#Line-By-Line-Profiling-with-%lprun)를 사용하여 시작합니다. 느린 점이있는 곳. – jakevdp

@jakevdp 의견을 보내 주셔서 감사합니다. 나는 라인 프로파일 러를 사용하기 위해 보았지만 Cython과 Python3에서 사용하는 방법을 먼저 배워야 할 것으로 보인다 ... 시간이 좀 걸린다. 도움이 될 경우 Cython 컴파일 결과에 주석 모드를 추가했습니다. – hopeso

나는 원래/비 -cython 코드에서 line profiler를 사용하여 어떤 작업이 느린 지 확인하는 것이 좋습니다. 그것들이 기본적인 numpy-primitives/vectorized 부분이라면, cython이 도움이되지 않는다는 것을 알 것입니다. – sascha

이 내가 프로파일을 갖고 있지 않기 때문에 매우 불완전한 대답은 나 :

참고로

, 다음은 컴파일에 사이 썬 주석입니다 무엇이든 시간을 정하거나 심지어 같은 대답을 준다는 것을 확인했습니다.

가 @cython.cdivision(True) 컴파일 지시문을 추가합니다 : 그러나 여기 사이 썬 생성 파이썬 코드의 양을 줄이는 몇 가지 제안 사항입니다. 즉, ZeroDivisionError은 부동 부서에서 발생하지 않으며 대신 NaN 값을 갖게됩니다. 오류를 발생시키지 않으려면이 작업을 수행하십시오.
변경 p_s = np.sqrt(...)에서 p_s = sqrt(...)으로 변경하십시오. 이렇게하면 단일 값에 대해서만 작동하는 감도없는 호출이 제거됩니다. 당신은 다른 곳에서 이런 짓을 한 것 같아요. 왜이 선을 놓쳤는 지 모르겠군요.
가능한 사용 고정 된 크기의 C 배열 대신 NumPy와 배열 : 크기는 컴파일 시간 (그리고 매우 작은)하고 반환하지 않으려는 알려진 경우
```
cdef double beta[3] 
# ... 
beta[0] = X[1]/X[5] 
beta[1] = X[3]/X[5] 
beta[2] = p_s/X[5] 
```
당신은이 작업을 수행 할 수 있습니다 그것. 이렇게하면 np.zeros에 대한 호출과 입력 된 numpy 배열에 할당하기위한 몇 가지 후속 유형 검사가 수행되지 않습니다. 나는 beta이 당신이 이것을 할 수있는 유일한 장소라고 생각합니다.
np.angle(complex(x/r, y/r))은 (을 libc.math에서 사용하여 atan2(y/r, x/r)으로 바꿀 수 있습니다.당신은 또한 (사이 썬은에서 자동 루프 변수의 유형을 따기 종종 좋은 그러나 여기 실패한 것으로 보인다) 나는 의심
getEMField에 빠르게 for 루프를 만드는 데 도움이 r
cdef int i에 의해 분할을 잃을 수 있습니다 list 및 01와 같은 특정 유형의 많은 가치가 없다
```
 E[0] = k2*x/0.045 + k10*rn**9*cos(9*arg) 
     E[1] = -k2*y/0.045 -k10*rn**9*sin(9*arg) 
```
: 그것은보다 전체 배열로 E 요소별로 요소 할당 빠릅니다이며 실제로 코드를 약간 느리게 만들 수 있습니다 (유형을 확인하는 데 시간이 낭비되기 때문).
더 큰 변화는 포인터보다는 np.zeros 할당 사용으로 GetEMField에 E 및 B을 통과하는 것입니다. 이렇게하면 equations (cdef double E[3])에 정적 C 배열로 할당 할 수 있습니다. 단점은 GetEMField이 cdef이어야하므로 더 이상 파이썬에서 호출 할 수 없다는 것입니다 (그러나 원하는 경우 Python 호출 가능 래퍼 함수도 만들 수 있습니다).

출처

2017-12-16 10:05:03 DavidW

답변

관련 문제