From 0464b832dfbf13c9860d92e70ebace08e81e9739 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 7 May 2024 12:27:44 +0200 Subject: [PATCH 1/2] 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, + ) From 41d19bb5d9dcbafae714260d8982da69223bba0c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 7 May 2024 13:17:49 +0200 Subject: [PATCH 2/2] Support event-assignments to compartments that aren't rate rule targets (#2425) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Event-assignments to compartments weren't supported yet, now they are. Compartments that are targets of event assignments are now processed as "species", the same way as it is already done for parameters that are targets of event assignments. Closes #2424. --------- Co-authored-by: Fabian Fröhlich --- documentation/python_interface.rst | 2 +- python/sdist/amici/sbml_import.py | 37 +++++++++++++++++++++++++++--- 2 files changed, 35 insertions(+), 4 deletions(-) diff --git a/documentation/python_interface.rst b/documentation/python_interface.rst index a919e925c4..2926abf963 100644 --- a/documentation/python_interface.rst +++ b/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 `_ (`current status `_). diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 608ceeed8f..5d71821f06 100644 --- a/python/sdist/amici/sbml_import.py +++ b/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.""" @@ -1633,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