### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 4, 2022, 3:08pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/1 "2022-10-04T15:08:45Z")
Hi,
I want to simulate a Jump process for biological applications quite precisely and fast. I developped a package [PiecewiseDeterministicMarkovProcesses.jl](https://github.com/rveltz/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](https://diffeq.sciml.ai/stable/basics/integrator/) 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
```julia
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**
>
> ```julia
> # 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) 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
>
> ```
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 4, 2022, 3:36pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/2 "2022-10-04T15:36:50Z")
Have you tried `FBDF` (`QNDF` would also be good to try)? In my experience, that is often very good for stiff systems.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 4, 2022, 3:49pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/3 "2022-10-04T15:49:34Z")
None of them return the result, it errors because of dtmin.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 4, 2022, 4:05pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/4 "2022-10-04T16:05:28Z")
That’s slightly odd. Another to try is `rodas5P` which isn’t documented yet, but should be better than `Rodas4P`.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 4, 2022, 4:17pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/5 "2022-10-04T16:17:45Z")
Yes it comes close (see above) and it much more precise. I tried to change the tolerance to decrease the computation time of `Rodas5` while having similar error to LSODA. I arrive at
```julia
--> norm difference = 0.008140039428781165 - solver = rodas5P
0.000922 seconds (1.40 k allocations: 107.453 KiB)
```
so roughly similar to LSODA but with much smaller allocations.
> However in my application, it is 10x slower maybe because of jacobian computation
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [October 4, 2022, 5:01pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/6 "2022-10-04T17:01:25Z")
You should read through [Solving Large Stiff Equations · DifferentialEquations.jl](https://diffeq.sciml.ai/stable/tutorials/advanced_ode_example/) and in particular, try jacobian free methods (Newton krylov).
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 4, 2022, 9:52pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/7 "2022-10-04T21:52:08Z")
> [@rveltz](#):
>
> ```julia
> --> norm difference = 0.01244811034121085 - solver = AutoTsit5-RS23
> 0.000439 seconds (329 allocations: 24.422 KiB)
>
> ```
Are you sure it’s stiff?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 5, 2022, 6:04am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/8 "2022-10-05T06:04:26Z")
We solve either
dx/dt = G(x) = 10x / x (= 10)
or
dx/dt = G(x) = -3x^2 / x (= x)
which is indeed not stiff but the last simplification is not done in the code and x can be small to zero.
Thank you for your remark, my example is maybe to simple for a MWE
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 5, 2022, 1:16pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/9 "2022-10-05T13:16:24Z")
LSODA switches between a stiff and non-stiff ODE solver. If the problem is non-stiff, then LSODA will switch over to the non-stiff method which will have a much lower overhead than the stuff for stiff equations, and it’ll beat every single one handedly. But of course, that’s just because stiff ODE solvers are not optimized for non-stiff ODEs!
I just had a collaborator make this mistake the other day. There’s a few benchmarks which directly show this behavior:
[https://benchmarks.sciml.ai/dev/Bio/egfr\_net/](https://benchmarks.sciml.ai/dev/Bio/egfr_net/)
[https://benchmarks.sciml.ai/dev/Bio/multisite2/#Implicit-Work-Precision-Diagrams](https://benchmarks.sciml.ai/dev/Bio/multisite2/#Implicit-Work-Precision-Diagrams)
etc.
The point is, if LSODA is lower than everything you’re trying, it may just be in the other class of methods. Or, auto-switch may just be better.
Note you probably want to try `AutoVern8(Rodas5P())` for this case. Tsit5 + Rosenbrock23 is only for very high tolerances.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 7, 2022, 10:02am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/10 "2022-10-07T10:02:35Z")
Thanks a lot for your input! In my use case, I was able to beat lsoda with `CVODE_BDF(linear_solver=:GMRES)`. I tried `AutoVern8(Rodas5P())` but it could not compute the first jump time.
```julia
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/xWByK/src/integrator_interface.jl:523
ERROR: AssertionError: Could not compute next jump time 1.
Return code = Unstable
0.0 < 0.0,
```
“0.0 \< 0.0” est un assert de PDMP
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 7, 2022, 12:27pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/11 "2022-10-07T12:27:31Z")
Try tweaking the switching parameters? If we can get this case for benchmarks that would be helpful
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 8, 2022, 6:18am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/12 "2022-10-08T06:18:50Z")
I tried with `AutoVern8(Rodas5P(), nonstifftol = 1//100)` but it did not succeed
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 8, 2022, 6:24am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/13 "2022-10-08T06:24:21Z")
Can we turn this into a SciMLBenchmarks work-precision diagram? That will make it more concrete for performance tracking.
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [October 8, 2022, 6:51am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/14 "2022-10-08T06:51:07Z")
Will do once the paper is out.
Just to mention, there is a part similar to this in the code,
```julia
dxdt = F(x, y)
dydt = G(x, y)
dzdt = [x, y] ∈ fixed_region
```
Hence, the vector field is only piecewise continuous although it should not have a big effect because `z` does not feedback on x,y
```julia
```
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 13, 2023, 8:55am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/15 "2023-04-13T08:55:22Z")
The paper is accepted! Hence, I can put the code online [here](https://github.com/rveltz/SynapseElife). Any improvement would be a major win for us
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 13, 2023, 9:06am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/16 "2023-04-13T09:06:55Z")
Can we get the ODEProblem definition out?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 13, 2023, 9:10am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/17 "2023-04-13T09:10:29Z")
Tough, it is a PDMP sampled with the CHV method. Maybe it is a benchamark for `JumpProcesses`?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [April 13, 2023, 9:12am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/18 "2023-04-13T09:12:39Z")
It would be good to test the new Coevolve against that. I don’t know if you’ve seen that paper but it tests against the CHV in a few problems. Can you share the MWE and we can translate it?
---
### Author:  [@rveltz](https://discourse.julialang.org/u/rveltz)
#### Post date: [April 13, 2023, 9:32am UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/19 "2023-04-13T09:32:18Z")
> I don’t know if you’ve seen that paper
Yes I have seen this paper but can’t find it again. I had a few comments about it but I forgot
> Can you share the MWE and we can translate it?
It is written for [PiecewiseDeterministicMarkovProcesses.jl](https://github.com/rveltz/PiecewiseDeterministicMarkovProcesses.jl), not for JumpProcesses. The drift is [here](https://github.com/rveltz/SynapseElife/blob/main/src/SynapseModel.jl#L2), the total rate is [here](https://github.com/rveltz/SynapseElife/blob/main/src/SynapseModel.jl#L179) and the [solve](https://github.com/rveltz/SynapseElife/blob/main/src/SynapseModel.jl#L426).
It would benefit from PeriodicCallback but this is not in PiecewiseDeterministicMarkovProcesses.jl
---
### Author:  [@gzagatti](https://discourse.julialang.org/u/gzagatti)
#### Post date: [April 13, 2023, 2:36pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/20 "2023-04-13T14:36:07Z")
I’m one of the authors of the Coevolve paper, Let me see if I can map this problem to Coevolve tomorrow and see how it does.
[Next page](https://discourse.julialang.org/t/help-me-beat-lsoda/88236.md?page=2)