Certain Continuous Callbacks does not work for JumpProcessess

[TorkelE]

A simple example where the continuous callback is never activated:

using Catalyst, DifferentialEquations

rn = @reaction_network begin
    (p,d), 0 <--> X
end

u0 = [:X => 100]
tspan = (0.0, 100.0)
p_vals = [:p => 10.0, :d => 1.0]

dprob = DiscreteProblem(rn, u0, tspan, p_vals)
jprob = JumpProblem(rn, dprob, Direct())

condition(u, t, integrator) = (u[1] - 25)
affect!(integrator) = println("Affect function activated.")
cb = ContinuousCallback(condition, affect!)

sol = solve(jprob, FunctionMap(); callback=cb)

That the callback should be activated can be checked through a modified affect function:

function condition(u, t, integrator)
    if (u[1] - 25) == 0
        println("Condition is true")
    end
    return u[1] - 25
end

or printing the solution:

using Plots
plot(sol)

E.g. using condition(u, t, integrator) = Float64(u[1] - 25) does not fix the problem.

Using something like

condition(u, t, integrator) = (t - 25)

and it works, which suggests that there is a problem with the discrete nature of these callbacks.

In practice, for JumpSimulations, a DiscreteCallback could be used here. However, this is basically a direct translation of a ContiniousCallback used for ODEs (it came up in the implementation of SBMLImports for JumpProblems). Is there a problem that these callbacks fail in the discrete case? It is rather niche, but it might still be worth pointing out in the docs, however, not sure if it is something we actually expect anyone to encounter again.

[isaacsas]

I don't think this is a JumpProcesses issue -- this seems an issue with using FunctionMap (which is not part of JumpProcesses).

[TorkelE]

Got it, didn' realise FunctionMap was outside JumpProcessess (not that I realized that FunctionMap was the issue either. I could use raise this at DiffEq instead.

[isaacsas]

Can you hand construct this example and get it to work with FunctionMap instead of going via Catalyst/MTK? That would help in seeing if FunctionMap is just not handling continuous callbacks correctly.

[isaacsas]

If that also doesn't work the issue should be moved to OrdinaryDiffEq (I can do that if it still doesn't work in that case).

[ChrisRackauckas]

I'm not what or how a continuouscallback would be defined for a FunctionMap/DiscreteProblem. There isn't even a continuous space for it to be evaluated on. We should just error if any are passed.

[isaacsas]

I think maybe the confusion is coming from the FunctionMap documentation here:

https://docs.sciml.ai/DiffEqDocs/stable/solvers/discrete_solve/

which gives the impression FunctionMap is for use with DiscreteProblems but supports continuous callbacks. Perhaps these should be split up, or if the intended use is FunctionMap with ODEProblem to support continuous callbacks that added?

[isaacsas]

We do, however, need to provide a nicer way to build continuous jump problems via Catalyst/MTK. There isn't a nice interface on it now for when one has VariableRateJumps or continuous callbacks. We should have an ODEProblem(::JumpSystem,...) in MTK that builds a dummy f!(du,u,p,t) and can then feed into JumpProblem(::JumpSystem, ::ODEProblem,...).

[isaacsas]

In fact, maybe we should a JumpProblem(::JumpSystem, aggregator, ...) in MTK, which would analyze the system and build a DiscreteProblem or ODEProblem as appropriate (i.e. depending on the presence of VariableRateJumps or ContinuousCallbacks)? Then we can drop the need to explicitly construct a DiscreteProblem or ODEProblem for users. Catalyst could then have a dispatch of JumpProblem(::ReactionSystem, aggregator) and drop the DiscreteProblem step.

edit: Of course, we'd have to figure out how to pass (u0, tspan, p) with this new interface.

[ChrisRackauckas]

In fact, maybe we should a JumpProblem(::JumpSystem, aggregator, ...) in MTK, which would analyze the system and build a DiscreteProblem or ODEProblem as appropriate (i.e. depending on the presence of VariableRateJumps or ContinuousCallbacks)? Then we can drop the need to explicitly construct a DiscreteProblem or ODEProblem for users. Catalyst could then have a dispatch of JumpProblem(::ReactionSystem, aggregator) and drop the DiscreteProblem step.

There is the growing DEProblem/System syntax which is fully automated, and yes I think we need to just build that out more and more.

But I think that's kind of different. condition(u,t,integrator) = u[1] - 0.5. The chemical species A goes from 0 to 1 in the interval. At what time t was there 1/2? That's just an ill-defined question. Either it doesn't have a sensible answer, or you say it happened at the time of the transition. But if it happened at the time of the transition, then it's always at a step time so its a DiscreteCallback.

[ChrisRackauckas]

It's not that it will only return an integer, it's that even if it does trigger, it triggers at exactly a step point, which is a discrete callback not a continuous callback. So it's either not sensical or it's a discrete callback.

[isaacsas]

Yeah, that’s why I removed those comments. I couldn’t think of a use case where a continuous callback with a pure jump process makes sense.

[isaacsas]

@TorkelE can we close this one? It sounds like whatever callbacks are being created are incorrectly being made as continuous when they should be discrete.

[TorkelE]

If you are happy with it I am happy, just figured I should bring it up when I encounter it.

[TorkelE]

Actually, one last concern/thing I would like to check.

Say that we want to define a cell division model, where we have a cell division even once a quantity X goes above a threshold. Tx. We can define this as a callback/even via MTK standard. Now, in MTK we need to define it as either a DiscreteEvent or a ContiniousEvent. Do we not encounter and issue that we here have to decide which to use depending on whether we want to simulate out model using ODE or Jumps?

The model is exactly the same, so what the event represent is no different. But due for simulation technical reasons we have to decide which type of simulation the model is intended for. This would basically also be the only case. Would it be possible/make sense to have continuous events treated as discrete events (with no loss of generality) for Jumps? I do think there might be something here we should think about.

[isaacsas]

I'm going to close this as this is a modeling choice a user has to make (or perhaps Catalyst/MTK can be made smarter), but I don't think this is something to handle at JumpProcesses level.

[isaacsas]

Feel free to open an issue on MTK or Catalyst to discuss / keep track of this there.