# Help me beat Lsoda

**URL:** <https://discourse.julialang.org/t/help-me-beat-lsoda/88236>\
**Category:** Performance\
**Created:** [October 4, 2022, 3:08pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236 "2022-10-04T15:08:45Z")\
**Posts on this page:** 7\
**Page:** 2

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [April 13, 2023, 6:43pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/21 "2023-04-13T18:43:28Z")

</div>

> [@gzagatti](#):
>
> I’m one of the authors of the Coevolve paper,

can you link your paper please?

---

<div class="post-metadata">

**Author:** ![gzagatti](https://avatars.discourse-cdn.com/v4/letter/g/97f17d/32.png) [@gzagatti](https://discourse.julialang.org/u/gzagatti)\
**Post date:** [April 14, 2023, 3:08pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/22 "2023-04-14T15:08:16Z")

</div>

> [@rveltz](#):
>
> can you link your paper please?

You can find a version of it in this PR: [JuliaCon 2023 Paper by gzagatti · Pull Request #304 · SciML/JumpProcesses.jl · GitHub](https://github.com/SciML/JumpProcesses.jl/pull/304)

> [@ChrisRackauckas](#):
>
> Can we get the ODEProblem definition out?

Next, I’ve started an attempt to map the PDMP problem to an ODEProblem in this fork: [https://github.com/gzagatti/SynapseElife/blob/coevolve/src/SynapseModelODEProblem.jl](https://github.com/gzagatti/SynapseElife/blob/coevolve/src/SynapseModelODEProblem.jl)

However, I think it will be difficult to use Coevolve here. The reason is that when the algorithm is searching for the next jump time it needs to compute the rate at the candidate time. But some rates depends on a continuous variable which needs to be brought forward when computing the rate to determine the next time. This is a known limitation of Coevolve as explained in the [documentation](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Defining-a-Variable-Rate-Jump).

> Note that it is ok if the `urate` bound would be violated within the `rateinterval` due to a change in `u` arising from another `ConstantRateJump` , `MassActionJump` or _bounded_ `VariableRateJump` being executed, as the chosen aggregator will then handle recalculating the rate bound and interval. _However, if the bound could be violated within the time interval due to a change in `u` arising from continuous dynamics such as a coupled ODE, SDE, or a general `VariableRateJump`, bounds should not be given._ This ensures the jump is classified as a general `VariableRateJump` and properly handled.

@ChrisRackauckas do you think there is a way around this?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 14, 2023, 3:26pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/23 "2023-04-14T15:26:48Z")

</div>

The ODE/SDE is bounded, so if one uses a reasonably conservative bound for the state variables then you get a somewhat conservative bound for the rates. Would that be fine enough?

---

<div class="post-metadata">

**Author:** ![gzagatti](https://avatars.discourse-cdn.com/v4/letter/g/97f17d/32.png) [@gzagatti](https://discourse.julialang.org/u/gzagatti)\
**Post date:** [April 14, 2023, 4:09pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/24 "2023-04-14T16:09:24Z")

</div>

The bound helps with computing the upper and lower bound which could help with fast rejecting candidates, but we still need to evaluate the actual rates when deciding whether to accept the proposed time or not. This is the part where we would need to bring the continuous part forward:

[https://github.com/SciML/JumpProcesses.jl/blob/fd47a8212711c85687e93b8bc9e90e82fb134d36/src/aggregators/coevolve.jl#L189](https://github.com/SciML/JumpProcesses.jl/blob/fd47a8212711c85687e93b8bc9e90e82fb134d36/src/aggregators/coevolve.jl#L189)

The good thing is that we could at least keep the same forwarded trajectory when updating all the dependent rates.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [April 14, 2023, 9:38pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/25 "2023-04-14T21:38:01Z")

</div>

Impressive that you could write that in such small time!  
If I follow well your arguments, you are implementing a thinning method and if the bound on the rate is not tight, you will spend time rejecting, is that right?

---

<div class="post-metadata">

**Author:** ![gzagatti](https://avatars.discourse-cdn.com/v4/letter/g/97f17d/32.png) [@gzagatti](https://discourse.julialang.org/u/gzagatti)\
**Post date:** [April 15, 2023, 3:21pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/26 "2023-04-15T15:21:01Z")

</div>

Thanks!

> [@rveltz](#):
>
> If I follow well your arguments, you are implementing a thinning method and if the bound on the rate is not tight, you will spend time rejecting, is that right?

Yes, that’s correct.

I’m trying to see if I can adapt the algorithm to work with continuous variable. I tried a few things today but it’s not quite there yet. Let’s see if I have more luck in a couple of days. I’ll keep you posted.

---

<div class="post-metadata">

**Author:** ![gzagatti](https://avatars.discourse-cdn.com/v4/letter/g/97f17d/32.png) [@gzagatti](https://discourse.julialang.org/u/gzagatti)\
**Post date:** [April 17, 2023, 4:54pm UTC](https://discourse.julialang.org/t/help-me-beat-lsoda/88236/27 "2023-04-17T16:54:29Z")

</div>

I have started an attempt at solving the problem using JumpProcesses.jl. However, the prototype is not yet working: [ODEProblem and CoevolveSynced by gzagatti · Pull Request #1 · rveltz/SynapseElife · GitHub](https://github.com/rveltz/SynapseElife/pull/1)

I hope we can get it working for additional benchmarks. I would appreciate some feedback.

[Previous page](https://discourse.julialang.org/t/help-me-beat-lsoda/88236.md?page=1)