fix: event-time sensitivity guard misses state-dependent triggers reached through SBML (#52)

CHANGELOG.md

@@ -264,6 +264,28 @@ in `CMakeLists.txt`) is derived from it.
   fields.
 
 ### Fixed
+- **The event-time sensitivity guard missed state-dependent triggers reached
+  through SBML, answering those models instead of refusing them (issue #52).**
+  The guard refuses forward sensitivities when an event's crossing time depends
+  on a requested parameter, and it decides that from the trigger's *bound
+  addresses*. It compared against species concentrations only. But ModelBuilder
+  registers a species as an ExprTk variable only when its name is still free, and
+  SBML models routinely give each species an observable of the same name — so the
+  species registration is skipped and a trigger token binds to the observable
+  total, never to `&sp.concentration`. Every SBML state-dependent trigger
+  therefore slipped the guard and was answered, with the event contributions
+  missing entirely. On AMICI's `neuron` fixture (Izhikevich, trigger `v > 30`,
+  which names no parameter but whose crossing time depends on `a` and `b` through
+  the trajectory) the returned sensitivities were 6x–135x off, uniformly in one
+  direction, across all four parameters.
+
+  The guard now tests against every address that carries live state: species
+  concentrations, observable totals, and rateOf accessors. An observable total is
+  a linear functional of the state and a rateOf accessor is dx/dt, so a trigger
+  reading either has a non-zero `dt*/dp` exactly as a concentration read does.
+  Refusing is unchanged as a policy — this is only the coverage of what counts as
+  state — and the message now says which of the three it saw and why that implies
+  a moving crossing time.
 - **The codegen cache did not invalidate on a codegen change, so a fix could be
   silently inert on a warm cache (issue #51).** The `.net` path keys its compiled
   `.so` on the model content plus the hand-maintained `_CODEGEN_VERSION` constant

python/bngsim/_simulator.py

@@ -976,7 +976,11 @@ def _raise_if_event_sensitivities(self, param_names: list[str] | None = None) ->
         (:func:`NetworkModel.event_sensitivity_unsupported_reason`), which knows
         each event's persistence/delay and — via the trigger's referenced
         variables — whether it is fixed-time and whether its crossing time
-        depends on a requested sensitivity parameter. ``param_names`` is the set
+        depends on a requested sensitivity parameter. "Fixed-time" is judged
+        against every address that carries live state, not just species
+        concentrations: an observable total or a rateOf accessor moves with the
+        trajectory too, so a trigger reading one has a parameter-dependent
+        crossing time even though it names no parameter (issue #52). ``param_names`` is the set
         of parameters whose sensitivities this call requests (defaults to
         ``self._sensitivity_params``); an IC-only request passes an empty list,
         which still exercises the persistence/delay/state-dependence checks.

python/tests/test_event_sensitivity.py

@@ -74,6 +74,19 @@
   </model>
 </sbml>"""
 
+
+def _sbml_event_with_state_trigger() -> str:
+    """``SBML_EVENT`` with its fixed-time trigger swapped for a state-dependent
+    one (``S < 5``) — the shape of AMICI's ``neuron`` fixture, whose ``v > 30``
+    reads a state variable while naming no parameter (issue #52)."""
+    fixed_time = """<apply><geq/>
+              <csymbol encoding="text"
+                definitionURL="http://www.sbml.org/sbml/symbols/time">t</csymbol>
+              <cn>1</cn></apply>"""
+    assert fixed_time in SBML_EVENT, "SBML_EVENT trigger changed; update this helper"
+    return SBML_EVENT.replace(fixed_time, "<apply><lt/><ci>S</ci><cn>5</cn></apply>")
+
+
 # ── Discontinuity-trigger model: a piecewise-time forcing pulse on parameter
 # `inp` drives production of X. n_discontinuity_triggers > 0 but n_events == 0 —
 # the pulse breaks the integrator step yet never jumps state, so forward
@@ -234,6 +247,40 @@ def test_state_dependent_trigger_raises(self):
         with pytest.raises(ValueError, match="state-dependent"):
             sim.run(t_span=(0, 10), n_points=11)
 
+    def test_sbml_state_dependent_trigger_raises(self):
+        """Issue #52: the same refusal, reached through SBML.
+
+        The guard tests the trigger's *bound addresses*, and ModelBuilder
+        registers a species as an ExprTk variable only when the name is free.
+        SBML models routinely give each species an observable of the same name,
+        so the species registration is skipped and the trigger's token binds to
+        the observable total instead of ``&sp.concentration``. Checking
+        concentrations alone therefore saw no state dependence here and answered
+        the model — on AMICI's ``neuron`` fixture (Izhikevich, trigger
+        ``v > 30``) the sensitivities came back 6x-135x off, uniformly in one
+        direction, rather than being refused.
+        """
+        m = bngsim.Model.from_sbml_string(_sbml_event_with_state_trigger())
+        assert m._core.n_events == 1
+        # Precondition for the bug: species and observable share the name, which
+        # is what pushed the trigger's binding onto the observable total.
+        assert "S" in list(m._core.species_names)
+        assert "S" in list(m._core.observable_names)
+
+        sim = bngsim.Simulator(m, method="ode", sensitivity_params=["k"])
+        with pytest.raises(ValueError, match="state-dependent"):
+            sim.run(t_span=(0, 10), n_points=11)
+
+    def test_sbml_state_dependent_trigger_refused_for_every_entry_point(self):
+        """``compute_all_sensitivities`` takes the same guard, so it must refuse
+        the same model rather than quietly returning a tensor missing the event
+        contributions."""
+        m = bngsim.Model.from_sbml_string(_sbml_event_with_state_trigger())
+        with pytest.raises(ValueError, match="state-dependent"):
+            bngsim.Simulator(m, method="ode").compute_all_sensitivities(
+                t_span=(0, 10), n_points=11
+            )
+
     def test_delayed_event_raises(self):
         m, _ = _decay_with_event("2.0", delay=1.0)
         sim = bngsim.Simulator(m, method="ode", sensitivity_params=["k"])

src/model.cpp

@@ -449,11 +449,35 @@ std::optional<std::string> NetworkModel::event_sensitivity_unsupported_reason(
         sens_param_addrs.insert(&match->value);
     }
 
-    // State-variable addresses (species concentrations the trigger may read).
-    std::unordered_set<const double *> species_addrs;
-    species_addrs.reserve(species.size());
+    // Addresses that carry live state, i.e. anything whose value moves with the
+    // trajectory. A trigger reading one of these has a crossing time that
+    // depends on the sensitivity parameters *through the state*, even when the
+    // trigger names no parameter at all, so dt*/dp is non-zero exactly as for an
+    // explicitly parameter-dependent trigger.
+    //
+    // Issue #52: this deliberately covers more than species concentrations.
+    // ModelBuilder registers a species as an ExprTk variable only when its name
+    // is not already taken (model_builder.cpp), and SBML models routinely give
+    // each species an observable of the same name — so in `v > 30` the token `v`
+    // binds to the *observable total*, not to &sp.concentration. Checking
+    // concentrations alone therefore missed every SBML state-dependent trigger
+    // and answered those models instead of refusing them: on AMICI's `neuron`
+    // (Izhikevich, trigger `v > 30`) the returned sensitivities were 6x-135x off,
+    // in one direction, across every parameter.
+    //
+    // An observable total is a linear functional of the state and a rateOf
+    // accessor is dx/dt, so both move with the trajectory just as a
+    // concentration does.
+    std::unordered_set<const double *> state_addrs;
+    state_addrs.reserve(species.size() + impl_->observables.size() + impl_->current_derivs.size());
     for (const Species &sp : species) {
-        species_addrs.insert(&sp.concentration);
+        state_addrs.insert(&sp.concentration);
+    }
+    for (const Observable &obs : impl_->observables) {
+        state_addrs.insert(&obs.total);
+    }
+    for (const double &deriv : impl_->current_derivs) {
+        state_addrs.insert(&deriv);
     }
 
     for (const Event &ev : events) {
@@ -471,10 +495,13 @@ std::optional<std::string> NetworkModel::event_sensitivity_unsupported_reason(
         const std::vector<const double *> refs =
             eval.referenced_variable_addresses(ev.trigger_expr_idx);
         for (const double *addr : refs) {
-            if (species_addrs.count(addr) != 0) {
+            if (state_addrs.count(addr) != 0) {
                 return "event '" + id +
-                       "' has a state-dependent trigger (it reads a species concentration); "
-                       "only fixed-time triggers are supported for forward sensitivity so far "
+                       "' has a state-dependent trigger (it reads a species concentration, an "
+                       "observable, or a rate); its crossing time therefore depends on the "
+                       "sensitivity parameters through the trajectory even though the trigger "
+                       "names none of them, so the event-time sensitivity dt*/dp is non-zero. "
+                       "Only fixed-time triggers are supported for forward sensitivity so far "
                        "(GH #212 Phase 2).";
             }
         }