Repository navigation
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #3439 +/- ##
==========================================
+ Coverage 80.11% 80.12% +0.01%
==========================================
Files 331 331
Lines 23325 23346 +21
Branches 1619 1624 +5
==========================================
+ Hits 18686 18707 +21
Misses 4633 4633
Partials 6 6
Flags with carried forward coverage won't be shown. Click here to find out more.
🚀 New features to boost your workflow:
|
ee59db3 to
6f0b8b4
Compare
288dd7a to
4db3e38
Compare
dweindl
added
events
74e54ed to
2277834
Compare
FFroehlich
left a comment
FFroehlich
left a comment
There was a problem hiding this comment.
I wonder whether we are taking the right approach here as I have the impression that we are just piling more and more complexity on this problem and it just keeps on failing in novel way.
The overall issue might be better solved by porting the root masking from the jax backend.
Piling up complexity - yes. And we won't be able to handle every adversarial case either way. However, I am not convinced that the JAX backend is currently on par in terms of handling simultaneous events and that it would be that much simpler there to implement. We could give the root masking a try, but I think most the recently introduced functionality/complexity would be required either way. I would prefer fixing all those other event-related bugs first, though. |
I would expect that the JAX hopefully isn't subject to as many issues due to it's restricted feature set, even though I wouldn't be suprised if there are many undiscovered bugs remaining. Specifically for this PR, root masking + tracking root derivative and event-induced root sign changes should be a substantially more robust solution. |
Yes, that was my point. Once the same feature set is supported there, I am not sure this will be substantially less complex.
I will push the alternative implementation shortly, then we can compare and decide. |
|
I opened two PRs with two alternative implementations:
Both PRs go beyond Diff: https://github.com/dweindl/AMICI/compare/root_filtering..root_masking |
|
The diff is tiny, but okay, let's go for the masking, which avoids additional messing with solver internals. |
Maybe I am not as up to date with the code base as I would need to be, but I think this would only be the first part of the refactoring I had in mind. The JAX backend works with pairs of trigger functions and turns off the rootfinding for the root that triggered until the other one triggers (gated by direction). This means that unless the root function is rapidly oscillation (I don't think we can address that case), any spurious re-triggering should be avoided. |
I might be missing the point. Where would it make a difference whether we add a positive and a negative root function (as we currently do for Heavisides), or track only one direction of a single one? Which part of the logic do you think would become obsolete?
Do you have an example in mind that would break #3490 / #3491 but work with the paired trigger functions? |
You can turn off event detection/root finding after the first triggering and only reactivate once the other direction triggers. I would hope that would render most of the re-triggering checks we have right now obsolete. The expected advantage would be that we don’t have to define any noise thresholds (either root function noise or trigger point closeness). Happy to give it a stab but I am already running pretty low on Claude tokens for the rest of the week. |
After a quick conversation with ChatGPT I concede my point. I think with #3490 in place there wouldn’t be much scope for improvement. |
`detect_secondary_events` skipped all roots that were crossed in the current discontinuity, as their values are within roundoff of zero. If an event assignment moved such a root back across zero, e.g., for `at x < thr: x = 1.5 * x`, the Heaviside variable was not updated. It was then incremented again at each subsequent crossing, giving wrong `piecewise` branches. For explicit trigger times, e.g., `at time >= s: s = s + 1`, the event did not fire again (since AMICI-dev#3407, which takes the crossing direction of explicit roots from the Heaviside variables). Instead of comparing the new root values to those at the trigger, only skip root functions that the preceding assignment did not change, and compare the others to the current Heaviside variables. Roots crossed at the current time are only considered crossed back if their value exceeds the change of the root function for a state perturbation by the integration tolerances, so that roundoff or the re-computation of algebraic states doesn't flip them. Exact zeros and NaNs don't count as crossings. Besides an assignment exactly onto a threshold, test one that lands one rounding error below a threshold computed at runtime (`0.1 + 0.2`), which reverted the switch without the noise threshold. Remove the unused `heaviside_differs`. Fixes AMICI-dev#3425. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2277834 to
750547a
Compare
…v#3480) If an event assignment crossed a `piecewise` condition in the negative direction without triggering another event, the jump in the right-hand side was missing from the forward and adjoint sensitivities. A typical case is a flag that is switched off by an event, e.g., `s = 1`, `z' = piecewise(c, s >= 1, 0)`, `E: at x >= 1: s = 0`. For each Heaviside function, there are root functions `r` and `-r`. Usually, the negative crossing of `r` is accompanied by a positive crossing of `-r`, which queues an (assignment-free) event whose sensitivity update includes the right-hand side jump. If the assigned variable was exactly at the threshold, `-r` goes from 0 to a positive value, which is no change of its Heaviside value, so no event is queued, and `update_heaviside` flips `h` after the last sensitivity update. Now, if no event is pending after the Heaviside update of a secondary discontinuity, the jump is applied separately: for FSA, `sx += (xdot_post - xdot_pre) * stau`; for ASA, the jump is stored in the `Discontinuity`, and `xB^T (xdot_post - xdot_pre)` is added to `stauB`, which is passed on to the triggering discontinuity via `stauB_secondary_`. For DAEs, the algebraic states are made consistent with the switched equations first (as after event applications); otherwise, a switch in an algebraic rule was still missing. The `SECONDARY_EVENT` warning no longer mentions this case. The new test checks FSA and ASA against the analytical gradient for the example above, the equivalent `piecewise(0, s < 1, c)`, the assignment in a secondary event, and the switch in an algebraic rule of a DAE. Not changed: for DAEs, if another event is pending after the Heaviside update, the algebraic states are only made consistent after that event's application, as before, so its assignments see algebraic states that are consistent with the old Heaviside values. Fixes AMICI-dev#3436. Possible textual conflicts with AMICI-dev#3439 (`detect_secondary_events`). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
Fix the Heaviside variables of roots that an event assignment crosses back at the same time. Fixes #3425.
detect_secondary_eventsskipped all roots that had been crossed in the current discontinuity, since their values are within roundoff of zero and the sign isn't reliable. If an event assignment moved such a root back across zero, e.g.,E1: at x < thr: x = 1.5 * x, the Heaviside variable kept the post-crossing value, and the next crossing incremented it again (h= 2, 3, …), sopiecewiseexpressions took the wrong branch. In the issue model, the simulation eventually failed withAMICI_TOO_MUCH_WORK. Similarly, an event moving its own explicit trigger time,E: at time >= s: s = s + 1, fired only once. This is a regression since #3407, which takes the crossing direction of explicit roots from the Heaviside variables. It affects, e.g., periodic dosing schedules written as one event.Now, root functions that the preceding event assignment didn't change are skipped, and the others are compared to the current Heaviside variables instead of to the root values at the trigger. A root crossed back gets
roots_found = -1(or+1for its negated counterpart, which is queued as an event without assignments, as for other secondary roots). Roots that were crossed at the current time are only considered crossed back if they are clearly on the other side: their value has to exceed the change of the root function for a state perturbation by the integration tolerances. This guards against spurious flips due to roundoff or the re-computation of algebraic states after an event (DAEs). For these roots, exact zeros and NaNs are never considered crossings.Tests:
test_event_crosses_own_trigger_back: the issue model against the analytical solution, the right-hand side at the final time point, the forward sensitivities against the analytical ones (includingdx/dthr = 0, which requires the event-time terms to cancel), and the forward and adjoint gradients.test_event_moves_own_trigger_time: the explicit-trigger-time model.test_event_assignment_onto_threshold: an assignment that sets a just-crossed root to zero, exactly (at x > xmax: x = xmax) or one rounding error below the threshold (xmax := 0.1 + 0.2computed at runtime,x = 0.3). The latter reverts the switch without the threshold.The
SECONDARY_EVENTwarning is now also logged for events that cross their own trigger back, e.g., at each execution ofE: at time >= s: s = s + T, x = x + Dwith forward or adjoint sensitivities. There, it is a false positive: the sensitivities are correct.🤖 Generated with Claude Code