From 7a9a16c37e7d0fcd888484b5a35a6227eead7693 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 May 2023 16:34:34 +0200 Subject: [PATCH 01/10] Update Doxygen (#2099) --- documentation/conf.py | 2 +- scripts/downloadAndBuildDoxygen.sh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/documentation/conf.py b/documentation/conf.py index 50834ae53e..96209e4c31 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -86,7 +86,7 @@ def install_mtocpp(): def install_doxygen(): """Get a more recent doxygen""" - version = "1.9.6" + version = "1.9.7" doxygen_exe = os.path.join( amici_dir, "ThirdParty", f"doxygen-{version}", "bin", "doxygen" ) diff --git a/scripts/downloadAndBuildDoxygen.sh b/scripts/downloadAndBuildDoxygen.sh index 503934cb31..c51c05c599 100755 --- a/scripts/downloadAndBuildDoxygen.sh +++ b/scripts/downloadAndBuildDoxygen.sh @@ -9,7 +9,7 @@ DOXYGEN_DIR="${AMICI_PATH}"/ThirdParty/doxygen cd "${AMICI_PATH}"/ThirdParty if [[ ! -d ${DOXYGEN_DIR} ]]; then # git clone --depth 1 https://github.com/doxygen/doxygen.git "${DOXYGEN_DIR}" - git clone --single-branch --branch Release_1_9_6 --depth 1 https://github.com/doxygen/doxygen.git "${DOXYGEN_DIR}" + git clone --single-branch --branch Release_1_9_7 --depth 1 https://github.com/doxygen/doxygen.git "${DOXYGEN_DIR}" fi cd "${DOXYGEN_DIR}" From 8d286c72c47630f7193f923eba354c4f18bf17ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Mon, 22 May 2023 15:37:52 +0100 Subject: [PATCH 02/10] update swig ignores (#2098) * add more swig ignores * Update amici.i * Update amici.i * Update amici.i * -s * 2304 * keep asserts * Update amici.i * Revert "2304" This reverts commit 29483825567ffa7db82b038c53b8c0beacc73682. * Revert "keep asserts" This reverts commit 3aa213b63c405b606939995b96f45b0793c97e1d. * Revert "-s" This reverts commit 731cfc2d353ada564c445c4698dbbffb75bbd822. --------- Co-authored-by: Daniel Weindl --- swig/amici.i | 3 +++ swig/model.i | 10 ++++++++++ swig/solver.i | 15 +++++++++++++++ 3 files changed, 28 insertions(+) diff --git a/swig/amici.i b/swig/amici.i index 6c796391e3..0015ca1bd0 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -154,12 +154,15 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData) %naturalvar amici::SimulationParameters::reinitialization_state_idxs_sim; %naturalvar amici::SimulationParameters::reinitialization_state_idxs_presim; +// DO NOT IGNORE amici::SimulationParameters, amici::ModelDimensions, amici::CpuTimer %ignore amici::ModelContext; %ignore amici::ContextManager; %ignore amici::ModelState; %ignore amici::ModelStateDerived; %ignore amici::unravel_index; %ignore amici::backtraceString; +%ignore amici::Logger; +%ignore amici::SimulationState; // Include before any other header which uses enums defined there %include "amici/defines.h" diff --git a/swig/model.i b/swig/model.i index d2c9f2eabd..3063590c21 100644 --- a/swig/model.i +++ b/swig/model.i @@ -84,6 +84,16 @@ using namespace amici; %ignore getObservableSigma; %ignore getObservableSigmaSensitivity; %ignore getUnobservedEventSensitivity; +%ignore fdsigmaydy; +%ignore fdspline_slopesdp; +%ignore fdspline_valuesdp; +%ignore fdtotal_cldp; +%ignore fdtotal_cldx_rdata; +%ignore fdx_rdatadp; +%ignore fdx_rdatadtcl; +%ignore fdx_rdatadx_solver; +%ignore fdsigmaydy; + diff --git a/swig/solver.i b/swig/solver.i index 319e654530..992842c409 100644 --- a/swig/solver.i +++ b/swig/solver.i @@ -39,6 +39,21 @@ using namespace amici; %ignore turnOffRootFinding; %ignore getRootInfo; %ignore updateAndReinitStatesAndSensitivities; +%ignore getCpuTime; +%ignore getCpuTimeB; +%ignore getLastOrder; +%ignore getNumErrTestFails; +%ignore getNumErrTestFailsB; +%ignore getNumNonlinSolvConvFails; +%ignore getNumNonlinSolvConvFailsB; +%ignore getNumRhsEvals; +%ignore getNumRhsEvalsB; +%ignore getNumSteps; +%ignore getNumStepsB; +%ignore gett; +%ignore startTimer; +%ignore switchForwardSensisOff; +%ignore timeExceeded; // Solver.__repr__ %pythoncode %{ From d0e8853e13e1a0b1b698cd2f190859b1c25f5c37 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 May 2023 19:57:10 +0200 Subject: [PATCH 03/10] PEtab-multilanguage proof of concept for PySB-PEtab (#1800) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Test new PEtab model abstraction * Update .github/workflows/test_petab_test_suite.yml * Run v2.0.0 test suite * pysb models * .. * .. * .. * .. * libpetab-python.git@model_pysb * generalize fixed parameters * .. * .. * .. * .. * cleanup * cleanup * add support for species reinitialization in pysb * fixup preequilibration * fixup preequilibration pysb * fixup pysb import * prevent warnings? * move resolve mapping to petab lib, adapt to updated api * fix error message * fixup condition checks * suppress model warnings * fixup export issues * cleanup * fixup * fixup compartment * simplify * Fix NameError: name 'sbml_model' is not defined * ? * make pattern matching exact so we only assign single species * Update petab_import_pysb.py * .. * Revert ".." This reverts commit 02a4f2402f0ef2f975a404720b0b0440fb19b085. * .. * deps * .. * Update test_benchmark_collection_models.yml * branch * Fix matplotlib conflict > The conflict is caused by: > The user requested matplotlib==3.5.3 > petab[vis] 0.2.0 depends on matplotlib>=3.6.0; extra == "vis" * merge? * e * .. * - libpetab-python.git@model_pysb * doc,cleanup * fix model name / dir * fixup, cleanup * cleanup; fixup merge * .. * pypi * pysb * .. * deps * pypi-petab * reqs --------- Co-authored-by: Fabian Fröhlich --- .github/workflows/test_petab_test_suite.yml | 7 +- python/sdist/amici/parameter_mapping.py | 3 + python/sdist/amici/petab_import.py | 134 +++--- python/sdist/amici/petab_import_pysb.py | 471 +++++++------------- python/sdist/amici/petab_objective.py | 260 ++++++----- python/sdist/amici/petab_util.py | 102 +++++ python/sdist/setup.cfg | 2 +- scripts/installAmiciSource.sh | 5 +- tests/petab_test_suite/conftest.py | 26 +- tests/petab_test_suite/test_petab_suite.py | 6 +- 10 files changed, 503 insertions(+), 513 deletions(-) create mode 100644 python/sdist/amici/petab_util.py diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml index 28a76678a7..70acd254df 100644 --- a/.github/workflows/test_petab_test_suite.yml +++ b/.github/workflows/test_petab_test_suite.yml @@ -43,8 +43,6 @@ jobs: libatlas-base-dev \ python3-venv - - run: pip3 install pysb petab - - name: Build BNGL run: | scripts/buildBNGL.sh @@ -59,6 +57,11 @@ jobs: run: | scripts/installAmiciSource.sh + - name: Install petab + run: | + source ./build/venv/bin/activate \ + && pip3 install wheel pytest shyaml pytest-cov pysb + # retrieve test models - name: Download and install PEtab test suite run: | diff --git a/python/sdist/amici/parameter_mapping.py b/python/sdist/amici/parameter_mapping.py index 716b61adc7..9f4d3b24dd 100644 --- a/python/sdist/amici/parameter_mapping.py +++ b/python/sdist/amici/parameter_mapping.py @@ -248,6 +248,9 @@ def _get_par(model_par, value, mapping): # condition table overrides must have been handled already, # e.g. by the PEtab parameter mapping, but parameters from # InitialAssignments may still be present. + if mapping[value] == model_par: + # prevent infinite recursion + raise return _get_par(value, mapping[value], mapping) if model_par in problem_parameters: # user-provided diff --git a/python/sdist/amici/petab_import.py b/python/sdist/amici/petab_import.py index a7a28dd328..909bf250ae 100644 --- a/python/sdist/amici/petab_import.py +++ b/python/sdist/amici/petab_import.py @@ -7,6 +7,7 @@ import argparse import importlib import logging +import math import os import re import shutil @@ -24,22 +25,20 @@ from _collections import OrderedDict from amici.logging import get_logger, log_execution_time, set_log_level from petab.C import * +from petab.models import MODEL_TYPE_PYSB, MODEL_TYPE_SBML from petab.parameters import get_valid_parameters_for_parameter_table from sympy.abc import _clash +from .petab_util import PREEQ_INDICATOR_ID, get_states_in_condition_table + try: - from amici.petab_import_pysb import PysbPetabProblem, import_model_pysb + from amici.petab_import_pysb import import_model_pysb except ModuleNotFoundError: # pysb not available - PysbPetabProblem = None import_model_pysb = None logger = get_logger(__name__, logging.WARNING) -# ID of model parameter that is to be added to SBML model to indicate -# preequilibration -PREEQ_INDICATOR_ID = "preequilibration_indicator" - def _add_global_parameter( sbml_model: libsbml.Model, @@ -97,25 +96,26 @@ def get_fixed_parameters( :return: List of IDs of parameters which are to be considered constant. """ - # initial concentrations for species or initial compartment sizes in - # condition table will need to be turned into fixed parameters - - # if there is no initial assignment for that species, we'd need - # to create one. to avoid any naming collision right away, we don't - # allow that for now - - # we can't handle them yet - compartments = [ - col - for col in petab_problem.condition_df - if petab_problem.sbml_model.getCompartment(col) is not None - ] - if compartments: - raise NotImplementedError( - "Can't handle initial compartment sizes " - "at the moment. Consider creating an " - f"initial assignment for {compartments}" - ) + if petab_problem.model.type_id == MODEL_TYPE_SBML: + # initial concentrations for species or initial compartment sizes in + # condition table will need to be turned into fixed parameters + + # if there is no initial assignment for that species, we'd need + # to create one. to avoid any naming collision right away, we don't + # allow that for now + + # we can't handle them yet + compartments = [ + col + for col in petab_problem.condition_df + if petab_problem.model.sbml_model.getCompartment(col) is not None + ] + if compartments: + raise NotImplementedError( + "Can't handle initial compartment sizes " + "at the moment. Consider creating an " + f"initial assignment for {compartments}" + ) # if we have a parameter table, all parameters that are allowed to be # listed in the parameter table, but are not marked as estimated, can be @@ -142,9 +142,6 @@ def get_fixed_parameters( estimated_parameters = petab_problem.parameter_df.index.values fixed_parameters = set(all_parameters) - set(estimated_parameters) - sbml_model = petab_problem.sbml_model - condition_df = petab_problem.condition_df - # Column names are model parameter IDs, compartment IDs or species IDs. # Thereof, all parameters except for any overridden ones should be made # constant. @@ -152,6 +149,7 @@ def get_fixed_parameters( # increase model reusability) # handle parameters in condition table + condition_df = petab_problem.condition_df if condition_df is not None: logger.debug(f"Condition table: {condition_df.shape}") @@ -167,30 +165,31 @@ def get_fixed_parameters( # of overriding if condition_df[p].dtype != "O" # p is a parameter - and sbml_model.getParameter(p) is not None - # but not a rule target - and sbml_model.getRuleByVariable(p) is None + and not petab_problem.model.is_state_variable(p) ) # Ensure mentioned parameters exist in the model. Remove additional ones # from list for fixed_parameter in fixed_parameters.copy(): # check global parameters - if not sbml_model.getParameter(fixed_parameter): + if not petab_problem.model.has_entity_with_id(fixed_parameter): + # TODO: could still exist as an output parameter? logger.warning( - f"Parameter or species '{fixed_parameter}'" - " provided in condition table but not present in" - " model. Ignoring." + f"Column '{fixed_parameter}' used in condition " + "table but not entity with the corresponding ID " + "exists. Ignoring." ) fixed_parameters.remove(fixed_parameter) - # exclude targets of rules or initial assignments - for fixed_parameter in fixed_parameters.copy(): - # check global parameters - if sbml_model.getInitialAssignmentBySymbol( - fixed_parameter - ) or sbml_model.getRuleByVariable(fixed_parameter): - fixed_parameters.remove(fixed_parameter) + if petab_problem.model.type_id == MODEL_TYPE_SBML: + # exclude targets of rules or initial assignments + sbml_model = petab_problem.model.sbml_model + for fixed_parameter in fixed_parameters.copy(): + # check global parameters + if sbml_model.getInitialAssignmentBySymbol( + fixed_parameter + ) or sbml_model.getRuleByVariable(fixed_parameter): + fixed_parameters.remove(fixed_parameter) return list(sorted(fixed_parameters)) @@ -300,17 +299,25 @@ def import_petab_problem( :return: The imported model. """ - # extract model name from pysb - if ( - PysbPetabProblem - and isinstance(petab_problem, PysbPetabProblem) - and model_name is None - ): + if petab_problem.model.type_id not in (MODEL_TYPE_SBML, MODEL_TYPE_PYSB): + raise NotImplementedError( + "Unsupported model type " + petab_problem.model.type_id + ) + + if petab_problem.mapping_df is not None: + # It's partially supported. Remove at your own risk... + raise NotImplementedError("PEtab v2.0.0 mapping tables are not yet supported.") + + model_name = model_name or petab_problem.model.model_id + + if petab_problem.model.type_id == MODEL_TYPE_PYSB and model_name is None: model_name = petab_problem.pysb_model.name + elif model_name is None and model_output_dir: + model_name = _create_model_name(model_output_dir) # generate folder and model name if necessary if model_output_dir is None: - if PysbPetabProblem and isinstance(petab_problem, PysbPetabProblem): + if petab_problem.model.type_id == MODEL_TYPE_PYSB: raise ValueError("Parameter `model_output_dir` is required.") model_output_dir = _create_model_output_dir_name( @@ -319,9 +326,6 @@ def import_petab_problem( else: model_output_dir = os.path.abspath(model_output_dir) - if model_name is None: - model_name = _create_model_name(model_output_dir) - # create folder if not os.path.exists(model_output_dir): os.makedirs(model_output_dir) @@ -341,7 +345,7 @@ def import_petab_problem( logger.info(f"Compiling model {model_name} to {model_output_dir}.") # compile the model - if PysbPetabProblem and isinstance(petab_problem, PysbPetabProblem): + if petab_problem.model.type_id == MODEL_TYPE_PYSB: import_model_pysb( petab_problem, model_name=model_name, @@ -664,13 +668,11 @@ def import_model_sbml( # TODO: to parameterize initial states or compartment sizes, we currently # need initial assignments. if they occur in the condition table, we - # create a new parameter initial_${startOrCompartmentID}. + # create a new parameter initial_${speciesOrCompartmentID}. # feels dirty and should be changed (see also #924) # - initial_states = [ - col for col in petab_problem.condition_df if element_is_state(sbml_model, col) - ] + initial_states = get_states_in_condition_table(petab_problem) fixed_parameters = [] if initial_states: # add preequilibration indicator variable @@ -691,7 +693,7 @@ def import_model_sbml( logger.debug( "Adding preequilibration indicator " f"constant {PREEQ_INDICATOR_ID}" ) - logger.debug(f"Adding initial assignments for {initial_states}") + logger.debug(f"Adding initial assignments for {initial_states.keys()}") for assignee_id in initial_states: init_par_id_preeq = f"initial_{assignee_id}_preeq" init_par_id_sim = f"initial_{assignee_id}_sim" @@ -782,7 +784,6 @@ def get_observation_model( :return: Tuple of dicts with observables, noise distributions, and sigmas. """ - if observable_df is None: return {}, {}, {} @@ -871,20 +872,6 @@ 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 @@ -892,7 +879,6 @@ def _parse_cli_args(): :return: Parsed CLI arguments from :mod:`argparse`. """ - parser = argparse.ArgumentParser( description="Import PEtab-format model into AMICI." ) diff --git a/python/sdist/amici/petab_import_pysb.py b/python/sdist/amici/petab_import_pysb.py index b3d7696dd1..63c1dd9681 100644 --- a/python/sdist/amici/petab_import_pysb.py +++ b/python/sdist/amici/petab_import_pysb.py @@ -6,323 +6,153 @@ """ import logging -import os -from itertools import chain +import re from pathlib import Path -from typing import Dict, Iterable, Optional, Union +from typing import Optional, Union import petab import pysb +import pysb.bng import sympy as sp -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 petab.models.sbml_model import SbmlModel +from petab.C import CONDITION_NAME, NOISE_FORMULA, OBSERVABLE_FORMULA +from petab.models.pysb_model import PySBModel from .logging import get_logger, log_execution_time, set_log_level +from .petab_util import PREEQ_INDICATOR_ID, get_states_in_condition_table logger = get_logger(__name__, logging.WARNING) -class PysbPetabProblem(petab.Problem): - """Representation of a PySB-model-based PEtab problem +def _add_observation_model(pysb_model: pysb.Model, petab_problem: petab.Problem): + """Extend PySB model by observation model as defined in the PEtab + observables table""" + + # add any required output parameters + local_syms = { + sp.Symbol.__str__(comp): comp + for comp in pysb_model.components + if isinstance(comp, sp.Symbol) + } + for formula in [ + *petab_problem.observable_df[OBSERVABLE_FORMULA], + *petab_problem.observable_df[NOISE_FORMULA], + ]: + sym = sp.sympify(formula, locals=local_syms) + for s in sym.free_symbols: + if not isinstance(s, pysb.Component): + p = pysb.Parameter(str(s), 1.0) + pysb_model.add_component(p) + local_syms[sp.Symbol.__str__(p)] = p + + # add observables and sigmas to pysb model + for observable_id, observable_formula, noise_formula in zip( + petab_problem.observable_df.index, + petab_problem.observable_df[OBSERVABLE_FORMULA], + petab_problem.observable_df[NOISE_FORMULA], + ): + obs_symbol = sp.sympify(observable_formula, locals=local_syms) + if observable_id in pysb_model.expressions.keys(): + obs_expr = pysb_model.expressions[observable_id] + else: + obs_expr = pysb.Expression(observable_id, obs_symbol) + pysb_model.add_component(obs_expr) + local_syms[observable_id] = obs_expr - This class extends :class:`petab.Problem` with a PySB model. - The model is augmented with the observation model based on the PEtab - observable table. - For now, a dummy SBML model is created which allows used the existing - SBML-PEtab API. + sigma_id = f"{observable_id}_sigma" + sigma_symbol = sp.sympify(noise_formula, locals=local_syms) + sigma_expr = pysb.Expression(sigma_id, sigma_symbol) + pysb_model.add_component(sigma_expr) + local_syms[sigma_id] = sigma_expr - :ivar pysb_model: - PySB model instance from of this PEtab problem. - """ +def _add_initialization_variables(pysb_model: pysb.Model, petab_problem: petab.Problem): + """Add initialization variables to the PySB model to support initial + conditions specified in the PEtab condition table. - def __init__(self, pysb_model: "pysb.Model" = None, *args, **kwargs): - """ - Constructor - - :param pysb_model: PySB model instance for this PEtab problem - :param args: See :meth:`petab.Problem.__init__` - :param kwargs: See :meth:`petab.Problem.__init__` - """ - flatten = kwargs.pop("flatten", False) - super().__init__(*args, **kwargs) - if flatten: - petab.flatten_timepoint_specific_output_overrides(self) - - self.pysb_model: "pysb.Model" = pysb_model - self._add_observation_model() - - if self.pysb_model is not None: - self.model = create_dummy_sbml( - self.pysb_model, - observable_ids=self.observable_df.index.values - if self.observable_df is not None - else None, - ) - - def _add_observation_model(self): - """Extend PySB model by observation model as defined in the PEtab - observables table""" - - # add any required output parameters - local_syms = { - sp.Symbol.__str__(comp): comp - for comp in self.pysb_model.components - if isinstance(comp, sp.Symbol) - } - for formula in [ - *self.observable_df[OBSERVABLE_FORMULA], - *self.observable_df[NOISE_FORMULA], - ]: - sym = sp.sympify(formula, locals=local_syms) - for s in sym.free_symbols: - if not isinstance(s, pysb.Component): - p = pysb.Parameter(str(s), 1.0, _export=False) - self.pysb_model.add_component(p) - local_syms[sp.Symbol.__str__(p)] = p - - # add observables and sigmas to pysb model - for observable_id, observable_formula, noise_formula in zip( - self.observable_df.index, - self.observable_df[OBSERVABLE_FORMULA], - self.observable_df[NOISE_FORMULA], - ): - obs_symbol = sp.sympify(observable_formula, locals=local_syms) - if observable_id in self.pysb_model.expressions.keys(): - obs_expr = self.pysb_model.expressions[observable_id] - else: - obs_expr = pysb.Expression(observable_id, obs_symbol, _export=False) - self.pysb_model.add_component(obs_expr) - local_syms[observable_id] = obs_expr - - sigma_id = f"{observable_id}_sigma" - sigma_symbol = sp.sympify(noise_formula, locals=local_syms) - sigma_expr = pysb.Expression(sigma_id, sigma_symbol, _export=False) - self.pysb_model.add_component(sigma_expr) - local_syms[sigma_id] = sigma_expr - - @staticmethod - 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. - - :param condition_file: - PEtab condition table - - :param measurement_file: - PEtab measurement table - - :param parameter_file: - PEtab parameter table - - :param visualization_files: - PEtab visualization tables - - :param observable_files: - PEtab observables tables - - :param pysb_model_file: - PySB model file - - :param flatten: - Flatten the petab problem - - :return: - Petab Problem - """ - - condition_df = measurement_df = parameter_df = visualization_df = None - observable_df = None - - if condition_file: - condition_df = petab.conditions.get_condition_df(condition_file) - - if measurement_file: - # If there are multiple tables, we will merge them - measurement_df = petab.core.concat_tables( - measurement_file, petab.measurements.get_measurement_df - ) - - if parameter_file: - parameter_df = petab.parameters.get_parameter_df(parameter_file) - - if visualization_files: - # If there are multiple tables, we will merge them - visualization_df = petab.core.concat_tables( - visualization_files, petab.core.get_visualization_df - ) + To parameterize initial states, we currently need initial assignments. + If they occur in the condition table, we create a new parameter + initial_${speciesID}. Feels dirty and should be changed (see also #924). + """ - if observable_files: - # If there are multiple tables, we will merge them - observable_df = petab.core.concat_tables( - observable_files, petab.observables.get_observable_df + initial_states = get_states_in_condition_table(petab_problem) + fixed_parameters = [] + 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 + if PREEQ_INDICATOR_ID in [c.name for c in pysb_model.components]: + raise AssertionError( + "Model already has a component with ID " + f"{PREEQ_INDICATOR_ID}. Cannot handle " + "species and compartments in condition table " + "then." ) - from amici.pysb_import import pysb_model_from_path - - return PysbPetabProblem( - pysb_model=pysb_model_from_path(pysb_model_file=pysb_model_file), - condition_df=condition_df, - measurement_df=measurement_df, - parameter_df=parameter_df, - observable_df=observable_df, - visualization_df=visualization_df, - flatten=flatten, + preeq_indicator = pysb.Parameter(PREEQ_INDICATOR_ID) + pysb_model.add_component(preeq_indicator) + # Can only reset parameters after preequilibration if they are fixed. + fixed_parameters.append(PREEQ_INDICATOR_ID) + logger.debug( + "Adding preequilibration indicator constant " f"{PREEQ_INDICATOR_ID}" ) - - @staticmethod - 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. - - NOTE: The PySB model is currently expected in the YAML file under - ``sbml_files``. - - :param yaml_config: - PEtab configuration as dictionary or YAML file name - - :param flatten: - Flatten the petab problem - - :return: - Petab Problem - """ - from petab.yaml import ( - assert_single_condition_and_sbml_file, - is_composite_problem, - load_yaml, + logger.debug(f"Adding initial assignments for {initial_states.keys()}") + + 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]: + if init_par_id in [c.name for c in pysb_model.components]: + raise ValueError( + "Cannot create parameter for initial assignment " + f"for {assignee_id} because an entity named " + f"{init_par_id} exists already in the model." + ) + p = pysb.Parameter(init_par_id) + pysb_model.add_component(p) + + species_idx = int(re.match(r"__s(\d+)$", assignee_id)[1]) + # use original model here since that's what was used to generate + # the ids in initial_states + species_pattern = petab_problem.model.model.species[species_idx] + + # species pattern comes from the _original_ model, but we only want + # to modify pysb_model, so we have to reconstitute the pattern using + # pysb_model + for c in pysb_model.components: + globals()[c.name] = c + species_pattern = pysb.as_complex_pattern(eval(str(species_pattern))) + + from pysb.pattern import match_complex_pattern + + formula = pysb.Expression( + f"initial_{assignee_id}_formula", + preeq_indicator * pysb_model.parameters[init_par_id_preeq] + + (1 - preeq_indicator) * pysb_model.parameters[init_par_id_sim], ) - - if isinstance(yaml_config, (str, Path)): - path_prefix = os.path.dirname(yaml_config) - yaml_config = load_yaml(yaml_config) - else: - path_prefix = "" - - if is_composite_problem(yaml_config): - raise ValueError( - "petab.Problem.from_yaml() can only be used for " - "yaml files comprising a single model. " - "Consider using " - "petab.CompositeProblem.from_yaml() instead." - ) - - if yaml_config[FORMAT_VERSION] != petab.__format_version__: - raise ValueError( - "Provided PEtab files are of unsupported version" - f"{yaml_config[FORMAT_VERSION]}. Expected " - f"{petab.__format_version__}." - ) - - problem0 = yaml_config["problems"][0] - - assert_single_condition_and_sbml_file(problem0) - - if isinstance(yaml_config[PARAMETER_FILE], list): - parameter_file = [ - os.path.join(path_prefix, f) for f in yaml_config[PARAMETER_FILE] - ] + pysb_model.add_component(formula) + + for initial in pysb_model.initials: + if match_complex_pattern(initial.pattern, species_pattern, exact=True): + logger.debug( + "The PySB model has an initial defined for species " + f"{assignee_id}, but this species also has an initial " + "value defined in the PEtab condition table. The SBML " + "initial assignment will be overwritten to handle " + "preequilibration and initial values specified by the " + "PEtab problem." + ) + initial.value = formula + break else: - parameter_file = os.path.join(path_prefix, yaml_config[PARAMETER_FILE]) - - return PysbPetabProblem.from_files( - pysb_model_file=os.path.join(path_prefix, problem0[SBML_FILES][0]), - measurement_file=[ - os.path.join(path_prefix, f) for f in problem0[MEASUREMENT_FILES] - ], - condition_file=os.path.join(path_prefix, problem0[CONDITION_FILES][0]), - parameter_file=parameter_file, - visualization_files=[ - os.path.join(path_prefix, f) - for f in problem0.get(VISUALIZATION_FILES, []) - ], - observable_files=[ - os.path.join(path_prefix, f) for f in problem0.get(OBSERVABLE_FILES, []) - ], - flatten=flatten, - ) - + # No initial in the pysb model, so add one + init = pysb.Initial(species_pattern, formula) + pysb_model.add_component(init) -def create_dummy_sbml( - pysb_model: "pysb.Model", observable_ids: Optional[Iterable[str]] = None -) -> SbmlModel: - """Create SBML dummy model for to use PySB models with PEtab. - - Model must at least contain PEtab problem parameter and noise parameters - for observables. - - :param pysb_model: PySB model - :param observable_ids: Observable IDs - :return: A dummy petab SBML model. - """ - import libsbml - - document = libsbml.SBMLDocument(3, 1) - dummy_sbml_model = document.createModel() - dummy_sbml_model.setTimeUnits("second") - dummy_sbml_model.setExtentUnits("mole") - dummy_sbml_model.setSubstanceUnits("mole") - - # mandatory if there are species - c = dummy_sbml_model.createCompartment() - c.setId("dummy_compartment") - c.setConstant(False) - - # parameters are required for parameter mapping - for parameter in pysb_model.parameters: - p = dummy_sbml_model.createParameter() - p.setId(parameter.name) - p.setConstant(True) - p.setValue(0.0) - - # noise parameters are required for every observable - for observable_id in observable_ids: - p = dummy_sbml_model.createParameter() - p.setId(f"noiseParameter1_{observable_id}") - p.setConstant(True) - p.setValue(0.0) - - # pysb observables and expressions are required in case they occur in - # the observableFormula or noiseFormula. - # as this code is only temporary and not performance-critical, we just add - # all of them. we just need an sbml entity with the same ID. sbml species - # seem to be the simplest, as parameters would interfere with parameter - # mapping later on - for component in chain(pysb_model.expressions, pysb_model.observables): - s = dummy_sbml_model.createSpecies() - s.setId(component.name) - s.setInitialAmount(0.0) - s.setHasOnlySubstanceUnits(False) - s.setBoundaryCondition(False) - s.setCompartment("dummy_compartment") - s.setConstant(False) - - return SbmlModel(sbml_model=dummy_sbml_model, sbml_document=document) + return fixed_parameters @log_execution_time("Importing PEtab model", logger) def import_model_pysb( - petab_problem: PysbPetabProblem, + petab_problem: petab.Problem, model_output_dir: Optional[Union[str, Path]] = None, verbose: Optional[Union[bool, int]] = True, model_name: Optional[str] = None, @@ -352,30 +182,55 @@ def import_model_pysb( logger.info("Importing model ...") - observable_table = petab_problem.observable_df - pysb_model = petab_problem.pysb_model + if not isinstance(petab_problem.model, PySBModel): + raise ValueError("Not a PySB model") + + # need to create a copy here as we don't want to modify the original + pysb.SelfExporter.cleanup() + og_export = pysb.SelfExporter.do_export + pysb.SelfExporter.do_export = False + pysb_model = pysb.Model( + base=petab_problem.model.model, + name=petab_problem.model.model_id, + ) + + _add_observation_model(pysb_model, petab_problem) + # generate species for the _original_ model + pysb.bng.generate_equations(petab_problem.model.model) + fixed_parameters = _add_initialization_variables(pysb_model, petab_problem) + pysb.SelfExporter.do_export = og_export - # For pysb, we only allow parameters in the condition table - # those must be pysb model parameters (either natively, or output - # parameters from measurement or condition table that have been added in - # PysbPetabProblem) + # check condition table for supported features, important to use pysb_model + # here, as we want to also cover output parameters model_parameters = [p.name for p in pysb_model.parameters] + condition_species_parameters = get_states_in_condition_table( + petab_problem, return_patterns=True + ) for x in petab_problem.condition_df.columns: if x == CONDITION_NAME: continue - if x not in model_parameters: - raise NotImplementedError( - "For PySB PEtab import, only model parameters, but no states " - "or compartments are allowed in the condition table." - f"Offending column: {x}" - ) + x = petab.mapping.resolve_mapping(petab_problem.mapping_df, x) + + # parameters + if x in model_parameters: + continue + + # species/pattern + if x in condition_species_parameters: + continue + + raise NotImplementedError( + "For PySB PEtab import, only model parameters and species, but " + "not compartments are allowed in the condition table. Offending " + f"column: {x}" + ) from .petab_import import get_fixed_parameters, petab_noise_distributions_to_amici - constant_parameters = get_fixed_parameters(petab_problem) + constant_parameters = get_fixed_parameters(petab_problem) + fixed_parameters - if observable_table is None: + if petab_problem.observable_df is None: observables = None sigmas = None noise_distrs = None @@ -383,12 +238,12 @@ def import_model_pysb( observables = [ expr.name for expr in pysb_model.expressions - if expr.name in observable_table.index + if expr.name in petab_problem.observable_df.index ] sigmas = {obs_id: f"{obs_id}_sigma" for obs_id in observables} - noise_distrs = petab_noise_distributions_to_amici(observable_table) + noise_distrs = petab_noise_distributions_to_amici(petab_problem.observable_df) from amici.pysb_import import pysb2amici diff --git a/python/sdist/amici/petab_objective.py b/python/sdist/amici/petab_objective.py index 2d0abc12da..f518724c82 100644 --- a/python/sdist/amici/petab_objective.py +++ b/python/sdist/amici/petab_objective.py @@ -8,6 +8,7 @@ import copy import logging import numbers +import re from typing import ( Any, Collection, @@ -28,6 +29,7 @@ import sympy as sp from amici.sbml_import import get_species_initial from petab.C import * # noqa: F403 +from petab.models import MODEL_TYPE_PYSB, MODEL_TYPE_SBML from sympy.abc import _clash from . import AmiciExpData, AmiciModel @@ -37,7 +39,13 @@ ParameterMappingForCondition, fill_in_parameters, ) -from .petab_import import PREEQ_INDICATOR_ID, element_is_state +from .petab_import import PREEQ_INDICATOR_ID +from .petab_util import get_states_in_condition_table + +try: + import pysb +except ImportError: + pysb = None logger = get_logger(__name__) @@ -466,15 +474,18 @@ def create_parameter_mapping( # Because AMICI globalizes all local parameters during model import, # we need to do that here as well to prevent parameter mapping errors # (PEtab does currently not care about SBML LocalParameters) - if petab_problem.sbml_document: - converter_config = libsbml.SBMLLocalParameterConverter().getDefaultProperties() - petab_problem.sbml_document.convert(converter_config) - else: - logger.debug( - "No petab_problem.sbml_document is set. Cannot convert " - "SBML LocalParameters. If the model contains " - "LocalParameters, parameter mapping will fail." - ) + if petab_problem.model.type_id == MODEL_TYPE_SBML: + if petab_problem.sbml_document: + converter_config = ( + libsbml.SBMLLocalParameterConverter().getDefaultProperties() + ) + petab_problem.sbml_document.convert(converter_config) + else: + logger.debug( + "No petab_problem.sbml_document is set. Cannot " + "convert SBML LocalParameters. If the model contains " + "LocalParameters, parameter mapping will fail." + ) default_parameter_mapping_kwargs = { "warn_unmapped": False, @@ -491,6 +502,7 @@ def create_parameter_mapping( measurement_df=petab_problem.measurement_df, parameter_df=petab_problem.parameter_df, observable_df=petab_problem.observable_df, + mapping_df=petab_problem.mapping_df, model=petab_problem.model, **dict(default_parameter_mapping_kwargs, **parameter_mapping_kwargs), ) @@ -507,6 +519,108 @@ def create_parameter_mapping( return parameter_mapping +def _get_initial_state_sbml( + petab_problem: petab.Problem, element_id: str +) -> Union[float, sp.Basic]: + 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()), locals=_clash + ) + 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." + ) + return value + + +def _get_initial_state_pysb( + petab_problem: petab.Problem, element_id: str +) -> Union[float, sp.Symbol]: + species_idx = int(re.match(r"__s(\d+)$", element_id)[1]) + species_pattern = petab_problem.model.model.species[species_idx] + from pysb.pattern import match_complex_pattern + + value = next( + ( + initial.value + for initial in petab_problem.model.model.initials + if match_complex_pattern(initial.pattern, species_pattern, exact=True) + ), + 0.0, + ) + if isinstance(value, pysb.Parameter): + if value.name in petab_problem.parameter_df.index: + value = value.name + else: + value = value.value + + return value + + +def _set_initial_state( + petab_problem, + condition_id, + element_id, + init_par_id, + par_map, + scale_map, + value, +): + value = petab.to_float_if_float(value) + if pd.isna(value): + if petab_problem.model.type_id == MODEL_TYPE_SBML: + value = _get_initial_state_sbml(petab_problem, element_id) + elif petab_problem.model.type_id == MODEL_TYPE_PYSB: + value = _get_initial_state_pysb(petab_problem, element_id) + + try: + value = float(value) + except (ValueError, TypeError): + if sp.nsimplify(value).is_Atom and ( + pysb is None or not isinstance(value, pysb.Component) + ): + # Get rid of multiplication with one + value = sp.nsimplify(value) + else: + raise NotImplementedError( + "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 {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 " + "defined in the SBML model." + ) + par_map[init_par_id] = value + if isinstance(value, float): + # numeric initial state + scale_map[init_par_id] = petab.LIN + else: + # parametric initial state + scale_map[init_par_id] = petab_problem.parameter_df[PARAMETER_SCALE].get( + value, petab.LIN + ) + + def create_parameter_mapping_for_condition( parameter_mapping_for_condition: petab.ParMappingDictQuadruple, condition: Union[pd.Series, Dict], @@ -563,12 +677,9 @@ def create_parameter_mapping_for_condition( # ExpData.x0, but in the case of preequilibration this would not allow for # resetting initial states. - states_in_condition_table = [ - col - for col in petab_problem.condition_df - if element_is_state(petab_problem.sbml_model, col) - ] - if states_in_condition_table: + if states_in_condition_table := get_states_in_condition_table( + petab_problem, condition + ): # 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...) @@ -579,88 +690,20 @@ 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_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, element_id] - ) - if pd.isna(value): - 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()), - locals=_clash, - ) - 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): - if sp.nsimplify(value).is_Atom: - # Get rid of multiplication with one - value = sp.nsimplify(value) - else: - raise NotImplementedError( - "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 {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 " - "defined in the SBML model." - ) - par_map[init_par_id] = value - if isinstance(value, float): - # numeric initial state - scale_map[init_par_id] = petab.LIN - else: - # parametric initial state - scale_map[init_par_id] = petab_problem.parameter_df[ - PARAMETER_SCALE - ].get(value, petab.LIN) - - for element_id in states_in_condition_table: + for element_id, (value, preeq_value) in states_in_condition_table.items(): # for preequilibration init_par_id = f"initial_{element_id}_preeq" - if condition.get(PREEQUILIBRATION_CONDITION_ID): - condition_id = condition[PREEQUILIBRATION_CONDITION_ID] + if ( + condition_id := condition.get(PREEQUILIBRATION_CONDITION_ID) + ) is not None: _set_initial_state( + petab_problem, condition_id, element_id, init_par_id, condition_map_preeq, condition_scale_map_preeq, + preeq_value, ) else: # need to set dummy value for preeq parameter anyways, as it @@ -673,11 +716,13 @@ def _set_initial_state( condition_id = condition[SIMULATION_CONDITION_ID] init_par_id = f"initial_{element_id}_sim" _set_initial_state( + petab_problem, condition_id, element_id, init_par_id, condition_map_sim, condition_scale_map_sim, + value, ) ########################################################################## @@ -757,26 +802,23 @@ def create_edatas( observable_ids = amici_model.getObservableIds() + measurement_groupvar = [SIMULATION_CONDITION_ID] + if PREEQUILIBRATION_CONDITION_ID in simulation_conditions: + measurement_groupvar.append(petab.PREEQUILIBRATION_CONDITION_ID) measurement_dfs = dict( - list( - petab_problem.measurement_df.groupby( - [petab.SIMULATION_CONDITION_ID, petab.PREEQUILIBRATION_CONDITION_ID] - if petab.PREEQUILIBRATION_CONDITION_ID in simulation_conditions - else petab.SIMULATION_CONDITION_ID - ) - ) + list(petab_problem.measurement_df.groupby(measurement_groupvar)) ) edatas = [] for _, condition in simulation_conditions.iterrows(): # Create amici.ExpData for each simulation - if petab.PREEQUILIBRATION_CONDITION_ID in condition: + if PREEQUILIBRATION_CONDITION_ID in condition: measurement_index = ( - condition.get(petab.SIMULATION_CONDITION_ID), - condition.get(petab.PREEQUILIBRATION_CONDITION_ID), + condition.get(SIMULATION_CONDITION_ID), + condition.get(PREEQUILIBRATION_CONDITION_ID), ) else: - measurement_index = condition.get(petab.SIMULATION_CONDITION_ID) + measurement_index = (condition.get(SIMULATION_CONDITION_ID),) edata = create_edata_for_condition( condition=condition, amici_model=amici_model, @@ -828,18 +870,16 @@ def create_edata_for_condition( edata.id += "+" + condition.get(PREEQUILIBRATION_CONDITION_ID) ########################################################################## # enable initial parameters reinitialization - 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 element_is_state(petab_problem.sbml_model, col) - ] + + states_in_condition_table = get_states_in_condition_table( + petab_problem, condition=condition + ) if condition.get(PREEQUILIBRATION_CONDITION_ID) and states_in_condition_table: state_ids = amici_model.getStateIds() state_idx_reinitalization = [ - state_ids.index(s) for s in states_in_condition_table + state_ids.index(s) + for s, (v, v_preeq) in states_in_condition_table.items() + if not np.isnan(v) ] edata.reinitialization_state_idxs_sim = state_idx_reinitalization logger.debug( diff --git a/python/sdist/amici/petab_util.py b/python/sdist/amici/petab_util.py new file mode 100644 index 0000000000..31f1ae1313 --- /dev/null +++ b/python/sdist/amici/petab_util.py @@ -0,0 +1,102 @@ +"""Various helper functions for working with PEtab problems.""" +import re +from typing import Dict, Tuple, Union + +import libsbml +import pandas as pd +import petab +from petab.C import PREEQUILIBRATION_CONDITION_ID, SIMULATION_CONDITION_ID +from petab.mapping import resolve_mapping +from petab.models import MODEL_TYPE_PYSB, MODEL_TYPE_SBML + +# ID of model parameter that is to be added to SBML model to indicate +# preequilibration +PREEQ_INDICATOR_ID = "preequilibration_indicator" + + +def get_states_in_condition_table( + petab_problem: petab.Problem, + condition: Union[Dict, pd.Series] = None, + return_patterns: bool = False, +) -> Dict[str, Tuple[Union[float, str, None], Union[float, str, None]]]: + """Get states and their initial condition as specified in the condition table. + + Returns: Dictionary: ``stateId -> (initial condition simulation, initial condition preequilibration)`` + """ + if petab_problem.model.type_id not in (MODEL_TYPE_SBML, MODEL_TYPE_PYSB): + raise NotImplementedError() + + species_check_funs = { + MODEL_TYPE_SBML: lambda x: _element_is_sbml_state( + petab_problem.sbml_model, x + ), + MODEL_TYPE_PYSB: lambda x: _element_is_pysb_pattern( + petab_problem.model.model, x + ), + } + states = { + resolve_mapping(petab_problem.mapping_df, col): (None, None) + if condition is None + else ( + petab_problem.condition_df.loc[ + condition[SIMULATION_CONDITION_ID], col + ], + petab_problem.condition_df.loc[ + condition[PREEQUILIBRATION_CONDITION_ID], col + ] + if PREEQUILIBRATION_CONDITION_ID in condition + else None, + ) + for col in petab_problem.condition_df.columns + if species_check_funs[petab_problem.model.type_id]( + resolve_mapping(petab_problem.mapping_df, col) + ) + } + + if petab_problem.model.type_id == MODEL_TYPE_PYSB: + if return_patterns: + return states + import pysb.pattern + + try: + spm = pysb.pattern.SpeciesPatternMatcher( + model=petab_problem.model.model + ) + except NotImplementedError as e: + raise NotImplementedError( + "Requires https://github.com/pysb/pysb/pull/570. " + "To use this functionality, update pysb via " + "`pip install git+https://github.com/FFroehlich/pysb@fix_pattern_matching`" + ) + + # expose model components as variables so we can evaluate patterns + for c in petab_problem.model.model.components: + globals()[c.name] = c + + states = { + f"__s{ix}": value + for pattern, value in states.items() + for ix in spm.match(eval(pattern), index=True, exact=True) + } + return states + + +def _element_is_pysb_pattern(model: "pysb.Model", element: str) -> bool: + """Check if element is a pysb pattern""" + if match := re.match(r"[a-zA-Z_][\w_]*\(", element): + return match[0][:-1] in [m.name for m in model.monomers] + return False + + +def _element_is_sbml_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 diff --git a/python/sdist/setup.cfg b/python/sdist/setup.cfg index 26255e77be..733a2a4938 100644 --- a/python/sdist/setup.cfg +++ b/python/sdist/setup.cfg @@ -44,7 +44,7 @@ include_package_data = True zip_safe = False [options.extras_require] -petab = petab>=0.2.0 +petab = petab>=0.2.1 pysb = pysb>=1.13.1 test = pytest diff --git a/scripts/installAmiciSource.sh b/scripts/installAmiciSource.sh index 0e4954acfa..aa330bef22 100755 --- a/scripts/installAmiciSource.sh +++ b/scripts/installAmiciSource.sh @@ -30,6 +30,7 @@ fi pip install -U "setuptools<64" pip install --upgrade pip wheel pip install --upgrade pip scipy matplotlib coverage pytest \ - pytest-cov cmake_build_extension numpy -pip install --verbose -e ${AMICI_PATH}/python/sdist[petab,test,pysb,vis] --no-build-isolation + pytest-cov cmake_build_extension numpy +pip install git+https://github.com/FFroehlich/pysb@fix_pattern_matching # pin to PR for SPM with compartments +pip install --verbose -e ${AMICI_PATH}/python/sdist[petab,test,vis] --no-build-isolation deactivate diff --git a/tests/petab_test_suite/conftest.py b/tests/petab_test_suite/conftest.py index 46b525da34..086740509c 100644 --- a/tests/petab_test_suite/conftest.py +++ b/tests/petab_test_suite/conftest.py @@ -32,9 +32,7 @@ def parse_selection(selection_str: str) -> List[int]: def pytest_addoption(parser): """Add pytest CLI options""" parser.addoption("--petab-cases", help="Test cases to run") - # TODO: re-enable in #1800 - # parser.addoption("--only-pysb", help="Run only PySB tests", - # action="store_true") + parser.addoption("--only-pysb", help="Run only PySB tests", action="store_true") parser.addoption( "--only-sbml", help="Run only SBML tests", @@ -59,22 +57,24 @@ def pytest_generate_tests(metafunc): if metafunc.config.getoption("--only-sbml"): argvalues = [ (case, "sbml", version) - for version in ("v1.0.0",) + for version in ("v1.0.0", "v2.0.0") for case in ( test_numbers if test_numbers else get_cases("sbml", version=version) ) ] - # TODO: re-enable in #1800 - # elif metafunc.config.getoption("--only-pysb"): - # argvalues = [ - # (case, 'pysb', "v1.0.0") - # for case in (test_numbers if test_numbers - # else get_cases("pysb", version="v1.0.0")) - # ] + elif metafunc.config.getoption("--only-pysb"): + argvalues = [ + (case, "pysb", "v2.0.0") + for case in ( + test_numbers + if test_numbers + else get_cases("pysb", version="v2.0.0") + ) + ] else: argvalues = [] - for version in ("v1.0.0",): - for format in ("sbml",): + for version in ("v1.0.0", "v2.0.0"): + for format in ("sbml", "pysb"): argvalues.extend( (case, format, version) for case in test_numbers or get_cases(format, version) diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 15218088ad..59e2ce7723 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -13,7 +13,7 @@ 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 PysbPetabProblem, import_petab_problem +from amici.petab_import import import_petab_problem from amici.petab_objective import ( create_parameterized_edatas, rdatas_to_measurement_df, @@ -55,7 +55,7 @@ def _test_case(case, model_type, version): problem = petab.Problem.from_yaml(yaml_file) # compile amici model - if case.startswith("0006") and model_type != "pysb": + if case.startswith("0006"): petab.flatten_timepoint_specific_output_overrides(problem) model_name = f"petab_{model_type}_test_case_{case}" f"_{version.replace('.', '_')}" model_output_dir = f"amici_models/{model_name}" @@ -174,7 +174,7 @@ def run(): n_success = 0 n_skipped = 0 n_total = 0 - for version in ("v1.0.0",): + for version in ("v1.0.0", "v2.0.0"): cases = petabtests.get_cases("sbml", version=version) n_total += len(cases) for case in cases: From 4d767e7a675210aeb387ae653b219dce8357e6ce Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 22 May 2023 20:23:19 +0200 Subject: [PATCH 04/10] CMake: Fix choosing SWIG via env variable (#2100) Only changing SWIG_EXECUTABLE will lead to inconsistent settings. Need to unset version and directory, so FindSWIG will try to set them according to SWIG_EXECUTABLE. --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index d7f180e8a4..0c0623f263 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -95,6 +95,8 @@ endif() if(DEFINED ENV{SWIG}) message(STATUS "Setting SWIG_EXECUTABLE to $ENV{SWIG} ($SWIG)") + unset(SWIG_VERSION CACHE) + unset(SWIG_DIR CACHE) set(SWIG_EXECUTABLE $ENV{SWIG}) endif() From 238635c53901600aaaa9f08dd41d73fd7a5675d0 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 23 May 2023 14:58:20 +0200 Subject: [PATCH 05/10] Fix unquoted filename (#2105) Fixes issues with paths containing blanks: `sh: line 1: C:/Program: No such file or directory` https://github.com/AMICI-dev/AMICI/actions/runs/5054632314/jobs/9069784824?pr=2104 --- cmake/version.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/version.cmake b/cmake/version.cmake index 3f7b054a03..8dfe45c3c6 100644 --- a/cmake/version.cmake +++ b/cmake/version.cmake @@ -3,7 +3,7 @@ if(Git_FOUND) execute_process( COMMAND sh -c - "${GIT_EXECUTABLE} describe --abbrev=4 --dirty=-dirty --always --tags | cut -c 2- | tr -d '\n' | sed s/-/./" + "'${GIT_EXECUTABLE}' describe --abbrev=4 --dirty=-dirty --always --tags | cut -c 2- | tr -d '\n' | sed s/-/./" WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} OUTPUT_VARIABLE PROJECT_VERSION_GIT) endif() From 451440c1d86108d29d071354e05bf8ad233b6f5c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 23 May 2023 13:06:01 +0200 Subject: [PATCH 06/10] CMake: Try FindBLAS if no other information was provided (#2104) Let CMake's FindBLAS try to find BLAS libraries in case the user didn't provide any other information via environment variables. Relates to #2103 --- CMakeLists.txt | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c0623f263..cdf0ff8fed 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -153,12 +153,18 @@ if(${BLAS} STREQUAL "MKL" OR DEFINED ENV{MKLROOT}) CACHE STRING "") endif() elseif(NOT DEFINED ENV{BLAS_LIBS} AND NOT DEFINED ENV{BLAS_CFLAGS}) - set(BLAS_INCLUDE_DIRS - "" - CACHE STRING "") - set(BLAS_LIBRARIES - -lcblas - CACHE STRING "") + # if nothing is specified via environment variables, let's try FindBLAS + find_package(BLAS) + if(NOT BLAS_FOUND) + # Nothing specified by the user and FindBLAS didn't find anything; let's try + # if cblas is available on the system paths. + set(BLAS_INCLUDE_DIRS + "" + CACHE STRING "") + set(BLAS_LIBRARIES + -lcblas + CACHE STRING "") + endif() endif() add_compile_definitions(AMICI_BLAS_${BLAS}) From aa87f4553cbb3aa04cfe35ce38bb3bb3ec441fd7 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 25 May 2023 13:02:15 +0200 Subject: [PATCH 07/10] CI: Increase number of starts in ExampleSplinesSwameye2003.ipynb (#2113) ... as this fails too often. --- .../ExampleSplinesSwameye2003.ipynb | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb b/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb index 32f7047a65..6ef92d4c49 100644 --- a/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb +++ b/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb @@ -105,7 +105,7 @@ "source": [ "# If running as a Github action, just do the minimal amount of work required to check whether the code is working\n", "if os.getenv('GITHUB_ACTIONS') is not None:\n", - " n_starts = 10\n", + " n_starts = 15\n", " pypesto_optimizer = pypesto.optimize.FidesOptimizer(verbose=logging.WARNING, options=dict(maxiter=10))\n", " pypesto_engine = pypesto.engine.SingleCoreEngine()" ] @@ -362,7 +362,7 @@ "### Maximum Likelihood estimation\n", "Using pyPESTO we can optimize for the parameter vector that maximizes the probability of observing the experimental data (maximum likelihood estimation).\n", "\n", - "A multistart method with local gradient-based optimization is used and the results of each multistart can be visualized in a waterfall plot. " + "A multistart method with local gradient-based optimization is used and the results of each multistart can be visualized in a waterfall plot." ] }, { @@ -908,7 +908,7 @@ " np.log10(regstrength) # parameter is specified as log10 scale in PEtab\n", " )\n", " regproblem = copy.deepcopy(pypesto_problem)\n", - " \n", + "\n", " # Load existing results if available\n", " if os.path.exists(f'{name}.h5'):\n", " regresult = pypesto.store.read_result(f'{name}.h5', problem=regproblem)\n", @@ -916,7 +916,7 @@ " regresult = None\n", " # Overwrite\n", " # regresult = None\n", - " \n", + "\n", " # Parallel multistart optimization with pyPESTO and FIDES\n", " if n_starts > 0:\n", " if regresult is None:\n", @@ -935,10 +935,10 @@ " regresult.optimize_result.sort()\n", " if regresult.optimize_result.x[0] is None:\n", " raise Exception(\"All multistarts failed (n_starts is probably too small)! If this error occurred during CI, just run the workflow again.\")\n", - " \n", + "\n", " # Save results to disk\n", " # pypesto.store.write_result(regresult, f'{name}.h5', overwrite=True)\n", - " \n", + "\n", " # Store result\n", " regproblems[regstrength] = regproblem\n", " regresults[regstrength] = regresult" @@ -1469,7 +1469,7 @@ " np.log10(regstrength) # parameter is specified as log10 scale in PEtab\n", " )\n", " regproblem = copy.deepcopy(pypesto_problem)\n", - " \n", + "\n", " # Load existing results if available\n", " if os.path.exists(f'{name}.h5'):\n", " regresult = pypesto.store.read_result(f'{name}.h5', problem=regproblem)\n", @@ -1477,7 +1477,7 @@ " regresult = None\n", " # Overwrite\n", " # regresult = None\n", - " \n", + "\n", " # Parallel multistart optimization with pyPESTO and FIDES\n", " if n_starts > 0:\n", " if regresult is None:\n", @@ -1496,10 +1496,10 @@ " regresult.optimize_result.sort()\n", " if regresult.optimize_result.x[0] is None:\n", " raise Exception(\"All multistarts failed (n_starts is probably too small)! If this error occurred during CI, just run the workflow again.\")\n", - " \n", + "\n", " # Save results to disk\n", " # pypesto.store.write_result(regresult, f'{name}.h5', overwrite=True)\n", - " \n", + "\n", " # Store result\n", " regproblems[regstrength] = regproblem\n", " regresults[regstrength] = regresult" @@ -1934,7 +1934,7 @@ "## Bibliography\n", "Schelker, M. et al. (2012). “Comprehensive estimation of input signals and dynamics in biochemical reaction networks”. In: Bioinformatics 28.18, pp. i529–i534. doi: [10.1093/bioinformatics/bts393](https://doi.org/10.1093/bioinformatics/bts393).\n", "\n", - "Swameye, I. et al. (2003). “Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling”. In: Proceedings of the National Academy of Sciences 100.3, pp. 1028–1033. doi: [10.1073/pnas.0237333100](https://doi.org/10.1073/pnas.0237333100)." + "Swameye, I. et al. (2003). “Identification of nucleocytoplasmic cycling as a remote sensor in cellular signaling by databased modeling”. In: Proceedings of the National Academy of Sciences 100.3, pp. 1028–1033. doi: [10.1073/pnas.0237333100](https://doi.org/10.1073/pnas.0237333100).\n" ] } ], From f4252b941969ebc9a6f1b147a3ae3e304385223c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 25 May 2023 14:08:12 +0200 Subject: [PATCH 08/10] Pin nbsphinx (#2114) Fixes: `/home/runner/work/AMICI/AMICI/documentation/example_errors.ipynb:2262: WARNING: undefined label: /example_errors.ipynb#unsuccessful_factorization` Link was working until recently. Not sure what changed. --- documentation/rtd_requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/rtd_requirements.txt b/documentation/rtd_requirements.txt index a4d68af842..64bc03e519 100644 --- a/documentation/rtd_requirements.txt +++ b/documentation/rtd_requirements.txt @@ -4,7 +4,7 @@ mock>=5.0.2 setuptools==67.7.2 pysb>=1.11.0 matplotlib==3.7.1 -nbsphinx>=0.9.1 +nbsphinx==0.9.1 nbformat==5.8.0 recommonmark>=0.7.1 sphinx_rtd_theme>=1.2.0 From 3b4e12daa755adc8e5601fe8f8b52d2401c1317c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 25 May 2023 14:08:54 +0200 Subject: [PATCH 09/10] Fix cblas error for models without solver states in combination with forward sensitivities (#2108) --- src/model.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/model.cpp b/src/model.cpp index d09f383b2a..2aa8ee72ee 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1069,15 +1069,17 @@ void Model::getObservableSensitivity( // dydx A[ny,nx_solver] * sx B[nx_solver,nplist] = sy C[ny,nplist] // M K K N M N // lda ldb ldc - setNaNtoZero(derived_state_.dydx_); - setNaNtoZero(derived_state_.sx_); - amici_dgemm( - BLASLayout::colMajor, BLASTranspose::noTrans, BLASTranspose::noTrans, - ny, nplist(), nx_solver, 1.0, derived_state_.dydx_.data(), ny, - derived_state_.sx_.data(), nx_solver, 1.0, derived_state_.dydp_.data(), - ny - ); + if (nx_solver) { + setNaNtoZero(derived_state_.dydx_); + setNaNtoZero(derived_state_.sx_); + amici_dgemm( + BLASLayout::colMajor, BLASTranspose::noTrans, + BLASTranspose::noTrans, ny, nplist(), nx_solver, 1.0, + derived_state_.dydx_.data(), ny, derived_state_.sx_.data(), + nx_solver, 1.0, derived_state_.dydp_.data(), ny + ); + } writeSlice(derived_state_.dydp_, sy); if (always_check_finite_) From fd1e22e98607ee8840099e4c5c3cec3e81a22816 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 25 May 2023 17:49:51 +0200 Subject: [PATCH 10/10] Fix compilation error for models with events and xdot=0 (#2111) Fixes #2110 --- python/sdist/amici/de_export.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 46d226b8be..025d0b2656 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -1999,9 +1999,10 @@ def _compute_equation(self, name: str) -> None: # need to check if equations are zero since we are using # symbols - if not smart_is_zero_matrix(self.eq("stau")[ie]): + if not smart_is_zero_matrix(self.eq("stau")[ie]) \ + and not smart_is_zero_matrix(self.eq("xdot")): tmp_eq += smart_multiply( - (self.sym("xdot_old") - self.sym("xdot")), + self.sym("xdot_old") - self.sym("xdot"), self.sym("stau").T, )