Hi,
I want to simulate a Jump process for biological applications quite precisely and fast. I developped a package PiecewiseDeterministicMarkovProcesses.jl based on DifferentialEquations.jl to simulate these jump processes.
This package has been optimized a lot to reduce allocations and be complient with AD. Note that for my application where the dimension is > 35 elements, I dont want to use StaticArrays.
In my application, I find lsoda to be the fastest and the most reliable despite the fact that I cannot use an ODE iterator with it (and hence it allocates). Thus, I want to be able to be as close as possible performance wise, to the lsoda solver but I dont manage to tweak the options to achieve this. If anybody has experience with this, I would not mind a little push.
Thank you.
PS: the MWE below gives the following timings. Note that we have an analytical expression of the jump process so we can attest the error. The results are
Solvers comparison
--> norm difference = 0.007230211093428807 - solver = lsoda
0.000935 seconds (9.61 k allocations: 502.719 KiB)
--> norm difference = 0.00043447196503620944 - solver = rodas5P
0.001447 seconds (2.08 k allocations: 160.797 KiB)
--> norm difference = 12.221574659298312 - solver = TRBDF2
0.002779 seconds (2.07 k allocations: 158.578 KiB)
--> norm difference = 0.010638187872245908 - solver = CVODEAdams
0.001777 seconds (17.38 k allocations: 480.062 KiB)
--> norm difference = 0.005070017916750658 - solver = CVODEBDF
0.002406 seconds (18.05 k allocations: 521.891 KiB)
--> norm difference = 0.0001112569666474883 - solver = rodas4P
0.002659 seconds (4.47 k allocations: 347.438 KiB)
--> norm difference = 0.005069408256076713 - solver = cvode
0.003115 seconds (24.58 k allocations: 820.281 KiB)
--> norm difference = 0.13379368165760752 - solver = Rosenbrock23
0.005789 seconds (14.07 k allocations: 1.071 MiB)
--> norm difference = 0.01244811034121085 - solver = AutoTsit5-RS23
0.000439 seconds (329 allocations: 24.422 KiB)
Code on PiecewiseDeterministicMarkovProcesses#master
# using Revise
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
const PDMP = PiecewiseDeterministicMarkovProcesses
function AnalyticalSampleCHV(xc0, xd0, ti, nj::Int64)
xch = [xc0[1]]
xdh = [xd0[1]]
th = [ti]
t = ti
while length(th)<nj
xc = xch[end]
xd = xdh[end]
S = -log(rand())
if mod(xd,2) == 0
t += 1/10*log(1+10*S/xc)
push!(xch,xc + 10 * S )
else
t += 1/(3xc)*(exp(3S)-1)
push!(xch,xc * exp(-3S) )
end
push!(xdh,xd + 1 )
push!(th,t)
S = -log(rand())
end
return th, xch, xdh
end
function F!(áş, xc, xd, parms, t)
if mod(xd[1], 2)==0
áş[1] = 10xc[1]
else
áş[1] = -3xc[1]^2
end
end
R(x) = x
function R!(rate, xc, xd, parms, t, issum::Bool)
# rate function
if issum == false
rate[1] = R(xc[1])
rate[2] = parms[1]
return 0., parms[1] + 50.
else
return R(xc[1]) + parms[1], parms[1] + 50.
end
end
xc0 = [1.0]
xd0 = [0, 0]
nu = [1 0;0 -1]
parms = [.0]
ti = 0.322156
tf = 100000.
nj = 50
errors = Float64[]
Random.seed!(8)
res_a_chv = AnalyticalSampleCHV(xc0,xd0,ti,nj)
problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
println("\n\nSolvers comparison")
for ode in [ (:lsoda,"lsoda"),
(Rodas5P(),"rodas5P"),
(TRBDF2(),"TRBDF2"),
(CVODE_Adams(),"CVODEAdams"),
(CVODE_BDF(),"CVODEBDF"),
(Rodas4P(),"rodas4P"),
(:cvode,"cvode"),
(Rosenbrock23(),"Rosenbrock23"),]
Random.seed!(8)
res = PDMP.solve(problem, CHV(ode[1]); n_jumps = nj, abstol = 1e-9, reltol = 1e-7)
printstyled(color=:green, "\n--> norm difference = ", norm(res.time - res_a_chv[1], Inf64), " - solver = ",ode[2],"\n")
Random.seed!(8)
res = @time PDMP.solve(problem, CHV(ode[1]); n_jumps = nj, abstol = 1e-9, reltol = 1e-7)
push!(errors,norm(res.time - res_a_chv[1], Inf64))
end