Fix handling of event assignments to compartments (#2428)

Previously, in some cases (see #2426), event assignments to compartments would result in incorrect concentrations because an incorrect post-event volume was used. This is fixed here.

Closes #2426.

python/sdist/amici/sbml_import.py

@@ -1613,7 +1613,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()
@@ -1644,23 +1646,23 @@ 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
                 ]:
-                    # If the species was not affected by an event assignment
+                    # If the species was not affected by an event assignment,
                     # then the old value should be updated.
                     if (
                         bolus[state_vector.index(species_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,
+        )