Wrappers for the SciPy differential equation solvers
Julia
Latest commit 411cf4c Nov 24, 2019
Type Name Latest commit message Commit time
Failed to load latest commit information.
.github/workflows Nov 24, 2019
src Nov 24, 2019
test
.codecov.yml Nov 23, 2019
.gitignore Nov 23, 2019
.travis.yml Nov 23, 2019
Project.toml Nov 24, 2019
appveyor.yml Nov 23, 2019

# SciPyDiffEq.jl

SciPyDiffEq.jl is a common interface binding for the SciPy solve_ivp module ordinary differential equation solvers. It uses the PyCall.jl interop in order to send the differential equation over to Python and solve it.

Note that this package isn't for production use and is mostly just for benchmarking and helping new users migrate models over to Julia. For more efficient solvers, see the DifferentialEquations.jl documentation.

## Installation

To install SciPyDiffEq.jl, use the following:

`Pkg.clone("https://github.com/JuliaDiffEq/SciPyDiffEq.jl")`

## Using SciPyDiffEq.jl

SciPyDiffEq.jl is simply a solver on the DiffEq common interface, so for details see the DifferentialEquations.jl documentation. The available algorithms are:

```SciPyDiffEq.RK45
SciPyDiffEq.RK23
SciPyDiffEq.BDF
SciPyDiffEq.LSODA```

## Example

```using SciPyDiffEq

function lorenz(u,p,t)
du1 = 10.0(u[2]-u[1])
du2 = u[1]*(28.0-u[3]) - u[2]
du3 = u[1]*u[2] - (8/3)*u[3]
[du1, du2, du3]
end
tspan = (0.0,10.0)
u0 = [1.0,0.0,0.0]
prob = ODEProblem(lorenz,u0,tspan)
sol = solve(prob,SciPyDiffEq.RK45())```

In the following we can measure the overhead and show that using SciPy from Julia is about 3x faster than using SciPy with Numba. Using SciPyDiffEq:

```using SciPyDiffEq, BenchmarkTools

function lorenz(u,p,t)
du1 = 10.0(u[2]-u[1])
du2 = u[1]*(28.0-u[3]) - u[2]
du3 = u[1]*u[2] - (8/3)*u[3]
[du1, du2, du3]
end
tspan = (0.0,100.0)
u0 = [1.0,0.0,0.0]
prob = ODEProblem(lorenz,u0,tspan)
@btime sol = solve(prob,SciPyDiffEq.RK45(),dense=false, abstol=1e-8,reltol=1e-8) # 2.760 s (4426860 allocations: 182.27 MiB)```

This gives 2.76s. Solving the equivalent problem with SciPy `odeint` is:

```import numpy as np
from scipy.integrate import odeint
import timeit
import numba
def f(u, t):
x, y, z = u
return [10.0 * (y - x), x * (28.0 - z) - y, x * y - 2.66 * z]

u0 = [1.0,0.0,0.0]
tspan = (0., 100.)
t = np.linspace(0, 100, 1001)
sol = odeint(f, u0, t)
def time_func():
odeint(f, u0, t, rtol = 1e-8, atol=1e-8)

_t = timeit.Timer(time_func).timeit(number=100)
print(_t) # 13.898981100000015 seconds```

which takes 13.89 seconds. Then using Numba JIT with nopython mode is:

```numba_f = numba.jit(f,nopython=True)
odeint(numba_f, u0, t,rtol = 1e-8, atol=1e-8)

def time_func():
odeint(numba_f, u0, t,rtol = 1e-8, atol=1e-8)

_t = timeit.Timer(time_func).timeit(number=100)
print(_t) # 8.05035870000006 seconds```

which takes 8 seconds. Solving it with SciPy `solve_ivp` is:

```import numpy as np
from scipy.integrate import solve_ivp
import timeit
import numba
def f(t,u):
x, y, z = u
return [10.0 * (y - x), x * (28.0 - z) - y, x * y - 2.66 * z]

u0 = [1.0,0.0,0.0]
tspan = (0.0, 100.0)
t = np.linspace(0, 100, 1001)
sol = solve_ivp(f,(0.0, 100.0),u0,t_eval=t)

def time_func():
solve_ivp(f,(0.0, 100.0),u0,t_eval=t)

_t = timeit.Timer(time_func).timeit(number=100)
print(_t) # 15.978812399999999 seconds```

and

```numba_f = numba.jit(f,nopython=True)
sol = solve_ivp(numba_f,(0.0, 100.0),u0,t_eval=t)

def time_func():
solve_ivp(numba_f,(0.0, 100.0),u0,t_eval=t)

_t = timeit.Timer(time_func).timeit(number=100)
print(_t) # 14.302745000000002 seconds```

which Numba seems to be unable to effectively accelerate. Together, this showcases a 3x speedup over the best SciPy+Numba setup by using the Julia based interface, (and 5x head-to-head via `solve_ivp`) so overhead concerns in future benchmarks are gone because any measurement here is accelerating SciPy more than standard accelerated use.

## Benchmarks

The following benchmarks demonstrate a 1000x-5,000,000x performance advantage for the pure-Julia methods over the Julia-accelerated (3x) SciPy ODE solvers across a range of stiff and non-stiff ODEs. These were ran with Julia 1.2, MATLAB 2019B, deSolve 1.2.5, and SciPy 1.3.1 after verifying negligible overhead on interop.

(Yes, it sounds ridiculous, so run the benchmark problems yourself in both Julia and Python to prove it.)

#### Non-Stiff Problem 1: Lotka-Volterra

```using ParameterizedFunctions, MATLABDiffEq, OrdinaryDiffEq, ODEInterface,
ODEInterfaceDiffEq, Plots, Sundials, SciPyDiffEq, deSolveDiffEq
using DiffEqDevTools
using LinearAlgebra

f = @ode_def_bare LotkaVolterra begin
dx = a*x - b*x*y
dy = -c*y + d*x*y
end a b c d
p = [1.5,1,3,1]
tspan = (0.0,10.0)
u0 = [1.0,1.0]
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
Dict(:alg=>dopri5())
Dict(:alg=>Tsit5())
Dict(:alg=>Vern7())
Dict(:alg=>MATLABDiffEq.ode45())
Dict(:alg=>MATLABDiffEq.ode113())
Dict(:alg=>SciPyDiffEq.RK45())
Dict(:alg=>SciPyDiffEq.LSODA())
Dict(:alg=>deSolveDiffEq.lsoda())
Dict(:alg=>deSolveDiffEq.ode45())
]

names = [
"Julia: DP5"
"Hairer: dopri5"
"Julia: Tsit5"
"Julia: Vern7"
"MATLAB: ode45"
"MATLAB: ode113"
"SciPy: RK45"
"SciPy: LSODA"
"deSolve: lsoda"
"deSolve: ode45"
]

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names = names,
appxsol=test_sol,dense=false,
save_everystep=false,numruns=100,maxiters=10000000,
timeseries_errors=false,verbose=false)
plot(wp,title="Non-stiff 1: Lotka-Volterra")
savefig("benchmark1.png")```

#### Non-Stiff Problem 2: Rigid Body

```f = @ode_def_bare RigidBodyBench begin
dy1  = -2*y2*y3
dy2  = 1.25*y1*y3
dy3  = -0.5*y1*y2 + 0.25*sin(t)^2
end
prob = ODEProblem(f,[1.0;0.0;0.9],(0.0,100.0))
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
Dict(:alg=>dopri5())
Dict(:alg=>Tsit5())
Dict(:alg=>Vern7())
Dict(:alg=>MATLABDiffEq.ode45())
Dict(:alg=>MATLABDiffEq.ode113())
Dict(:alg=>SciPyDiffEq.RK45())
Dict(:alg=>SciPyDiffEq.LSODA())
Dict(:alg=>deSolveDiffEq.lsoda())
Dict(:alg=>deSolveDiffEq.ode45())
]

names = [
"Julia: DP5"
"Hairer: dopri5"
"Julia: Tsit5"
"Julia: Vern7"
"MATLAB: ode45"
"MATLAB: ode113"
"SciPy: RK45"
"SciPy: LSODA"
"deSolve: lsoda"
"deSolve: ode45"
]

abstols = 1.0 ./ 10.0 .^ (6:13)
reltols = 1.0 ./ 10.0 .^ (3:10)
wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names = names,
appxsol=test_sol,dense=false,
save_everystep=false,numruns=100,maxiters=10000000,
timeseries_errors=false,verbose=false)
plot(wp,title="Non-stiff 2: Rigid-Body")
savefig("benchmark2.png")```

#### Stiff Problem 1: ROBER Shorter and Simpler for SciPy

```rober = @ode_def begin
dy₁ = -k₁*y₁+k₃*y₂*y₃
dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e3),[0.04,3e7,1e4])
sol = solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (7:8)
reltols = 1.0 ./ 10.0 .^ (3:4);

setups = [Dict(:alg=>Rosenbrock23())
Dict(:alg=>TRBDF2())
Dict(:alg=>rodas())
Dict(:alg=>MATLABDiffEq.ode23s())
Dict(:alg=>MATLABDiffEq.ode15s())
Dict(:alg=>SciPyDiffEq.LSODA())
Dict(:alg=>SciPyDiffEq.BDF())
Dict(:alg=>deSolveDiffEq.lsoda())
]

names = [
"Julia: Rosenbrock23"
"Julia: TRBDF2"
"Hairer: rodas"
"MATLAB: ode23s"
"MATLAB: ode15s"
"SciPy: LSODA"
"SciPy: BDF"
"deSolve: lsoda"
]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names = names,print_names = true,
dense=false,verbose = false,
save_everystep=false,appxsol=test_sol,
maxiters=Int(1e5))
plot(wp,title="Stiff 1: ROBER Short")
savefig("benchmark31.png")```

#### Rober Standard length

SciPy Omitted due to failures at higher tolerances and because it's too slow to finish in a day! See note below.

```prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
sol = solve(prob,CVODE_BDF(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);

setups = [Dict(:alg=>Rosenbrock23())
Dict(:alg=>TRBDF2())
Dict(:alg=>rodas())
Dict(:alg=>MATLABDiffEq.ode23s())
Dict(:alg=>MATLABDiffEq.ode15s())
#Dict(:alg=>SciPyDiffEq.LSODA())
#Dict(:alg=>SciPyDiffEq.BDF())
Dict(:alg=>deSolveDiffEq.lsoda())
]

names = [
"Julia: Rosenbrock23"
"Julia: TRBDF2"
"Hairer: rodas"
"MATLAB: ode23s"
"MATLAB: ode15s"
#"SciPy: LSODA"
#"SciPy: BDF"
"deSolve: lsoda"
]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names = names,print_names = true,
dense=false,verbose = false,
save_everystep=false,appxsol=test_sol,
maxiters=Int(1e5))
plot(wp,title="Stiff 1: ROBER Standard Length")
savefig("benchmark32.png")```

#### Stiff Problem 2: HIRES

```f = @ode_def Hires begin
dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
dy2 = 1.71*y1 - 8.75*y2
dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
0.43*y6 + 0.69*y7
dy7 = 280.0*y6*y8 - 1.81*y7
dy8 = -280.0*y6*y8 + 1.81*y7
end

u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(f,u0,(0.0,321.8122))

sol = solve(prob,Rodas5(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

abstols = 1.0 ./ 10.0 .^ (5:8)
reltols = 1.0 ./ 10.0 .^ (1:4);
setups = [Dict(:alg=>Rosenbrock23())
Dict(:alg=>TRBDF2())
Dict(:alg=>rodas())
Dict(:alg=>MATLABDiffEq.ode23s())
Dict(:alg=>MATLABDiffEq.ode15s())
Dict(:alg=>SciPyDiffEq.LSODA())
Dict(:alg=>SciPyDiffEq.BDF())
Dict(:alg=>deSolveDiffEq.lsoda())
]

names = [
"Julia: Rosenbrock23"
"Julia: TRBDF2"
"Hairer: rodas"
"MATLAB: ode23s"
"MATLAB: ode15s"
"SciPy: LSODA"
"SciPy: BDF"
"deSolve: lsoda"
]

wp = WorkPrecisionSet(prob,abstols,reltols,setups;
names = names,print_names = true,
save_everystep=false,appxsol=test_sol,
maxiters=Int(1e5),numruns=100)
plot(wp,title="Stiff 2: Hires")
savefig("benchmark4.png")```

## Test Problem Failures

Note that these benchmarks also showcase that the SciPy integrators seem to fail on standard stiff ODE benchmark problems. This is probably even more worrisome than the timing results...

A simple testing script is:

```using ParameterizedFunctions, SciPyDiffEq
rober = @ode_def begin
dy₁ = -k₁*y₁+k₃*y₂*y₃
dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
dy₃ =  k₂*y₂^2
end k₁ k₂ k₃
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e3),[0.04,3e7,1e4])
@show solve(prob,SciPyDiffEq.LSODA(),reltol=1e-1,abstol=1e-5).retcode
@show solve(prob,SciPyDiffEq.LSODA(),reltol=1e-2,abstol=1e-6).retcode
@show solve(prob,SciPyDiffEq.LSODA(),reltol=1e-3,abstol=1e-7).retcode

@show solve(prob,SciPyDiffEq.BDF(),reltol=1e-1,abstol=1e-5).retcode
@show solve(prob,SciPyDiffEq.BDF(),reltol=1e-2,abstol=1e-6).retcode
@show solve(prob,SciPyDiffEq.BDF(),reltol=1e-3,abstol=1e-7).retcode```

For timings, notice that the standard case seems to take more than a day in Python. This is what the standard test case is (takes <1ms in Julia/C/Fortran). If anyone runs the solve long enough to let it finish, please let me know how long this takes (I only tried running it for a day). The following Python code runs it directly:

```from scipy.integrate import solve_ivp
def rober(t,u):
k1 = 0.04
k2 = 3e7
k3 = 1e4
y1, y2, y3 = u
dy1 = -k1*y1+k3*y2*y3
dy2 =  k1*y1-k2*y2*y2-k3*y2*y3
dy3 =  k2*y2*y2
return [dy1,dy2,dy3]

u0  = [1.0,0.0,0.0]
tspan = (0.0,1e5)
solve_ivp(rober,tspan,u0,t_eval = [0.0,1e5])```
You can’t perform that action at this time.