Как я могу использовать Cython, чтобы быстрее решить дифференциальное уравнение?

Question

Я хотел бы уменьшить время, затрачиваемое Сипи на решение дифференциального уравнения.

На практике я использовал пример, описанный в Python, в научных вычислениях в качестве шаблона. Потому что odeint берет на себя функцию f в качестве аргумента я написал эту функцию как статически типизированную версию Cython и надеялся, что время выполнения odeint значительно уменьшится.

Функция f содержится в файле с именем ode.pyx следующее:

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

def f(y, t, params):
  cdef double theta = y[0], omega = y[1]
  cdef double Q = params[0], d = params[1], Omega = params[2]
  cdef double derivs[2]
  derivs[0] = omega
  derivs[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)
  return derivs

def fCMath(y, double t, params):
  cdef double theta = y[0], omega = y[1]
  cdef double Q = params[0], d = params[1], Omega = params[2]
  cdef double derivs[2]
  derivs[0] = omega
  derivs[1] = -omega/Q + sin(theta) + d*cos(Omega*t)
  return derivs

Затем я создаю файл setup.py чтобы выполнить функцию:

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

setup(ext_modules=cythonize('ode.pyx'))

Скрипт, решающий дифференциальное уравнение (также содержащий версию Python f) называется solveODE.py и выглядит как:

import ode
import numpy as np
from scipy.integrate import odeint
import time

def f(y, t, params):
    theta, omega = y
    Q, d, Omega = params
    derivs = [omega,
             -omega/Q + np.sin(theta) + d*np.cos(Omega*t)]
    return derivs

params = np.array([2.0, 1.5, 0.65])
y0 = np.array([0.0, 0.0])
t = np.arange(0., 200., 0.05)

start_time = time.time()
odeint(f, y0, t, args=(params,))
print("The Python Code took: %.6s seconds" % (time.time() - start_time))

start_time = time.time()
odeint(ode.f, y0, t, args=(params,))
print("The Cython Code took: %.6s seconds ---" % (time.time() - start_time))

start_time = time.time()
odeint(ode.fCMath, y0, t, args=(params,))
print("The Cython Code incorpoarting two of DavidW_s suggestions took: %.6s seconds ---" % (time.time() - start_time))

Затем я запускаю:

python setup.py build_ext --inplace
python solveODE.py 

в терминале.

Время для версии Python составляет примерно 0,055 секунды, в то время как версия Cython занимает примерно 0,04 секунды.

Есть ли у кого-нибудь рекомендации по улучшению моей попытки решить дифференциальное уравнение, желательно без использования самой процедуры odeint с Cython?

редактировать

Я включил предложение DavidW в два файла ode.pyx а также solveODE.py Выполнение кода с этими предложениями заняло всего примерно 0,015 секунды.

Answer 1

Самое простое изменение (которое, вероятно, принесет вам много пользы) - это использование математической библиотеки C sin а также cos для операций над отдельными числами вместо числа. Призыв к numpy и время, потраченное на разработку, что это не массив, довольно дорого.

from libc.math cimport sin, cos

    # later
    -omega/Q + sin(theta) + d*cos(Omega*t)

Я был бы соблазн назначить тип для ввода d (ни один из других вводов не может быть легко набран без изменения интерфейса):

def f(y, double t, params):

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

Answer 2

tldr; используйте numba.jit для ускорения в 3 раза...

У меня нет большого опыта работы с Cython, но моя машина, похоже, получает такое же время вычислений для вашей строго Python-версии, поэтому мы должны быть в состоянии примерно сравнить яблоки с яблоками. я использовал numba скомпилировать функцию f (который я немного переписал, чтобы сделать его более приятным с компилятором).

def f(y, t, params):
    return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])

numba_f = numba.jit(f)

падение в numba_f вместо вашего ode.f дает мне этот вывод...

The Python Code took: 0.0468 seconds
The Numba Code took: 0.0155 seconds

Я тогда задавался вопросом, могу ли я дублировать odeint а также компилировать с Numba, чтобы ускорить процесс еще дальше... (я не мог)

Вот мой интегратор числовых дифференциальных уравнений Рунге-Кутты:

#function f is provided inline (not as an arg)
def runge_kutta(y0, steps, dt, args=()): #improvement on euler's method. *note: time steps given in number of steps and dt
    Y = np.empty([steps,y0.shape[0]])
    Y[0] = y0
    t = 0
    n = 0
    for n in range(steps-1):
        #calculate coeficients
        k1 = f(Y[n], t, args) #(euler's method coeficient) beginning of interval
        k2 = f(Y[n] + (dt * k1 / 2), t + (dt/2), args) #interval midpoint A
        k3 = f(Y[n] + (dt * k2 / 2), t + (dt/2), args) #interval midpoint B
        k4 = f(Y[n] + dt * k3, t + dt, args) #interval end point

        Y[n + 1] = Y[n] + (dt/6) * (k1 + 2*k2 + 2*k3 + k4) #calculate Y(n+1)
        t += dt #calculate t(n+1)
    return Y

Наивные циклические функции обычно бывают самыми быстрыми после компиляции, хотя, возможно, их можно было бы реструктурировать для большей скорости. Я должен отметить, что это дает другой ответ, чем odeint, отклоняясь на целых 0,001 примерно после 2000 шагов, и совершенно иным после 3000. Для версии функции numba я просто заменил f с numba_fи добавил компиляцию с @numba.jit в качестве декоратора. В этом случае, как и ожидалось, чисто Python-версия очень медленная, но версия Numba не быстрее, чем Numba с odeint (опять же, мммм).

using custom integrator
The Python Code took: 0.2340 seconds
The Numba Code took: 0.0156 seconds

Вот пример компиляции заранее. У меня нет необходимого набора инструментов на этом компьютере для компиляции, и у меня нет администратора для его установки, так что это выдает ошибку, что у меня нет требуемого компилятора, но он должен работать иначе.

import numpy as np
from numba.pycc import CC

cc = CC('diffeq')

@cc.export('func', 'f8[:](f8[:], f8, f8[:])')
def func(y, t, params):
    return np.array([y[1], -y[1]/params[0] + np.sin(y[0]) + params[1]*np.cos(params[2]*t)])

cc.compile()

Answer 3

Если другие ответят на этот вопрос с помощью других модулей, я мог бы также вмешаться:

Я являюсь автором JiTCODE, который принимает ODE, записанную в символах SymPy, а затем преобразует этот ODE в код C для модуля Python, компилирует этот код C, загружает результат и использует его как производную для ODE SciPy. Ваш пример, переведенный на JiTCODE, выглядит так:

from jitcode import jitcode, provide_basic_symbols
import numpy as np
from sympy import sin, cos
import time

Q = 2.0
d = 1.5
Ω = 0.65

t, y = provide_basic_symbols()

f = [
    y(1),
    -y(1)/Q + sin(y(0)) + d*cos(Ω*t)
    ]

initial_state = np.array([0.0,0.0])

ODE = jitcode(f)
ODE.set_integrator("lsoda")
ODE.set_initial_value(initial_state,0.0)

start_time = time.time()
data = np.vstack(ODE.integrate(T) for T in np.arange(0.05, 200., 0.05))
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))

Это занимает 0,11 секунды, что ужасно медленно по сравнению с решениями, основанными на odeint, но это не из-за фактической интеграции, а из-за способа обработки результатов: odeint непосредственно создает массив внутренне, это делается через Python здесь. В зависимости от того, что вы делаете, это может быть решающим недостатком, но это быстро становится неактуальным для более грубой выборки или больших дифференциальных уравнений.

Итак, давайте удалим сбор данных и просто посмотрим на интеграцию, заменив последние строки следующим:

ODE = jitcode(f)
ODE.set_integrator("lsoda", max_step=0.05, nsteps=1e10)
ODE.set_initial_value(initial_state,0.0)

start_time = time.time()
ODE.integrate(200.0)
end_time = time.time()
print("JiTCODE took: %.6s seconds" % (end_time - start_time))

Обратите внимание, что я установил max_step=0.05 заставить интегратора сделать как минимум столько же шагов, сколько в вашем примере, и убедиться, что единственное отличие состоит в том, что результаты интеграции не сохраняются в каком-либо массиве. Это выполняется за 0,010 секунд.

Answer 4

NumbaLSODA занимает 0,00088 секунд (в 17 раз быстрее, чем Cython).

      from NumbaLSODA import lsoda_sig, lsoda
import numba as nb
import numpy as np
import time

@nb.cfunc(lsoda_sig)
def f(t, y_, dy, p_):
    p = nb.carray(p_, (3,))
    y = nb.carray(y_, (2,))
    theta, omega = y
    Q, d, Omega = p
    dy[0] = omega
    dy[1] = -omega/Q + np.sin(theta) + d*np.cos(Omega*t)

funcptr = f.address # address to ODE function
y0 = np.array([0.0, 0.0])
data = np.array([2.0, 1.5, 0.65])
t = np.arange(0., 200., 0.05)

start_time = time.time()
usol, success = lsoda(funcptr, y0, t, data = data)
print("NumbaLSODA took: %.8s seconds ---" % (time.time() - start_time))

результат

      NumbaLSODA took: 0.000880 seconds ---