fix h, doc, rename

Author: dweindl

include/amici/model_dae.h

@@ -446,7 +446,7 @@ class Model_DAE : public Model {
      * @param x Vector with the states
      * @param p parameter vector
      * @param k constants vector
-     * @param h heavyside vector
+     * @param h Heaviside vector
      * @param dx Vector with the derivative states
      * @param w vector with helper variables
      */

models/model_nested_events_py/h.h

@@ -1,3 +1,3 @@
-#define injection h[0]
-#define Heaviside_1 h[1]
-#define Heaviside_2 h[2]
+#define Heaviside_1 h[0]
+#define Heaviside_2 h[1]
+#define injection h[2]

python/sdist/amici/de_model.py

@@ -7,6 +7,7 @@
 import logging
 import re
 from collections.abc import Callable, Sequence
+from contextlib import suppress
 from itertools import chain
 from operator import itemgetter
 from typing import TYPE_CHECKING
@@ -651,9 +652,12 @@ def num_events_solver(self) -> int:
         :return:
             number of event symbols (length of the root vector in AMICI)
         """
-        constant_syms = self._static_symbols(["k", "p", "w"])
+        # TODO(performance): we could include constant `x` here as well
+        #   (dx/dt = 0 AND not target of event assignments)
+        #   this will require passing `x` to `fexplicit_roots`
+        static_syms = self._static_symbols(["k", "p", "w"])
         return sum(
-            not event.has_explicit_trigger_times(constant_syms)
+            not event.has_explicit_trigger_times(static_syms)
             for event in self.events()
         )
 
@@ -1137,6 +1141,7 @@ def generate_basic_variables(self) -> None:
         Generates the symbolic identifiers for all variables in
         ``DEModel._variable_prototype``
         """
+        self.parse_events()
         self._reorder_events()
 
         for var in self._variable_prototype:
@@ -1163,11 +1168,6 @@ def parse_events(self) -> None:
         for expr in self._expressions:
             expr.set_val(self._process_heavisides(expr.get_val(), roots))
 
-        # remove all possible Heavisides from roots, which may arise from
-        # the substitution of `'w'` in `_collect_heaviside_roots`
-        for root in roots:
-            root.set_val(self._process_heavisides(root.get_val(), roots))
-
         # Now add the found roots to the model components
         for root in roots:
             # skip roots of SBML events, as these have already been added
@@ -1181,20 +1181,57 @@ def _reorder_events(self) -> None:
         Re-order events - first those that require root tracking,
         then the others.
         """
-        constant_syms = self._static_symbols(["k", "p", "w"])
+        # Currently, the C++ simulations relies on the order of events:
+        #  those that require numerical root-finding must come first, then
+        #  those with explicit trigger times that don't depend on dynamic
+        #  variables.
+        #  TODO: This re-ordering here is a bit ugly, because we already need
+        #   to generate certain model equations to perform this ordering.
+        #   Ideally, we'd split froot into explicit and implicit parts during
+        #   code generation instead (as already done for jax models).
+        static_syms = self._static_symbols(["k", "p", "w"])
+
+        # ensure that we don't have computed any root-related symbols/equations
+        #  yet, because the re-ordering might invalidate them
+        # check after `self._static_symbols` which itself generates certain
+        #  equations
+        if (
+            generated := set(self._syms)
+            | set(self._eqs)
+            | set(self._sparsesyms)
+            | set(self._sparseeqs)
+        ) and (
+            "root" in generated
+            or any(
+                name.startswith("droot") or name.endswith("droot")
+                for name in generated
+            )
+        ):
+            raise AssertionError(
+                "This function must be called before computing any "
+                "root-related symbols/equations. "
+                "The following symbols/equations are already "
+                f"generated: {generated}"
+            )
+
         self._events = list(
             chain(
                 itertools.filterfalse(
-                    lambda e: e.has_explicit_trigger_times(constant_syms),
+                    lambda e: e.has_explicit_trigger_times(static_syms),
                     self._events,
                 ),
                 filter(
-                    lambda e: e.has_explicit_trigger_times(constant_syms),
+                    lambda e: e.has_explicit_trigger_times(static_syms),
                     self._events,
                 ),
             )
         )
 
+        # regenerate after re-ordering
+        with suppress(KeyError):
+            del self._syms["h"]
+        self.sym("h")
+
     def get_appearance_counts(self, idxs: list[int]) -> list[int]:
         """
         Counts how often a state appears in the time derivative of

python/sdist/amici/exporters/sundials/de_export.py

@@ -914,14 +914,14 @@ def _get_create_splines_body(self):
     def _get_explicit_roots_body(self) -> list[str]:
         events = self.model.events()
         lines = []
-        constant_syms = self.model._static_symbols(["k", "p", "w"])
+        static_syms = self.model._static_symbols(["k", "p", "w"])
 
         for event_idx, event in enumerate(events):
             if not (
                 tigger_times := {
                     tt
                     for tt in event.get_trigger_times()
-                    if tt.free_symbols.issubset(constant_syms)
+                    if tt.free_symbols.issubset(static_syms)
                 }
             ):
                 continue

python/sdist/amici/importers/pysb/__init__.py

@@ -398,7 +398,6 @@ def ode_model_from_pysb_importer(
 
     _process_stoichiometric_matrix(model, ode, fixed_parameters)
 
-    ode.parse_events()
     ode.generate_basic_variables()
 
     return ode

python/sdist/amici/importers/sbml/__init__.py

@@ -700,9 +700,7 @@ def _build_ode_model(
         if hybridization:
             ode_model._process_hybridization(hybridization)
 
-        ode_model.parse_events()
         # substitute SBML-rateOf constructs
-        #  must be done after parse_events, but before generate_basic_variables
         self._process_sbml_rate_of(ode_model)
 
         # fill in 'self._sym' based on prototypes and components in ode_model

python/tests/test_pysb.py

@@ -328,7 +328,7 @@ def test_names_and_ids(pysb_example_presimulation_module):
 
 
 @skip_on_valgrind
-def test_heavyside_and_special_symbols():
+def test_heaviside_and_special_symbols():
     pysb.SelfExporter.cleanup()  # reset pysb
     pysb.SelfExporter.do_export = True
 

src/model_dae.cpp

@@ -95,7 +95,8 @@ void Model_DAE::froot(
 ) {
     std::ranges::fill(root, 0.0);
     auto const x_pos = compute_x_pos(x);
-    // TODO too costly? only dynamic expressions?
+    // TODO(performance) only dynamic expressions? consider storing a flag
+    //  for whether froot actually depends on `w`.
     fw(t, N_VGetArrayPointerConst(x_pos));
     froot(
         root.data(), t, N_VGetArrayPointerConst(x_pos),

src/model_ode.cpp

@@ -84,7 +84,8 @@ void Model_ODE::froot(
     realtype const t, const_N_Vector x, gsl::span<realtype> root
 ) {
     auto const x_pos = compute_x_pos(x);
-    // TODO too costly? only dynamic expressions?
+    // TODO(performance) only dynamic expressions? consider storing a flag
+    //  for whether froot actually depends on `w`.
     fw(t, N_VGetArrayPointerConst(x_pos));
     std::ranges::fill(root, 0.0);
     froot(