Addition of a Runge-Kutta version

.gitignore

@@ -8,3 +8,6 @@ Photo.run
 build
 *.a
 *.dylib
+_skbuild
+numbalsoda.egg-info
+.vscode

CMakeLists.txt

@@ -11,7 +11,7 @@ endif()
 
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON)
-set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_CXX_STANDARD 17)
 
 # lsoda
 add_library(lsoda SHARED 
@@ -41,12 +41,21 @@ add_library(solve_ivp SHARED
 set_target_properties(solve_ivp PROPERTIES PREFIX "lib")
 target_link_libraries(solve_ivp dop853_mod)
 
+# RKF 45
+add_library(rk45 SHARED
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/RK45.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/src/RK45_wrapper.cpp
+)
+set_target_properties(rk45 PROPERTIES PREFIX "lib")
+
 if (SKBUILD)
   install(TARGETS lsoda DESTINATION numbalsoda)
   install(TARGETS dop853 DESTINATION numbalsoda)
   install(TARGETS solve_ivp DESTINATION numbalsoda)
+  install(TARGETS rk45 DESTINATION numbalsoda)
 else()
   install(TARGETS lsoda DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/)
   install(TARGETS dop853 DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/)
   install(TARGETS solve_ivp DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/)
+  install(TARGETS rk45 DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/numbalsoda/)
 endif()

README.md

@@ -4,7 +4,7 @@
 
 `numbalsoda` also wraps the `dop853` explicit Runge-Kutta method from this repository: https://github.com/jacobwilliams/dop853
 
-This package is very similar to `scipy.integrate.solve_ivp` ([see here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)), when you set `method = 'LSODA'` or `method = DOP853`. But, `scipy.integrate.solve_ivp` invokes the python interpreter every time step which can be slow. Also, `scipy.integrate.solve_ivp` can not be used within numba jit-compiled python functions. In contrast, `numbalsoda` never invokes the python interpreter during integration and can be used within a numba compiled function which makes `numbalsoda` a lot faster than scipy for most problems, and achieves similar performance to Julia's DifferentialEquations.jl in some cases (see `benchmark` folder).
+This package is very similar to `scipy.integrate.solve_ivp` ([see here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)), when you set `method = 'LSODA'`, `method = DOP853` or `method = RK45`. But, `scipy.integrate.solve_ivp` invokes the python interpreter every time step which can be slow. Also, `scipy.integrate.solve_ivp` can not be used within numba jit-compiled python functions. In contrast, `numbalsoda` never invokes the python interpreter during integration and can be used within a numba compiled function which makes `numbalsoda` a lot faster than scipy for most problems, and achieves similar performance to Julia's DifferentialEquations.jl in some cases (see `benchmark` folder).
 
 ## Installation
 Conda:
@@ -19,7 +19,7 @@ python -m pip install numbalsoda
 ## Basic usage
 
 ```python
-from numbalsoda import lsoda_sig, lsoda, dop853
+from numbalsoda import lsoda_sig, lsoda, dop853, rk45
 from numba import njit, cfunc
 import numpy as np
 
@@ -39,6 +39,9 @@ usol, success = lsoda(funcptr, u0, t_eval, data = data)
 # integrate with dop853 method
 usol1, success1 = dop853(funcptr, u0, t_eval, data = data)
 
+# integrate with rk45 method (automatic adaptive time stepping)
+tsol2, usol2, tf2, success2 = rk45(funcptr, u0, -1.0, np.amin(t_eval), np.amax(t_eval), 100000, data = data)
+
 # usol = solution
 # success = True/False
 ```

benchmark/README.md

@@ -12,6 +12,7 @@ Also, Scipy might be faster than numbalsoda for stiff problems with > ~500 ODEs
 | -------------------------- | ---------- | ---------------------- |
 | numbalsoda lsoda           | 2.594 ms   | 1x                     |
 | numbalsoda dop853          | 0.700 ms   | 0.27x                  |
+| numbalsoda rk45            | 5.785 ms   | 2.23x
 | Scipy LSODA                | 109.349 ms | 42x                    |
 | Scipy RK45                 | 292.826 ms | 113x                   |
 | Scipy DOP853               | 180.997 ms | 70x                    |

benchmark/benchmark_lorenz.py

@@ -1,5 +1,5 @@
 import numpy as np
-from numbalsoda import lsoda_sig, lsoda, dop853
+from numbalsoda import lsoda_sig, lsoda, dop853, rk45
 from scipy.integrate import solve_ivp
 import timeit
 import numba as nb
@@ -24,14 +24,16 @@ def f_sp(t, u, sigma, rho, beta):
 
 u0 = np.array([1.0,0.0,0.0])
 tspan = np.array([0., 100.])
-t_eval = np.linspace(tspan[0], tspan[1], 1001)
 args = np.array([10.0,28.0,8./3.])
 args_tuple = tuple(args)
 rtol = 1.0e-8
 atol = 1.0e-8
 
+tsol_nb_rk, usol_nb_rk, tf_rk, success = rk45(funcptr, u0, -1.0, tspan[0], tspan[1], 30000, data=args, rtol=rtol, atol=atol)
+assert np.isclose(tf_rk, tspan[1])
+t_eval = tsol_nb_rk
+usol_nb, success = lsoda(funcptr, u0, t_eval, args, rtol=rtol, atol=atol)
 usol_sp = solve_ivp(f_sp, tspan, u0, t_eval = t_eval, args=args_tuple, rtol=rtol, atol=atol, method='LSODA')
-usol_nb, success = lsoda(funcptr, u0, t_eval, args, rtol=rtol, atol=atol) 
 
 @nb.njit(boundscheck=False)
 def time_nb():
@@ -41,6 +43,10 @@ def time_nb():
 def time_dop853():
     usol_nb, success = dop853(funcptr, u0, t_eval, args, rtol=rtol, atol=atol)  
 
+@nb.njit(boundscheck=False)
+def time_rk45():
+    _, _, _, _ = rk45(funcptr, u0, -1.0, tspan[0], tspan[1], 30000, data=args, rtol=rtol, atol=atol)
+
 def time_sp_LSODA():
     usol_sp = solve_ivp(f_sp, tspan, u0, t_eval = t_eval, args=args_tuple, rtol=rtol, atol=atol, method='LSODA')
 
@@ -53,17 +59,20 @@ def time_sp_DOP853():
 # time
 time_nb()
 time_dop853()
+time_rk45()
 time_sp_LSODA()
 time_sp_RK45()
 time_sp_DOP853()
-iters = 10
+iters = 100
 t_nb = timeit.Timer(time_nb).timeit(number=iters)/iters*1e3
 t_dop853 = timeit.Timer(time_dop853).timeit(number=iters)/iters*1e3
+t_rk45 = timeit.Timer(time_rk45).timeit(number=iters)/iters*1e3
 t_sp_LSODA = timeit.Timer(time_sp_LSODA).timeit(number=iters)/iters*1e3
 t_sp_RK45 = timeit.Timer(time_sp_RK45).timeit(number=iters)/iters*1e3
 t_sp_DOP853 = timeit.Timer(time_sp_DOP853).timeit(number=iters)/iters*1e3
 print("numbalsoda lsoda time =",'%.3f'%t_nb,'ms')
 print("numbalsoda dop853 time =",'%.3f'%t_dop853,'ms')
+print("numbalsoda rk45 time =",'%.3f'%t_rk45,'ms')
 print("Scipy LSODA time =",'%.3f'%t_sp_LSODA,'ms')
 print("Scipy RK45 time =",'%.3f'%t_sp_RK45,'ms')
 print("Scipy DOP853 time =",'%.3f'%t_sp_DOP853,'ms')

numbalsoda/__init__.py

@@ -1,3 +1,4 @@
 from .driver import lsoda_sig, lsoda, address_as_void_pointer
 from .driver_dop853 import dop853
-from .driver_solve_ivp import solve_ivp, ODEResult
+from .driver_solve_ivp import solve_ivp, ODEResult
+from .driver_rk45 import rk45_sig, rk45

numbalsoda/driver_rk45.py

@@ -0,0 +1,53 @@
+import ctypes as ct
+import numba as nb
+from numba import njit, types
+import numpy as np
+import os
+import platform
+
+rk45_sig = types.void(types.double,
+                      types.CPointer(types.double),
+                      types.CPointer(types.double),
+                      types.CPointer(types.double))
+
+rootdir = os.path.dirname(os.path.realpath(__file__))+'/'
+
+if platform.uname()[0] == "Windows":
+    name = "librk45.dll"
+elif platform.uname()[0] == "Linux":
+    name = "librk45.so"
+else:
+    name = "librk45.dylib"
+librk45 = ct.CDLL(rootdir+name)
+rk45_wrapper = librk45.rk45_wrapper
+rk45_wrapper.argtypes = [ct.c_void_p, ct.c_int, ct.c_void_p, ct.c_void_p,\
+                         ct.c_double, ct.c_double, ct.c_double, ct.c_int,\
+                         ct.c_void_p, ct.c_void_p, ct.c_double, ct.c_double, ct.c_double,\
+                         ct.c_void_p, ct.c_void_p, ct.c_void_p]
+rk45_wrapper.restype = None
+
+@njit
+def rk45(funcptr, u0, dt0, t0, tf, itf, data = np.array([0.0], np.float64), rtol = 1.0e-3, atol = 1.0e-6, mxstep = 10000.0):
+    neq = len(u0)
+    usol = np.empty((itf + 1, neq),dtype=np.float64)
+    tsol = np.empty(itf + 1, dtype=np.float64)
+    success = np.array([999], np.int32)
+    actual_final_time = np.array([-1.0], dtype=np.float64)
+    actual_final_iteration = np.array([-1], dtype=np.int32)
+    rk45_wrapper(funcptr, neq, u0.ctypes.data, data.ctypes.data, dt0, t0, tf, itf,
+                 usol.ctypes.data, tsol.ctypes.data, rtol, atol, mxstep, success.ctypes.data,
+                 actual_final_time.ctypes.data, actual_final_iteration.ctypes.data)
+    success_ = True
+    if success != 1:
+        success_ = False
+    return tsol[:actual_final_iteration[0]], usol[:actual_final_iteration[0],:], actual_final_time[0], success_
+
+@nb.extending.intrinsic
+def address_as_void_pointer(typingctx, src):
+    """ returns a void pointer from a given memory address """
+    from numba import types 
+    from numba.core import cgutils
+    sig = types.voidptr(src)
+    def codegen(cgctx, builder, sig, args):
+        return builder.inttoptr(args[0], cgutils.voidptr_t)
+    return sig, codegen

setup.py

@@ -8,7 +8,7 @@
 setup(
     name="numbalsoda",
     packages=['numbalsoda'],
-    version='0.3.5',
+    version='0.3.6.dev0',
     license='MIT',
     install_requires=['numpy','numba'],
     author = 'Nicholas Wogan',