Continious Callback does not work when `condition` is not properly written as a function of time

[TorkelE]

This came up in SBMLImporter, think the bug only appeared this week, so was likely introduced relative recently.

Here the callback should be triggered, but never is.

using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

# Create ODEProblem
@parameters k1 k2 C
@variables S1(t) S2(t) S3(t)
eqs = [
    D(S1) ~ - (k1 / C) * S1 * S2 + k2 * S3,
    D(S2) ~ - (k1 / C) * S1 * S2 + k2 * S3,
    D(S3) ~ (k1 / C) * S1 * S2 - k2 * S3,
]
@mtkbuild model = ODESystem(eqs, t)
u0 = [S1 => 1.0, S2 => 2.0, S3 => 1.0]
ps = [k1 => 0.75, k2 => 0.25, C => 1.0]
ode = ODEProblem(model, u0, (0.0, 2.0), ps)

# Create callabck.
function myaffect!(integrator)
    println("Activated")
    integrator[S2] = 1 * integrator.ps[C]
    return nothing
end
function mycond(u, t, integrator)
    return integrator[S1] / integrator.ps[C] - 0.75
end
cb = ContinuousCallback(mycond, myaffect!)

# Solve problem.
sol = solve(ode; callback = cb)
plot(sol; idxs = S1/C - 0.75)

[ChrisRackauckas]

This is a user error because the condition is only written to evaluate the end point.

function mycond(u, t, integrator)
    return u[1] / integrator.ps[C] - 0.75
end
cb = ContinuousCallback(mycond, myaffect!)

# Solve problem.
sol = solve(ode, Tsit5(); callback = cb)
plot(sol; idxs = S1/C - 0.75)

If you properly use the u (or integrator(t)[1]) to evaluate at the t point then it's fine. "stopped working" working is not correct here: this would have never worked because it would never split signs because it always evaluates to a single condition.

[ChrisRackauckas]

Note that we purposely do not document because you have to be careful to evaluate w.r.t. t which isn't natural in the lowered form when u is a vector since it loses the symbolic interface.

[TorkelE]

Wait, I think I remember this from back in a while. Continious callback's integrator field does not work properly with MTK, that is the issue, right? It was designed pre-MTK where one have to use u. However, since MTK does not support index-based indexing, one cannot use u either.

This was something that we wanted to support long-term though, as an alternative to adding symbolic events directly to the model?

In the meantime e, wouldn't it make sense to throw an error?

[ChrisRackauckas]

Wait, I think I remember this from back in a while. Continious callback's integrator field does not work properly with MTK, that is the issue, right?

It's not that it doesn't work properly, it's that the integrator is the same in every call. This is always true for the ODE interface. It's only passed there in case you need to access properties like integrator.p. But integrator.u is the same in every call of the rootfind. Thus you cannot rootfind using integrator.u because it will never hit span over zero. You would have to use integrator(t). t is what is changed between calls to the rootfind, and u = integrator(t) is given as the thing to use with it.

t was designed pre-MTK where one have to use u. However, since MTK does not support index-based indexing, one cannot use u either.

Yes, the callback system was pre-MTK so u = integrator(t) "just use u and you're good" is a good advice, until MTK comes along and u is a vector that doesn't have symbolic indexing. What we should do is make integrator(t; idxs = :S1) work well with the integrator, but that's a separate issue to open. Even then though, there are major other reasons to use the MTK symbolic events.

This was something that we wanted to support long-term though, as an alternative to adding symbolic events directly to the model?

I don't see how we would? The ODE solver cannot do without doubling the full memory which would break PDEs. Also, you really need to use symbolic events for other reasons. The MTK symbolic events are doing a lot more than just being a nice interface, they are kind of required in order to do DAEs correctly. So you're getting a worse interface if you avoid it and then you also have to do some other numerical tricks which I assume most users (even you) don't even realize are required in order to make this robust in terms of reinitialization.

There's a lot that can be done and continue to improve here. For example, Shampine has a whole discussion about how to improve the robustness of event handling which requires modifying and adding ODEs based on the potential event functions: https://dl.acm.org/doi/10.1145/103147.103149 . It's a good read that's a simple introduction into the kinds of things we can do with symbolic-numeric algorithms that we cannot do if people black-box the condition function. So more and more, like for other things in MTK, if you use the symbolic interfaces you get more back: most robustness, speed, etc. because if you black-box we cannot do the fancier transformations. For this reason, we are trying as much as we can to push people towards the interfaces that are easier, safer, and give us more information about what's being solved: it's just better for everyone.

The most fallback case we document is ImperativeEffect, which still has the symbolic conditions, and that's on purpose because it both solves the problem here while giving a full imperative affect! function, which is a good enough middle ground for things like complex clocking and real-world hardware drivers.

In the meantime e, wouldn't it make sense to throw an error?

Since this has been the standard documented behavior of the event handling system since 2016, no. It's an oddity that the integrator is integrator.t == final t always and that you have to integrator(t) but it has always been the documented behavior and I don't see how that can change. It's pretty fundamental to the design because otherwise you'd need to copy the integrator, which in many circumstances would just blow the memory and you'd go OOM. So it's not entirely possible to make it so that integrator is correct time, and there would be a lot of adverse consequences to it.

If there was any change, it could be potentially to do condition(t, integrator) instead and force the user to integrator(t). But then that's such a minor optimization in most cases but leads to an interface where it's much easier to do things wrong. Also, I don't like doing breaking changes unless they are absolutely necessary, and this doesn't seem necessary.

I just can't see how we could make that integrator automatically update without all sorts of unintended consequences.

[TorkelE]

Sorry, I am unclear (and wrong!).

u is what should be used, agreed. However, what we would want is to make MTK-style symbolic indexing work on u? (i.e. right now that requires using a integer index, which we strongly should advice against) That should be technically possible, right? (although maybe not at the top of our list of priorities)

[ChrisRackauckas]

That's technically possible but not easy, since it's just built with integrator(t), so it uses the same u type as the integrator, which we want to be a vector since canonicalization ensures we have top performance and minimal recompilation. We could add an interface to the event handler so that the u is then wrapped with some other type information, but that would need something new. So... I'd hope people just use the callback interfaces? That would also solve other problems as I mentioned above about robustness and DAE reinitialization. So we could do something along these lines but we'd also not want users to really use that 😅 so it would be lower priority. But yes, it is possible.

[TorkelE]

Sounds good.

Yes, I agree that the native MTK event interface should cover most practical applications, so this is not really something we need very soon. It might be good to have a FAQ somewhere how to find the correct index for people who need it, and then at some point when more urgent things are sorted out it would be good to have MTK interfacing for u in place.

[ChrisRackauckas]

Yes, should probably bring this back up with the v10, as then describing what you have to do manually to make "as good of a callback" is a bit more concrete (and also will be a better pitch for, "please don't" especially in the DAE sense 😅)