Merge branch 'develop' into fix_2430_simultaneous_ea

documentation/python_interface.rst

@@ -26,7 +26,7 @@ AMICI can import :term:`SBML` models via the
 Status of SBML support in Python-AMICI
 ++++++++++++++++++++++++++++++++++++++
 
-Python-AMICI currently **passes 1247 out of the 1821 (~68%) test cases** from
+Python-AMICI currently **passes 1252 out of the 1821 (~68%) test cases** from
 the semantic
 `SBML Test Suite <https://github.com/sbmlteam/sbml-test-suite/>`_
 (`current status <https://github.com/AMICI-dev/AMICI/actions>`_).

python/sdist/amici/sbml_import.py

@@ -561,9 +561,9 @@ def _build_ode_model(
         assert dxdt.shape[0] - len(self.symbols[SymbolId.SPECIES]) == len(
             self.symbols.get(SymbolId.ALGEBRAIC_STATE, [])
         ), (
-            self.symbols[SymbolId.SPECIES],
+            self.symbols.get(SymbolId.SPECIES),
             dxdt,
-            self.symbols[SymbolId.SPECIES],
+            self.symbols.get(SymbolId.ALGEBRAIC_STATE),
         )
 
         # correct time derivatives for compartment changes
@@ -954,6 +954,7 @@ def _process_species(self) -> None:
             }
 
         self._convert_event_assignment_parameter_targets_to_species()
+        self._convert_event_assignment_compartment_targets_to_species()
         self._process_species_initial()
         self._process_rate_rules()
 
@@ -1564,6 +1565,36 @@ def _convert_event_assignment_parameter_targets_to_species(self):
                 "dt": sp.Float(0),
             }
 
+    def _convert_event_assignment_compartment_targets_to_species(self):
+        """Find compartments that are event assignment targets and convert
+        those compartments to species."""
+        for event in self.sbml.getListOfEvents():
+            for event_assignment in event.getListOfEventAssignments():
+                if event_assignment.getMath() is None:
+                    # Ignore event assignments with no change in value.
+                    continue
+                variable = symbol_with_assumptions(
+                    event_assignment.getVariable()
+                )
+                if variable not in self.compartments:
+                    continue
+                if variable in self.symbols[SymbolId.SPECIES]:
+                    # Compartments with rate rules are already present as
+                    # species
+                    continue
+
+                self.symbols[SymbolId.SPECIES][variable] = {
+                    "name": str(variable),
+                    "init": self.compartments[variable],
+                    # 'compartment': None,  # can ignore for amounts
+                    "constant": False,
+                    "amount": True,
+                    # 'conversion_factor': 1.0,  # can be ignored
+                    "index": len(self.symbols[SymbolId.SPECIES]),
+                    "dt": sp.Float(0),
+                }
+                del self.compartments[variable]
+
     @log_execution_time("processing SBML events", logger)
     def _process_events(self) -> None:
         """Process SBML events."""
@@ -1613,7 +1644,9 @@ def get_empty_bolus_value() -> sp.Float:
             # parse the boluses / event assignments
             bolus = [get_empty_bolus_value() for _ in state_vector]
             event_assignments = event.getListOfEventAssignments()
-            compartment_event_assignments = set()
+            compartment_event_assignments: set[tuple[sp.Symbol, sp.Expr]] = (
+                set()
+            )
             for event_assignment in event_assignments:
                 variable_sym = symbol_with_assumptions(
                     event_assignment.getVariable()
@@ -1631,7 +1664,7 @@ def get_empty_bolus_value() -> sp.Float:
                         "Could not process event assignment for "
                         f"{str(variable_sym)}. AMICI currently only allows "
                         "event assignments to species; parameters; or, "
-                        "compartments with rate rules, at the moment."
+                        "compartments."
                     )
                 try:
                     # Try working with the formula now to detect errors
@@ -1644,19 +1677,19 @@ def get_empty_bolus_value() -> sp.Float:
                         "expressions as event assignments."
                     )
                 if variable_sym in concentration_species_by_compartment:
-                    compartment_event_assignments.add(variable_sym)
+                    compartment_event_assignments.add((variable_sym, formula))
 
                 for (
                     comp,
                     assignment,
                 ) in self.compartment_assignment_rules.items():
                     if variable_sym not in assignment.free_symbols:
                         continue
-                    compartment_event_assignments.add(comp)
+                    compartment_event_assignments.add((comp, formula))
 
             # Update the concentration of species with concentration units
             # in compartments that were affected by the event assignments.
-            for compartment_sym in compartment_event_assignments:
+            for compartment_sym, formula in compartment_event_assignments:
                 for species_sym in concentration_species_by_compartment[
                     compartment_sym
                 ]:

python/tests/test_events.py

@@ -15,6 +15,7 @@
     create_amici_model,
     create_sbml_model,
 )
+from numpy.testing import assert_allclose
 
 
 @pytest.fixture(
@@ -746,3 +747,85 @@ def test_handling_of_fixed_time_point_event_triggers():
         assert (rdata.x[(rdata.ts >= 3)] == 3).all()
 
         check_derivatives(amici_model, amici_solver, edata=None)
+
+
+def test_multiple_event_assignment_with_compartment():
+    """see https://github.com/AMICI-dev/AMICI/issues/2426"""
+    ant_model = """
+    model test_events_multiple_assignments
+        compartment event_target = 1
+        event_target' = 0
+        species species_in_event_target in event_target = 1
+        unrelated = 2
+
+        # use different order of event assignments for the two events
+        at (time > 5): unrelated = 4, event_target = 10
+        at (time > 10): event_target = 1, unrelated = 2
+    end
+    """
+    # watch out for too long path names on windows ...
+    module_name = "tst_mltple_ea_w_cmprtmnt"
+    with TemporaryDirectory(prefix=module_name, delete=False) as outdir:
+        antimony2amici(
+            ant_model,
+            model_name=module_name,
+            output_dir=outdir,
+            verbose=True,
+        )
+        model_module = amici.import_model_module(
+            module_name=module_name, module_path=outdir
+        )
+        amici_model = model_module.getModel()
+        assert amici_model.ne == 2
+        assert amici_model.ne_solver == 0
+        assert amici_model.nx_rdata == 3
+        amici_model.setTimepoints(np.linspace(0, 15, 16))
+        amici_solver = amici_model.getSolver()
+        rdata = amici.runAmiciSimulation(amici_model, amici_solver)
+        assert rdata.status == amici.AMICI_SUCCESS
+        idx_event_target = amici_model.getStateIds().index("event_target")
+        idx_unrelated = amici_model.getStateIds().index("unrelated")
+        idx_species_in_event_target = amici_model.getStateIds().index(
+            "species_in_event_target"
+        )
+
+        assert_allclose(
+            rdata.x[(rdata.ts < 5) & (rdata.ts > 10), idx_event_target],
+            1,
+            rtol=0,
+            atol=1e-15,
+        )
+        assert_allclose(
+            rdata.x[(5 < rdata.ts) & (rdata.ts < 10), idx_event_target],
+            10,
+            rtol=0,
+            atol=1e-15,
+        )
+        assert_allclose(
+            rdata.x[(rdata.ts < 5) & (rdata.ts > 10), idx_unrelated],
+            2,
+            rtol=0,
+            atol=1e-15,
+        )
+        assert_allclose(
+            rdata.x[(5 < rdata.ts) & (rdata.ts < 10), idx_unrelated],
+            4,
+            rtol=0,
+            atol=1e-15,
+        )
+        assert_allclose(
+            rdata.x[
+                (rdata.ts < 5) & (rdata.ts > 10), idx_species_in_event_target
+            ],
+            1,
+            rtol=0,
+            atol=1e-15,
+        )
+        assert_allclose(
+            rdata.x[
+                (5 < rdata.ts) & (rdata.ts < 10), idx_species_in_event_target
+            ],
+            0.1,
+            rtol=0,
+            atol=1e-15,
+        )