Help me beat Lsoda

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

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

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.

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?

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?

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

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

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?

Thanks!

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.

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

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