¿Cómo puedo usar Cython bien para resolver una ecuación diferencial más rápido?

Me gustaría reducir el tiempo que Scipy’s odeint toma para resolver una ecuación diferencial.

Para practicar, utilicé el ejemplo cubierto en Python en cálculos científicos como plantilla. Como odeint toma una función f como argumento, escribí esta función como una versión de Cython tipificada estáticamente y esperaba que el tiempo de ejecución de odeint disminuyera significativamente.

La función f está contenida en un archivo llamado ode.pyx siguiente manera:

 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 

Luego creo un archivo setup.py para completar la función:

 from distutils.core import setup from Cython.Build import cythonize setup(ext_modules=cythonize('ode.pyx')) 

El script que resuelve la ecuación diferencial (que también contiene la versión de Python de f ) se denomina solveODE.py y se ve así:

 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)) 

Entonces corro

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

en la terminal

El tiempo para la versión de python es de aproximadamente 0.055 segundos, mientras que la versión de Cython toma aproximadamente 0.04 segundos.

¿Alguien tiene una recomendación para mejorar mi bash de resolver la ecuación diferencial, preferiblemente sin retoques con la rutina odeint, con Cython?

Editar

Incorporé la sugerencia de DavidW en los dos archivos ode.pyx y solveODE.py Tomó solo 0.015 segundos ejecutar el código con estas sugerencias.

El cambio más fácil de hacer (que probablemente te hará ganar mucho) es usar la biblioteca matemática de C para los operaciones con números únicos en lugar de números. La llamada a numpy y el tiempo dedicado a resolver que no es una matriz es bastante costoso.

 from libc.math cimport sin, cos # later -omega/Q + sin(theta) + d*cos(Omega*t) 

Estaría tentado de asignar un tipo a la entrada d (ninguna de las otras entradas se escribe fácilmente sin cambiar la interfaz):

 def f(y, double t, params): 

Creo que también me gustaría devolver una lista como lo haces en tu versión de Python. No creo que gane mucho usando una matriz de C.

tldr; usar numba.jit para acelerar 3x …

No tengo mucha experiencia con cython, pero mi máquina parece tener tiempos de cómputo similares para su versión estrictamente python, por lo que deberíamos poder comparar aproximadamente manzanas con manzanas. Utilicé numba para comstackr la función f (que reescribí ligeramente para hacer que suene mejor con el comstackdor).

 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) 

soltar numba_f en lugar de tu ode.f me da esta salida …

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

Luego me pregunté si podría duplicar odeint y también comstackr con numba para acelerar aún más las cosas … (no podría)

Aquí está mi integrador de ecuación diferencial numérica Runge-Kutta:

 #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 

Las funciones de bucle ingenuas son típicamente las más rápidas una vez comstackdas, aunque esto probablemente podría ser reestructurado para una velocidad un poco mejor. Debo tener en cuenta que esto da una respuesta diferente a odeint , que se desvía tanto como .001 después de unos 2000 pasos, y es completamente diferente después de 3000. Para la versión numba de la función, simplemente reemplacé f con numba_f , y agregué la comstackción con @numba.jit como decorador. En este caso, como se esperaba, la versión de python pura es muy lenta, pero la versión numba no es más rápida que la numba con odeint (de nuevo, ymmv).

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

Aquí hay un ejemplo de comstackción anticipada. No tengo la cadena de herramientas necesaria en esta computadora para comstackr, y no tengo admin para instalarla, así que esto me da un error de que no tengo el comstackdor requerido, pero debería funcionar de otra manera.

 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() 

Si otros responden a esta pregunta usando otros módulos, yo también podría hacer el timbre en:

Soy el autor de JiTCODE , que acepta una ODE escrita en símbolos SymPy y luego convierte esta ODE en código C para un módulo Python, comstack este código C, carga el resultado y lo utiliza como un derivado para la ODE de SciPy . Su ejemplo traducido a JiTCODE se ve así:

 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)) 

Esto toma 0.11 segundos, lo que es terriblemente lento en comparación con las soluciones basadas en odeint , pero no se debe a la integración real, sino a la forma en que se manejan los resultados: mientras que odeint crea una matriz de manera interna, esto se hace a través de Python. Dependiendo de lo que haga, esto puede ser una desventaja crucial, pero esto rápidamente se vuelve irrelevante para un muestreo más grueso o ecuaciones diferenciales más grandes.

Entonces, eliminemos la recostackción de datos y observemos la integración, reemplazando las últimas líneas con lo siguiente:

 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)) 

Tenga en cuenta que configuro max_step=0.05 para forzar al integrador a realizar al menos tantos pasos como en su ejemplo y asegurar que la única diferencia es que los resultados de la integración no se almacenan en algún arreglo. Esto funciona en 0.010 segundos.