Detection of simultaneous events with VectorContinuousCallback

[KalelR]

Hi

the VectorContinuousCallback approach cannot currently deal with events that occur at the same time. If 2+ simultaneous events occur, the callback tells the solver to keep trying to step between them with dt = 0.0 and so the code doesn't halt. As I understood, this is because one of the events is detected, and the integrator steps to that event time. But then it detects the other simultaneous event, and tries to step to that time, which is the same (so dt = 0). When it does so, it detects the first even again, and repeats this cycle.

A simple example for this is for two cosine oscillations with the same frequency. Events occur when the oscillation up-crosses a threshold, with the effect of that event time being pushed to an array or event times (not efficient, but just for the example).
If oscillations' initial conditions are different, they cross the threshold at different times, and the code works perfectly. Even if the ics are very similar and the crossing times are very similar.
Making the ics identical, however, leads to the error, which is silent (if you try to run, just make sure to cancel after some point so the memory usage doesn't explode!).

using OrdinaryDiffEq
function oscilate!(du, u, p, t)
    ω1 = p[1]; ω2 = p[2];
    du[1] = cos(ω1*t)
    du[2] = cos(ω2*t)
end

function condition!(out, u, t, integrator)
    uth = integrator.p[3]
    for i in eachindex(out)
        out[i] = u[i] - uth
    end
end

function affect!(integrator, idx_event)
    crosstimes = integrator.p[4]
    push!(crosstimes[idx_event], integrator.t)
end

# uth = 0.1; ω1 = pi/6; ω2 = pi/6; u0s = [-0.8, 0]; case = "cos-example-faraway-$ω1-$ω2-$u0s"
# uth = 0.1; ω1 = pi/6; ω2 = pi/6; u0s = [-0.8, -0.7]; case = "cos-example-closer-not-sim-$ω1-$ω2-$u0s-uth_$uth"
uth = 0.1; ω1 = pi/6; ω2 = pi/6; u0s = [-0.8, -0.8]; case = "cos-example-sim-$ω1-$ω2-$u0s-uth_$uth"

cb = VectorContinuousCallback(condition!, affect!, 2; affect_neg! = nothing);
T = 50
crosstimes = [Float64[] for i=1:2];
prob = ODEProblem(oscilate!, u0s,  (0, T), [ω1, ω2, uth, crosstimes])
sol = solve(prob, Tsit5(); callback=cb)

using CairoMakie
ts = 0:0.01:T
sol_ds = reduce(hcat, sol.(ts))
fig = Figure()
ax = Axis(fig[1,1])
colors = [:blue, :red]
for i=1:2
    lines!(ax, ts, sol_ds[i,:], color=colors[i])
    length(crosstimes[i]) > 0 && scatter!(ax, crosstimes[i], [uth for i in crosstimes[i]], color=colors[i])
end
fig

This issue comes in a similar vein to Issue #519, by @freemin7, where @ChrisRackauckas suggested a solution.

I think the solution here may be simpler though. And also maybe it should be implemented for all cases since the user might not know whether events can occur simultaneously. For instance, synchronization of events is an important research topic, e.g. in neuroscience, and its occurrence is not necessarily known a priori. This is what brought me to this issue, and I asked this question on discourse for my specific problem.

I'd really appreciate suggestions for how to solve this. I have one proposal which I'm unsure about. This would be to define a BitArray vector_events_idxs stored in CallbackCache and to modify the find_callback_time function to

1. detect when events occur simultaneously and store that in the array and
2. "nudge" all events that were found.

I think this solves the problem. But then what to do with affect!, since currently it receives an Int idx?
A way to sidestep the issue would be to keep affect! as it is, and just tell the user to use integrator.callback_cache.vector_event_idxs. This does not sound very clean, though. Ideas?

A related point is: for large systems (such as a typical use-case in neuroscience) maybe several events will be triggered either simultaneously or very close in time (e.g. when neurons spike synchronizedly). The code works for "very close in time", but it will have several very small steps between these events, as it jumps from one event to another just after. Would it be interesting to consider a faster approach, such as detecting an event, searching for other events within a window centered on the event, and detecting them all on one go and skip them? Please correct me if I'm wrong, but this would be similar to the proposal by @ChrisRackauckas, and it might speed up code in those circumstances.

[ChrisRackauckas]

Would it be interesting to consider a faster approach, such as detecting an event, searching for other events within a window centered on the event, and detecting them all on one go and skip them? Please correct me if I'm wrong, but this would be similar to the proposal by @ChrisRackauckas, and it might speed up code in those circumstances.

yup that's my proposal, so I'll close this as a duplicate but will work this example in as one of the tests when I get around to this, hopefully soon.