From 6f778dd0c29646d67c4c90e4dc257d32ee5f6ce1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Thu, 7 Apr 2022 08:30:18 -0400 Subject: [PATCH 01/11] Speedup for models with conservation laws (#1765) * implement speedup dtotal_cldx_rdata, dx_rdatadx_solver * fixup flattening * fixup * speedup dwdx * simplify flattening * cleanup sbml import * Apply suggestions from code review Co-authored-by: Daniel Weindl * fixup merge Co-authored-by: Daniel Weindl --- python/amici/ode_export.py | 92 +++++++++++-------- python/amici/ode_model.py | 119 ++++++++++++++++++++----- python/amici/pysb_import.py | 170 +++++++++++------------------------- python/amici/sbml_import.py | 41 ++------- python/tests/test_pysb.py | 3 +- 5 files changed, 211 insertions(+), 214 deletions(-) diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index e3da48b416..ff25b184c2 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -886,8 +886,7 @@ def add_component(self, component: ModelQuantity, def add_conservation_law(self, state: sp.Symbol, total_abundance: sp.Symbol, - state_expr: sp.Expr, - abundance_expr: sp.Expr) -> None: + coefficients: Dict[sp.Symbol, sp.Expr]) -> None: r""" Adds a new conservation law to the model. A conservation law is defined by the conserved quantity :math:`T = \sum_i(a_i * x_i)`, where @@ -901,16 +900,8 @@ def add_conservation_law(self, :param total_abundance: symbolic identifier of the total abundance (:math:`T/a_j`) - :param state_expr: - symbolic algebraic formula that replaces the state. This is - used to compute the numeric value of ``state`` during simulations. - :math:`x_j = T/a_j - \sum_{i≠j}(a_i * x_i)/a_j` - - :param abundance_expr: - symbolic algebraic formula that computes the value of the - conserved quantity. This is used to update the numeric value for - ``total_abundance`` after (re-)initialization. - :math:`T/a_j = \sum_{i≠j}(a_i * x_i)/a_j + x_j` + :param coefficients: + Dictionary of coefficients {x_i: a_i} """ try: ix = [ @@ -923,20 +914,29 @@ def add_conservation_law(self, state_id = self._states[ix].get_id() + # \sum_{i≠j}(a_i * x_i)/a_j + target_expression = sp.Add(*( + c_i*x_i for x_i, c_i in coefficients.items() if x_i != state + )) / coefficients[state] + + # x_j = T/a_j - \sum_{i≠j}(a_i * x_i)/a_j + state_expr = total_abundance - target_expression + + # T/a_j = \sum_{i≠j}(a_i * x_i)/a_j + x_j + abundance_expr = target_expression + state_id + self.add_component( Expression(state_id, str(state_id), state_expr), insert_first=True ) - self.add_component( - ConservationLaw( - total_abundance, - f'total_{state_id}', - abundance_expr - ) + cl = ConservationLaw( + total_abundance, f'total_{state_id}', abundance_expr, + coefficients, state_id ) - self._states[ix].set_conservation_law(state_expr) + self.add_component(cl) + self._states[ix].set_conservation_law(cl) def get_observable_transformations(self) -> List[ObservableTransformation]: """ @@ -1219,14 +1219,14 @@ def _generate_symbol(self, name: str, *, from_sbml: bool = False) -> None: self._syms[name] = sp.Matrix([ state.get_id() for state in self._states - if state._conservation_law is None + if not state.has_conservation_law() ]) return elif name == 'sx0': self._syms[name] = sp.Matrix([ f's{state.get_id()}_0' for state in self._states - if state._conservation_law is None + if not state.has_conservation_law() ]) return elif name == 'sx_rdata': @@ -1422,15 +1422,13 @@ def _compute_equation(self, name: str) -> None: elif name == 'xdot': self._eqs[name] = sp.Matrix([ - s.get_dt() for s in self._states - if s._conservation_law is None + state.get_dt() for state in self._states + if not state.has_conservation_law() ]) elif name == 'x_rdata': self._eqs[name] = sp.Matrix([ - state.get_id() - if state._conservation_law is None - else state._conservation_law + state.get_x_rdata() for state in self._states ]) @@ -1438,14 +1436,14 @@ def _compute_equation(self, name: str) -> None: self._eqs[name] = sp.Matrix([ state.get_id() for state in self._states - if state._conservation_law is None + if not state.has_conservation_law() ]) elif name == 'sx_solver': self._eqs[name] = sp.Matrix([ self.sym('sx_rdata')[ix] for ix, state in enumerate(self._states) - if state._conservation_law is None + if not state.has_conservation_law() ]) elif name == 'sx0': @@ -1484,8 +1482,13 @@ def _compute_equation(self, name: str) -> None: self._x0_fixedParameters_idx]) elif name == 'dtotal_cldx_rdata': - # not correctly parsed in regex - self._derivative('total_cl', 'x_rdata') + x_rdata = self.sym('x_rdata') + self._eqs[name] = sp.Matrix( + [ + [cl.get_ncoeff(xr) for xr in x_rdata] + for cl in self._conservationlaws + ] + ) elif name == 'dtcldx': # this is always zero @@ -1498,8 +1501,13 @@ def _compute_equation(self, name: str) -> None: elif name == 'dx_rdatadx_solver': if self.num_cons_law(): - self._eqs[name] = smart_jacobian(self.eq('x_rdata'), - self.sym('x')) + x_solver = self.sym('x') + self._eqs[name] = sp.Matrix( + [ + [state.get_dx_rdata_dx_solver(xs) for xs in x_solver] + for state in self._states + ] + ) else: # so far, dx_rdatadx_solver is only required for sx_rdata # in case of no conservation laws, C++ code will directly use @@ -1617,6 +1625,16 @@ def _compute_equation(self, name: str) -> None: # force symbols self._eqs[name] = self.sym(name) + elif name == 'dwdx': + x = self.sym('x') + self._eqs[name] = sp.Matrix([ + [-cl.get_ncoeff(xs) for xs in x] + # the insert first in ode_model._add_conservation_law() means + # that we need to reverse the order here + for cl in reversed(self._conservationlaws) + ]) .col_join(smart_jacobian(self.eq('w')[self.num_cons_law():,:], + x)) + elif match_deriv: self._derivative(match_deriv.group(1), match_deriv.group(2), name) @@ -1881,16 +1899,16 @@ def _equation_from_component(self, name: str, component: str) -> None: [comp.get_val() for comp in getattr(self, component)] ) - def get_conservation_laws(self) -> List[Tuple[sp.Symbol, sp.Basic]]: + def get_conservation_laws(self) -> List[Tuple[sp.Symbol, sp.Expr]]: """Returns a list of states with conservation law set :return: list of state identifiers """ return [ - (state.get_id(), state._conservation_law) + (state.get_id(), state.get_x_rdata()) for state in self._states - if state._conservation_law is not None + if state.has_conservation_law() ] def _generate_value(self, name: str) -> None: @@ -1958,7 +1976,7 @@ def state_has_conservation_law(self, ix: int) -> bool: :return: boolean indicating if conservation_law is not None """ - return self._states[ix]._conservation_law is not None + return self._states[ix].has_conservation_law() def get_solver_indices(self) -> Dict[int, int]: """ @@ -2904,7 +2922,7 @@ def _write_model_header_cpp(self) -> None: [ str(idx) for idx, state in enumerate(self.model._states) - if state._conservation_law is None + if not state.has_conservation_law() ] ), 'REINIT_FIXPAR_INITCOND': diff --git a/python/amici/ode_model.py b/python/amici/ode_model.py index 60d44147d2..a5885216f3 100644 --- a/python/amici/ode_model.py +++ b/python/amici/ode_model.py @@ -131,6 +131,73 @@ def set_val(self, val: sp.Expr): self._value = cast_to_sym(val, 'value') +class ConservationLaw(ModelQuantity): + """ + A conservation law defines the absolute the total amount of a + (weighted) sum of states + + """ + def __init__(self, + identifier: sp.Symbol, + name: str, + value: sp.Expr, + coefficients: Dict[sp.Symbol, sp.Expr], + state_id: sp.Symbol): + """ + Create a new ConservationLaw instance. + + :param identifier: + unique identifier of the ConservationLaw + + :param name: + individual name of the ConservationLaw (does not need to be + unique) + + :param value: formula (sum of states) + + :param coefficients: + coefficients of the states in the sum + + :param state_id: + identifier of the state that this conservation law replaces + """ + self._state_expr: sp.Symbol = identifier - (value - state_id) + self._coefficients: Dict[sp.Symbol, sp.Expr] = coefficients + self._ncoeff: sp.Expr = coefficients[state_id] + super(ConservationLaw, self).__init__(identifier, name, value) + + def get_state(self) -> sp.Symbol: + """ + Get the identifier of the state that this conservation law replaces + + :return: identifier of the state + """ + return self._state_id + + def get_ncoeff(self, state_id) -> Union[sp.Expr, int, float]: + """ + Computes the normalized coefficient a_i/a_j where i is the index of + the provided state_id and j is the index of the state that is + replaced by this conservation law. This can be used to compute both + dtotal_cl/dx_rdata (=ncoeff) and dx_rdata/dx_solver (=-ncoeff). + + :param state_id: + identifier of the state + + :return: normalized coefficent of the state + """ + return self._coefficients.get(state_id, 0.0) / self._ncoeff + + def get_x_rdata(self): + """ + Returns the expression that allows computation of x_rdata for the state + that this conservation law replaces. + + :return: x_rdata expression + """ + return self._state_expr + + class State(ModelQuantity): """ A State variable defines an entity that evolves with time according to @@ -171,10 +238,9 @@ def __init__(self, """ super(State, self).__init__(identifier, name, init) self._dt = cast_to_sym(dt, 'dt') - self._conservation_law = None + self._conservation_law: Union[ConservationLaw, None] = None - def set_conservation_law(self, - law: sp.Expr) -> None: + def set_conservation_law(self, law: ConservationLaw) -> None: """ Sets the conservation law of a state. @@ -185,9 +251,9 @@ def set_conservation_law(self, linear sum of states that if added to this state remain constant over time """ - if not isinstance(law, sp.Expr): - raise TypeError(f'conservation law must have type sympy.Expr, ' - f'was {type(law)}') + if not isinstance(law, ConservationLaw): + raise TypeError(f'conservation law must have type ConservationLaw' + f', was {type(law)}') self._conservation_law = law @@ -219,30 +285,37 @@ def get_free_symbols(self) -> Set[sp.Symbol]: """ return self._dt.free_symbols.union(self._value.free_symbols) + def has_conservation_law(self): + """ + Checks whether this state has a conservation law assigned. -class ConservationLaw(ModelQuantity): - """ - A conservation law defines the absolute the total amount of a - (weighted) sum of states + :return: True if assigned, False otherwise + """ + return self._conservation_law is not None - """ - def __init__(self, - identifier: sp.Symbol, - name: str, - value: sp.Expr): + def get_x_rdata(self): """ - Create a new ConservationLaw instance. + Returns the expression that allows computation of x_rdata for this + state, accounting for conservation laws. - :param identifier: - unique identifier of the ConservationLaw + :return: x_rdata expression + """ + if self._conservation_law is None: + return self.get_id() + else: + return self._conservation_law.get_x_rdata() - :param name: - individual name of the ConservationLaw (does not need to be - unique) + def get_dx_rdata_dx_solver(self, state_id): + """ + Returns the expression that allows computation of ``dx_rdata_dx_solver`` for this + state, accounting for conservation laws. - :param value: formula (sum of states) + :return: dx_rdata_dx_solver expression """ - super(ConservationLaw, self).__init__(identifier, name, value) + if self._conservation_law is None: + return sp.Integer(self._identifier == state_id) + else: + return -self._conservation_law.get_ncoeff(state_id) class Observable(ModelQuantity): diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py index f666d1bb4a..90e5302f3e 100644 --- a/python/amici/pysb_import.py +++ b/python/amici/pysb_import.py @@ -1109,37 +1109,17 @@ def _construct_conservation_from_prototypes( conservation_laws = [] for monomer_name in cl_prototypes: target_index = cl_prototypes[monomer_name]['target_index'] + coefficients = dict() - # T = sum_i(a_i * x_i) - # x_j = (T - sum_i≠j(a_i * x_i))/a_j - # law: sum_i≠j(a_i * x_i))/a_j - # state: x_j - target_expression = sp.Add(*( - sp.Symbol(f'__s{ix}') - * extract_monomers(specie).count(monomer_name) - for ix, specie in enumerate(pysb_model.species) - if ix != target_index - )) / extract_monomers(pysb_model.species[ - target_index - ]).count(monomer_name) - # normalize by the stoichiometry of the target species - target_state = sp.Symbol(f'__s{target_index}') - # = x_j - - total_abundance = sp.Symbol(f'tcl__s{target_index}') - # = T/a_j - - state_expr = total_abundance - target_expression - # x_j = T/a_j - sum_i≠j(a_i * x_i)/a_j - - abundance_expr = target_expression + target_state - # T/a_j = sum_i≠j(a_i * x_i)/a_j + x_j + for ix, specie in enumerate(pysb_model.species): + count = extract_monomers(specie).count(monomer_name) + if count > 0: + coefficients[sp.Symbol(f'__s{ix}')] = count conservation_laws.append({ - 'state': target_state, - 'total_abundance': total_abundance, - 'state_expr': state_expr, - 'abundance_expr': abundance_expr, + 'state': sp.Symbol(f'__s{target_index}'), + 'total_abundance': sp.Symbol(f'tcl__s{target_index}'), + 'coefficients': coefficients, }) return conservation_laws @@ -1163,14 +1143,10 @@ def _add_conservation_for_constant_species( for ix in range(ode_model.num_states_rdata()): if ode_model.state_is_constant(ix): - target_state = sp.Symbol(f'__s{ix}') - total_abundance = sp.Symbol(f'tcl__s{ix}') - conservation_laws.append({ - 'state': target_state, - 'total_abundance': total_abundance, - 'state_expr': total_abundance, - 'abundance_expr': target_state, + 'state': sp.Symbol(f'__s{ix}'), + 'total_abundance': sp.Symbol(f'tcl__s{ix}'), + 'coefficients': {sp.Symbol(f'__s{ix}'): 1.0} }) @@ -1186,83 +1162,58 @@ def _flatten_conservation_laws( conservation_law_subs = \ _get_conservation_law_subs(conservation_laws) - while len(conservation_law_subs): + while conservation_law_subs: for cl in conservation_laws: - if _sub_matches_cl( - conservation_law_subs, - cl['state_expr'], - cl['state'] + # only update if we changed something + if any( + _apply_conseration_law_sub(cl, sub) + for sub in conservation_law_subs ): - # this optimization is done by subs anyways, but we dont - # want to recompute the subs if we did not change anything - valid_subs = _select_valid_cls( - conservation_law_subs, - cl['state'] - ) - if len(valid_subs) > 0: - cl['state_expr'] = cl['state_expr'].subs(valid_subs) - conservation_law_subs = \ - _get_conservation_law_subs(conservation_laws) - - -def _select_valid_cls(subs: Iterable[Tuple[sp.Symbol, sp.Basic]], - state: sp.Symbol) -> List[Tuple[sp.Symbol, sp.Basic]]: - """ - Subselect substitutions such that we do not end up with conservation - laws that are self-referential - - :param subs: - substitutions in tuple format - - :param state: - target symbolic state to which substitutions will be applied - - :return: - list of valid substitutions - """ - return [ - sub - for sub in subs - if str(state) not in [str(symbol) for symbol in sub[1].free_symbols] - ] + conservation_law_subs = \ + _get_conservation_law_subs(conservation_laws) -def _sub_matches_cl(subs: Iterable[Tuple[sp.Symbol, sp.Basic]], - state_expr: sp.Basic, - state: sp.Basic) -> bool: +def _apply_conseration_law_sub(cl: ConservationLaw, + sub: Tuple[sp.Symbol, ConservationLaw]) -> bool: """ - Checks whether any of the substitutions in subs will be applied to - state_expr - - :param subs: - substitutions in tuple format + Applies a substitution to a conservation law by replacing the + coefficient of the state of the - :param state_expr: - target symbolic expressions in which substitutions will be applied + :param cl: + conservation law - :param state: target symbolic state to which substitutions will - be applied + :param sub: + substitution to apply, tuple of (state to be replaced, conservation + law) - :return: - boolean indicating positive match + :return: boolean flag indicating whether the substitution was applied """ + coeff = cl['coefficients'].get(sub[0], 0.0) + if coeff == 0.0 or cl['state'] == sub[0]: + return False - sub_symbols = set( - sub[0] - for sub in subs - if str(state) not in [ - str(symbol) for symbol in sub[1].free_symbols - ] - ) + del cl['coefficients'][sub[0]] + # x_j = T/b_j - sum_{i≠j}(x_i * b_i) / b_j + # don't need to account for totals here as we can simply + # absorb that into the new total + for k, v in sub[1].items(): + if k == sub[0]: + continue + update = - coeff * v / sub[1][sub[0]] - return len(sub_symbols.intersection(state_expr.free_symbols)) > 0 + if k in cl['coefficients']: + cl['coefficients'][k] += update + else: + cl['coefficients'][k] = update + + return True def _get_conservation_law_subs( conservation_laws: List[ConservationLaw] -) -> List[Tuple[sp.Symbol, sp.Basic]]: +) -> List[Tuple[sp.Symbol, Dict[sp.Symbol, sp.Expr]]]: """ - Computes a list of (state, law) tuples for conservation laws that still + Computes a list of (state, coeffs) tuples for conservation laws that still appear in other conservation laws :param conservation_laws: @@ -1272,30 +1223,15 @@ def _get_conservation_law_subs( list of tuples containing substitution rules to be used with sympy subs """ - free_symbols_cl = _conservation_law_variables(conservation_laws) return [ - (cl['state'], cl['state_expr']) for cl in conservation_laws - if cl['state'] in free_symbols_cl + (cl['state'], cl['coefficients']) for cl in conservation_laws + if any( + cl['state'] in other_cl['coefficients'] + for other_cl in conservation_laws + if other_cl != cl + ) ] - -def _conservation_law_variables( - conservation_laws: List[ConservationLaw]) -> Set[sp.Symbol]: - """ - Construct the set of all free variables from a list of conservation laws - - :param conservation_laws: - list of conservation laws - - :return: - free variables in conservation laws - """ - variables = set() - for cl in conservation_laws: - variables |= cl['state_expr'].free_symbols - return variables - - def has_fixed_parameter_ic(specie: pysb.core.ComplexPattern, pysb_model: pysb.Model, ode_model: ODEModel) -> bool: diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index 6387708bd4..df00522698 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -1631,47 +1631,19 @@ def _add_conservation_for_non_constant_species( compartment_sizes = [all_compartment_sizes[i] for i in state_idxs] target_state_id = all_state_ids[target_state_model_idx] - target_compartment = all_compartment_sizes[target_state_model_idx] - target_state_coeff = coefficients[0] total_abundance = symbol_with_assumptions(f'tcl_{target_state_id}') - # \sum coeff * state * volume - abundance_expr = sp.Add(*[ - state_id * coeff * compartment - for state_id, coeff, compartment - in zip(state_ids, coefficients, compartment_sizes) - ]) - new_conservation_laws.append({ 'state': target_state_id, 'total_abundance': total_abundance, - 'state_expr': - (total_abundance - (abundance_expr - - target_state_id * target_compartment - * target_state_coeff)) - / target_state_coeff / target_compartment, - 'abundance_expr': abundance_expr + 'coefficients': { + state_id: coeff * compartment + for state_id, coeff, compartment + in zip(state_ids, coefficients, compartment_sizes) + }, }) species_to_be_removed.add(target_state_model_idx) - # replace eliminated states by their state expressions, taking care of - # any (non-cyclic) dependencies - state_exprs = { - cl['state']: cl['state_expr'] - for cl in itt.chain(conservation_laws, new_conservation_laws) - } - try: - sorted_state_exprs = toposort_symbols(state_exprs) - except CircularDependencyError as e: - raise AssertionError( - "Circular dependency detected in conservation laws. " - "This should not have happened." - ) from e - - for cl in new_conservation_laws: - cl['state_expr'] = smart_subs_dict(cl['state_expr'], - sorted_state_exprs) - conservation_laws.extend(new_conservation_laws) # list of species that are not determined by conservation laws @@ -2090,8 +2062,7 @@ def _add_conservation_for_constant_species( conservation_laws.append({ 'state': target_state, 'total_abundance': total_abundance, - 'state_expr': total_abundance, - 'abundance_expr': target_state, + 'coefficients': {target_state: 1.0}, }) # mark species to delete from stoichiometric matrix species_solver.pop(ix) diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py index df8558fadd..9d373c6d65 100644 --- a/python/tests/test_pysb.py +++ b/python/tests/test_pysb.py @@ -165,8 +165,7 @@ def test_compare_to_pysb_simulation(example): solver.setRelativeTolerance(rtol) rdata = amici.runAmiciSimulation(model_pysb, solver) - # check agreement of species simulation - + # check agreement of species simulations assert np.isclose(rdata['x'], pysb_simres.species, 1e-4, 1e-4).all() From f6ee397dc9051f39cb0ea3caec4d715142222789 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 7 Apr 2022 17:14:28 +0200 Subject: [PATCH 02/11] Update petab test suite code (#1768) Adapt to separate sbml and pysb test suite (https://github.com/PEtab-dev/petab_test_suite/pull/47) --- python/amici/petab_import_pysb.py | 17 +++++++------- python/tests/test_petab_simulate.py | 2 +- tests/petab_test_suite/conftest.py | 20 +++++++++++------ tests/petab_test_suite/test_petab_suite.py | 26 +++++++++++----------- 4 files changed, 36 insertions(+), 29 deletions(-) diff --git a/python/amici/petab_import_pysb.py b/python/amici/petab_import_pysb.py index 3d33906e5b..9f71ad11b8 100644 --- a/python/amici/petab_import_pysb.py +++ b/python/amici/petab_import_pysb.py @@ -8,16 +8,17 @@ import logging import os from itertools import chain -from typing import List, Dict, Union, Optional, Tuple, Iterable +from pathlib import Path +from typing import Dict, Iterable, List, Optional, Tuple, Union import libsbml import petab import pysb import sympy as sp -from petab.C import (CONDITION_NAME, OBSERVABLE_TRANSFORMATION, LIN, - OBSERVABLE_FORMULA, NOISE_FORMULA, FORMAT_VERSION, - PARAMETER_FILE, SBML_FILES, CONDITION_FILES, - MEASUREMENT_FILES, VISUALIZATION_FILES, OBSERVABLE_FILES) +from petab.C import (CONDITION_FILES, CONDITION_NAME, FORMAT_VERSION, + MEASUREMENT_FILES, NOISE_FORMULA, OBSERVABLE_FILES, + OBSERVABLE_FORMULA, PARAMETER_FILE, SBML_FILES, + VISUALIZATION_FILES) from . import petab_import from .logging import get_logger, log_execution_time, set_log_level @@ -194,7 +195,7 @@ def from_yaml(yaml_config: Union[Dict, str], """ from petab.yaml import (load_yaml, is_composite_problem, assert_single_condition_and_sbml_file) - if isinstance(yaml_config, str): + if isinstance(yaml_config, (str, Path)): path_prefix = os.path.dirname(yaml_config) yaml_config = load_yaml(yaml_config) else: @@ -328,7 +329,7 @@ def import_model_pysb( set_log_level(logger, verbose) - logger.info(f"Importing model ...") + logger.info("Importing model ...") observable_table = petab_problem.observable_df pysb_model = petab_problem.pysb_model @@ -355,6 +356,7 @@ def import_model_pysb( if observable_table is None: observables = None sigmas = None + noise_distrs = None else: observables = [expr.name for expr in pysb_model.expressions if expr.name in observable_table.index] @@ -364,7 +366,6 @@ def import_model_pysb( noise_distrs = petab_import.petab_noise_distributions_to_amici( observable_table) - from amici.pysb_import import pysb2amici pysb2amici(pysb_model, model_output_dir, verbose=True, observables=observables, diff --git a/python/tests/test_petab_simulate.py b/python/tests/test_petab_simulate.py index ecb0a2b168..a992673702 100644 --- a/python/tests/test_petab_simulate.py +++ b/python/tests/test_petab_simulate.py @@ -13,7 +13,7 @@ def petab_problem() -> petab.Problem: """Create a PEtab problem for use in tests.""" test_case = '0001' - test_case_dir = Path(petabtests.SBML_DIR) / petabtests.CASES_LIST[0] + test_case_dir = Path(petabtests.SBML_DIR) / test_case petab_yaml_path = test_case_dir / petabtests.problem_yaml_name(test_case) return petab.Problem.from_yaml(str(petab_yaml_path)) diff --git a/tests/petab_test_suite/conftest.py b/tests/petab_test_suite/conftest.py index 2a4c9cdd69..abcd9d821d 100644 --- a/tests/petab_test_suite/conftest.py +++ b/tests/petab_test_suite/conftest.py @@ -3,7 +3,7 @@ from typing import List import re import sys -import petabtests +from petabtests.core import get_cases def parse_selection(selection_str: str) -> List[int]: @@ -43,6 +43,7 @@ def pytest_generate_tests(metafunc): # Run for all PEtab test suite cases if "case" in metafunc.fixturenames \ and "model_type" in metafunc.fixturenames: + # Get CLI option cases = metafunc.config.getoption("--petab-cases") if cases: @@ -50,15 +51,20 @@ def pytest_generate_tests(metafunc): test_numbers = parse_selection(cases) else: # Run all tests - test_numbers = petabtests.CASES_LIST + test_numbers = None if metafunc.config.getoption("--only-sbml"): - model_types = ['sbml'] + test_numbers = test_numbers if test_numbers else get_cases("sbml") + argvalues = [(case, 'sbml') for case in test_numbers] + elif metafunc.config.getoption("--only-pysb"): - model_types = ['pysb'] + test_numbers = test_numbers if test_numbers else get_cases("pysb") + argvalues = [(case, 'pysb') for case in test_numbers] else: - model_types = ['sbml', 'pysb'] + argvalues = [] + for format in ('sbml', 'pysb'): + test_numbers = test_numbers if test_numbers else get_cases( + format) + argvalues += [(case, format) for case in test_numbers] - argvalues = [(case, model_type) for model_type in model_types - for case in test_numbers] metafunc.parametrize("case,model_type", argvalues) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 97d4bc6866..33a6c47e18 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -3,7 +3,6 @@ """Run PEtab test suite (https://github.com/PEtab-dev/petab_test_suite)""" import logging -import os import sys import amici @@ -11,12 +10,12 @@ import petabtests import pytest from _pytest.outcomes import Skipped +from amici import SteadyStateSensitivityMode from amici.gradient_check import check_derivatives as amici_check_derivatives from amici.logging import get_logger, set_log_level -from amici.petab_import import import_petab_problem, PysbPetabProblem -from amici.petab_objective import ( - simulate_petab, rdatas_to_measurement_df, create_parameterized_edatas) -from amici import SteadyStateSensitivityMode +from amici.petab_import import PysbPetabProblem, import_petab_problem +from amici.petab_objective import (create_parameterized_edatas, + rdatas_to_measurement_df, simulate_petab) logger = get_logger(__name__, logging.DEBUG) set_log_level(get_logger("amici.petab_import"), logging.DEBUG) @@ -46,17 +45,17 @@ def _test_case(case, model_type): # load if model_type == "sbml": - case_dir = os.path.join(petabtests.SBML_DIR, case) + case_dir = petabtests.SBML_DIR / case # import petab problem - yaml_file = os.path.join(case_dir, petabtests.problem_yaml_name(case)) + yaml_file = case_dir / petabtests.problem_yaml_name(case) problem = petab.Problem.from_yaml(yaml_file) elif model_type == "pysb": import pysb pysb.SelfExporter.cleanup() pysb.SelfExporter.do_export = True - case_dir = os.path.join(petabtests.PYSB_DIR, case) + case_dir = petabtests.PYSB_DIR / case # import petab problem - yaml_file = os.path.join(case_dir, petabtests.problem_yaml_name(case)) + yaml_file = case_dir / petabtests.problem_yaml_name(case) problem = PysbPetabProblem.from_yaml(yaml_file, flatten=case.startswith('0006')) else: @@ -152,9 +151,10 @@ def run(): n_success = 0 n_skipped = 0 - for case in petabtests.CASES_LIST: + cases = petabtests.get_cases('sbml') + for case in cases: try: - test_case(case) + test_case(case, 'sbml') n_success += 1 except Skipped: n_skipped += 1 @@ -163,9 +163,9 @@ def run(): logger.error(f"Case {case} failed.") logger.error(e) - logger.info(f"{n_success} / {len(petabtests.CASES_LIST)} successful, " + logger.info(f"{n_success} / {len(cases)} successful, " f"{n_skipped} skipped") - if n_success != len(petabtests.CASES_LIST): + if n_success != len(cases): sys.exit(1) From 6879ccd803589aeb29c22b61eb86d3ace11c4f78 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Thu, 7 Apr 2022 23:08:46 +0200 Subject: [PATCH 03/11] Steady state tolerance factor (#1758) Add option to specify steady state (sensitivity) tolerances relative to integration tolerances based on a multiplier. **NOTE: Relaxes default steady state (sensitivity) tolerances by a factor of 100.** Co-authored-by: Daniel Weindl --- include/amici/serialization.h | 2 + include/amici/solver.h | 64 +++++++++++++++--- src/hdf5.cpp | 19 +++++- src/solver.cpp | 36 ++++++++-- src/steadystateproblem.cpp | 2 +- tests/cpp/testOptions.h5 | Bin 339328 -> 336344 bytes tests/cpp/unittests/testMisc.cpp | 53 +++++++++++++++ .../generateTestConfig/example_steadystate.py | 29 ++++++++ tests/petab_test_suite/conftest.py | 7 +- tests/petab_test_suite/test_petab_suite.py | 9 ++- 10 files changed, 198 insertions(+), 23 deletions(-) diff --git a/include/amici/serialization.h b/include/amici/serialization.h index 6cd6e7865b..3ba9ee7ffd 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -62,8 +62,10 @@ void serialize(Archive &ar, amici::Solver &s, const unsigned int /*version*/) { ar &s.rtol_fsa_; ar &s.quad_atol_; ar &s.quad_rtol_; + ar &s.ss_tol_factor_; ar &s.ss_atol_; ar &s.ss_rtol_; + ar &s.ss_tol_sensi_factor_; ar &s.ss_atol_sensi_; ar &s.ss_rtol_sensi_; ar &s.maxsteps_; diff --git a/include/amici/solver.h b/include/amici/solver.h index edb761a82f..42fa67f5c2 100644 --- a/include/amici/solver.h +++ b/include/amici/solver.h @@ -44,7 +44,7 @@ namespace amici { * * NOTE: Any changes in data members here must be propagated to copy ctor, * equality operator, serialization functions in serialization.h, and - * amici::hdf5::readSolverSettingsFromHDF5 in hdf5.cpp. + * amici::hdf5::(read/write)SolverSettings(From/To)HDF5 in hdf5.cpp. */ class Solver { public: @@ -384,6 +384,26 @@ class Solver { */ void setAbsoluteToleranceQuadratures(double atol); + /** + * @brief returns the steady state simulation tolerance factor. + * + * Steady state simulation tolerances are the product of the simulation + * tolerances and this factor, unless manually set with + * `set(Absolute/Relative)ToleranceSteadyState()`. + * @return steady state simulation tolerance factor + */ + double getSteadyStateToleranceFactor() const; + + /** + * @brief set the steady state simulation tolerance factor. + * + * Steady state simulation tolerances are the product of the simulation + * tolerances and this factor, unless manually set with + * `set(Absolute/Relative)ToleranceSteadyState()`. + * @param factor tolerance factor (non-negative number) + */ + void setSteadyStateToleranceFactor(double factor); + /** * @brief returns the relative tolerance for the steady state problem * @return relative tolerance @@ -408,6 +428,26 @@ class Solver { */ void setAbsoluteToleranceSteadyState(double atol); + /** + * @brief returns the steady state sensitivity simulation tolerance factor. + * + * Steady state sensitivity simulation tolerances are the product of the + * sensitivity simulation tolerances and this factor, unless manually set + * with `set(Absolute/Relative)ToleranceSteadyStateSensi()`. + * @return steady state simulation tolerance factor + */ + double getSteadyStateSensiToleranceFactor() const; + + /** + * @brief set the steady state sensitivity simulation tolerance factor. + * + * Steady state sensitivity simulation tolerances are the product of the + * sensitivity simulation tolerances and this factor, unless manually set + * with `set(Absolute/Relative)ToleranceSteadyStateSensi()`. + * @param factor tolerance factor (non-negative number) + */ + void setSteadyStateSensiToleranceFactor(double factor); + /** * @brief returns the relative tolerance for the sensitivities of the * steady state problem @@ -855,7 +895,7 @@ class Solver { std::vector const& getLastOrder() const { return order_; } - + /** * @brief Returns how convergence checks for steadystate computation are performed. * @return boolean flag indicating newton step (true) or the right hand side (false) @@ -863,7 +903,7 @@ class Solver { bool getNewtonStepSteadyStateCheck() const { return newton_step_steadystate_conv_; } - + /** * @brief Returns how convergence checks for steadystate computation are performed. * @return boolean flag indicating state and sensitivity equations (true) or only state variables (false). @@ -871,7 +911,7 @@ class Solver { bool getSensiSteadyStateCheck() const { return check_sensi_steadystate_conv_; } - + /** * @brief Sets how convergence checks for steadystate computation are performed. * @param flag boolean flag to pick newton step (true) or the right hand side (false, default) @@ -879,7 +919,7 @@ class Solver { void setNewtonStepSteadyStateCheck(bool flag) { newton_step_steadystate_conv_ = flag; } - + /** * @brief Sets for which variables convergence checks for steadystate computation are performed. * @param flag boolean flag to pick state and sensitivity equations (true, default) or only state variables (false). @@ -1753,23 +1793,29 @@ class Solver { /** relative tolerances for backward quadratures */ realtype quad_rtol_ {1e-8}; + /** steady state simulation tolerance factor */ + realtype ss_tol_factor_ {1e2}; + /** absolute tolerances for steadystate computation */ realtype ss_atol_ {NAN}; /** relative tolerances for steadystate computation */ realtype ss_rtol_ {NAN}; - /** absolute tolerances for steadystate computation */ + /** steady state sensitivity simulation tolerance factor */ + realtype ss_tol_sensi_factor_ {1e2}; + + /** absolute tolerances for steadystate sensitivity computation */ realtype ss_atol_sensi_ {NAN}; - /** relative tolerances for steadystate computation */ + /** relative tolerances for steadystate sensitivity computation */ realtype ss_rtol_sensi_ {NAN}; RDataReporting rdata_mode_ {RDataReporting::full}; - + /** whether newton step should be used for convergence steps */ bool newton_step_steadystate_conv_ {false}; - + /** whether sensitivities should be checked for convergence to steadystate */ bool check_sensi_steadystate_conv_ {true}; diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 9c006dfa1f..168841e317 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -657,6 +657,10 @@ void writeSolverSettingsToHDF5(Solver const& solver, H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), "quad_rtol", &dbuffer, 1); + dbuffer = solver.getSteadyStateToleranceFactor(); + H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), + "ss_tol_factor", &dbuffer, 1); + dbuffer = solver.getAbsoluteToleranceSteadyState(); H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), "ss_atol", &dbuffer, 1); @@ -665,6 +669,10 @@ void writeSolverSettingsToHDF5(Solver const& solver, H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), "ss_rtol", &dbuffer, 1); + dbuffer = solver.getSteadyStateSensiToleranceFactor(); + H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), + "ss_tol_sensi_factor", &dbuffer, 1); + dbuffer = solver.getAbsoluteToleranceSteadyStateSensi(); H5LTset_attribute_double(file.getId(), hdf5Location.c_str(), "ss_atol_sensi", &dbuffer, 1); @@ -773,7 +781,6 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, getDoubleScalarAttribute(file, datasetPath, "rtol_fsa")); } - if(attributeExists(file, datasetPath, "atolB")) { solver.setAbsoluteToleranceB( getDoubleScalarAttribute(file, datasetPath, "atolB")); @@ -794,6 +801,11 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, getDoubleScalarAttribute(file, datasetPath, "quad_rtol")); } + if(attributeExists(file, datasetPath, "ss_tol_factor")) { + solver.setSteadyStateToleranceFactor( + getDoubleScalarAttribute(file, datasetPath, "ss_tol_factor")); + } + if(attributeExists(file, datasetPath, "ss_atol")) { solver.setAbsoluteToleranceSteadyState( getDoubleScalarAttribute(file, datasetPath, "ss_atol")); @@ -804,6 +816,11 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, getDoubleScalarAttribute(file, datasetPath, "ss_rtol")); } + if(attributeExists(file, datasetPath, "ss_tol_sensi_factor")) { + solver.setSteadyStateSensiToleranceFactor( + getDoubleScalarAttribute(file, datasetPath, "ss_tol_sensi_factor")); + } + if(attributeExists(file, datasetPath, "ss_atol_sensi")) { solver.setAbsoluteToleranceSteadyStateSensi( getDoubleScalarAttribute(file, datasetPath, diff --git a/src/solver.cpp b/src/solver.cpp index 05e8d8d8b2..091c8bb7f2 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -29,8 +29,10 @@ Solver::Solver(const Solver &other) linsol_(other.linsol_), atol_(other.atol_), rtol_(other.rtol_), atol_fsa_(other.atol_fsa_), rtol_fsa_(other.rtol_fsa_), atolB_(other.atolB_), rtolB_(other.rtolB_), quad_atol_(other.quad_atol_), - quad_rtol_(other.quad_rtol_), ss_atol_(other.ss_atol_), - ss_rtol_(other.ss_rtol_), ss_atol_sensi_(other.ss_atol_sensi_), + quad_rtol_(other.quad_rtol_), ss_tol_factor_(other.ss_tol_factor_), + ss_atol_(other.ss_atol_), ss_rtol_(other.ss_rtol_), + ss_tol_sensi_factor_(other.ss_tol_sensi_factor_), + ss_atol_sensi_(other.ss_atol_sensi_), ss_rtol_sensi_(other.ss_rtol_sensi_), rdata_mode_(other.rdata_mode_), newton_step_steadystate_conv_(other.newton_step_steadystate_conv_), check_sensi_steadystate_conv_(other.check_sensi_steadystate_conv_), @@ -781,8 +783,19 @@ void Solver::setAbsoluteToleranceQuadratures(const double atol) { applyTolerancesASA(iMem); } +double Solver::getSteadyStateToleranceFactor() const { + return static_cast(ss_tol_factor_); +} + +void Solver::setSteadyStateToleranceFactor(const double ss_tol_factor) { + if (ss_tol_factor < 0) + throw AmiException("ss_tol_factor must be a non-negative number"); + + ss_tol_factor_ = static_cast(ss_tol_factor); +} + double Solver::getRelativeToleranceSteadyState() const { - return static_cast(isNaN(ss_rtol_) ? rtol_ : ss_rtol_); + return static_cast(isNaN(ss_rtol_) ? rtol_ * ss_tol_factor_ : ss_rtol_); } void Solver::setRelativeToleranceSteadyState(const double rtol) { @@ -793,7 +806,7 @@ void Solver::setRelativeToleranceSteadyState(const double rtol) { } double Solver::getAbsoluteToleranceSteadyState() const { - return static_cast(isNaN(ss_atol_) ? atol_ : ss_atol_); + return static_cast(isNaN(ss_atol_) ? atol_ * ss_tol_factor_ : ss_atol_); } void Solver::setAbsoluteToleranceSteadyState(const double atol) { @@ -803,8 +816,19 @@ void Solver::setAbsoluteToleranceSteadyState(const double atol) { ss_atol_ = static_cast(atol); } +double Solver::getSteadyStateSensiToleranceFactor() const { + return static_cast(ss_tol_sensi_factor_); +} + +void Solver::setSteadyStateSensiToleranceFactor(const double ss_tol_sensi_factor) { + if (ss_tol_sensi_factor < 0) + throw AmiException("ss_tol_sensi_factor must be a non-negative number"); + + ss_tol_sensi_factor_ = static_cast(ss_tol_sensi_factor); +} + double Solver::getRelativeToleranceSteadyStateSensi() const { - return static_cast(isNaN(ss_rtol_sensi_) ? rtol_ : ss_rtol_sensi_); + return static_cast(isNaN(ss_rtol_sensi_) ? rtol_ * ss_tol_sensi_factor_ : ss_rtol_sensi_); } void Solver::setRelativeToleranceSteadyStateSensi(const double rtol) { @@ -815,7 +839,7 @@ void Solver::setRelativeToleranceSteadyStateSensi(const double rtol) { } double Solver::getAbsoluteToleranceSteadyStateSensi() const { - return static_cast(isNaN(ss_atol_sensi_) ? atol_ : ss_atol_sensi_); + return static_cast(isNaN(ss_atol_sensi_) ? atol_ * ss_tol_sensi_factor_ : ss_atol_sensi_); } void Solver::setAbsoluteToleranceSteadyStateSensi(const double atol) { diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index be16d56981..336d98e1ea 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -445,7 +445,7 @@ realtype SteadystateProblem::getWrmsNorm(const AmiVector &x, N_VAddConst(ewt.getNVector(), atol, ewt.getNVector()); /* ewt = 1/ewt (ewt = 1/(rtol*x+atol)) */ N_VInv(ewt.getNVector(), ewt.getNVector()); - /* wrms = sqrt(sum((xdot/ewt)**2)) */ + /* wrms = sqrt(sum((xdot/ewt)**2)/n) where n = size of state vector */ return N_VWrmsNorm(const_cast(xdot.getNVector()), ewt.getNVector()); } diff --git a/tests/cpp/testOptions.h5 b/tests/cpp/testOptions.h5 index 407b6ff4fe09b3419c0a9c54a6cbb967be67c7a1..1794ea6e63da3850a0bffc2e755ac37fc555d383 100644 GIT binary patch delta 16075 zcmcgzf1FdrwZAiS7jBY8vKBG?+@P{Vpc;7k3_to_fai1NhsG8(Qi0f7U01-UwOtiT zkow+MvCVps*)Cv%f`%XM%Nq2q2s9SdNI`?ObyaGk1zmoqL8`8bzBe;xCiiA7{Kv!3 z$DDKSxpU7s-*e8)ojaR&t6%5Ne#rpt9M(Q^B__prxaIpp#w24{0fM2AYr6I=q zgjCbK>_7f5|6{lw&DM~=;epYrDk0{o{9i)ysoti$*W4w>P7tfhLCHs+2YDENv?x(e z1wlgcsVez2spp?lsaG%ghIt@wka~f41wT=^z-18-Xz>G}*d!HO>xrP0fjP++rxF&G zY`FKwr1qZ&A~Ohvh-R{9Q|eax{#>&>l0cMzXnv!*R1pP-gp=E>R7+H5Jj z5u>xrRa1tbsCkuVM7B7^sKW7UTH1%2DX)?BUOu^ew<{5au?uF)Ws|Zaqkmz;G(S1J z`@8{x6kJjEMkF6MYnFV@#`h%D_>n|29*}#uQoT#^M)ugfH(^{f)#DTQ$;W^3^bS&5 zA&l>X-DUk&+EyYXj9x^SB-HyU zaiLZGTCUB&#v`1H>@biod#vsCDwl2lzSIbkUawgb?$gi0IQgh~_=DIu;Tu>ZtTPFZr&mAUEJg>w2_O&-)ng!R*q2XHU#FM%!TeIPO7jD>)J@LpH2sdvqr)M|2?19jh^RY+^7Fk!p zg7F2g=ElQg@J7v=vk9BXZJTIL?Oo@JZv#z7sIl>^xv{kqh0G({*5i)t+@r(IKG+q8A1(P^Wd#U2 zB;Pp`{vH)6<5UUy&}P95MTum|JJL7tZx1G&(< zTvn_%P|idLYgf9M=zT6?Eddw1TG|)VlndRBn_smCyU|dBjn}|V-)rXn!@JzR!}7UD zd6AI#M zJ!}OYxE0Yavgj^!+^QyM1@bt|q;Gvd>Wp4suCr?tbb7cK>cpe}!C@hZ-rO2)14$Z! zWLZK2Rf1&5&<1WC>m8I`_H$rpkrBLkJak|--u#I~X7p!(%*4+)GM?^gc(SiAGN<(( zP9kR_XyqitBu74&camQ2^}~IlOh~#PGG`2%#+x~>6VhFS%v%^d!$ArDuW50~Cr5z9 z08q4K-%7m>$#>Q;?>-TgLP}>LbNk0xQW_#Wj@;+M=W=RgJ9g3M`XNy(M5|@8?Hrs& zgEaL=+26iVg-COMWDfZ{>v4>m9W`@Axn6M7c^sc6lJAeKJ_HRyl$qG%aoB6#va-%N zAI!8-a-k6dBvO4gGUpt448`cM(bAOoC}C*4&4mMkVjz*kXhG)MTy)VeYu!2^-`*TXwKjz&i(T#V}Q#R@6k*i#`@nhaJ)&6dXk_>xMX zIxpftEtnFhuS6gs`S?_jE2QFMEA?6>-*FN1K_XQM0aAr9_m#%nRWAIW!tG+(h5uCu zgPpnyPhzCPDmqTO+b|BWa*SMwRrE4}6;6+{_FvA6$zdy9<*lkSrh$z%Oa@L)1)Pqu z(HTx8JVO}7t0{7YAR`4Ao+>DgQ|wAX@^T7YDa(n;R5_9P6BJHICrZOlOaz=96Bj+g zCP*ArICalF6;2(;!R#nq8%4FGbrv$eymSurUyFjIG`j*e6sFHb(L6KwxOYlxkNMRx zKa1)rBiBPuyBZ$nk1|gsa-r)0xzG*eKxA$fB=0rFYX#{)=)yl@ltd2auCo^P#=Q|{ zH}fKM3|t9Zt`UIq&aw6^K&z1}bn7>PDl7n?X3hbi5^?cR-g=YtsykT8h_<&-NI6oe z+XcyxBAGh`$(N})ikPl!e&n^e3ayzjiPjA6Mzrab>|CV!#_dmQy z$apqu&5Yb`rSBCDaEbJ`d5OTrZ4Yd{xg6OG2HfOr=RN4e)L-R{jqSU!YHox}D@JA#c&w^aR6kk=T*DCprHOw8={MO3mw}HDT znjgJtJxg2O`+g&2qW5iNWV!b}D!pxZ6kxIk?g3LYBYWSYs`uHZLhsuGCfYEu-j|iV zkBu&R-_zU-v66^bg5%R5Z;4}#<+mD|h6?PC2+LeEQo>z{{i=NP%7+mX

z!YZ-WILMI;gjA5=JXGRXtu^?R-L@1bdENkJ|U zL0!`v3q~=&q~sbdE(cV4jt2-0*oI>7e}~`(rDwkyfII;|S`re;-U{llGfSvMq+&sK za5bQ#qqnp#rvC2`KZ$ge^NGg$p_<9vpsnuleGXLUOZYA4(X2|n&c6fSG>ZjU=27b* z$fMR{0t1eW_cKBAJPN-rNS;UOj~OZR$bVSk5j+fe)vI?8V`PG~C|V!O<^8BH{R9I(igZkH6*?(9$TD-~U^-EAyq0imcZ#Pbr>9A;d|V zKA@Y&py%E8h9A9WEUGgM*r)S?Js!UVDOZK*J@296X8xXt zC)^8^Vg?rP9K~uZ-#HRYGU25knjgB{lYCMO9Mm+&U1Cs`P~ovKn$l~Kd^2qsiP>`D z@i$NM3X>WLH$uw^^~KIt7~EVtQruf~Jdi9cTsAggOgeEtT zMx5}{`Tg(~lJJRuvrRWVL{MT$bA}*!q9g|}awnbJ9}h&ITklN5qlmR`0KT9XPJ8LZ zX?U~S9^Ja%k6|9Chp)gph*67q+hKcvwx#ro&3J@yE*L0V4%>6XE_4oUV;jmi|1G+I zCSFHM=Sjm@3SY=cvA8&=Y?yWdCuxK)If6pUJ5F#6BNx!fD10~8>a0DZ@CDV_e?DEk z18;XLFfPI(U-2ah>hwkU6w9l`rC@E84*LW1Z+rzdq~Q|G-&!VhOsx!hQS#ahkSlPs zi!1d4lCP)ljKsIHJzBsgHUaR7O%ZPLJvv-3+~f$)zhgJs36Lxez&8r~2Wl+g_&-(p0m;`dq5~J=p`Q2+a6Ne=8|o7?by*`xH)69- zeU_we7t!H?^fw8T7h$9Uq{m(O4n}RalDFX>7&zZ(4SopsN4U@k^=A)kjeJ7Y-`1b8 z`U~DA>(65wBf6Zcuo#SS!rQP2;4Lixcztwf9M{rFki8xMqHI-e5hqu2VY3i&c{{|d^z8PLohK9y(&{^w(+AOZC#Ku?;eBcTARka!TRi$1;@^O0XK3tnkJwBQ3 z-#iSUHLnrwbI`QM1j#`Qnj#Xp`&hX?!pQAbEYDF(KP(;2JPfFnhi>VCn#(JVh5h}J z(I76gTGwsDYuT+B5C3)jjkSh?4t(AH~$w%lr83kfP4x5XAz6IyG$F#7MxbR-#Ji`!q z>C$#SS+lp=mzYpN*qoEzCUXE`dEjk5AuPQNCX^v0nVs026Mc(Dox?Y;5J53q55bgcT~`Om*=#-_maC>`?wzS848 z2m=J)$0BWy0Bgu~V!~zM$oq8nM|k0qF8OU$by3HyqPwbvR{CogL+KSHU((>oA)&2l ztrsdwL6DPtaSX`)QqO<3;Lj|#?|Kku4YbBfo;c|jn)v5n2V}1Cp`E#W>p6M^vwgD# zgGybM*7^M|3(Oo89&U_|$(cCXL-;GdOgr_^#afm6{eDQG@AL zhp0G^Aznnx2~$qC_1;wjcYyziUc|X>B?Ur+{|yMV6k18iX9!XL<9u&Ae-hb3@_h-* zX?UqTOOU)KtaqUAp*j&P7!ugxoE?{r+)OSc9M7nwtIs=hS9>>Vm1Acr?h3BqF zf3TxKkA=x3`X)wWkjvV36P!u7EHNK&&1uOQv@~8nkFK6b){*W@ zdU$>@c%JN0*~Y&}&o6d(K2CVfwv=z1Okic%@t!S!x_JZ-aW!lL)hK;# z7P-=smWuFK>xdlA&7iZRWZ{yiR8@uHxK%tXATMpkbc_?dVH?pJkW{Kr4!hCc4lQh<*$N#ie2d(NzMk)$egq>^>IeGV*z=w2ACT*{5m-_u?yb+%?PE z)4&=lF|L+X<#!%so9e2$8@vZac-FW>E ze|+oh8l8I%3It03Uai$jt$t3eJASnAj+>RP2HI7nf7dj_oUn4Z#}QqTVYk4e9T`T9 z2HvagbQPN_g1MkaFdJns8+rs&^65rrFc+#|{=^REB0HG#?O^^?1(Pj(yB4Y#-l}3~ zkKQAOi&YHi?jt0QQ*mqhR&w2G4(89`pb}*AFn4QQ9uORpL!d?HZLfWiCYCoe#NLNaek$L!< zI7$72WamV`f0`p|*EIcdGp;k%{F9KjK}b`Rp^-%3OJkae`IpAjdJ-gBsvKZ(EDNuJ z0Bb$Ergs&LA6(-IaQ#|B=d|lvDB2{alHBh^O!;eC=@AiWe)I}&RN-a08eb~|n`C$M z1%=xl+Zy<^noD@S7yd^D7IHQ4*y+txIw+*qx!d8P@hpMJKHW32?1c^j+pr5RMZ`xz z?U^2-d>TU8l?5%8tKKcZ2cLF^a;pmEE;TE^Xou3WL%Gu#%B?Dt>~08Ud2DOoE5wh@ zwyVWKM|VCC4~-%j%5A4Nl%j{$?}mpq>|v2$>%CWib(VcIyTp5F?5}VwV)B*UzH%hR z@sj#iAd;*pXiX)QYKgp^Oza;|57x|f@Ao<#^ASki6+s{Ku*>vBTl6*cd||!kf(4TP|1mArMk}_hAAt&AkJUS20bp+fQpN_3Wy1{~ewg zJ4oQPK^IYb=?4khV7=27{S-7HDNBbdLmBQO;sA*E1Zb&Hh9n>T8szFtzw=BXw_-gUL#G3X9JM>DLd+ROj zb5>rh-%SoMs1Ld>CxHI zrCyA_ln;qt`@p)cT4!hKc713q>dPX84}1ur!um(n>Kc7M@;KI<+&0|N#~Od1-XD{c zn=O}F5;P5xEnjpduUA*^j5K^E9)SDkA_;=SKuaY-RPs&df?Op*w_MsPJ--L8hx;(A zWfDYdbiNGSf(LM~#JgpFs33U~Xni?JlV`C7;Xp>pEFg&*-Fncg_dz&SLzfTNH+vjo zTOOp>+G+5DyR#OatMZ^!V^#k_Kh0|wM^LEew|X5KZcTlJ$2vWlIU>>k>u1u57wYRs zeVuSq%y}a?DNe7H&J`p-!WIw|PlTrPLC6lTG*dK4ucaapT1PRlNCfhaGj(fDwa(6l z9r%`RE$+*rMD%_#UgSJRt*Oy(L7v37Au6d0b&(9kiy#?r?n3MS@%qbXNkC3%3crXw z#ff*22Kgr9m;GAlIr_KI@JqdvT+S)IyuG{=8eVxJu}fHHzXGK}pygixlM`}4pX_RC zz_8urghoc{MWw4s&*8ez$VZJero4f;_$L8l9cz^ul+B&?V)<5QnRx_0Id) zpGinQH5KHlAYyHmT`@Nnt~a;YT`>_71(79*ut1v6L#0av$#+`mVcqjD57u?!q@LX9C#5Yd2W^5&=Z|zB^v5&rPJ1SyJ3r&|5QTQHn z1S(?mIH-vD4B;`KZxdH=k|wTXADgEMDGy2VN+B(9X+%ibt}4pVNt0DcF(?GpSW^2LEuvK`m?J7_&Pfc_?}|= zTAjZ2F5PI@n9nGd3iFZM-gFVoF}{Rn!r2%vf{8NtAOtKT4xQtr zm5Jb6A|}d7T@O_P>@|T8J`BGq5i)XXM4yc$!{KieB%dmSs|3mWLYSfl zRbPm#VPZ&z%|BPku=+lg4BHq5$*`YM&=(wAYvQXAS||HL(^}{YktZQ;DY0(qdb~0n zLhE#T&n0N#jAFPCsTBh@q!O{=cZl!mTg5B`m=pjnmKkWZi9Y=Bw94hqZ?+te(vllcMJusT< zH(Chx@29)p)0;@>WjTe^DA&DmXFSb6@!W-{zif>@qz^^-hL`E|!}^NTwJl3GvA>*g z%V)D)dGVJt=pa-u@xRL(`n=kIj#l@_5-ncWuH^gQM5}lJ&V^Oc$!1o+@!N4kQ09+@u@(LCv)g?LGoo`f}G`5zok?T~`ge+b}oJih0QfTN!YIC1onyU|; zRy4&1kt42fyzy%pXGa-DztZ{Nt1VtxtZ*is+}qB+QlIo2 zcxo_=3EzN2HC5_;UhI!GpHaoX0T1w3gO(b-Uh)mJ#a$(vOtx{w8OY4kdr9I<5jcL_ zp|%{tfLD-wzDW*ZnY_kl#R64r Uh=2H`yGK^Bqlv?+;BTw{4=uR{w*UYD delta 17182 zcmd5^d301o)_+yi3B9m12u=rqv<12{IDCk8kMgrNnL))8%Bh>eQ2Akr*q z3m#2)kcuqYO2CM5>`a3%%3>7Ih@!1BW3%XV#w~2xB1=Y5zpA=b>DSCT-*=o7&o}-- zZryiZy}G~Ot$XWLy*9a{d-y>2iXt4oq`qGjCQY~F%E`SQp*yg~|Er#k=?|r(w=Fgg z6sdSu@(FWjQ8fN-pM3n+NxuAe+Pe{p=jyTHR+Q4m zItt^d;np17I^~U(5iGA1rS})3rMfuRVYw%w(3V&Pmt_%JJ zBXft<$jZVIoL5o>7V3_mY*)m;U33S=ffZ`t;fLC9;5X8;-Z;$B!=aV&-FjziYkRKV zyt2^dS_PvjHfGy#2TY@d$QD1pdI7@z)qFsy5-Sd`WfBcujBMyhuJc@H#lbaP*T{7G z&{c&v`gCP|x$fUGH$HCl*A{yl8qeyb9SVqNgJvwd9QVM)wN9h+z94gK={P($MJq2y zMe%{_euHsxy*^`6`LNh4*+3;_yLAibt?)++uFJNomhFc1mG!6e^oAKY@x0!#WD`p0 zU26+*!r&uLUB$e3PI-~~Yi0fR<)zNJ@A*kql73#Jm)uNt=+#e@+MOGr_o${na6^qf z(|Nq#_#)h1foPfh9qkg8kMsMBf=FI8tZEijnaBwU0xT*uO>CbP?U;epOFJTw_tE||%L-L0G zFR(jTFbdteba{wB~Nd)XHu5lZ%091_1kcx z$~393AmGdD>%T&J*)A6ee~BPfMyA81lwP?j-%)ZWv;wmceJ39m*x!;Pcw3`eE+Yka zj|?L0eOvF>RK6f7bEKM3TIL3Erx)(m8Yh5-vTKoTs}D3DY3I8<+b%5I(JMf&JW8(- z?e1-J-?bJDil0W6_3Nqsdg8MA4?wReih7bF^u9je#=Z58{g9p->B1#@M4-Ey!3OTu zv=Z44dGU9^V>dT^q5tA3Wg8Mzer_ zPX&D0eira?Phy4@H+>LalgvdV$?kh`8j9O9-L`9pAHroRk$LPH|5#Qiw0;g=m|_`6 zKGP33CyDDr7G`JAI#jGTf0j>!lbH)&aTj9pWO&iWd`jO_3$xb#v&L(ZlMeeFS78^a3i}Ll@XExkHxag_1(XJg3qtb>u z3Js%@q&hKi7H~P@G;q1I6@1>?3S2r|UF2iCB7Ibcz9iTk8H49q&c(PZfHhy%zfO3dkBdYQQ|(IiQ5>RGp*Q#;`4vVMaV@veup;OGabkK zk1uSWdeddUJ3auds5KWDiYqMn&3Y)jwu@6<+S1GP*1f?8Bey9N=|Q#?m`g&5j*N`8I! zaxI^=f&0XDPH!m0L`kpZR!Q%YKS`%o{0Zod-4D-C%tLgUpV7NO(VJ!(vF2}jOw!9= z*NA^%8=<1E;s=4z$a1iPe>E*8+ldvRQH+*6kZl)~?TQE4a5Ie7E&@htmvV=ovQ`M+ zIHP4@q2;9&FJt01BmE#^GuN|Ny60UCoSJ&q`;VWa$9&!5rM;r)qIP=2JwzHRGu)E! zVzl`f(pj0!crVIq6*!%lz&5-BNZz>^=5KmAR9RdtykupTKC_(h%c?AXlB;-`T_vuG zGP|0IlHCyxNp^QWEZyJxFtBS{;zFH$R@K?3CApp$j11~*O&hm-!Ao4v$~wC_n={wT z(wzqQtGJCloo$zt?bc`5U^Co0H^@5s7w!(!*-gSTUT4<{iDfdD#$J!2{L02Cl-dol z)NX)Mo9Vm#nH42>bB0*lvO%V1lj~4&Uy)FlOKt;{+~=U=ZV_Jcl8d(r-&h2r7+tuT z&1q3{w{fkgx!Z-QsJX8(QP$jARdY8YS#vi+&25BdN8f;&yGu>o#y2Fp)gOU{thr6u z^oDoIn!7KXwy)MfR=ap^$E;iN-y`i^d&HaFT12y+!A&Q#&mJm6=5&Fyy z;FrJI5N)@NUgvIEb~F9AZ$ieAr1AYJRE2c^%8zo3z|G8LW!3fo%+|8qVz!dt zr^*Ww{s$LD#l;7OR+QY7(2A1VEVQEJ{*#Gv)+XMOv$kxH^xMBjmfX)I?q9;{xa5c; zwurrhZ@J8}$e7TKMxh1i&e?>Rmerz)M&TVsYm^;=wWhUX9w027(k>!jv6Nlrl?x2Q z(#{8gm)V94w?!h{TBbUb)upXN*=E(IvVK0Tn2FpZ$`Ynr+*eJqbRG5+7^_>-vh75$ z!>NKXEVluaUC|+sYJIHhOfa^54Z+xQN|3;*X*nZA9*h-7g~%6xt*4o|O}t~1;5U+y z@HY^Nrq0_r5?g=dk=PkKOn%bApOFG927k>_=s}pjn&HN8APQ3i^BKm6pSEz3rf0Kt z-|QdxYGDPdjk!YUy^Ae}ItWu9gl;PnWe|G5l|k5eT)Lb%4nb(zs@As*FR*ro;W*`n zV=M<-akYcKnuG6kxH`(QcE)tF)lz`FVR)TGo#|XFE+3wj``wBx(+wh0Q~*Ycp(AD2 zB-<^6*pN@c7EQg{C7V~-XK%;JLVp8T>ME&{y@qt-6jc5WxMz1_qU z(@&8Z#CI;Hv$|L|;XVRG9z(v1gvet^>&``5;bb7^GjRuf*vVd#T8(F%_(mHJpG(i= z;W~$@W8;Ue!a3} z4;u?cZ9yU1EjNH(QRyF?Z5NU4+N;@UPf|f;<4{a*dYZ&Z<8|N{TU}kxC7Kw|_HM2L zC37}W%EUG*7Y^oBwhRR-O&$MgRN4|{;7x)yKShSFd=+`Hb)+%mI<{_fG-VJFi;&<1 zur;%59S`jG8VT$+{89K8-D<;rkXY>c#|tsdlvBCH$!)ns=)Bu(_*N#dM8Uajnkcz# znFQRLj;QE)4oiV`28=-@h%T4}Z^yKMGwNr&I299w)3bnk+%TQ)z7;<>*>nOYQh-~Q z4DV~%ECysd^1AF_+9maDJ8cED-BoPRd!V$G1hWIGdKM&y;KJ;Ct7W@kG3Zr9rJ>oj zW9*v=iqI_fjFbGVf>g}LdZ&^&skmRD#&QFj63z25-&wKYJB5-@p6FaAZZm3U;9IRY zHrwcUH`~<04YNf=*)ts*`|n{>uDxDQZlt06a0Q8mWm+&fN=H4xHWnTNhs=SKTnvGe zyaxg&@-TNozJUpUc-HeQ5!yyxWuD+751oceah(^LrsYhOu~YpBkDaZ=C|iujv9rGW zjG74UjBss;Zh91-!oK?z97`d#OkDD}C-kJu8{R6JH(bxbxS#4QF{1t`YtIvd`=q-{oI zCFj_+RyyTg3mj*<{Z%-&)?t0}#QTUhiZcgM4wKLNsAz#T2ChoJspkT=H~kXuT^0qt z)AigRcq`k=#i}0c&E+NOB&&B@#Pza^81vhNme&CfUHBSi?`ma`F0^L`KOubNMc@B-;W#nb7 zC)(l9F~fEUfiuU@82FwcMe9Gu*E=eXfk&}#*b0H2A9zk+v4_IAa~4V`Z0=?N%G;i)Fcq{r)UF*=`*SdbhH3i+1LeF$oH-^#am%%9uuQ z6P%*645*f{V%d}Y4zvB!S?2}6%?{FP*>3m_*Vp$({Mdq4xoePt$i>*+SLy7{0n` zr$Y?&H3&)N5?G>|kQP(k)5$nIjj$simz^$+lY4FHo>?e$A=xN(B!H3~tQ3LOvfa>; zWTAu;PI@qhY$6F4!vMA=5`DM`Aaxd60aBg@8TzGs7Vu>@;i2CSmzTp6pZ&VvdzOh0SxpOEN zA$N{8ffO9BKne{9kQ}Bbu)WO9G|?Md4bKQNqaF9D1LI66IY&yAOy7m+{de@uYw&@u zD(#gF;U9u&S`+^agr-V~3%3lS-GiARZylZne~5c%-tk9)WhzG5A|8tTip&E!)m@4z!DVz~7XzTO!&E+E($d zhos>n@V8YO7Hx#KdT^O+2kStuN`^O>?MOhjBX@ybVds6;Y}ZNC0#XsFVc~EI>}G^# zfRS9+Ol=EE@x@7C3N2kr?(8|q+>@pi$#%&U5s~867W@5l%F`^r`2!*Z#e$_=gdktA zP&9$C!vwxwna-5EjIw9QK@{KIwO6kA<|^Vh2HZ?KBU0;^L1}tOyY6{GnhV}3%YX?z zh?_YwZffs>%objpg;3}$gysN3cb|ojY?sV2BXqw+XkiASM@$HnD})v)gxESha=!xR zeg#aX@Ae1A7L0%=@ZL`crSL#Jj~G{9huN`G;8D03Stcjh;|i>0609{@U{#(4*5d%y zV`qUS+a-^i!CEfBLVh~%N(&q#`=chX9#LSeR$wKO5q*>_&Bc-B3bN%2vdk0OA2Jer zhNrgsuy=*AW+j;-KEZRCCYy+>l3E2-4L}tf1JgPUr`9FwdG^1lYBlkY=*BGENcuUb z4AHYvSZuxI>Nyi|EiXWnwZu7E;`pQG6(O>!!1v3d61mn*3ON>L&N_*irw%Z)XF50L zzR=#w3|s4~a_C}viLSa34a6_1FiUTNi_Yy3W+ksn>^$26JL~&kkQ7J4(H&>O)d+A! zwu6;Qs4>}2G@9YsDdBooU{uBq=|421^@>7kj~T693awoVt+u|6xjTuweHUf6^hf`W z)|(2gmc4M%`Y{okDrmn%tLS4uOKUNsRedlUt(JX&R`}y=v}$F$X`czL)F*&e>I=Y# zZQLE79I=oFcLcr~!tOMk`vrRkXcZr# zvtlm-t>R=hbgAP+Z1y+hW>@aCY&&hBSFin{bT-6M*={_}(1IOeoIC}2FS&B4r4Y5O zthD4)XA3)3K}T35I#=27&m3CAdz9SDHbYJmfgKJgmv4zgp7ZkEdWn`q)EG8Zp7O-KJsYDmh# zbn@t_BNrh8Z9)`hPo0Ix&z|UPC#xWQzrgJj$_joU)s=}OMW^JRR)*pj3;X=&wto?$ z(^!t@5+(t%HB4F?viup@Gb=oi(t|Ez&|C}KB`wp;{<0f3mhjnF9OH-2PyrysIyNQ==t1l zJ{{eCnNG$GF79PKu*}khv-8+vmKW`revO%@Nd7NNUB6Dic=`htDqyZf15qzyK)z)P zdkMH7|K6U#j5A?NTn6?KGqF+&`~2uC<}w5g*~(xV^D|jTTmuD6We+l#E*CKIH>Ji) zxmd-OzvU9&WGcBrSVp<>N+G((uwDbdWAezVYRf#-+jyYbl21?{W6f%gW*?(&4g4K& z0}B>MSuu<$ILN{eF!41OnOiCb8LwVuc@HgUl4VXYY$|(V3I)W!n${>2ys|STd|~01 zc3!zUQ+B4+ny_%n3$gdb$TDvvNl<95AHa~qvcm?MA618`OuiHrF==p*=C2(kB9hR9l74JHSdv7q`#NJ|LJ2CD2KpI#Y0-Qo%2Ai;-=`i=`wCo^94L?7NInKeF1}qe-(n~aV>ln5*DKPEM%du z6}5Sh(2Cl;n29oQE}9_&r~6D9I9JS+*p8oRtO;3WA$w*+1?uurS(hJCb-5Y_4t4pl zws?V$442Bf{8V-v6v=i;6!a=yVshI?#Z^*m<2=d+t?KgQsxDW7Q>-r2O3Ee<|H5z? z*FjyTOjI#b{J5;kk3(I~^!uOG`j>_^gQ!em1Tw! zl;w2-C0>?ka}8t~(ir7RQI*$oNmS(x!cbJ@=Y$AVxt7ULm7kDRxk6Uu3R#saWR5Xy zs__F!{E|vB6?LW*6W9g=_gZxOxNy$YBXxqqkZ zoG-$4&r597B4B3XxQLfccndhie#b-6HVH&J#+7=QeJZk*Nh(}|aXAet;*erY{r)$l zn1(mN-AszXsW=1;|E;3K1e^F&`7JOv#fL&Sr9DqmOgQw&;pZNcKuj0DGoq{ui+4CV)_eJzRaFZJu8M z1lFC91nxYMA<%>|9h}R)3bH{@#hKYdLe?Di>FNBcyDgK^Njl&sOARS9VHHY0Q+DDKuNuD2LKkDYE=-yQI|&0$$%)_Vw8VADX^v}+2&&Bga3P0|Dq&kg zJ#yp;f2tdOwL{Lm4*xkY%5aYy-Lb9!_f!Va9(1lfr~JHPX`n*HbVEI#G?MPk@?<$_sok)+BI7zfOcU5fqM*9nq|B&PL zt~vQ4Qx-7QHgRwt?UqxDV||RH-EunHaHOo!ew}S=dAi=E7kw) U|HwcS`8iG!8<4{e;ZLIf0lAp@ vec1{ 1, 2, 4, 3 }; diff --git a/tests/generateTestConfig/example_steadystate.py b/tests/generateTestConfig/example_steadystate.py index 517508f9b2..1142c70732 100755 --- a/tests/generateTestConfig/example_steadystate.py +++ b/tests/generateTestConfig/example_steadystate.py @@ -21,6 +21,9 @@ def writeNoSensi(filename): ex = ExampleSteadystate() ex.modelOptions['ts'] = np.append(np.linspace(0, 100, 50), np.inf) + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 + ex.writeToFile(filename, '/model_steadystate/nosensi/') @@ -30,6 +33,8 @@ def writeSensiForward(filename): ex.modelOptions['ts'] = np.append(np.linspace(0, 100, 50), np.inf) ex.solverOptions['sens_ind'] = np.arange(0, ex.numP) ex.solverOptions['sensi'] = 1 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiforward/') @@ -40,6 +45,8 @@ def writeSensiForwardPlist(filename): ex.modelOptions['ts'] = np.append(np.linspace(0, 100, 50), np.inf) ex.solverOptions['sens_ind'] = [3, 1, 2, 4] ex.solverOptions['sensi'] = 1 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiforwardplist/') @@ -51,6 +58,8 @@ def writeSensiForwardDense(filename): ex.solverOptions['sens_ind'] = np.arange(0, ex.numP) ex.solverOptions['sensi'] = 1 ex.solverOptions['linsol'] = 1 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiforwarddense/') @@ -62,6 +71,8 @@ def writeSensiForwardErrorInt(filename): ex.solverOptions['sensi'] = 1 ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 100 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiforwarderrorint/') @@ -75,6 +86,8 @@ def writeSensiForwardErrorNewt(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 100 ex.solverOptions['newton_maxsteps'] = 2 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiforwarderrornewt/') @@ -110,6 +123,8 @@ def writeSensiFwdNewtonPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 20 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensifwdnewtonpreeq/') @@ -145,6 +160,8 @@ def writeSensiAdjNewtonPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 20 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiadjnewtonpreeq/') @@ -180,6 +197,8 @@ def writeSensiFwdSimPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensifwdsimpreeq/') @@ -216,6 +235,8 @@ def writeSensiFwdSimPreeqFSA(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensifwdsimpreeqFSA/') @@ -251,6 +272,8 @@ def writeSensiAdjSimPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiadjsimpreeq/') @@ -287,6 +310,8 @@ def writeSensiAdjSimPreeqFSA(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiadjsimpreeqFSA/') @@ -308,6 +333,8 @@ def writeSensiFwdByhandPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensifwdbyhandpreeq/') @@ -347,6 +374,8 @@ def writeSensiAdjByhandPreeq(filename): ex.solverOptions['linsol'] = 9 ex.solverOptions['maxsteps'] = 10000 ex.solverOptions['newton_maxsteps'] = 0 + ex.solverOptions['ss_tol_factor'] = 1.0 + ex.solverOptions['ss_tol_sensi_factor'] = 1.0 ex.writeToFile(filename, '/model_steadystate/sensiadjbyhandpreeq/') diff --git a/tests/petab_test_suite/conftest.py b/tests/petab_test_suite/conftest.py index abcd9d821d..84aca11026 100644 --- a/tests/petab_test_suite/conftest.py +++ b/tests/petab_test_suite/conftest.py @@ -56,15 +56,12 @@ def pytest_generate_tests(metafunc): if metafunc.config.getoption("--only-sbml"): test_numbers = test_numbers if test_numbers else get_cases("sbml") argvalues = [(case, 'sbml') for case in test_numbers] - elif metafunc.config.getoption("--only-pysb"): test_numbers = test_numbers if test_numbers else get_cases("pysb") argvalues = [(case, 'pysb') for case in test_numbers] else: argvalues = [] for format in ('sbml', 'pysb'): - test_numbers = test_numbers if test_numbers else get_cases( - format) - argvalues += [(case, format) for case in test_numbers] - + argvalues.extend((case, format) + for case in test_numbers or get_cases(format)) metafunc.parametrize("case,model_type", argvalues) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 33a6c47e18..6c31ebbfa5 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -68,9 +68,16 @@ def _test_case(case, model_type): model = import_petab_problem( problem, model_output_dir=model_output_dir, force_compile=True) + solver = model.getSolver() + solver.setSteadyStateToleranceFactor(1.0) # simulate - ret = simulate_petab(problem, model, log_level=logging.DEBUG) + ret = simulate_petab( + problem, + model, + solver=solver, + log_level=logging.DEBUG, + ) rdatas = ret['rdatas'] chi2 = sum(rdata['chi2'] for rdata in rdatas) From 5e986f3c9c5094dbdb0cc9fa715023ab36bb856b Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 02:07:53 +0200 Subject: [PATCH 04/11] Add support for pathlib.Path (#1769) Closes #1721 --- python/amici/__init__.py | 7 +++--- python/amici/ode_export.py | 5 ++-- python/amici/petab_import.py | 40 +++++++++++++++---------------- python/amici/petab_import_pysb.py | 27 +++++++++++++-------- python/amici/pysb_import.py | 8 +++---- python/amici/sbml_import.py | 11 +++++---- 6 files changed, 53 insertions(+), 45 deletions(-) diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 9bf1f28855..53502f9750 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -24,8 +24,9 @@ import os import re import sys +from pathlib import Path from types import ModuleType as ModelModule -from typing import Optional +from typing import Optional, Union def _get_amici_path(): @@ -127,8 +128,8 @@ def getModel(self) -> amici.Model: class add_path: """Context manager for temporarily changing PYTHONPATH""" - def __init__(self, path: str): - self.path: str = path + def __init__(self, path: Union[str, Path]): + self.path: str = str(path) def __enter__(self): if self.path: diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index ff25b184c2..a1f40409b6 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -19,6 +19,7 @@ import sys from dataclasses import dataclass from itertools import chain, starmap +from pathlib import Path from string import Template from typing import (Any, Callable, Dict, List, Optional, Sequence, Set, Tuple, Union) @@ -2208,7 +2209,7 @@ class ODEExporter: def __init__( self, ode_model: ODEModel, - outdir: Optional[str] = None, + outdir: Optional[Union[Path, str]] = None, verbose: Optional[Union[bool, int]] = False, assume_pow_positivity: Optional[bool] = False, compiler: Optional[str] = None, @@ -3090,7 +3091,7 @@ def _write_module_setup(self) -> None: template_data ) - def set_paths(self, output_dir: Optional[str] = None) -> None: + def set_paths(self, output_dir: Optional[Union[str, Path]] = None) -> None: """ Set output paths for the model and create if necessary diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index e2b69e424a..c227595a11 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -16,6 +16,7 @@ from _collections import OrderedDict from itertools import chain from typing import List, Dict, Union, Optional, Tuple +from pathlib import Path import amici import libsbml @@ -229,7 +230,7 @@ def constant_species_to_parameters(sbml_model: 'libsbml.Model') -> List[str]: def import_petab_problem( petab_problem: petab.Problem, - model_output_dir: str = None, + model_output_dir: Union[str, Path, None] = None, model_name: str = None, force_compile: bool = False, **kwargs) -> 'amici.Model': @@ -321,31 +322,25 @@ def import_petab_problem( return model -def _create_model_output_dir_name(sbml_model: 'libsbml.Model') -> str: +def _create_model_output_dir_name(sbml_model: 'libsbml.Model') -> Path: """ Find a folder for storing the compiled amici model. If possible, use the sbml model id, otherwise create a random folder. The folder will be located in the `amici_models` subfolder of the current folder. """ - BASE_DIR = os.path.abspath("amici_models") - - # create base directory - if not os.path.exists(BASE_DIR): - os.makedirs(BASE_DIR) - + BASE_DIR = Path("amici_models").absolute() + BASE_DIR.mkdir(exist_ok=True) # try sbml model id sbml_model_id = sbml_model.getId() if sbml_model_id: - model_output_dir = os.path.join(BASE_DIR, sbml_model_id) - else: - # create random folder name - model_output_dir = tempfile.mkdtemp(dir=BASE_DIR) + return BASE_DIR / sbml_model_id - return model_output_dir + # create random folder name + return Path(tempfile.mkdtemp(dir=BASE_DIR)) -def _create_model_name(folder: str) -> str: +def _create_model_name(folder: Union[str, Path]) -> str: """ Create a name for the model. Just re-use the last part of the folder. @@ -353,7 +348,10 @@ def _create_model_name(folder: str) -> str: return os.path.split(os.path.normpath(folder))[-1] -def _can_import_model(model_name: str, model_output_dir: str) -> bool: +def _can_import_model( + model_name: str, + model_output_dir: Union[str, Path] +) -> bool: """ Check whether a module of that name can already be imported. """ @@ -370,12 +368,12 @@ def _can_import_model(model_name: str, model_output_dir: str) -> bool: @log_execution_time('Importing PEtab model', logger) def import_model_sbml( - sbml_model: Union[str, 'libsbml.Model'], - condition_table: Optional[Union[str, pd.DataFrame]] = None, - observable_table: Optional[Union[str, pd.DataFrame]] = None, - measurement_table: Optional[Union[str, pd.DataFrame]] = None, + sbml_model: Union[str, Path, 'libsbml.Model'], + condition_table: Optional[Union[str, Path, pd.DataFrame]] = None, + observable_table: Optional[Union[str, Path, pd.DataFrame]] = None, + measurement_table: Optional[Union[str, Path, pd.DataFrame]] = None, model_name: Optional[str] = None, - model_output_dir: Optional[str] = None, + model_output_dir: Optional[Union[str, Path]] = None, verbose: Optional[Union[bool, int]] = True, allow_reinit_fixpar_initcond: bool = True, **kwargs) -> amici.SbmlImporter: @@ -484,7 +482,7 @@ def import_model_sbml( logger.info(f'Observables: {len(observables)}') logger.info(f'Sigmas: {len(sigmas)}') - if not len(sigmas) == len(observables): + if len(sigmas) != len(observables): raise AssertionError( f'Number of provided observables ({len(observables)}) and sigmas ' f'({len(sigmas)}) do not match.') diff --git a/python/amici/petab_import_pysb.py b/python/amici/petab_import_pysb.py index 9f71ad11b8..170a9ed45a 100644 --- a/python/amici/petab_import_pysb.py +++ b/python/amici/petab_import_pysb.py @@ -9,7 +9,7 @@ import os from itertools import chain from pathlib import Path -from typing import Dict, Iterable, List, Optional, Tuple, Union +from typing import Dict, Iterable, Optional, Tuple, Union import libsbml import petab @@ -105,13 +105,20 @@ def _add_observation_model(self): local_syms[sigma_id] = sigma_expr @staticmethod - def from_files(condition_file: str = None, - measurement_file: Union[str, Iterable[str]] = None, - parameter_file: Union[str, List[str]] = None, - visualization_files: Union[str, Iterable[str]] = None, - observable_files: Union[str, Iterable[str]] = None, - pysb_model_file: str = None, - flatten: bool = False) -> 'PysbPetabProblem': + def from_files( + condition_file: + Union[str, Path, Iterable[Union[str, Path]]] = None, + measurement_file: + Union[str, Path, Iterable[Union[str, Path]]] = None, + parameter_file: + Union[str, Path, Iterable[Union[str, Path]]] = None, + visualization_files: + Union[str, Path, Iterable[Union[str, Path]]] = None, + observable_files: + Union[str, Path, Iterable[Union[str, Path]]] = None, + pysb_model_file: Union[str, Path] = None, + flatten: bool = False + ) -> 'PysbPetabProblem': """ Factory method to load model and tables from files. @@ -176,7 +183,7 @@ def from_files(condition_file: str = None, ) @staticmethod - def from_yaml(yaml_config: Union[Dict, str], + def from_yaml(yaml_config: Union[Dict, Path, str], flatten: bool = False) -> 'PysbPetabProblem': """ Factory method to load model and tables as specified by YAML file. @@ -305,7 +312,7 @@ def create_dummy_sbml( @log_execution_time('Importing PEtab model', logger) def import_model_pysb( petab_problem: PysbPetabProblem, - model_output_dir: Optional[str] = None, + model_output_dir: Optional[Union[str, Path]] = None, verbose: Optional[Union[bool, int]] = True, **kwargs ) -> None: diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py index 90e5302f3e..e2185f0c10 100644 --- a/python/amici/pysb_import.py +++ b/python/amici/pysb_import.py @@ -9,6 +9,7 @@ import logging import os import sys +from pathlib import Path from typing import (Any, Callable, Dict, Iterable, List, Optional, Set, Tuple, Union) @@ -22,8 +23,7 @@ _parse_special_functions, generate_measurement_symbol, noise_distribution_to_cost_function, - noise_distribution_to_observable_transformation, - RESERVED_SYMBOLS) + noise_distribution_to_observable_transformation) from .logging import get_logger, log_execution_time, set_log_level from .ode_export import (Constant, Expression, LogLikelihood, ODEExporter, ODEModel, Observable, Parameter, SigmaY, State) @@ -36,7 +36,7 @@ def pysb2amici( model: pysb.Model, - output_dir: str = None, + output_dir: Optional[Union[str, Path]] = None, observables: List[str] = None, constant_parameters: List[str] = None, sigmas: Dict[str, str] = None, @@ -1369,7 +1369,7 @@ def _get_changed_stoichiometries( return changed_stoichiometries -def pysb_model_from_path(pysb_model_file: str) -> pysb.Model: +def pysb_model_from_path(pysb_model_file: Union[str, Path]) -> pysb.Model: """Load a pysb model module and return the :class:`pysb.Model` instance :param pysb_model_file: Full or relative path to the PySB model module diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index df00522698..a5fd812a12 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -11,6 +11,7 @@ import os import re import warnings +from pathlib import Path from typing import (Any, Callable, Dict, Iterable, List, Optional, Tuple, Union) @@ -19,7 +20,7 @@ from . import has_clibs from .constants import SymbolId -from .import_utils import (CircularDependencyError, RESERVED_SYMBOLS, +from .import_utils import (RESERVED_SYMBOLS, _check_unsupported_functions, _get_str_symbol_identifiers, _parse_special_functions, @@ -115,7 +116,7 @@ class SbmlImporter: """ def __init__(self, - sbml_source: Union[str, sbml.Model], + sbml_source: Union[str, Path, sbml.Model], show_sbml_warnings: bool = False, from_file: bool = True) -> None: """ @@ -138,7 +139,7 @@ def __init__(self, else: self.sbml_reader: sbml.SBMLReader = sbml.SBMLReader() if from_file: - sbml_doc = self.sbml_reader.readSBMLFromFile(sbml_source) + sbml_doc = self.sbml_reader.readSBMLFromFile(str(sbml_source)) else: sbml_doc = self.sbml_reader.readSBMLFromString(sbml_source) self.sbml_doc = sbml_doc @@ -202,11 +203,11 @@ def _reset_symbols(self) -> None: Reset the symbols attribute to default values """ self.symbols = copy.deepcopy(default_symbols) - self._local_symbols = dict() + self._local_symbols = {} def sbml2amici(self, model_name: str = None, - output_dir: str = None, + output_dir: Union[str, Path] = None, observables: Dict[str, Dict[str, str]] = None, constant_parameters: Iterable[str] = None, sigmas: Dict[str, Union[str, float]] = None, From 93c7304a03746fb01081cad01ee8d31ed724beff Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 02:09:47 +0200 Subject: [PATCH 05/11] Fix deprecation warning for pandas.DataFrame.append in rdatas_to_measurement_df (#1770) ... and make more efficient --- python/amici/petab_objective.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py index 0b57d2739e..7ed62bb509 100644 --- a/python/amici/petab_objective.py +++ b/python/amici/petab_objective.py @@ -711,14 +711,11 @@ def rdatas_to_measurement_df( :return: A dataframe built from the rdatas in the format of `measurement_df`. """ - - df = pd.DataFrame(columns=list(measurement_df.columns)) - simulation_conditions = petab.get_simulation_conditions( measurement_df) observable_ids = model.getObservableIds() - + rows = [] # iterate over conditions for (_, condition), rdata in zip(simulation_conditions.iterrows(), rdatas): # current simulation matrix @@ -747,10 +744,9 @@ def rdatas_to_measurement_df( # change measurement entry row_sim[MEASUREMENT] = measurement_sim - # append to dataframe - df = df.append(row_sim, ignore_index=True) + rows.append(row_sim) - return df + return pd.DataFrame(rows) def rdatas_to_simulation_df( From 7a5c1f4d9106821b128cae5f5302ecb0b3ffcc48 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 15:10:18 +0200 Subject: [PATCH 06/11] CI: Automatically detect available SBML semantic test cases (#1773) ... and do some cleanup along the way. Also fixes a `KeyError` when raising an exception using for reserved symbols as model entity identifiers. Closes #1771 --- .../test_sbml_semantic_test_suite.yml | 5 +- documentation/python_interface.rst | 3 +- python/amici/__init__.py | 7 +- python/amici/sbml_import.py | 18 +-- scripts/run-SBMLTestsuite.sh | 8 +- tests/conftest.py | 61 ++++++++-- tests/testSBMLSuite.py | 108 +++++++----------- 7 files changed, 117 insertions(+), 93 deletions(-) diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml index 04406d44d7..458f010fd0 100644 --- a/.github/workflows/test_sbml_semantic_test_suite.yml +++ b/.github/workflows/test_sbml_semantic_test_suite.yml @@ -13,6 +13,7 @@ on: - python/amici/import_utils.py - scripts/run-SBMLTestsuite.sh - tests/testSBMLSuite.py + - tests/conftest.py check_suite: types: [requested] workflow_dispatch: @@ -25,8 +26,8 @@ jobs: strategy: fail-fast: false matrix: - cases: ["1 - 250", "251 - 500", "501 - 750", "751 - 1000", - "1000-1250", "1251-1780"] + cases: ["1-250", "251-500", "501-750", "751-1000", + "1000-1250", "1251-"] python-version: [ 3.8 ] steps: diff --git a/documentation/python_interface.rst b/documentation/python_interface.rst index bea321b4ab..b4137c27e3 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 996 out of the 1780 (~56%) test cases** from +Python-AMICI currently **passes 1014 out of the 1821 (~56%) test cases** from the semantic `SBML Test Suite `_ (`current status `_). @@ -60,6 +60,7 @@ The following SBML test suite tags are currently supported * Concentration * ConstantSpecies * ConversionFactors +* DefaultValue * EventT0Firing * HasOnlySubstanceUnits * InitialValueReassigned diff --git a/python/amici/__init__.py b/python/amici/__init__.py index 53502f9750..97de3323b7 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -142,8 +142,10 @@ def __exit__(self, exc_type, exc_value, traceback): pass -def import_model_module(module_name: str, - module_path: Optional[str] = None) -> ModelModule: +def import_model_module( + module_name: str, + module_path: Optional[Union[Path, str]] = None +) -> ModelModule: """ Import Python module of an AMICI model @@ -154,6 +156,7 @@ def import_model_module(module_name: str, :return: The model module """ + module_path = str(module_path) # ensure we will find the newly created module importlib.invalidate_caches() diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index a5fd812a12..ead75d2f40 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -621,16 +621,20 @@ def add_local_symbol(self, key: str, value: sp.Expr): :param value: local symbol value """ - if key in itt.chain(self._local_symbols.keys(), - ['True', 'False', 'pi']): + if key in self._local_symbols.keys(): raise SBMLException( f'AMICI tried to add a local symbol {key} with value {value}, ' f'but {key} was already instantiated with ' - f'{self._local_symbols[key]}. This means either that there ' - f'are multiple SBML element with SId {key}, which is ' - f'invalid SBML, or that the employed SId {key} is a special ' - 'reserved symbol in AMICI. This can be fixed by renaming ' - f'the element with SId {key}.' + f'{self._local_symbols[key]}. This means that there ' + f'are multiple SBML elements with SId {key}, which is ' + f'invalid SBML. This can be fixed by renaming ' + f'the elements with SId {key}.' + ) + if key in {'True', 'False', 'true', 'false', 'pi'}: + raise SBMLException( + f'AMICI tried to add a local symbol {key} with value {value}, ' + f'but {key} is a reserved symbol in AMICI. This can be fixed ' + f'by renaming the element with SId {key}.' ) self._local_symbols[key] = value diff --git a/scripts/run-SBMLTestsuite.sh b/scripts/run-SBMLTestsuite.sh index 1767a16dfa..97c2033751 100755 --- a/scripts/run-SBMLTestsuite.sh +++ b/scripts/run-SBMLTestsuite.sh @@ -14,9 +14,11 @@ pip show pytest-xdist > /dev/null 2>&1 || pip install pytest-xdist pip install coverage pytest-cov if [[ -z "$*" ]]; then - args="1-1780" # run all tests + # run all tests + cases="" else - args="$@" # use user selection + # use user selection + cases="--cases=$@" fi # delete old result directory and recreate @@ -26,5 +28,5 @@ if [[ -d "${RESULT_DIR}" ]]; then fi mkdir "${RESULT_DIR}" -pytest ./tests/testSBMLSuite.py --cases="${args}" -rfsE -n auto \ +pytest ./tests/testSBMLSuite.py $cases -rfsE -n auto \ --cov=amici --cov-report=xml:"coverage_SBMLSuite.xml" --cov-append diff --git a/tests/conftest.py b/tests/conftest.py index 89eb1342da..a213d1901c 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -2,7 +2,8 @@ import re import sys -from typing import List +from pathlib import Path +from typing import List, Tuple, Set import pytest @@ -10,8 +11,17 @@ # stores passed SBML semantic test suite IDs passed_ids = [] +SBML_SEMANTIC_CASES_DIR = \ + Path(__file__).parent / 'sbml-test-suite' / 'cases' / 'semantic' -def parse_selection(selection_str: str) -> List[int]: + +@pytest.fixture +def sbml_semantic_cases_dir() -> Path: + """directory with sbml semantic test cases""" + return SBML_SEMANTIC_CASES_DIR + + +def parse_selection(selection_str: str, last: int) -> List[int]: """ Parse comma-separated list of integer ranges, return selected indices as integer list @@ -20,7 +30,7 @@ def parse_selection(selection_str: str) -> List[int]: """ indices = [] for group in selection_str.split(','): - if not re.match(r'^(?:-?\d+)|(?:\d+(?:-\d+))$', group): + if not re.match(r'^(?:-?\d+|\d+-\d*)$', group): print("Invalid selection", group) sys.exit() spl = group.split('-') @@ -28,11 +38,19 @@ def parse_selection(selection_str: str) -> List[int]: indices.append(int(spl[0])) elif len(spl) == 2: begin = int(spl[0]) if spl[0] else 0 - end = int(spl[1]) + end = int(spl[1]) if spl[1] else last indices.extend(range(begin, end + 1)) return indices +def get_all_semantic_case_ids(): + """Get iterator over test sorted IDs of all cases in the SBML semantic + suite""" + pattern = re.compile(r'\d{5}') + return sorted(str(x.name) for x in SBML_SEMANTIC_CASES_DIR.iterdir() + if pattern.match(x.name)) + + def pytest_addoption(parser): """Add pytest CLI options""" parser.addoption("--cases", help="Test cases to run") @@ -47,13 +65,11 @@ def pytest_generate_tests(metafunc): cases = metafunc.config.getoption("cases") if cases: # Run selected tests - test_numbers = set(parse_selection(cases)) + last_id = int(list(get_all_semantic_case_ids())[-1]) + test_numbers = sorted(set(parse_selection(cases, last_id))) else: # Run all tests - test_numbers = set(range(1, 1781)) - - # We skip this test due to NaNs in the Jacobian - test_numbers -= {1395} + test_numbers = get_all_semantic_case_ids() metafunc.parametrize("test_number", test_numbers) @@ -78,7 +94,6 @@ def write_passed_tags(passed_ids, out=sys.stdout): passed_component_tags = set() passed_test_tags = set() - from testSBMLSuite import get_tags_for_test for test_id in passed_ids: cur_component_tags, cur_test_tags = get_tags_for_test(test_id) passed_component_tags |= cur_component_tags @@ -100,3 +115,29 @@ def pytest_runtest_logreport(report: "TestReport") -> None: test_case_id = re.sub(r'^.*::test_sbml_testsuite_case\[(\d+)].*$', r'\1', report.nodeid) passed_ids.append(test_case_id) + + +def get_tags_for_test(test_id: str) -> Tuple[Set[str], Set[str]]: + """Get sbml test suite tags for the given test ID + + Returns: + Tuple of set of strings for componentTags and testTags + """ + current_test_path = SBML_SEMANTIC_CASES_DIR / test_id + info_file = current_test_path / f'{test_id}-model.m' + with open(info_file) as f: + component_tags = set() + test_tags = set() + for line in f: + if line.startswith('testTags:'): + test_tags = set( + re.split(r'[ ,:]', line[len('testTags:'):].strip())) + test_tags.discard('') + if line.startswith('componentTags:'): + component_tags = set( + re.split(r'[ ,:]', line[len('componentTags:'):].strip())) + component_tags.discard('') + if test_tags and component_tags: + return component_tags, test_tags + print(f"No componentTags or testTags found for test case {test_id}.") + return component_tags, test_tags diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index 187637a15f..3187720a8f 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -4,21 +4,17 @@ [https://github.com/sbmlteam/sbml-test-suite/releases] Usage: - python tests/testSBMLSuite.py SELECTION - SELECTION can be e.g.: `1`, `1,3`, or `-3,4,6-7` to select specific - test cases or 1-1780 to run all. - - pytest tests.testSBMLSuite -n CORES --cases SELECTION + pytest tests.testSBMLSuite -n CORES --cases=SELECTION CORES can be an integer or `auto` for all available cores. - SELECTION same as above. + SELECTION can be e.g.: `1`, `1,3`, `-3,4,6-7`, or `100-` to select + specific test cases. If `--cases` is omitted, all cases are run. """ import copy import os -import re import shutil import sys -from typing import Set, Tuple +from pathlib import Path import libsbml as sbml import numpy as np @@ -29,14 +25,10 @@ from amici.constants import SymbolId from amici.gradient_check import check_derivatives -# directory with sbml semantic test cases -TEST_PATH = os.path.join(os.path.dirname(__file__), 'sbml-test-suite', 'cases', - 'semantic') - @pytest.fixture(scope="session") -def result_path(): - return os.path.join(os.path.dirname(__file__), 'amici-semantic-results') +def result_path() -> Path: + return Path(__file__).parent / 'amici-semantic-results' @pytest.fixture(scope="function", autouse=True) @@ -52,10 +44,17 @@ def sbml_test_dir(): sys.path = old_path -def test_sbml_testsuite_case(test_number, result_path): +def test_sbml_testsuite_case( + test_number, + result_path, + sbml_semantic_cases_dir +): test_id = format_test_id(test_number) model_dir = None + if test_id == "01395": + pytest.skip("NaNs in the Jacobian") + # test cases for which sensitivities are to be checked # key: case ID; value: epsilon for finite differences sensitivity_check_cases = { @@ -64,19 +63,17 @@ def test_sbml_testsuite_case(test_number, result_path): } try: - current_test_path = os.path.join(TEST_PATH, test_id) + current_test_path = sbml_semantic_cases_dir / test_id # parse expected results - results_file = os.path.join(current_test_path, - f'{test_id}-results.csv') + results_file = current_test_path / f'{test_id}-results.csv' results = pd.read_csv(results_file, delimiter=',') results.rename(columns={c: c.replace(' ', '') for c in results.columns}, inplace=True) # setup model - model_dir = os.path.join(os.path.dirname(__file__), 'SBMLTestModels', - test_id) + model_dir = Path(__file__).parent / 'SBMLTestModels' / test_id model, solver, wrapper = compile_model( current_test_path, test_id, model_dir, generate_sensitivity_code=test_id in sensitivity_check_cases) @@ -151,7 +148,7 @@ def verify_results( requested_concentrations) # simulated may contain `object` dtype columns and `expected` may - # contain `np.int64` columns so we cast everything to `np.float64`. + # contain `np.int64` columns, so we cast everything to `np.float64`. for variable in variables: assert np.isclose( simulated[variable].astype(np.float64).values, @@ -170,6 +167,7 @@ def amounts_to_concentrations( ): """ Convert AMICI simulated amounts to concentrations + Convert from concentration to amount: C=n/V n=CV (multiply by V) @@ -221,7 +219,7 @@ def concentrations_to_amounts( def write_result_file( simulated: pd.DataFrame, test_id: str, - result_path: str + result_path: Path ): """ Create test result file for upload to @@ -231,7 +229,7 @@ def write_result_file( """ # TODO: only states are reported here, not compartments or parameters - filename = os.path.join(result_path, f'{test_id}.csv') + filename = result_path / f'{test_id}.csv' simulated.to_csv(filename, index=False) @@ -271,19 +269,19 @@ def apply_settings(settings, solver, model): return atol, rtol -def compile_model(path, test_id, model_dir, +def compile_model(sbml_dir: Path, test_id: str, model_dir: Path, generate_sensitivity_code: bool = False): """Import the given test model to AMICI""" - sbml_file = find_model_file(path, test_id) - - wrapper = amici.SbmlImporter(sbml_file) + model_dir.mkdir(parents=True, exist_ok=True) - if not os.path.exists(model_dir): - os.makedirs(model_dir) + sbml_file = find_model_file(sbml_dir, test_id) + sbml_importer = amici.SbmlImporter(sbml_file) model_name = f'SBMLTest{test_id}' - wrapper.sbml2amici(model_name, output_dir=model_dir, - generate_sensitivity_code=generate_sensitivity_code) + sbml_importer.sbml2amici( + model_name, output_dir=model_dir, + generate_sensitivity_code=generate_sensitivity_code + ) # settings model_module = amici.import_model_module(model_name, model_dir) @@ -291,28 +289,28 @@ def compile_model(path, test_id, model_dir, model = model_module.getModel() solver = model.getSolver() - return model, solver, wrapper + return model, solver, sbml_importer -def find_model_file(current_test_path: str, test_id: str): +def find_model_file(current_test_path: Path, test_id: str) -> Path: """Find model file for the given test (guess filename extension)""" - sbml_file = os.path.join(current_test_path, f'{test_id}-sbml-l3v2.xml') + sbml_file = current_test_path / f'{test_id}-sbml-l3v2.xml' - # fallback l3v1 - if not os.path.isfile(sbml_file): - sbml_file = os.path.join(current_test_path, f'{test_id}-sbml-l3v1.xml') + if not sbml_file.is_file(): + # fallback l3v1 + sbml_file = current_test_path / f'{test_id}-sbml-l3v1.xml' - # fallback l2v5 - if not os.path.isfile(sbml_file): - sbml_file = os.path.join(current_test_path, f'{test_id}-sbml-l2v5.xml') + if not sbml_file.is_file(): + # fallback l2v5 + sbml_file = current_test_path / f'{test_id}-sbml-l2v5.xml' return sbml_file -def read_settings_file(current_test_path: str, test_id: str): +def read_settings_file(current_test_path: Path, test_id: str): """Read settings for the given test""" - settings_file = os.path.join(current_test_path, f'{test_id}-settings.txt') + settings_file = current_test_path / f'{test_id}-settings.txt' settings = {} with open(settings_file) as f: for line in f: @@ -325,29 +323,3 @@ def read_settings_file(current_test_path: str, test_id: str): def format_test_id(test_id) -> str: """Format numeric to 0-padded string""" return f"{test_id:0>5}" - - -def get_tags_for_test(test_id) -> Tuple[Set[str], Set[str]]: - """Get sbml test suite tags for the given test ID - - Returns: - Tuple of set of strings for componentTags and testTags - """ - current_test_path = os.path.join(TEST_PATH, test_id) - info_file = os.path.join(current_test_path, f'{test_id}-model.m') - with open(info_file) as f: - component_tags = set() - test_tags = set() - for line in f: - if line.startswith('testTags:'): - test_tags = set( - re.split(r'[ ,:]', line[len('testTags:'):].strip())) - test_tags.discard('') - if line.startswith('componentTags:'): - component_tags = set( - re.split(r'[ ,:]', line[len('componentTags:'):].strip())) - component_tags.discard('') - if test_tags and component_tags: - return component_tags, test_tags - print(f"No componentTags or testTags found for test case {test_id}.") - return component_tags, test_tags From a2e543796d9f478673afd8cbd60fd2286912dc12 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 16:32:36 +0200 Subject: [PATCH 07/11] Remove long deprecated sbml2amici arguments modelName and constantParameters (#1774) * Remove long deprecated sbml2amici arguments modelName and constantParameters (and update notebooks) * Remove duplicated test `python/tests/test_misc.py::test_sbml2amici_no_observables` also present in `python/tests/test_sbml_import.py` --- python/amici/sbml_import.py | 73 ++++++------------- .../ExampleEquilibrationLogic.ipynb | 6 +- python/tests/test_misc.py | 26 ------- 3 files changed, 24 insertions(+), 81 deletions(-) diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index ead75d2f40..616f82e102 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -205,24 +205,25 @@ def _reset_symbols(self) -> None: self.symbols = copy.deepcopy(default_symbols) self._local_symbols = {} - def sbml2amici(self, - model_name: str = None, - output_dir: Union[str, Path] = None, - observables: Dict[str, Dict[str, str]] = None, - constant_parameters: Iterable[str] = None, - sigmas: Dict[str, Union[str, float]] = None, - noise_distributions: Dict[str, Union[str, Callable]] = None, - verbose: Union[int, bool] = logging.ERROR, - assume_pow_positivity: bool = False, - compiler: str = None, - allow_reinit_fixpar_initcond: bool = True, - compile: bool = True, - compute_conservation_laws: bool = True, - simplify: Callable = lambda x: sp.powsimp(x, deep=True), - cache_simplify: bool = False, - log_as_log10: bool = True, - generate_sensitivity_code: bool = True, - **kwargs) -> None: + def sbml2amici( + self, + model_name: str, + output_dir: Union[str, Path] = None, + observables: Dict[str, Dict[str, str]] = None, + constant_parameters: Iterable[str] = None, + sigmas: Dict[str, Union[str, float]] = None, + noise_distributions: Dict[str, Union[str, Callable]] = None, + verbose: Union[int, bool] = logging.ERROR, + assume_pow_positivity: bool = False, + compiler: str = None, + allow_reinit_fixpar_initcond: bool = True, + compile: bool = True, + compute_conservation_laws: bool = True, + simplify: Callable = lambda x: sp.powsimp(x, deep=True), + cache_simplify: bool = False, + log_as_log10: bool = True, + generate_sensitivity_code: bool = True, + ) -> None: """ Generate and compile AMICI C++ files for the model provided to the constructor. @@ -315,22 +316,8 @@ def sbml2amici(self, """ set_log_level(logger, verbose) - if 'constantParameters' in kwargs: - logger.warning('Use of `constantParameters` as argument name ' - 'is deprecated and will be removed in a future ' - 'version. Please use `constant_parameters` as ' - 'argument name.') - - if constant_parameters is not None: - raise ValueError('Cannot specify constant parameters using ' - 'both `constantParameters` and ' - '`constant_parameters` as argument names.') - - constant_parameters = kwargs.pop('constantParameters', []) - - elif constant_parameters is None: - constant_parameters = [] - constant_parameters = list(constant_parameters) + constant_parameters = list(constant_parameters) \ + if constant_parameters else [] if sigmas is None: sigmas = {} @@ -338,24 +325,6 @@ def sbml2amici(self, if noise_distributions is None: noise_distributions = {} - if model_name is None: - model_name = kwargs.pop('modelName', None) - if model_name is None: - raise ValueError('Missing argument: `model_name`') - else: - logger.warning('Use of `modelName` as argument name is ' - 'deprecated and will be removed in a future' - ' version. Please use `model_name` as ' - 'argument name.') - else: - if 'modelName' in kwargs: - raise ValueError('Cannot specify model name using both ' - '`modelName` and `model_name` as argument ' - 'names.') - - if len(kwargs): - raise ValueError(f'Unknown arguments {kwargs.keys()}.') - self._reset_symbols() self.sbml_parser_settings.setParseLog( sbml.L3P_PARSE_LOG_AS_LOG10 if log_as_log10 else diff --git a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb index 1e24e62a56..40ecb9f7aa 100644 --- a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb +++ b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb @@ -153,12 +153,12 @@ "sbml_importer.sbml2amici(model_reduced_name,\n", " model_reduced_output_dir,\n", " observables=observables,\n", - " constantParameters=constantParameters,\n", + " constant_parameters=constantParameters,\n", " sigmas=sigmas)\n", "sbml_importer.sbml2amici(model_name,\n", " model_output_dir,\n", " observables=observables,\n", - " constantParameters=constantParameters,\n", + " constant_parameters=constantParameters,\n", " sigmas=sigmas,\n", " compute_conservation_laws=False)" ] @@ -1204,4 +1204,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/python/tests/test_misc.py b/python/tests/test_misc.py index 9bb9ebede3..c644aacade 100644 --- a/python/tests/test_misc.py +++ b/python/tests/test_misc.py @@ -28,32 +28,6 @@ def test_parameter_scaling_from_int_vector(): assert scale_vector[2] == amici.ParameterScaling.none -def test_sbml2amici_no_observables(): - """Test model generation works for model without observables""" - - # test model - document = libsbml.SBMLDocument(3, 1) - model = document.createModel() - model.setTimeUnits("second") - model.setExtentUnits("mole") - model.setSubstanceUnits('mole') - c1 = model.createCompartment() - c1.setId('C1') - model.addCompartment(c1) - s1 = model.createSpecies() - s1.setId('S1') - s1.setCompartment('C1') - model.addSpecies(s1) - - sbml_importer = amici.sbml_import.SbmlImporter(sbml_source=model, - from_file=False) - tmpdir = TemporaryDirectory() - sbml_importer.sbml2amici(modelName="test", - output_dir=tmpdir.name, - observables=None, - compute_conservation_laws=False) - - def test_hill_function_dwdx(): """Kinetic laws with Hill functions, may lead to NaNs in the Jacobian if involved states are zero if not properly arranged symbolically. From e9fa6461f94236ecb152b94904ea996f9d4a7bb0 Mon Sep 17 00:00:00 2001 From: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> Date: Fri, 8 Apr 2022 21:04:44 +0200 Subject: [PATCH 08/11] Specify initial timepoint with ExpData (#1776) --- include/amici/edata.h | 1 + src/edata.cpp | 3 +++ 2 files changed, 4 insertions(+) diff --git a/include/amici/edata.h b/include/amici/edata.h index 8c1e821aea..4c12596a08 100644 --- a/include/amici/edata.h +++ b/include/amici/edata.h @@ -519,6 +519,7 @@ class ConditionContext : public ContextManager { std::vector original_sx0_; std::vector original_parameters_; std::vector original_fixed_parameters_; + realtype original_tstart_; std::vector original_timepoints_; std::vector original_parameter_list_; std::vector original_scaling_; diff --git a/src/edata.cpp b/src/edata.cpp index 1842cf669f..cc37a36455 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -336,6 +336,7 @@ ConditionContext::ConditionContext(Model *model, const ExpData *edata, : model_(model), original_parameters_(model->getParameters()), original_fixed_parameters_(model->getFixedParameters()), + original_tstart_(model->t0()), original_timepoints_(model->getTimepoints()), original_parameter_list_(model->getParameterList()), original_scaling_(model->getParameterScale()), @@ -454,6 +455,7 @@ void ConditionContext::applyCondition(const ExpData *edata, break; } + model_->setT0(edata->tstart_); if(edata->nt()) { // fixed parameter in model are superseded by those provided in edata model_->setTimepoints(edata->getTimepoints()); @@ -475,6 +477,7 @@ void ConditionContext::restore() model_->setParameters(original_parameters_); model_->setFixedParameters(original_fixed_parameters_); + model_->setT0(original_tstart_); model_->setTimepoints(original_timepoints_); model_->setReinitializeFixedParameterInitialStates( original_reinitialize_fixed_parameter_initial_states_); From 8d19766d672f69c49254a1030135dbe17875ff9a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 21:56:19 +0200 Subject: [PATCH 09/11] Fix Rule-handling in PEtab import (#1753) RateRules-targets in the condition table were previously ignored. Closes #1750 Also fixes a bug where InitialAssignments were not honored for elements that occurred as column in the condition table and were set to `NaN` (i.e. "use model value"). --- .github/workflows/test_petab_test_suite.yml | 5 +- python/amici/parameter_mapping.py | 20 +++--- python/amici/petab_import.py | 28 ++++++-- python/amici/petab_objective.py | 76 +++++++++++++-------- tests/petab_test_suite/test_petab_suite.py | 33 ++++++--- 5 files changed, 108 insertions(+), 54 deletions(-) diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml index b9a6d1bd6e..f089b544b9 100644 --- a/.github/workflows/test_petab_test_suite.yml +++ b/.github/workflows/test_petab_test_suite.yml @@ -63,9 +63,10 @@ jobs: # retrieve test models - name: Download and install PEtab test suite run: | - git clone --depth 1 --branch develop https://github.com/PEtab-dev/petab_test_suite \ + git clone --depth 1 --branch develop \ + https://github.com/PEtab-dev/petab_test_suite \ && source ./build/venv/bin/activate \ - && cd petab_test_suite && pip3 install -e . && cd .. + && cd petab_test_suite && pip3 install -e . - name: Run PEtab-related unit tests run: | diff --git a/python/amici/parameter_mapping.py b/python/amici/parameter_mapping.py index 23eb190dc1..930cd02a81 100644 --- a/python/amici/parameter_mapping.py +++ b/python/amici/parameter_mapping.py @@ -191,14 +191,18 @@ def fill_in_parameters_for_condition( # Parameter mapping may contain parameter_ids as values, these *must* # be replaced - def _get_par(model_par, value): + def _get_par(model_par, value, mapping): """Replace parameter IDs in mapping dicts by values from problem_parameters where necessary""" if isinstance(value, str): - # estimated parameter - # (condition table overrides must have been handled already, - # e.g. by the PEtab parameter mapping) - return problem_parameters[value] + try: + # estimated parameter + return problem_parameters[value] + except KeyError: + # condition table overrides must have been handled already, + # e.g. by the PEtab parameter mapping, but parameters from + # InitialAssignments may still be present. + return _get_par(value, mapping[value], mapping) if model_par in problem_parameters: # user-provided return problem_parameters[model_par] @@ -208,11 +212,11 @@ def _get_par(model_par, value): # constant value return value - map_preeq_fix = {key: _get_par(key, val) + map_preeq_fix = {key: _get_par(key, val, map_preeq_fix) for key, val in map_preeq_fix.items()} - map_sim_fix = {key: _get_par(key, val) + map_sim_fix = {key: _get_par(key, val, map_sim_fix) for key, val in map_sim_fix.items()} - map_sim_var = {key: _get_par(key, val) + map_sim_var = {key: _get_par(key, val, map_sim_var) for key, val in map_sim_var.items()} # If necessary, (un)scale parameters diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index c227595a11..f1183b451f 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -86,7 +86,8 @@ def get_fixed_parameters( # remove overridden parameters (`object`-type columns) fixed_parameters = [p for p in fixed_parameters if condition_df[p].dtype != 'O' - and sbml_model.getParameter(p) is not None] + and sbml_model.getParameter(p) is not None + and sbml_model.getRuleByVariable(p) is None] # must be unique if len(fixed_parameters) != len(set(fixed_parameters)): raise AssertionError( @@ -512,12 +513,11 @@ def import_model_sbml( # create a new parameter initial_${startOrCompartmentID}. # feels dirty and should be changed (see also #924) # + initial_states = [col for col in condition_df - if sbml_model.getSpecies(col) is not None] - initial_sizes = [col for col in condition_df - if sbml_model.getCompartment(col) is not None] + if element_is_state(sbml_model, col)] fixed_parameters = [] - if len(initial_states) or len(initial_sizes): + if initial_states: # add preequilibration indicator variable # NOTE: would only be required if we actually have preequilibration # adding it anyways. can be optimized-out later @@ -534,8 +534,8 @@ def import_model_sbml( logger.debug("Adding preequilibration indicator " f"constant {PREEQ_INDICATOR_ID}") logger.debug("Adding initial assignments for " - f"{initial_sizes + initial_states}") - for assignee_id in chain(initial_sizes, initial_states): + f"{initial_states}") + for assignee_id in initial_states: init_par_id_preeq = f"initial_{assignee_id}_preeq" init_par_id_sim = f"initial_{assignee_id}_sim" for init_par_id in [init_par_id_preeq, init_par_id_sim]: @@ -689,6 +689,20 @@ def show_model_info(sbml_model: 'libsbml.Model'): logger.info(f'Reactions: {len(sbml_model.getListOfReactions())}') +def element_is_state(sbml_model: libsbml.Model, sbml_id: str) -> bool: + """Does the element with ID `sbml_id` correspond to a state variable? + """ + if sbml_model.getCompartment(sbml_id) is not None: + return True + if sbml_model.getSpecies(sbml_id) is not None: + return True + if (rule := sbml_model.getRuleByVariable(sbml_id)) is not None \ + and rule.getTypeCode() == libsbml.SBML_RATE_RULE: + return True + + return False + + def _parse_cli_args(): """ Parse command line arguments diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py index 7ed62bb509..1e6dfb0e7f 100644 --- a/python/amici/petab_objective.py +++ b/python/amici/petab_objective.py @@ -22,7 +22,7 @@ from . import AmiciModel, AmiciExpData from .logging import get_logger, log_execution_time -from .petab_import import PREEQ_INDICATOR_ID +from .petab_import import PREEQ_INDICATOR_ID, element_is_state from .parameter_mapping import ( fill_in_parameters, ParameterMappingForCondition, ParameterMapping) @@ -337,11 +337,11 @@ def create_parameter_mapping_for_condition( # ExpData.x0, but in the case of preequilibration this would not allow for # resetting initial states. - species_in_condition_table = [ + states_in_condition_table = [ col for col in petab_problem.condition_df - if petab_problem.sbml_model.getSpecies(col) is not None] - - if species_in_condition_table: + if element_is_state(petab_problem.sbml_model, col) + ] + if states_in_condition_table: # set indicator fixed parameter for preeq # (we expect here, that this parameter was added during import and # that it was not added by the user with a different meaning...) @@ -352,14 +352,33 @@ def create_parameter_mapping_for_condition( condition_map_sim[PREEQ_INDICATOR_ID] = 0.0 condition_scale_map_sim[PREEQ_INDICATOR_ID] = LIN - def _set_initial_concentration(condition_id, species_id, init_par_id, - par_map, scale_map): + def _set_initial_state(condition_id, element_id, init_par_id, + par_map, scale_map): value = petab.to_float_if_float( - petab_problem.condition_df.loc[condition_id, species_id]) + petab_problem.condition_df.loc[condition_id, element_id]) if pd.isna(value): - value = get_species_initial( - petab_problem.sbml_model.getSpecies(species_id) - ) + element = petab_problem.sbml_model.getElementBySId(element_id) + type_code = element.getTypeCode() + initial_assignment = petab_problem.sbml_model\ + .getInitialAssignmentBySymbol(element_id) + if initial_assignment: + initial_assignment = sp.sympify( + libsbml.formulaToL3String(initial_assignment.getMath()) + ) + if type_code == libsbml.SBML_SPECIES: + value = get_species_initial(element) \ + if initial_assignment is None else initial_assignment + elif type_code == libsbml.SBML_PARAMETER: + value = element.getValue()\ + if initial_assignment is None else initial_assignment + elif type_code == libsbml.SBML_COMPARTMENT: + value = element.getSize()\ + if initial_assignment is None else initial_assignment + else: + raise NotImplementedError( + f"Don't know what how to handle {element_id} in " + "condition table.") + try: value = float(value) except (ValueError, TypeError): @@ -368,11 +387,11 @@ def _set_initial_concentration(condition_id, species_id, init_par_id, value = sp.nsimplify(value) else: raise NotImplementedError( - "Cannot handle non-trivial expressions for " - f"species initial for {species_id}: {value}") + "Cannot handle non-trivial initial state " + f"expression for {element_id}: {value}") # this should be a parameter ID value = str(value) - logger.debug(f'The species {species_id} has no initial value ' + logger.debug(f'The species {element_id} has no initial value ' f'defined for the condition {condition_id} in ' 'the PEtab conditions table. The initial value is ' f'now set to {value}, which is the initial value ' @@ -383,16 +402,17 @@ def _set_initial_concentration(condition_id, species_id, init_par_id, scale_map[init_par_id] = petab.LIN else: # parametric initial state - scale_map[init_par_id] = petab_problem.parameter_df.loc[ - value, PARAMETER_SCALE] + scale_map[init_par_id] = \ + petab_problem.parameter_df[PARAMETER_SCALE]\ + .get(value, petab.LIN) - for species_id in species_in_condition_table: + for element_id in states_in_condition_table: # for preequilibration - init_par_id = f'initial_{species_id}_preeq' + init_par_id = f'initial_{element_id}_preeq' if condition.get(PREEQUILIBRATION_CONDITION_ID): condition_id = condition[PREEQUILIBRATION_CONDITION_ID] - _set_initial_concentration( - condition_id, species_id, init_par_id, condition_map_preeq, + _set_initial_state( + condition_id, element_id, init_par_id, condition_map_preeq, condition_scale_map_preeq) else: # need to set dummy value for preeq parameter anyways, as it @@ -403,9 +423,9 @@ def _set_initial_concentration(condition_id, species_id, init_par_id, # for simulation condition_id = condition[SIMULATION_CONDITION_ID] - init_par_id = f'initial_{species_id}_sim' - _set_initial_concentration( - condition_id, species_id, init_par_id, condition_map_sim, + init_par_id = f'initial_{element_id}_sim' + _set_initial_state( + condition_id, element_id, init_par_id, condition_map_sim, condition_scale_map_sim) ########################################################################## @@ -547,22 +567,22 @@ def create_edata_for_condition( edata.id += "+" + condition.get(PREEQUILIBRATION_CONDITION_ID) ########################################################################## # enable initial parameters reinitialization - species_in_condition_table = [ + states_in_condition_table = [ col for col in petab_problem.condition_df if not pd.isna(petab_problem.condition_df.loc[ condition[SIMULATION_CONDITION_ID], col]) - and petab_problem.sbml_model.getSpecies(col) is not None + and element_is_state(petab_problem.sbml_model, col) ] if condition.get(PREEQUILIBRATION_CONDITION_ID) \ - and species_in_condition_table: + and states_in_condition_table: state_ids = amici_model.getStateIds() state_idx_reinitalization = [state_ids.index(s) - for s in species_in_condition_table] + for s in states_in_condition_table] edata.reinitialization_state_idxs_sim = state_idx_reinitalization logger.debug("Enabling state reinitialization for condition " f"{condition.get(PREEQUILIBRATION_CONDITION_ID, '')} - " f"{condition.get(SIMULATION_CONDITION_ID)} " - f"{species_in_condition_table}") + f"{states_in_condition_table}") ########################################################################## # timepoints diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 6c31ebbfa5..14f319d77b 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -1,15 +1,16 @@ #!/usr/bin/env python3 - """Run PEtab test suite (https://github.com/PEtab-dev/petab_test_suite)""" import logging import sys -import amici +import pandas as pd import petab import petabtests import pytest from _pytest.outcomes import Skipped + +import amici from amici import SteadyStateSensitivityMode from amici.gradient_check import check_derivatives as amici_check_derivatives from amici.logging import get_logger, set_log_level @@ -105,16 +106,26 @@ def _test_case(case, model_type): simulations_match = petabtests.evaluate_simulations( [simulation_df], gt_simulation_dfs, tol_simulations) + logger.log(logging.DEBUG if simulations_match else logging.ERROR, + f"Simulations: match = {simulations_match}") + if not simulations_match: + with pd.option_context('display.max_rows', None, + 'display.max_columns', None, + 'display.width', 200): + logger.log(logging.DEBUG, f"x_ss: {model.getStateIds()} " + f"{[rdata.x_ss for rdata in rdatas]}") + logger.log(logging.ERROR, + f"Expected simulations:\n{gt_simulation_dfs}") + logger.log(logging.ERROR, + f"Actual simulations:\n{simulation_df}") logger.log(logging.DEBUG if chi2s_match else logging.ERROR, f"CHI2: simulated: {chi2}, expected: {gt_chi2}," f" match = {chi2s_match}") logger.log(logging.DEBUG if simulations_match else logging.ERROR, f"LLH: simulated: {llh}, expected: {gt_llh}, " f"match = {llhs_match}") - logger.log(logging.DEBUG if simulations_match else logging.ERROR, - f"Simulations: match = {simulations_match}") - check_derivatives(problem, model) + check_derivatives(problem, model, solver) if not all([llhs_match, simulations_match]) or not chi2s_match: logger.error(f"Case {case} failed.") @@ -124,19 +135,23 @@ def _test_case(case, model_type): logger.info(f"Case {case} passed.") -def check_derivatives(problem: petab.Problem, model: amici.Model) -> None: +def check_derivatives( + problem: petab.Problem, + model: amici.Model, + solver: amici.Solver +) -> None: """Check derivatives using finite differences for all experimental conditions Arguments: problem: PEtab problem model: AMICI model matching ``problem`` + solver: AMICI solver """ problem_parameters = {t.Index: getattr(t, petab.NOMINAL_VALUE) for t in problem.parameter_df.itertuples()} - solver = model.getSolver() - solver.setSensitivityMethod(amici.SensitivityMethod_forward) - solver.setSensitivityOrder(amici.SensitivityOrder_first) + solver.setSensitivityMethod(amici.SensitivityMethod.forward) + solver.setSensitivityOrder(amici.SensitivityOrder.first) # Required for case 9 to not fail in # amici::NewtonSolver::computeNewtonSensis model.setSteadyStateSensitivityMode( From 342e74377a338d252a019b0c4459a6e8d058547f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Fri, 8 Apr 2022 16:35:37 -0400 Subject: [PATCH 10/11] Improve efficiency of newton step conv check (#1775) * improve efficiency of newton step conv check * Update src/steadystateproblem.cpp Co-authored-by: Daniel Weindl * fix test by increasing tolerance Co-authored-by: Daniel Weindl --- include/amici/solver.h | 3 ++- python/tests/test_preequilibration.py | 1 + src/steadystateproblem.cpp | 10 ++++++++-- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/include/amici/solver.h b/include/amici/solver.h index 42fa67f5c2..bf705338af 100644 --- a/include/amici/solver.h +++ b/include/amici/solver.h @@ -897,7 +897,8 @@ class Solver { } /** - * @brief Returns how convergence checks for steadystate computation are performed. + * @brief Returns how convergence checks for steadystate computation are performed. If activated, + * convergence checks are limited to every 25 steps in the simulation solver to limit performance impact. * @return boolean flag indicating newton step (true) or the right hand side (false) */ bool getNewtonStepSteadyStateCheck() const { diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index d87f10b8e1..0c390e0456 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -338,6 +338,7 @@ def test_newton_solver_equilibration(preeq_fixture): amici.SteadyStateSensitivityMode.newtonOnly] solver.setNewtonStepSteadyStateCheck(True) + solver.setRelativeToleranceSteadyState(1e-12) for equil_meth in settings: # set sensi method diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index 336d98e1ea..df5104dd9e 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -630,7 +630,12 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, /* If run after Newton's method checks again if it converged */ wrms_ = getWrms(model, sensitivityFlag); int sim_steps = 0; - + + int convergence_check_frequency = 1; + + if (newton_step_conv_) + convergence_check_frequency = 25; + while (true) { /* One step of ODE integration reason for tout specification: @@ -654,7 +659,8 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, /* getWrms needs to be called before getWrmsFSA such that the linear system is prepared for newton type convergence check */ if (wrms_ < conv_thresh && check_sensi_conv_ && - sensitivityFlag == SensitivityMethod::forward) { + sensitivityFlag == SensitivityMethod::forward && + sim_steps % convergence_check_frequency == 0) { updateSensiSimulation(solver); wrms_ = getWrmsFSA(model); } From 244ebbb4fcdea16438fb6b1ddd8af1f8aed5bc3b Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Fri, 8 Apr 2022 22:57:40 +0200 Subject: [PATCH 11/11] Bump version number, update changelog --- CHANGELOG.md | 36 +++++++++++++++++++++++++++++++++++- version.txt | 2 +- 2 files changed, 36 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9823652959..4bde64a0b1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,41 @@ ## v0.X Series +### v0.11.28 (2022-04-08) + +New features: +* Added `Solver.setSteadyStateToleranceFactor` and + `Solver.setSteadyStateSensiToleranceFactor` to specify a steady state + tolerance factor by @dilpath in https://github.com/AMICI-dev/AMICI/pull/1758 + + **NOTE: This also relaxed the default steady state (sensitivity)** + **tolerances by a factor of 100.** +* Added support for `pathlib.Path` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1769 +* Allow specifying initial timepoint with `ExpData` by @dilpath in + https://github.com/AMICI-dev/AMICI/pull/1776 + +Performance: +* Speedup for models with conservation laws by @FFroehlich in + https://github.com/AMICI-dev/AMICI/pull/1765 +* Improved efficiency of newton step convergence check by @FFroehlich in + https://github.com/AMICI-dev/AMICI/pull/1775 + +Fixes: +* Fixed deprecation warning for pandas.DataFrame.append in + `rdatas_to_measurement_df` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1770 +* Fixed Rule-target handling in PEtab import by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1753 + +Removed functionality: +* Removed long deprecated `sbml2amici` arguments `modelName` + and `constantParameters` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1774 + +**Full Changelog**: +https://github.com/AMICI-dev/AMICI/compare/v0.11.27...v0.11.28 + ### v0.11.27 (2022-04-04) New features: @@ -24,7 +59,6 @@ Performance: * Optional parallel computation of derivatives during model import by @dweindl in https://github.com/AMICI-dev/AMICI/pull/1740 * Sparsify jacobian by @FFroehlich in https://github.com/AMICI-dev/AMICI/pull/1766 -* Speedup for models with conservation laws by @FFroehlich in https://github.com/AMICI-dev/AMICI/pull/1765 * Speedup conservation law computation by @FFroehlich in https://github.com/AMICI-dev/AMICI/pull/1754 * Exploit stoichiometric matrix in pysb import by @FFroehlich in diff --git a/version.txt b/version.txt index 7d7d637254..2b71711fca 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.11.27 +0.11.28