Hello,

I have a question about the behaviour of callbacks using DifferentialEquations.jl, because I find the behaviour of DiscreteCallbacks counterintuitive.

Consider the following code, i.e. an adaptation of this:

```
using DifferentialEquations
using JumpProcesses
using Plots
using Random
using StatsBase
u_thr = 2.
u_base = 1.
f(u, p, t) = 1.01 * u
function continuous_condition(u, t, integrator)
u - u_thr
end
function continuous_affect!(integrator)
integrator.u = u_base
push!(continuous_crosses, integrator.t)
end
function discrete_condition(u, t, integrator)
u .> u_thr
end
function discrete_affect!(integrator)
integrator.u = u_base
push!(discrete_crosses, integrator.t)
end
#jump process
rate(u, p, t) = 1. #in Hz
function affect!(integrator)
integrator.u += sample([0., 1.])
end
function find_greater_equal_thresh(sol)
spikes = findall(sol.u .>= u_thr)
spiketimes = sol.t[spikes]
return spikes, spiketimes
end
begin
continuous_crosses = Float64[]
discrete_crosses = Float64[]
Random.seed!(2)
spike_cb = ContinuousCallback(continuous_condition, continuous_affect!; rootfind = SciMLBase.RightRootFind);
discrete_condition_cb = DiscreteCallback(discrete_condition, discrete_affect!);
cbs = CallbackSet(spike_cb, discrete_condition_cb);
spike_cb = ContinuousCallback(continuous_condition, continuous_affect!; rootfind = SciMLBase.RightRootFind);
discrete_condition_cb = DiscreteCallback(discrete_condition, discrete_affect!)
cbs = CallbackSet(spike_cb, discrete_condition_cb)
u0 = 3 / 2
tspan = (0.0, 1.5)
prob = ODEProblem(f, u0, tspan)
jump = ConstantRateJump(rate, affect!)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, ImplicitEuler(), callback = cbs, save_everystep = false);
spikes, spiketimes= find_greater_equal_thresh(sol)
#plot solution
p = plot(sol.t, sol.u, linewidth = 2, title = "Continuous vs discrete callbacks",
xlabel = "Time (s)", alpha = 0.3)
hline!([u_thr], label = "Threshold", linestyle = :dash, linewidth = 0.7)
vline!(discrete_crosses, label = "Discrete callback", linestyle = :dash, linewidth = 0.7)
vline!(continuous_crosses, label = "Continuous callback", linestyle = :dash, linewidth = 0.7)
display(p)
savefig(p, "08042024/outputs/callback-differences.pdf")
end
begin
println("number of discrete callbacks: ", length(discrete_crosses))#1
println("number of continuous callbacks: ", length(continuous_crosses))#2
println("number of time frames spent at or above threshold: ", length(findall(sol.u .>= u_thr)))#3
println("times when u is above threshold: ", sol.t[findall(sol.u .>= u_thr)])#4
println("timesteps when u is above threshold: ", findall(sol.u .>= u_thr))#5
println("times of discrete callbacks: ", discrete_crosses)#6
println("times of continuous callbacks: ", continuous_crosses)#7
end
```

The code produces the following graph:

**The graph suggests that there are 3 crossings of the threshold, and therefore there should be 3 timesteps at which u > u_thr holds. However, this is not the case because there are 5!**

Here are the outputs of the print statements:

number of discrete callbacks: 2 #1

number of continuous callbacks: 1 #2

number of time frames spent at or above threshold: 5 #3

times when u is above threshold: [0.007120977006082045, 0.007120977006082045, 0.5984251769043971, 0.5984251769043971, 1.2748637298022953] #4

timesteps when u is above threshold: [3, 4, 7, 8, 10] #5

times of discrete callbacks: [0.007120977006082045, 0.5984251769043971] #6

times of continuous callbacks: [1.2748637298022953] #7

The print statements say that despite there being only three callback calls (three crossings) the dynamical variable is not reset immediately after a `DiscreteCallback`

call. Instead, you see that a `DiscreteCallback`

call always has two time steps at which `u`

is above threshold (i.e. fulfils the callback condition to be reset).

Instead, when `ContinuousCallback`

is called, the dynamical variable `u`

is immediately reset to `u_base`

. **After a DiscreteCallback the dynamical variable stays above threshold instead of being reset at the next timestep. Why is this?**

On the one hand I feel that this is fine detail, but this plays a large role when analysing the solution of an `ODEProblem`

after the fact. If I want to find the timepoints where `u`

crossed threshold, I have to consider if the crossing was called by a `ContinuousCallback`

or a `DiscreteCallback`

, which I cannot know after the code is run.

Note: I need both a continuous and a discrete callback here because there are two types of dynamics of `u`

(I think): if `u`

crosses threshold because of the normal evolution of the differential equation, a `ContinuousCallback`

will detect its crossing the threshold. Instead, if the crossing is caused by a jump from the JumpProcess, I need a `DiscreteCallback`

.