From 0464b832dfbf13c9860d92e70ebace08e81e9739 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 7 May 2024 12:27:44 +0200 Subject: [PATCH] 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 | 12 +++-- python/tests/test_events.py | 83 +++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 5 deletions(-) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index a26c6e7d4a..608ceeed8f 100644 --- a/python/sdist/amici/sbml_import.py +++ b/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,7 +1646,7 @@ 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, @@ -1652,15 +1654,15 @@ def get_empty_bolus_value() -> sp.Float: ) 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)] diff --git a/python/tests/test_events.py b/python/tests/test_events.py index 05eddab59f..d16877fd2e 100644 --- a/python/tests/test_events.py +++ b/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, + )