Cythonize интегратор уравнения в частных производных

Я пытаюсь ускорить интегратор конечных разностей для уравнения в частных производных, используя Cython. Я не уверен, что мне нужно сделать, чтобы Cython правильно работал с массивами numpy.

Функция члена диффузии, которую я использую,

def laplacian(var, dh2):
    """ (1D array, dx^2) -> laplacian(1D array)
    periodic_laplacian_1D_4th_order
    Implementing the 4th order 1D laplacian with periodic condition
    """
    lap = numpy.zeros_like(var)
    lap[1:]    = (4.0/3.0)*var[:-1]
    lap[0]     = (4.0/3.0)*var[1]
    lap[:-1]  += (4.0/3.0)*var[1:]
    lap[-1]   += (4.0/3.0)*var[0]
    lap       += (-5.0/2.0)*var

    lap[2:]   += (-1.0/12.0)*var[:-2]
    lap[:2]   += (-1.0/12.0)*var[-2:]
    lap[:-2]  += (-1.0/12.0)*var[2:]
    lap[-2:]  += (-1.0/12.0)*var[:2]

    return lap / dh2

А правые части уравнений модели равны

from derivatives import laplacian

def dbdt(b,w,p,m,d,dx2):
    """ db/dt of Modified Klausmeier """
    return w*b**2 - m*b + laplacian(b,dx2)

def dwdt(b,w,p,m,d,dx2):
    """ dw/dt of Modified Klausmeier """
    return p - w - w*b**2 + d*laplacian(b,dx2)

Как я могу оптимизировать эти функции с помощью Cython?

У меня есть репозиторий на Github для моего рабочего кода, который интегрирует модель Грея-Скотта — интегратор моделей Грея-Скотта.


python numpy cython pde
person Ohm    schedule 09.09.2015    source источник
comment
Что вы пробовали? Вы получите дополнительную помощь, если у вас возникнут определенные проблемы с реализацией версии cython.   -  person hpaulj    schedule 10.09.2015

Ответы (2)


arrow_upward
1
arrow_downward

Чтобы эффективно использовать Cython, вы должны сделать все циклы явными и убедиться, что cython -a показывает как можно меньше вызовов Python. Первая попытка будет:

import numpy as np
cimport numpy as np
cimport cython
@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
def laplacian(double [::1] var, double dh2):
    """ (1D array, dx^2) -> laplacian(1D array)
    periodic_laplacian_1D_4th_order
    Implementing the 4th order 1D laplacian with periodic condition
    """
    cdef int n = var.shape[0]
    cdef double[::1] lap = np.zeros(n)
    cdef int i
    for i in range(0, n-1):
        lap[1+i] = (4.0/3.0)*var[i]
    lap[0]     = (4.0/3.0)*var[1]
    for i in range(0, n-1):
        lap[i]  += (4.0/3.0)*var[1+i]
    lap[n-1]   += (4.0/3.0)*var[0]
    for i in range(0, n):
        lap[i]       += (-5.0/2.0)*var[i]

    for i in range(0, n-2):
        lap[2+i]   += (-1.0/12.0)*var[i]
    for i in range(0, 2):
        lap[i]   += (-1.0/12.0)*var[n - 2 + i]
    for i in range(0, n-2):
        lap[i]   += (-1.0/12.0)*var[i+2]
    for i in range(0, 2):
        lap[n-2+i]  += (-1.0/12.0)*var[i]
    for i in range(0, n):
        lap[i]  /= dh2
    return lap

Теперь это дает вам:

$ python -m timeit -s 'import numpy as np; from lap import laplacian; var = np.random.rand(1000000); dh2 = .01' 'laplacian(var, dh2)'
100 loops, best of 3: 11.5 msec per loop

в то время как код NumPy дал:

100 loops, best of 3: 18.5 msec per loop

Обратите внимание, что Cython можно дополнительно оптимизировать, объединяя циклы и т. д.

Я также попытался использовать настроенную (т. е. не зафиксированную в мастере) версию Pythran без изменения исходного кода Python, У меня было такое же ускорение, как и в версии Cython, без хлопот с преобразованием кода:

#pythran export laplacian(float [], float)
import numpy
def laplacian(var, dh2):
    """ (1D array, dx^2) -> laplacian(1D array)
    periodic_laplacian_1D_4th_order
    Implementing the 4th order 1D laplacian with periodic condition
    """
    lap = numpy.zeros_like(var)
    lap[1:]    = (4.0/3.0)*var[:-1]
    lap[0]     = (4.0/3.0)*var[1]
    lap[:-1]  += (4.0/3.0)*var[1:]
    lap[-1]   += (4.0/3.0)*var[0]
    lap       += (-5.0/2.0)*var

    lap[2:]   += (-1.0/12.0)*var[:-2]
    lap[:2]   += (-1.0/12.0)*var[-2:]
    lap[:-2]  += (-1.0/12.0)*var[2:]
    lap[-2:]  += (-1.0/12.0)*var[:2]

    return lap / dh2

Конвертировано с:

$ pythran lap.py -O3

И я получаю:

100 loops, best of 3: 11.6 msec per loop
person serge-sans-paille    schedule 11.09.2015

arrow_upward
0
arrow_downward

Итак, я думаю, что понял это, хотя я не уверен, что это самый оптимизированный способ сделать это:

import numpy as np
cimport numpy as np
cdef laplacian(np.ndarray[np.float64_t, ndim=1] var,np.float64_t dh2):
    """ (1D array, dx^2) -> laplacian(1D array)
    periodic_laplacian_1D_4th_order
    Implementing the 4th order 1D laplacian with periodic condition
    """
    lap = np.zeros_like(var)
    lap[1:]    = (4.0/3.0)*var[:-1]
    lap[0]     = (4.0/3.0)*var[1]
    lap[:-1]  += (4.0/3.0)*var[1:]
    lap[-1]   += (4.0/3.0)*var[0]
    lap       += (-5.0/2.0)*var

    lap[2:]   += (-1.0/12.0)*var[:-2]
    lap[:2]   += (-1.0/12.0)*var[-2:]
    lap[:-2]  += (-1.0/12.0)*var[2:]
    lap[-2:]  += (-1.0/12.0)*var[:2]

    return lap / dh2

Я использовал следующий setup.py

from distutils.core import setup
from Cython.Build import cythonize

setup(
    ext_modules = cythonize("derivatives_c.pyx")
)

Приветствуются любые советы по улучшению..

person Ohm    schedule 11.09.2015
comment
Это не даст вам никакого хорошего ускорения, потому что большинство вычислений все еще выполняется на Python. Попробуйте cython -a в своем коде! Чтобы получить максимальную отдачу от cython, вы должны сделать цикл явным... - person serge-sans-paille; 11.09.2015