Как я могу оптимизировать спектр Jonswap в Cython - PullRequest
Новая
       6

Как я могу оптимизировать спектр Jonswap в Cython

0 голосов
/ 07 ноября 2018

Я пытаюсь улучшить производительность вычисления спектра Jonswap с использованием Cython. Но это кажется намного медленнее, чем оригинальный код. Как я могу улучшить это?

код Cython: 

from libc.math cimport exp
from libc.stdlib cimport malloc
import numpy as np
cimport numpy as np

DTYPE_float = np.float64
ctypedef np.float64_t DTYPE_float_t

def jonswap(np.ndarray[DTYPE_float_t, ndim=1, mode ='c'] w, double Hs, double Tp, double gamma = 3.7):
    '''
    get Jonswap spectra
    :param w: np.ndarray Angular Frequency
    '''
    cdef:
        int n = w.shape[0]
        double *sigma = <double*>malloc(n * sizeof(double)) 
        double *a = <double*>malloc(n * sizeof(double)) 
        int i 
    cdef double wp
    cdef np.ndarray[DTYPE_float_t, ndim=1, mode='c'] sj = np.ones(n, dtype=DTYPE_float)

    wp = 2 * np.pi / Tp
    for i in range(n):
        sigma[i] = 0.07 if w[i] < wp else 0.09
        a[i] = exp(-0.5 * pow((w[i] - wp) / (sigma[i] * w[i]), 2.0))
        sj[i] = 320 * pow(Hs, 2) * pow(w[i], -5.0) / pow(Tp, 4) * exp(-1950 * pow(w[i], -4) / pow(Tp, 4)) * pow(gamma, a[i])

    return sj

оригинальный код: 

def jonswap(w: np.ndarray, Hs: float, Tp: float, gamma: float = 3.7) -> np.ndarray:
    '''
    get Jonswap spectra
    :param w: np.ndarray Angular Frequency
    '''
    omega = w
    wp = 2 * np.pi / Tp
    sigma = np.where(omega < wp, 0.07, 0.09)
    a = np.exp(-0.5 * np.power((omega - wp) / (sigma * omega), 2.0))
    sj = 320 * np.power(Hs, 2) * np.power(omega, -5.0) / np.power(Tp, 4) * \
          np.exp(-1950 * np.power(omega, -4) / np.power(Tp, 4)) * np.power(gamma, a)

    return sj

1 Ответ

0 голосов
/ 08 ноября 2018

Ваш оригинальный код - все векторизованные nupy ops, поэтому возможности для улучшения ограничены. Запуск Cython с флагом аннотации (-a) указывает на следующие возможные улучшения.

  • используйте libc pow вместо встроенного Python
  • семантика обхода / обхода elide
  • используйте деление c, чтобы отключить проверку деления на 0 (если это безопасно!)

Новая версия Cython 

from libc.math cimport exp, pow
from libc.stdlib cimport malloc
import numpy as np
cimport numpy as np
cimport cython

DTYPE_float = np.float64
ctypedef np.float64_t DTYPE_float_t

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def cy_jonswap(np.ndarray[DTYPE_float_t, ndim=1, mode ='c'] w, double Hs, double Tp, double gamma = 3.7):
    '''
    get Jonswap spectra
    :param w: np.ndarray Angular Frequency
    '''
    cdef:
        int n = w.shape[0]
        double *sigma = <double*>malloc(n * sizeof(double)) 
        double *a = <double*>malloc(n * sizeof(double)) 
        int i 
    cdef double wp
    cdef np.ndarray[DTYPE_float_t, ndim=1, mode='c'] sj = np.ones(n, dtype=DTYPE_float)

    wp = 2 * np.pi / Tp
    with nogil:
        for i in range(n):
            sigma[i] = 0.07 if w[i] < wp else 0.09
            a[i] = exp(-0.5 * pow((w[i] - wp) / (sigma[i] * w[i]), 2.0))
            sj[i] = 320 * pow(Hs, 2) * pow(w[i], -5.0) / pow(Tp, 4) * exp(-1950 * pow(w[i], -4) / pow(Tp, 4)) * pow(gamma, a[i])

    return sj

Задержка 

w = np.random.randn(1_000_000)

%timeit cy_jonswap(w, .5, .5)
289 ms ± 7.34 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit jonswap(w, .5, .5)
411 ms ± 26.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Также обратите внимание, что в вашей версии Cython вы теряете память для sigma и a

...