diff --git a/.coveragerc b/.coveragerc index cf7406e8ba..5c8d382e6f 100644 --- a/.coveragerc +++ b/.coveragerc @@ -6,6 +6,7 @@ omit = */amici_without_hdf5.py */ThirdParty/* *.template.* + /tmp/* parallel = true [report] diff --git a/.github/workflows/deploy_branch.yml b/.github/workflows/deploy_branch.yml index 0e0ee6e61c..3e706300f3 100644 --- a/.github/workflows/deploy_branch.yml +++ b/.github/workflows/deploy_branch.yml @@ -5,7 +5,7 @@ jobs: sdist: name: Deploy Python Source Distribution - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -33,7 +33,7 @@ jobs: scripts/buildSdist.sh - name: "Upload artifact: sdist" - uses: actions/upload-artifact@v1 + uses: actions/upload-artifact@v3 with: name: sdist - path: python/sdist/dist + path: python/sdist/dist/amici-*.gz diff --git a/.github/workflows/deploy_protected.yml b/.github/workflows/deploy_protected.yml index b5f02fc38d..ac8b858a16 100644 --- a/.github/workflows/deploy_protected.yml +++ b/.github/workflows/deploy_protected.yml @@ -14,7 +14,7 @@ jobs: # https://github.com/marketplace/actions/publish-docker name: Deploy Dockerhub - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -28,7 +28,7 @@ jobs: - uses: actions/checkout@master - run: git archive -o docker/amici.tar.gz --format=tar.gz HEAD - name: Publish to Registry - uses: elgohr/Publish-Docker-Github-Action@2.8 + uses: elgohr/Publish-Docker-Github-Action@v4 with: name: dweindl/amici username: ${{ secrets.DOCKER_USERNAME }} diff --git a/.github/workflows/deploy_release.yml b/.github/workflows/deploy_release.yml index 48fcd50df5..77fe30a9be 100644 --- a/.github/workflows/deploy_release.yml +++ b/.github/workflows/deploy_release.yml @@ -8,7 +8,7 @@ jobs: pypi: name: Deploy PyPI - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -45,7 +45,7 @@ jobs: bioSimulatorsUpdateCliAndDockerImage: name: Release to BioSimulators needs: pypi - runs-on: ubuntu-latest + runs-on: ubuntu-22.04 env: # Owner/repository-id for the GitHub repository for the downstream command-line interface and Docker image DOWNSTREAM_REPOSITORY: biosimulators/Biosimulators_AMICI diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml index 9aa9b13f3f..b691231860 100644 --- a/.github/workflows/test_benchmark_collection_models.yml +++ b/.github/workflows/test_benchmark_collection_models.yml @@ -17,7 +17,7 @@ jobs: build: name: Benchmark Collection - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: diff --git a/.github/workflows/test_doc.yml b/.github/workflows/test_doc.yml index 4bcc6beea3..9742de92eb 100644 --- a/.github/workflows/test_doc.yml +++ b/.github/workflows/test_doc.yml @@ -16,7 +16,7 @@ jobs: doxygen: name: Test Doxygen - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -51,7 +51,7 @@ jobs: sphinx: name: Test Sphinx - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: diff --git a/.github/workflows/test_install.yml b/.github/workflows/test_install.yml index 58ad6693a5..dd75ff576e 100644 --- a/.github/workflows/test_install.yml +++ b/.github/workflows/test_install.yml @@ -5,7 +5,7 @@ jobs: archive: name: Archive Install - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -53,7 +53,7 @@ jobs: sdist: name: sdist Install - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: diff --git a/.github/workflows/test_performance.yml b/.github/workflows/test_performance.yml index 1a7e12c53b..4fd5ba0294 100644 --- a/.github/workflows/test_performance.yml +++ b/.github/workflows/test_performance.yml @@ -19,7 +19,7 @@ jobs: build: name: Performance Test - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: @@ -66,7 +66,7 @@ jobs: AMICI_IMPORT_NPROCS=2 check_time.sh petab_import python tests/performance/test.py import - name: "Upload artifact: CS_Signalling_ERBB_RAS_AKT_petab" - uses: actions/upload-artifact@v1 + uses: actions/upload-artifact@v3 with: name: CS_Signalling_ERBB_RAS_AKT path: CS_Signalling_ERBB_RAS_AKT/CS_Signalling_ERBB_RAS_AKT_petab diff --git a/.github/workflows/test_petab_test_suite.yml b/.github/workflows/test_petab_test_suite.yml index 4594317582..98f5e32271 100644 --- a/.github/workflows/test_petab_test_suite.yml +++ b/.github/workflows/test_petab_test_suite.yml @@ -16,7 +16,7 @@ jobs: build: name: PEtab Testsuite - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 env: ENABLE_GCOV_COVERAGE: TRUE diff --git a/.github/workflows/test_pypi.yml b/.github/workflows/test_pypi.yml index d3c4341135..be090ec029 100644 --- a/.github/workflows/test_pypi.yml +++ b/.github/workflows/test_pypi.yml @@ -12,7 +12,7 @@ jobs: fail-fast: false matrix: python-version: [3.8, 3.9, '3.10'] - os: [ubuntu-20.04, macos-latest] + os: [ubuntu-22.04, macos-latest] runs-on: ${{ matrix.os }} diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index 69b09e9dde..aceb352d85 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -6,7 +6,7 @@ jobs: name: Tests Ubuntu # TODO: prepare image with more deps preinstalled - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 env: ENABLE_GCOV_COVERAGE: "TRUE" @@ -87,7 +87,7 @@ jobs: - name: Python tests run: | source build/venv/bin/activate \ - && pip3 install coverage pytest pytest-cov \ + && pip3 install coverage pytest pytest-cov pytest-rerunfailures \ && pytest \ --ignore-glob=*petab* \ --cov=amici \ diff --git a/.github/workflows/test_python_ver_matrix.yml b/.github/workflows/test_python_ver_matrix.yml index de5b6fb373..e896e79681 100644 --- a/.github/workflows/test_python_ver_matrix.yml +++ b/.github/workflows/test_python_ver_matrix.yml @@ -14,7 +14,7 @@ jobs: build: name: Python Version Matrix - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 continue-on-error: ${{ matrix.experimental }} env: AMICI_SKIP_CMAKE_TESTS: "TRUE" diff --git a/.github/workflows/test_sbml_semantic_test_suite.yml b/.github/workflows/test_sbml_semantic_test_suite.yml index 7901ef8640..31da1cdc9d 100644 --- a/.github/workflows/test_sbml_semantic_test_suite.yml +++ b/.github/workflows/test_sbml_semantic_test_suite.yml @@ -21,7 +21,7 @@ on: jobs: build: name: SBML Semantic Test Suite - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: fail-fast: false @@ -47,7 +47,7 @@ jobs: - run: AMICI_PARALLEL_COMPILE=2 ./scripts/run-SBMLTestsuite.sh ${{ matrix.cases }} - name: "Upload artifact: SBML semantic test suite results" - uses: actions/upload-artifact@v1 + uses: actions/upload-artifact@v3 with: name: amici-semantic-results path: tests/amici-semantic-results diff --git a/.github/workflows/test_valgrind.yml b/.github/workflows/test_valgrind.yml index 5e63d84ee9..ef0d592056 100644 --- a/.github/workflows/test_valgrind.yml +++ b/.github/workflows/test_valgrind.yml @@ -14,7 +14,7 @@ jobs: name: Tests Valgrind # TODO: prepare image with more deps preinstalled - runs-on: ubuntu-20.04 + runs-on: ubuntu-22.04 strategy: matrix: diff --git a/CHANGELOG.md b/CHANGELOG.md index 07e7c3c6f8..b3ff0593a2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,44 @@ ## v0.X Series +### v0.11.30 (2022-07-07) + +Features: +* Allow overriding model name during PySB import by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1801 +* Added __repr__ for parameter mapping classes by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1799 +* More informative warning messages for NaNs/Infs by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1798 +* Moved `sim_steps` increment by @plakrisenko in + https://github.com/AMICI-dev/AMICI/pull/1806 +* Re-arranged application of parameters from `ExpData` to avoid initial + sensitivities error by @dilpath in + https://github.com/AMICI-dev/AMICI/pull/1805 +* Checking for unused parameters in `simulate_petab` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1816 +* Add `create_parameter_mapping` kwarg forwarding by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1820 + +Other +* Remove `constant_species_to_parameters` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1809 + +Fixes +* Fixed handling of SBML models given as `pathlib.Path` in + `petab_import.import_model_sbml by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1808 +* Fixed missing CPU time reset by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1814 +* Raise in `simulate_petab` with `scaled_parameters=True` + `problem_parameters=None` by @dweindl in + https://github.com/AMICI-dev/AMICI/pull/1819 + +... + +**Full Changelog**: +https://github.com/AMICI-dev/AMICI/compare/v0.11.29...v0.11.30 + ### v0.11.29 (2022-05-06) ## What's Changed diff --git a/documentation/amici_refs.bib b/documentation/amici_refs.bib index 9eaf4c9db7..30a68b159e 100644 --- a/documentation/amici_refs.bib +++ b/documentation/amici_refs.bib @@ -1025,6 +1025,20 @@ @article {Froehlich2022.02.17.480899 journal = {bioRxiv} } +@Article{SluijsMaa2022, + author = {van Sluijs, Bob and Maas, Roel J. M. and van der Linden, Ardjan J. and de Greef, Tom F. A. and Huck, Wilhelm T. S.}, + journal = {Nature Communications}, + title = {A microfluidic optimal experimental design platform for forward design of cell-free genetic networks}, + year = {2022}, + issn = {2041-1723}, + number = {1}, + pages = {3626}, + volume = {13}, + abstract = {Cell-free protein synthesis has been widely used as a “breadboard” for design of synthetic genetic networks. However, due to a severe lack of modularity, forward engineering of genetic networks remains challenging. Here, we demonstrate how a combination of optimal experimental design and microfluidics allows us to devise dynamic cell-free gene expression experiments providing maximum information content for subsequent non-linear model identification. Importantly, we reveal that applying this methodology to a library of genetic circuits, that share common elements, further increases the information content of the data resulting in higher accuracy of model parameters. To show modularity of model parameters, we design a pulse decoder and bistable switch, and predict their behaviour both qualitatively and quantitatively. Finally, we update the parameter database and indicate that network topology affects parameter estimation accuracy. Utilizing our methodology provides us with more accurate model parameters, a necessity for forward engineering of complex genetic networks.}, + doi = {10.1038/s41467-022-31306-3}, + refid = {van Sluijs2022}, +} + @Comment{jabref-meta: databaseType:bibtex;} @Comment{jabref-meta: grouping: diff --git a/documentation/conf.py b/documentation/conf.py index f9dbc3fb94..2560b70cd5 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -236,7 +236,10 @@ def install_doxygen(): intersphinx_mapping = { 'pysb': ('https://pysb.readthedocs.io/en/stable/', None), - 'petab': ('https://petab.readthedocs.io/en/stable/', None), + 'petab': ( + 'https://petab.readthedocs.io/projects/libpetab-python/en/latest/', + None + ), 'pandas': ('https://pandas.pydata.org/docs/', None), 'numpy': ('https://numpy.org/devdocs/', None), 'sympy': ('https://docs.sympy.org/latest/', None), diff --git a/documentation/references.md b/documentation/references.md index e21c044115..6c786737e7 100644 --- a/documentation/references.md +++ b/documentation/references.md @@ -1,6 +1,6 @@ # References -List of publications using AMICI. Total number is 66. +List of publications using AMICI. Total number is 67. If you applied AMICI in your work and your publication is missing, please let us know via a new Github issue. @@ -12,6 +12,9 @@ If you applied AMICI in your work and your publication is missing, please let us

Schmucker, Robin, Gabriele Farina, James Faeder, Fabian Fröhlich, Ali Sinan Saglam, and Tuomas Sandholm. 2022. “Combination Treatment Optimization Using a Pan-Cancer Pathway Model.” PLOS Computational Biology 17 (12): 1–22. https://doi.org/10.1371/journal.pcbi.1009689.

+
+

Sluijs, Bob van, Roel J. M. Maas, Ardjan J. van der Linden, Tom F. A. de Greef, and Wilhelm T. S. Huck. 2022. “A Microfluidic Optimal Experimental Design Platform for Forward Design of Cell-Free Genetic Networks.” Nature Communications 13 (1): 3626. https://doi.org/10.1038/s41467-022-31306-3.

+

Stapor, Paul, Leonard Schmiester, Christoph Wierling, Simon Merkt, Dilan Pathirana, Bodo M. H. Lange, Daniel Weindl, and Jan Hasenauer. 2022. “Mini-batch optimization enables training of ODE models on large-scale datasets.” Nature Communications 13 (1): 34. https://doi.org/10.1038/s41467-021-27374-6.

diff --git a/include/amici/amici.h b/include/amici/amici.h index b3fe221de9..9e4a9f632f 100644 --- a/include/amici/amici.h +++ b/include/amici/amici.h @@ -96,16 +96,6 @@ class AmiciApplication { * @param ... arguments to be formatted */ void errorF(const char *identifier, const char *format, ...) const; - - /** - * @brief Checks the values in an array for NaNs and Infs - * - * @param array array - * @param fun name of calling function - * @return AMICI_RECOVERABLE_ERROR if a NaN/Inf value was found, - * AMICI_SUCCESS otherwise - */ - int checkFinite(gsl::span array, const char *fun); }; /** diff --git a/include/amici/model.h b/include/amici/model.h index 1897f9b295..31e1f8b25e 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -65,6 +65,28 @@ enum class ModelQuantity { qBdot_ss, xBdot_ss, JSparseB_ss, + deltax, + deltasx, + deltaxB, + k, + p, + ts, + dJydy, + dJydy_matlab, + deltaqB, + dsigmaydp, + dsigmaydy, + dJydsigma, + dJydx, + dzdx, + dzdp, + dJrzdsigma, + dJrzdz, + dJzdsigma, + dJzdz, + drzdp, + drzdx, + dsigmazdp, }; extern const std::map model_quantity_to_str; @@ -1869,7 +1891,7 @@ class Model : public AbstractModel, public ModelDimensions { /** vector of bools indicating whether state variables are to be assumed to * be positive */ std::vector state_is_non_negative_; - + /** Vector of booleans indicating the initial boolean value for every event trigger function. Events at t0 * can only trigger if the initial value is set to `false`. Must be specified during model compilation by * setting the `initialValue` attribute of an event trigger. */ diff --git a/python/amici/parameter_mapping.py b/python/amici/parameter_mapping.py index 930cd02a81..85969ed386 100644 --- a/python/amici/parameter_mapping.py +++ b/python/amici/parameter_mapping.py @@ -16,8 +16,10 @@ import numbers -from typing import Any, Dict, List, Union +import warnings +from typing import Any, Dict, List, Union, Set from collections.abc import Sequence +from itertools import chain import amici import numpy as np @@ -91,6 +93,27 @@ def __init__( scale_map_sim_fix = {key: LIN for key in map_sim_fix} self.scale_map_sim_fix = scale_map_sim_fix + def __repr__(self): + return (f"{self.__class__.__name__}(" + f"map_sim_var={repr(self.map_sim_var)}," + f"scale_map_sim_var={repr(self.scale_map_sim_var)}," + f"map_preeq_fix={repr(self.map_preeq_fix)}," + f"scale_map_preeq_fix={repr(self.scale_map_preeq_fix)}," + f"map_sim_fix={repr(self.map_sim_fix)}," + f"scale_map_sim_fix={repr(self.scale_map_sim_fix)})") + + @property + def free_symbols(self) -> Set[str]: + """Get IDs of all (symbolic) parameters present in this mapping""" + return { + p for p in chain( + self.map_sim_var.values(), + self.map_preeq_fix.values(), + self.map_sim_fix.values() + ) + if isinstance(p, str) + } + class ParameterMapping(Sequence): r"""Parameter mapping for multiple conditions. @@ -111,8 +134,7 @@ def __init__( self.parameter_mappings = parameter_mappings def __iter__(self): - for mapping in self.parameter_mappings: - yield mapping + yield from self.parameter_mappings def __getitem__(self, item): return self.parameter_mappings[item] @@ -122,17 +144,27 @@ def __len__(self): def append( self, - parameter_mapping_for_condition: ParameterMappingForCondition): + parameter_mapping_for_condition: ParameterMappingForCondition + ): """Append a condition specific parameter mapping.""" self.parameter_mappings.append(parameter_mapping_for_condition) + def __repr__(self): + return f"{self.__class__.__name__}({repr(self.parameter_mappings)})" + + @property + def free_symbols(self) -> Set[str]: + """Get IDs of all (symbolic) parameters present in this mapping""" + return set.union(*(mapping.free_symbols for mapping in self)) + def fill_in_parameters( edatas: List[amici.ExpData], problem_parameters: Dict[str, numbers.Number], scaled_parameters: bool, parameter_mapping: ParameterMapping, - amici_model: AmiciModel) -> None: + amici_model: AmiciModel +) -> None: """Fill fixed and dynamic parameters into the edatas (in-place). :param edatas: @@ -151,6 +183,11 @@ def fill_in_parameters( :param amici_model: AMICI model. """ + if unused_parameters := (set(problem_parameters.keys()) + - parameter_mapping.free_symbols): + warnings.warn("The following problem parameters were not used: " + + str(unused_parameters), RuntimeWarning) + for edata, mapping_for_condition in zip(edatas, parameter_mapping): fill_in_parameters_for_condition( edata, problem_parameters, scaled_parameters, @@ -355,7 +392,7 @@ def scale_parameters_dict( :param petab_scale_dict: Target scales of ``values`` """ - if not value_dict.keys() == petab_scale_dict.keys(): + if value_dict.keys() != petab_scale_dict.keys(): raise AssertionError("Keys don't match.") for key, value in value_dict.items(): @@ -378,7 +415,7 @@ def unscale_parameters_dict( :param petab_scale_dict: Target scales of ``values`` """ - if not value_dict.keys() == petab_scale_dict.keys(): + if value_dict.keys() != petab_scale_dict.keys(): raise AssertionError("Keys don't match.") for key, value in value_dict.items(): diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index bc520f9d20..0396e6515a 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -75,7 +75,7 @@ def _add_global_parameter(sbml_model: libsbml.Model, def get_fixed_parameters( sbml_model: 'libsbml.Model', condition_df: Optional[pd.DataFrame] = None, - const_species_to_parameters: bool = False) -> List[str]: +) -> List[str]: """ Determine, set and return fixed model parameters. @@ -89,11 +89,6 @@ def get_fixed_parameters( :param sbml_model: libsbml.Model instance - :param const_species_to_parameters: - If `True`, species which are marked constant within the SBML model - will be turned into constant parameters *within* the given - `sbml_model`. - :return: List of IDs of parameters which are to be considered constant. """ @@ -127,22 +122,6 @@ def get_fixed_parameters( else: fixed_parameters = [] - # Others are optional - if const_species_to_parameters: - # Turn species which are marked constant in the SBML model into - # parameters - constant_species = constant_species_to_parameters(sbml_model) - - logger.debug("Constant species converted to parameters: " - + str(len(constant_species))) - logger.info("Non-constant species " - + str(len(sbml_model.getListOfSpecies()))) - - # ... and append them to the list of fixed_parameters - for species in constant_species: - if species not in fixed_parameters: - fixed_parameters.append(species) - # Ensure mentioned parameters exist in the model. Remove additional ones # from list for fixed_parameter in fixed_parameters[:]: @@ -237,30 +216,6 @@ def species_to_parameters(species_ids: List[str], return transformables -def constant_species_to_parameters(sbml_model: 'libsbml.Model') -> List[str]: - """ - Convert constant species in the SBML model to constant parameters. - - This can be used e.g. for setting up models with condition-specific - constant species for PEtab, since there it is not possible to specify - constant species in the condition table. - - :param sbml_model: - SBML Model - - :return: - List of IDs of SBML species that have been turned into constants - """ - transformables = [] - for species in sbml_model.getListOfSpecies(): - if not species.getConstant() and not species.getBoundaryCondition(): - continue - - transformables.append(species.getId()) - - return species_to_parameters(transformables, sbml_model) - - def import_petab_problem( petab_problem: petab.Problem, model_output_dir: Union[str, Path, None] = None, @@ -280,7 +235,7 @@ def import_petab_problem( :param model_name: Name of the generated model. If model file name was provided, this defaults to the file name without extension, otherwise - the SBML model ID will be used. + the model ID will be used. :param force_compile: Whether to compile the model even if the target folder is not empty, @@ -303,12 +258,9 @@ def import_petab_problem( else: model_output_dir = os.path.abspath(model_output_dir) - if PysbPetabProblem and isinstance(petab_problem, PysbPetabProblem): - if model_name is None: - model_name = petab_problem.pysb_model.name - else: - raise ValueError( - "Argument model_name currently not allowed for pysb models") + if PysbPetabProblem and isinstance(petab_problem, PysbPetabProblem) \ + and model_name is None: + model_name = petab_problem.pysb_model.name elif model_name is None: model_name = _create_model_name(model_output_dir) @@ -333,6 +285,7 @@ def import_petab_problem( if PysbPetabProblem and isinstance(petab_problem, PysbPetabProblem): import_model_pysb( petab_problem, + model_name=model_name, model_output_dir=model_output_dir, **kwargs) else: @@ -453,7 +406,7 @@ def import_model_sbml( set_log_level(logger, verbose) - logger.info(f"Importing model ...") + logger.info("Importing model ...") # Get PEtab tables observable_df = petab.get_observable_df(observable_table) @@ -477,16 +430,19 @@ def import_model_sbml( logger.info(f"Model name is '{model_name}'.\n" f"Writing model code to '{model_output_dir}'.") + if isinstance(sbml_model, Path): + sbml_model = str(sbml_model) + # Load model if isinstance(sbml_model, str): # from file sbml_reader = libsbml.SBMLReader() sbml_doc = sbml_reader.readSBMLFromFile(sbml_model) - sbml_model = sbml_doc.getModel() else: # Create a copy, because it will be modified by SbmlImporter sbml_doc = sbml_model.getSBMLDocument().clone() - sbml_model = sbml_doc.getModel() + + sbml_model = sbml_doc.getModel() show_model_info(sbml_model) @@ -511,6 +467,8 @@ def import_model_sbml( if observable_df is not None: observables, noise_distrs, sigmas = \ get_observation_model(observable_df) + else: + observables = noise_distrs = sigmas = None logger.info(f'Observables: {len(observables)}') logger.info(f'Sigmas: {len(sigmas)}') diff --git a/python/amici/petab_import_pysb.py b/python/amici/petab_import_pysb.py index 170a9ed45a..90c12e88d8 100644 --- a/python/amici/petab_import_pysb.py +++ b/python/amici/petab_import_pysb.py @@ -314,6 +314,7 @@ def import_model_pysb( petab_problem: PysbPetabProblem, model_output_dir: Optional[Union[str, Path]] = None, verbose: Optional[Union[bool, int]] = True, + model_name: Optional[str] = None, **kwargs ) -> None: """ @@ -329,11 +330,13 @@ def import_model_pysb( :param verbose: Print/log extra information. + :param model_name: + Name of the generated model module + :param kwargs: Additional keyword arguments to be passed to :meth:`amici.pysb_import.pysb2amici`. """ - set_log_level(logger, verbose) logger.info("Importing model ...") @@ -374,7 +377,10 @@ def import_model_pysb( observable_table) from amici.pysb_import import pysb2amici - pysb2amici(pysb_model, model_output_dir, verbose=True, + pysb2amici(model=pysb_model, + output_dir=model_output_dir, + model_name=model_name, + verbose=True, observables=observables, sigmas=sigmas, constant_parameters=constant_parameters, diff --git a/python/amici/petab_objective.py b/python/amici/petab_objective.py index 0cb8b4d9d2..95c25a87f0 100644 --- a/python/amici/petab_objective.py +++ b/python/amici/petab_objective.py @@ -66,18 +66,18 @@ def simulate_petab( will be used). To be provided as dict, mapping PEtab problem parameters to SBML IDs. :param simulation_conditions: - Result of `petab.get_simulation_conditions`. Can be provided to save - time if this has be obtained before. - Not required if `edatas` and `parameter_mapping` are provided. + Result of :py:func:`petab.get_simulation_conditions`. Can be provided + to save time if this has be obtained before. + Not required if ``edatas`` and ``parameter_mapping`` are provided. :param edatas: Experimental data. Parameters are inserted in-place for simulation. :param parameter_mapping: Optional precomputed PEtab parameter mapping for efficiency, as - generated by `create_parameter_mapping`. + generated by :py:func:`create_parameter_mapping`. :param scaled_parameters: - If True, problem_parameters are assumed to be on the scale provided - in the PEtab parameter table and will be unscaled. If False, they - are assumed to be in linear scale. + If ``True``, ``problem_parameters`` are assumed to be on the scale + provided in the PEtab parameter table and will be unscaled. + If ``False``, they are assumed to be in linear scale. :param log_level: Log level, see :mod:`amici.logging` module. :param num_threads: @@ -90,10 +90,8 @@ def simulate_petab( :return: Dictionary of - * cost function value (LLH), - * const function sensitivity w.r.t. parameters (SLLH), - (**NOTE**: Sensitivities are computed for the scaled parameters) - * list of `ReturnData` (RDATAS), + * cost function value (``LLH``), + * list of :class:`amici.amici.ReturnData` (``RDATAS``), corresponding to the different simulation conditions. For ordering of simulation conditions, see @@ -109,7 +107,10 @@ def simulate_petab( # Use PEtab nominal values as default problem_parameters = {t.Index: getattr(t, NOMINAL_VALUE) for t in petab_problem.parameter_df.itertuples()} - scaled_parameters = False + if scaled_parameters: + raise NotImplementedError( + "scaled_parameters=True in combination with " + "problem_parameters=None is currently not supported.") # number of amici simulations will be number of unique # (preequilibrationConditionId, simulationConditionId) pairs. @@ -185,15 +186,15 @@ def create_parameterized_edatas( will be used). To be provided as dict, mapping PEtab problem parameters to SBML IDs. :param scaled_parameters: - If True, problem_parameters are assumed to be on the scale provided - in the PEtab parameter table and will be unscaled. If False, they - are assumed to be in linear scale. + If ``True``, ``problem_parameters`` are assumed to be on the scale + provided in the PEtab parameter table and will be unscaled. + If ``False``, they are assumed to be in linear scale. :param parameter_mapping: Optional precomputed PEtab parameter mapping for efficiency, as - generated by `create_parameter_mapping`. + generated by :func:`create_parameter_mapping`. :param simulation_conditions: - Result of `petab.get_simulation_conditions`. Can be provided to save - time if this has been obtained before. + Result of :func:`petab.get_simulation_conditions`. Can be provided to + save time if this has been obtained before. :return: List with one :class:`amici.amici.ExpData` per simulation condition, @@ -236,21 +237,26 @@ def create_parameter_mapping( simulation_conditions: Union[pd.DataFrame, List[Dict]], scaled_parameters: bool, amici_model: AmiciModel, + **parameter_mapping_kwargs, ) -> ParameterMapping: """Generate AMICI specific parameter mapping. :param petab_problem: PEtab problem :param simulation_conditions: - Result of `petab.get_simulation_conditions`. Can be provided to save - time if this has been obtained before. + Result of :func:`petab.get_simulation_conditions`. Can be provided to + save time if this has been obtained before. :param scaled_parameters: - If True, problem_parameters are assumed to be on the scale provided - in the PEtab parameter table and will be unscaled. If False, they + If ``True``, problem_parameters are assumed to be on the scale provided + in the PEtab parameter table and will be unscaled. If ``False``, they are assumed to be in linear scale. :param amici_model: AMICI model. - + :param parameter_mapping_kwargs: + Optional keyword arguments passed to + :func:`petab.get_optimization_to_simulation_parameter_mapping`. + To allow changing fixed PEtab problem parameters (``estimate=0``), + use ``fill_fixed_parameters=False``. :return: List of the parameter mappings. """ @@ -272,13 +278,25 @@ def create_parameter_mapping( "SBML LocalParameters. If the model contains " "LocalParameters, parameter mapping will fail.") - prelim_parameter_mapping = \ - petab_problem.get_optimization_to_simulation_parameter_mapping( - warn_unmapped=False, scaled_parameters=scaled_parameters, - allow_timepoint_specific_numeric_noise_parameters= + default_parameter_mapping_kwargs = { + "warn_unmapped": False, + "scaled_parameters": scaled_parameters, + "allow_timepoint_specific_numeric_noise_parameters": not petab.lint.observable_table_has_nontrivial_noise_formula( - petab_problem.observable_df - ) + petab_problem.observable_df), + } + if parameter_mapping_kwargs is None: + parameter_mapping_kwargs = {} + + prelim_parameter_mapping = \ + petab.get_optimization_to_simulation_parameter_mapping( + condition_df=petab_problem.condition_df, + measurement_df=petab_problem.measurement_df, + parameter_df=petab_problem.parameter_df, + observable_df=petab_problem.observable_df, + sbml_model=petab_problem.sbml_model, + **dict(default_parameter_mapping_kwargs, + **parameter_mapping_kwargs) ) parameter_mapping = ParameterMapping() @@ -303,8 +321,8 @@ def create_parameter_mapping_for_condition( :param parameter_mapping_for_condition: Preliminary parameter mapping for condition. :param condition: - pandas.DataFrame row with preequilibrationConditionId and - simulationConditionId. + :class:`pandas.DataFrame` row with ``preequilibrationConditionId`` and + ``simulationConditionId``. :param petab_problem: Underlying PEtab problem. :param amici_model: @@ -493,8 +511,8 @@ def create_edatas( :param petab_problem: Underlying PEtab problem. :param simulation_conditions: - Result of `petab.get_simulation_conditions`. Can be provided to save - time if this has be obtained before. + Result of :func:`petab.get_simulation_conditions`. Can be provided to + save time if this has be obtained before. :return: List with one :class:`amici.amici.ExpData` per simulation condition, @@ -547,10 +565,10 @@ def create_edata_for_condition( Sets timepoints, observed data and sigmas. :param condition: - pandas.DataFrame row with preequilibrationConditionId and - simulationConditionId. + :class:`pandas.DataFrame` row with ``preequilibrationConditionId`` and + ``simulationConditionId``. :param measurement_df: - pandas.DataFrame with measurements for the given condition. + :class:`pandas.DataFrame` with measurements for the given condition. :param amici_model: AMICI model :param petab_problem: @@ -721,20 +739,20 @@ def rdatas_to_measurement_df( measurement_df: pd.DataFrame) -> pd.DataFrame: """ Create a measurement dataframe in the PEtab format from the passed - `rdatas` and own information. + ``rdatas`` and own information. :param rdatas: A sequence of rdatas with the ordering of - `petab.get_simulation_conditions`. + :func:`petab.get_simulation_conditions`. :param model: - AMICI model used to generate `rdatas`. + AMICI model used to generate ``rdatas``. :param measurement_df: - PEtab measurement table used to generate `rdatas`. + PEtab measurement table used to generate ``rdatas``. :return: - A dataframe built from the rdatas in the format of `measurement_df`. + A dataframe built from the rdatas in the format of ``measurement_df``. """ simulation_conditions = petab.get_simulation_conditions( measurement_df) @@ -778,10 +796,11 @@ def rdatas_to_simulation_df( rdatas: Sequence[amici.ReturnData], model: AmiciModel, measurement_df: pd.DataFrame) -> pd.DataFrame: - """Create a PEtab simulation dataframe from amici.ReturnDatas. + """Create a PEtab simulation dataframe from + :class:`amici.amici.ReturnData` s. - See ``rdatas_to_measurement_df`` for details, only that model outputs - will appear in column "simulation" instead of "measurement".""" + See :func:`rdatas_to_measurement_df` for details, only that model outputs + will appear in column ``simulation`` instead of ``measurement``.""" df = rdatas_to_measurement_df(rdatas=rdatas, model=model, measurement_df=measurement_df) diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py index e2185f0c10..bc7ac80bd3 100644 --- a/python/amici/pysb_import.py +++ b/python/amici/pysb_import.py @@ -51,6 +51,7 @@ def pysb2amici( # See https://github.com/AMICI-dev/AMICI/pull/1672 cache_simplify: bool = False, generate_sensitivity_code: bool = True, + model_name: Optional[str] = None, ): r""" Generate AMICI C++ files for the provided model. @@ -124,6 +125,10 @@ def pysb2amici( :param generate_sensitivity_code: if set to ``False``, code for sensitivity computation will not be generated + + :param model_name: + Name for the generated model module. If None, :attr:`pysb.Model.name` + will be used. """ if observables is None: observables = [] @@ -133,6 +138,8 @@ def pysb2amici( if sigmas is None: sigmas = {} + model_name = model_name or model.name + set_log_level(logger, verbose) ode_model = ode_model_from_pysb_importer( model, constant_parameters=constant_parameters, @@ -146,7 +153,7 @@ def pysb2amici( exporter = ODEExporter( ode_model, outdir=output_dir, - model_name=model.name, + model_name=model_name, verbose=verbose, assume_pow_positivity=assume_pow_positivity, compiler=compiler, diff --git a/python/tests/test_edata.py b/python/tests/test_edata.py new file mode 100644 index 0000000000..184f5e2ce1 --- /dev/null +++ b/python/tests/test_edata.py @@ -0,0 +1,47 @@ +"""Tests related to amici.ExpData via Python""" +import numpy as np + +import amici + +from test_sbml_import import model_units_module + + +def test_edata_sensi_unscaling(model_units_module): + """ + ExpData parameters should be used for unscaling initial state + sensitivities. + """ + parameters0 = (5, 5) + parameters1 = (2, 2) + + sx0 = (3, 3, 3, 3) + + parameter_scales_log10 = \ + [amici.ParameterScaling.log10.value]*len(parameters0) + amici_parameter_scales_log10 = \ + amici.parameterScalingFromIntVector(parameter_scales_log10) + + model = model_units_module.getModel() + model.setTimepoints(np.linspace(0, 1, 3)) + model.setParameterScale(parameter_scales_log10) + model.setParameters(parameters0) + + solver = model.getSolver() + solver.setSensitivityOrder(amici.SensitivityOrder.first) + + edata0 = amici.ExpData(model) + edata0.pscale = amici_parameter_scales_log10 + edata0.parameters = parameters0 + edata0.sx0 = sx0 + + edata1 = amici.ExpData(model) + edata1.pscale = amici_parameter_scales_log10 + edata1.parameters = parameters1 + edata1.sx0 = sx0 + + rdata0 = amici.runAmiciSimulation(model, solver, edata0) + rdata1 = amici.runAmiciSimulation(model, solver, edata1) + + # The initial state sensitivities are as specified. + assert np.isclose(rdata0.sx0.flatten(), sx0).all() + assert np.isclose(rdata1.sx0.flatten(), sx0).all() diff --git a/python/tests/test_petab_import.py b/python/tests/test_petab_import.py index 5de3f182df..796bbd6232 100644 --- a/python/tests/test_petab_import.py +++ b/python/tests/test_petab_import.py @@ -38,19 +38,3 @@ def simple_sbml_model(): species_ref.setSpecies(name) return document, model - - -def test_constant_species_to_parameters(simple_sbml_model): - """Test conversion from species to constant parameters""" - - document, model = simple_sbml_model - - amici_petab_import.constant_species_to_parameters(model) - - assert len(list(model.getListOfParameters())) == 1 - assert len(list(model.getListOfSpecies())) == 0 - - r = model.getReaction(0) - assert len(list(r.getListOfReactants())) == 0 - assert len(list(r.getListOfProducts())) == 0 - assert len(list(r.getListOfModifiers())) == 0 diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 64bf2610f1..5c20eb709b 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -2,6 +2,7 @@ import os import re import shutil +from pathlib import Path from urllib.request import urlopen import libsbml @@ -182,23 +183,20 @@ def model_steadystate_module(): shutil.rmtree(outdir, ignore_errors=True) -@pytest.fixture +@pytest.fixture(scope='session') def model_units_module(): - sbml_file = os.path.join(os.path.dirname(__file__), '..', - 'examples', 'example_units', - 'model_units.xml') - sbml_importer = amici.SbmlImporter(sbml_file) - - outdir = 'test_model_units' + sbml_file = Path(__file__).parent / '..' / 'examples' \ + / 'example_units' / 'model_units.xml' module_name = 'test_model_units' - sbml_importer.sbml2amici(model_name=module_name, - output_dir=outdir) - model_module = amici.import_model_module(module_name=module_name, - module_path=outdir) - yield model_module + sbml_importer = amici.SbmlImporter(sbml_file) - shutil.rmtree(outdir, ignore_errors=True) + with TemporaryDirectory() as outdir: + sbml_importer.sbml2amici(model_name=module_name, + output_dir=outdir) + + yield amici.import_model_module(module_name=module_name, + module_path=outdir) def test_presimulation(sbml_example_presimulation_module): diff --git a/python/tests/test_sbml_import_special_functions.py b/python/tests/test_sbml_import_special_functions.py index 67e1368f30..7b7cb7a8f7 100644 --- a/python/tests/test_sbml_import_special_functions.py +++ b/python/tests/test_sbml_import_special_functions.py @@ -50,6 +50,8 @@ def model_special_likelihoods(): shutil.rmtree(outdir, ignore_errors=True) +# FD check fails occasionally, so give some extra tries +@pytest.mark.flaky(reruns=5) def test_special_likelihoods(model_special_likelihoods): """Test special likelihood functions.""" model = model_special_likelihoods.getModel() @@ -65,7 +67,7 @@ def test_special_likelihoods(model_special_likelihoods): # make sure measurements are smaller for non-degenerate probability y = edata.getObservedData() - y = tuple([val * np.random.uniform(0, 1) for val in y]) + y = tuple(val * np.random.uniform(0, 1) for val in y) edata.setObservedData(y) # set sigmas @@ -100,7 +102,7 @@ def test_special_likelihoods(model_special_likelihoods): solver.setSensitivityMethod(sensi_method) solver.setSensitivityOrder(amici.SensitivityOrder.first) check_derivatives( - model, solver, edata, atol=1e-1, rtol=1e-1, + model, solver, edata, atol=1e-4, rtol=1e-3, check_least_squares=False) # Test for m > y, i.e. in region with 0 density @@ -110,7 +112,7 @@ def test_special_likelihoods(model_special_likelihoods): # make sure measurements are smaller for non-degenerate probability y = edata.getObservedData() - y = tuple([val * np.random.uniform(0.5, 3) for val in y]) + y = tuple(val * np.random.uniform(0.5, 3) for val in y) edata.setObservedData(y) edata.setObservedDataStdDev(sigmas) diff --git a/python/tests/util.py b/python/tests/util.py index cac811b115..dfc867f3b8 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -11,6 +11,7 @@ SensitivityMethod, SensitivityOrder ) +from amici.gradient_check import _check_close def create_amici_model(sbml_model, model_name) -> AmiciModel: @@ -131,14 +132,16 @@ def check_trajectories_without_sensitivities( solver = amici_model.getSolver() solver.setAbsoluteTolerance(1e-15) rdata = runAmiciSimulation(amici_model, solver=solver) - np.testing.assert_almost_equal(rdata['x'], result_expected_x, decimal=5) + _check_close(rdata['x'], result_expected_x, field="x", + rtol=5e-5, atol=1e-13) # Show that we can do arbitrary precision here (test 8 digits) solver = amici_model.getSolver() solver.setAbsoluteTolerance(1e-15) solver.setRelativeTolerance(1e-12) rdata = runAmiciSimulation(amici_model, solver=solver) - np.testing.assert_almost_equal(rdata['x'], result_expected_x, decimal=8) + _check_close(rdata['x'], result_expected_x, field="x", + rtol=5e-9, atol=1e-13) def check_trajectories_with_forward_sensitivities( @@ -157,8 +160,10 @@ def check_trajectories_with_forward_sensitivities( solver.setSensitivityOrder(SensitivityOrder.first) solver.setSensitivityMethod(SensitivityMethod.forward) rdata = runAmiciSimulation(amici_model, solver=solver) - np.testing.assert_almost_equal(rdata['x'], result_expected_x, decimal=5) - np.testing.assert_almost_equal(rdata['sx'], result_expected_sx, decimal=5) + _check_close(rdata['x'], result_expected_x, field="x", + rtol=1e-5, atol=1e-13) + _check_close(rdata['sx'], result_expected_sx, field="sx", + rtol=1e-5, atol=1e-7) # Show that we can do arbitrary precision here (test 8 digits) solver = amici_model.getSolver() @@ -169,5 +174,7 @@ def check_trajectories_with_forward_sensitivities( solver.setAbsoluteToleranceFSA(1e-15) solver.setRelativeToleranceFSA(1e-13) rdata = runAmiciSimulation(amici_model, solver=solver) - np.testing.assert_almost_equal(rdata['x'], result_expected_x, decimal=8) - np.testing.assert_almost_equal(rdata['sx'], result_expected_sx, decimal=8) + _check_close(rdata['x'], result_expected_x, field="x", + rtol=1e-10, atol=1e-12) + _check_close(rdata['sx'], result_expected_sx, field="sx", + rtol=1e-10, atol=1e-9) diff --git a/src/amici.cpp b/src/amici.cpp index cc20bd5927..bce248e94f 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -238,6 +238,14 @@ AmiciApplication::runAmiciSimulation(Solver& solver, rdata->cpu_time_total = static_cast(clock() - start_time_total) * 1000.0 / CLOCKS_PER_SEC; + // verify that reported CPU times are plausible + gsl_EnsuresDebug(rdata->cpu_time <= rdata->cpu_time_total); + gsl_EnsuresDebug(rdata->cpu_timeB <= rdata->cpu_time_total); + gsl_EnsuresDebug(rdata->preeq_cpu_time <= rdata->cpu_time_total); + gsl_EnsuresDebug(rdata->preeq_cpu_timeB <= rdata->cpu_time_total); + gsl_EnsuresDebug(rdata->posteq_cpu_time <= rdata->cpu_time_total); + gsl_EnsuresDebug(rdata->posteq_cpu_timeB <= rdata->cpu_time_total); + return rdata; } @@ -302,27 +310,4 @@ AmiciApplication::errorF(const char* identifier, const char* format, ...) const error(identifier, str); } -int -AmiciApplication::checkFinite(gsl::span array, const char* fun) -{ - Expects(array.size() - <= static_cast(std::numeric_limits::max())); - - for (int idx = 0; idx < (int)array.size(); idx++) { - if (isNaN(array[idx])) { - warningF("AMICI:NaN", - "AMICI encountered a NaN value for %s[%i]!", - fun, idx); - return AMICI_RECOVERABLE_ERROR; - } - if (isInf(array[idx])) { - warningF("AMICI:Inf", - "AMICI encountered an Inf value for %s[%i]!", - fun, idx); - return AMICI_RECOVERABLE_ERROR; - } - } - return AMICI_SUCCESS; -} - } // namespace amici diff --git a/src/edata.cpp b/src/edata.cpp index cc37a36455..e9ced1b018 100644 --- a/src/edata.cpp +++ b/src/edata.cpp @@ -380,6 +380,16 @@ void ConditionContext::applyCondition(const ExpData *edata, model_->setParameterScale(edata->pscale); } + // this needs to be set in the model before handling initial state + // sensitivities, which may be unscaled using model parameter values + if(!edata->parameters.empty()) { + if(edata->parameters.size() != (unsigned) model_->np()) + throw AmiException("Number of parameters (%d) in model does not" + " match ExpData (%zd).", + model_->np(), edata->parameters.size()); + model_->setParameters(edata->parameters); + } + if(!edata->x0.empty()) { if(edata->x0.size() != (unsigned) model_->nx_rdata) throw AmiException("Number of initial conditions (%d) in model does" @@ -397,14 +407,6 @@ void ConditionContext::applyCondition(const ExpData *edata, model_->setInitialStateSensitivities(edata->sx0); } - if(!edata->parameters.empty()) { - if(edata->parameters.size() != (unsigned) model_->np()) - throw AmiException("Number of parameters (%d) in model does not" - " match ExpData (%zd).", - model_->np(), edata->parameters.size()); - model_->setParameters(edata->parameters); - } - model_->setReinitializeFixedParameterInitialStates( edata->reinitializeFixedParameterInitialStates); diff --git a/src/model.cpp b/src/model.cpp index e6097fa88f..6e914bb320 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -46,6 +46,29 @@ const std::map model_quantity_to_str { {ModelQuantity::qBdot_ss, "qBdot_ss"}, {ModelQuantity::xBdot_ss, "xBdot_ss"}, {ModelQuantity::JSparseB_ss, "JSparseB_ss"}, + {ModelQuantity::deltax, "deltax"}, + {ModelQuantity::deltasx, "deltasx"}, + {ModelQuantity::deltaxB, "deltaxB"}, + {ModelQuantity::k, "k"}, + {ModelQuantity::p, "p"}, + {ModelQuantity::ts, "ts"}, + {ModelQuantity::dJydy, "dJydy"}, + {ModelQuantity::dJydy_matlab, "dJydy"}, + {ModelQuantity::deltaqB, "deltaqB"}, + {ModelQuantity::dsigmaydp, "dsigmaydp"}, + {ModelQuantity::dsigmaydy, "dsigmaydy"}, + {ModelQuantity::dJydsigma, "dJydsigma"}, + {ModelQuantity::dJydx, "dJydx"}, + {ModelQuantity::dzdp, "dzdp"}, + {ModelQuantity::dzdx, "dzdx"}, + {ModelQuantity::dJrzdsigma, "dJrzdsigma"}, + {ModelQuantity::dJrzdz, "dJrzdz"}, + {ModelQuantity::dJzdsigma, "dJzdsigma"}, + {ModelQuantity::dJzdz, "dJzdz"}, + {ModelQuantity::drzdp, "drzdp"}, + {ModelQuantity::drzdx, "drzdx"}, + {ModelQuantity::dsigmazdp, "dsigmazdp"}, + }; /** @@ -154,7 +177,7 @@ Model::Model(ModelDimensions const & model_dimensions, simulation_parameters_.pscale, state_.unscaledParameters); state_.fixedParameters = simulation_parameters_.fixedParameters; state_.plist = simulation_parameters_.plist; - + root_initial_values_.resize(ne, true); /* If Matlab wrapped: dxdotdp is a full AmiVector, @@ -1166,7 +1189,7 @@ void Model::addStateEventUpdate(AmiVector &x, const int ie, const realtype t, state_.h.data(), ie, xdot.data(), xdot_old.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.deltax_, "deltax"); + checkFinite(derived_state_.deltax_, ModelQuantity::deltax); } // update @@ -1195,7 +1218,8 @@ void Model::addStateSensitivityEventUpdate(AmiVectorArray &sx, const int ie, state_.total_cl.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.deltasx_, "deltasx"); + checkFinite(derived_state_.deltasx_, ModelQuantity::deltasx, + nplist()); } amici_daxpy(nx_solver, 1.0, derived_state_.deltasx_.data(), 1, @@ -1217,7 +1241,7 @@ void Model::addAdjointStateEventUpdate(AmiVector &xB, const int ie, xdot_old.data(), xB.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.deltaxB_, "deltaxB"); + checkFinite(derived_state_.deltaxB_, ModelQuantity::deltaxB); } // apply update @@ -1243,7 +1267,7 @@ void Model::addAdjointQuadratureEventUpdate( } if (always_check_finite_) { - app->checkFinite(derived_state_.deltaqB_, "deltaqB"); + checkFinite(derived_state_.deltaqB_, ModelQuantity::deltaqB, nplist()); } } @@ -1292,6 +1316,8 @@ int Model::checkFinite(gsl::span array, case ModelQuantity::Jv: case ModelQuantity::JvB: case ModelQuantity::JDiag: + case ModelQuantity::deltax: + case ModelQuantity::deltaxB: if(hasStateIds()) { element_id = getStateIdsSolver()[flat_index]; } @@ -1306,6 +1332,16 @@ int Model::checkFinite(gsl::span array, element_id = getExpressionIds()[flat_index]; } break; + case ModelQuantity::k: + if(hasFixedParameterIds()) { + element_id = getFixedParameterIds()[flat_index]; + } + break; + case ModelQuantity::p: + if(hasParameterIds()) { + element_id = getParameterIds()[flat_index]; + } + break; default: break; } @@ -1327,15 +1363,19 @@ int Model::checkFinite(gsl::span array, element_id.c_str() ); - // check upstream - app->checkFinite(state_.fixedParameters, "k"); - app->checkFinite(state_.unscaledParameters, "p"); - app->checkFinite(simulation_parameters_.ts_, "t"); - if(!always_check_finite_) { - // don't check twice if always_check_finite_ is true - app->checkFinite(derived_state_.w_, "w"); + // check upstream, without infinite recursion + if(model_quantity != ModelQuantity::k + && model_quantity != ModelQuantity::p + && model_quantity != ModelQuantity::ts) + { + checkFinite(state_.fixedParameters, ModelQuantity::k); + checkFinite(state_.unscaledParameters, ModelQuantity::p); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); + if(!always_check_finite_ && model_quantity != ModelQuantity::w) { + // don't check twice if always_check_finite_ is true + checkFinite(derived_state_.w_, ModelQuantity::w); + } } - return AMICI_RECOVERABLE_ERROR; } @@ -1368,12 +1408,49 @@ int Model::checkFinite(gsl::span array, case ModelQuantity::sy: case ModelQuantity::ssigmay: case ModelQuantity::dydp: - row_id += " " + getObservableIds()[row]; - col_id += " " + getParameterIds()[plist(col)]; + case ModelQuantity::dsigmaydp: + if(hasObservableIds()) + row_id += " " + getObservableIds()[row]; + if(hasParameterIds()) + col_id += " " + getParameterIds()[plist(col)]; break; case ModelQuantity::dydx: - row_id += " " + getObservableIds()[row]; - col_id += " " + getStateIdsSolver()[col]; + if(hasObservableIds()) + row_id += " " + getObservableIds()[row]; + if(hasStateIds()) + col_id += " " + getStateIdsSolver()[col]; + break; + case ModelQuantity::deltasx: + if(hasStateIds()) + row_id += " " + getStateIdsSolver()[row]; + if(hasParameterIds()) + col_id += " " + getParameterIds()[plist(col)]; + break; + case ModelQuantity::dJydy: + case ModelQuantity::dJydy_matlab: + case ModelQuantity::dJydsigma: + if(hasObservableIds()) + col_id += " " + getObservableIds()[col]; + break; + case ModelQuantity::dJydx: + case ModelQuantity::dzdx: + case ModelQuantity::drzdx: + if(hasStateIds()) + col_id += " " + getStateIdsSolver()[col]; + break; + case ModelQuantity::deltaqB: + case ModelQuantity::dzdp: + case ModelQuantity::drzdp: + case ModelQuantity::dsigmazdp: + if(hasParameterIds()) + col_id += " " + getParameterIds()[plist(col)]; + break; + case ModelQuantity::dsigmaydy: + if(hasObservableIds()) { + auto obs_ids = getObservableIds(); + row_id += " " + obs_ids[row]; + col_id += " " + obs_ids[col]; + } break; default: break; @@ -1399,10 +1476,10 @@ int Model::checkFinite(gsl::span array, ); // check upstream - app->checkFinite(state_.fixedParameters, "k"); - app->checkFinite(state_.unscaledParameters, "p"); - app->checkFinite(simulation_parameters_.ts_, "t"); - app->checkFinite(derived_state_.w_, "w"); + checkFinite(state_.fixedParameters, ModelQuantity::k); + checkFinite(state_.unscaledParameters, ModelQuantity::p); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); + checkFinite(derived_state_.w_, ModelQuantity::w); return AMICI_RECOVERABLE_ERROR; } @@ -1489,10 +1566,10 @@ int Model::checkFinite(SUNMatrix m, ModelQuantity model_quantity, realtype t) co ); // check upstream - app->checkFinite(state_.fixedParameters, "k"); - app->checkFinite(state_.unscaledParameters, "p"); - app->checkFinite(simulation_parameters_.ts_, "t"); - app->checkFinite(derived_state_.w_, "w"); + checkFinite(state_.fixedParameters, ModelQuantity::k); + checkFinite(state_.unscaledParameters, ModelQuantity::p); + checkFinite(simulation_parameters_.ts_, ModelQuantity::ts); + checkFinite(derived_state_.w_, ModelQuantity::w); return AMICI_RECOVERABLE_ERROR; } @@ -1789,7 +1866,8 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) { } if (always_check_finite_) { - app->checkFinite(derived_state_.dsigmaydp_, "dsigmaydp"); + checkFinite(derived_state_.dsigmaydp_, ModelQuantity::dsigmaydp, + nplist()); } } @@ -1818,7 +1896,7 @@ void Model::fdsigmaydy(const int it, const ExpData *edata) { } if (always_check_finite_) { - app->checkFinite(derived_state_.dsigmaydy_, "dsigmaydy"); + checkFinite(derived_state_.dsigmaydy_, ModelQuantity::dsigmaydy, ny); } } @@ -1871,9 +1949,9 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { derived_state_.dJydy_.at(iyt).refresh(); if (always_check_finite_) { - app->checkFinite( - gsl::make_span(derived_state_.dJydy_.at(iyt).get()), - "dJydy"); + checkFinite( + gsl::make_span(derived_state_.dJydy_.at(iyt).get()), + ModelQuantity::dJydy, ny); } } } else { @@ -1887,10 +1965,13 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { state_.fixedParameters.data(), derived_state_.y_.data(), derived_state_.sigmay_.data(), edata.getObservedDataPtr(it)); - } - if (always_check_finite_) { - // get dJydy slice (ny) for current timepoint and observable - app->checkFinite(derived_state_.dJydy_matlab_, "dJydy"); + if (always_check_finite_) { + // get dJydy slice (ny) for current timepoint and observable + checkFinite( + gsl::span( + &derived_state_.dJydy_matlab_[iyt * ny * nJ], ny * nJ), + ModelQuantity::dJydy, ny); + } } } } @@ -1905,17 +1986,20 @@ void Model::fdJydsigma(const int it, const AmiVector &x, const ExpData &edata) { fsigmay(it, &edata); for (int iyt = 0; iyt < nytrue; iyt++) { - if (edata.isSetObservedData(it, iyt)) + if (edata.isSetObservedData(it, iyt)) { // get dJydsigma slice (ny) for current timepoint and observable fdJydsigma(&derived_state_.dJydsigma_.at(iyt * ny * nJ), iyt, state_.unscaledParameters.data(), state_.fixedParameters.data(), derived_state_.y_.data(), derived_state_.sigmay_.data(), edata.getObservedDataPtr(it)); - } - - if (always_check_finite_) { - app->checkFinite(derived_state_.dJydsigma_, "dJydsigma"); + if (always_check_finite_) { + checkFinite( + gsl::span( + &derived_state_.dJydsigma_.at(iyt * ny * nJ), ny * nJ), + ModelQuantity::dJydsigma, ny); + } + } } } @@ -1997,7 +2081,7 @@ void Model::fdJydx(const int it, const AmiVector &x, const ExpData &edata) { } if (always_check_finite_) { - app->checkFinite(derived_state_.dJydx_, "dJydx"); + checkFinite(derived_state_.dJydx_, ModelQuantity::dJydx, nx_solver); } } @@ -2023,7 +2107,7 @@ void Model::fdzdp(const int ie, const realtype t, const AmiVector &x) { } if (always_check_finite_) { - app->checkFinite(derived_state_.dzdp_, "dzdp"); + checkFinite(derived_state_.dzdp_, ModelQuantity::dzdp, nplist()); } } @@ -2038,7 +2122,7 @@ void Model::fdzdx(const int ie, const realtype t, const AmiVector &x) { state_.h.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.dzdx_, "dzdx"); + checkFinite(derived_state_.dzdx_, ModelQuantity::dzdx, nx_solver); } } @@ -2064,7 +2148,7 @@ void Model::fdrzdp(const int ie, const realtype t, const AmiVector &x) { } if (always_check_finite_) { - app->checkFinite(derived_state_.drzdp_, "drzdp"); + checkFinite(derived_state_.drzdp_, ModelQuantity::drzdp, nplist()); } } @@ -2079,7 +2163,7 @@ void Model::fdrzdx(const int ie, const realtype t, const AmiVector &x) { state_.h.data()); if (always_check_finite_) { - app->checkFinite(derived_state_.drzdx_, "drzdx"); + checkFinite(derived_state_.drzdx_, ModelQuantity::drzdx, nx_solver); } } @@ -2141,7 +2225,8 @@ void Model::fdsigmazdp(const int ie, const int nroots, const realtype t, } if (always_check_finite_) { - app->checkFinite(derived_state_.dsigmazdp_, "dsigmazdp"); + checkFinite(derived_state_.dsigmazdp_, ModelQuantity::dsigmazdp, + nplist()); } } @@ -2162,12 +2247,15 @@ void Model::fdJzdz(const int ie, const int nroots, const realtype t, state_.fixedParameters.data(), derived_state_.z_.data(), derived_state_.sigmaz_.data(), edata.getObservedEventsPtr(nroots)); + if (always_check_finite_) { + checkFinite( + gsl::span( + &derived_state_.dJzdz_.at(iztrue * nz * nJ), + nz * nJ), + ModelQuantity::dJzdz, nz); + } } } - - if (always_check_finite_) { - app->checkFinite(derived_state_.dJzdz_, "dJzdz"); - } } void Model::fdJzdsigma(const int ie, const int nroots, const realtype t, @@ -2187,12 +2275,15 @@ void Model::fdJzdsigma(const int ie, const int nroots, const realtype t, state_.fixedParameters.data(), derived_state_.z_.data(), derived_state_.sigmaz_.data(), edata.getObservedEventsPtr(nroots)); + if (always_check_finite_) { + checkFinite( + gsl::span( + &derived_state_.dJzdsigma_.at(iztrue * nz * nJ), + nz * nJ), + ModelQuantity::dJzdsigma, nz); + } } } - - if (always_check_finite_) { - app->checkFinite(derived_state_.dJzdsigma_, "dJzdsigma"); - } } void Model::fdJzdp(const int ie, const int nroots, realtype t, @@ -2303,12 +2394,15 @@ void Model::fdJrzdz(const int ie, const int nroots, const realtype t, state_.unscaledParameters.data(), state_.fixedParameters.data(), derived_state_.rz_.data(), derived_state_.sigmaz_.data()); + if (always_check_finite_) { + checkFinite( + gsl::span( + &derived_state_.dJrzdz_.at(iztrue * nz * nJ), + nz * nJ), + ModelQuantity::dJrzdz, nz); + } } } - - if (always_check_finite_) { - app->checkFinite(derived_state_.dJrzdz_, "dJrzdz"); - } } void Model::fdJrzdsigma(const int ie, const int nroots, const realtype t, @@ -2327,12 +2421,15 @@ void Model::fdJrzdsigma(const int ie, const int nroots, const realtype t, state_.unscaledParameters.data(), state_.fixedParameters.data(), derived_state_.rz_.data(), derived_state_.sigmaz_.data()); + if (always_check_finite_) { + checkFinite( + gsl::span( + &derived_state_.dJrzdsigma_.at(iztrue * nz * nJ), + nz * nJ), + ModelQuantity::dJrzdsigma, nz); + } } } - - if (always_check_finite_) { - app->checkFinite(derived_state_.dJrzdsigma_, "dJrzdsigma"); - } } void Model::fw(const realtype t, const realtype *x) { diff --git a/src/solver.cpp b/src/solver.cpp index 091c8bb7f2..37752261c9 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -172,6 +172,8 @@ void Solver::setup(const realtype t0, Model *model, const AmiVector &x0, /* calculate consistent DAE initial conditions (no effect for ODE) */ if (model->nt() > 1) calcIC(model->getTimepoint(1)); + + cpu_time_ = 0.0; } void Solver::setupB(int *which, const realtype tf, Model *model, @@ -203,6 +205,8 @@ void Solver::setupB(int *which, const realtype tf, Model *model, applyQuadTolerancesASA(*which); setStabLimDetB(*which, stldet_); + + cpu_timeB_ = 0.0; } void Solver::setupSteadystate(const realtype t0, Model *model, const AmiVector &x0, diff --git a/src/solver_idas.cpp b/src/solver_idas.cpp index 0039f6de0a..46d3097b18 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -1050,7 +1050,8 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, return AMICI_MAX_TIME_EXCEEDED; } - if (t > 1e200 && !model->app->checkFinite(gsl::make_span(x), "fxdot")) { + if (t > 1e200 + && !model->checkFinite(gsl::make_span(x), ModelQuantity::xdot)) { /* when t is large (typically ~1e300), CVODES may pass all NaN x to fxdot from which we typically cannot recover. To save time on normal execution, we do not always want to check finiteness diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index a2acec9596..1da5fc31bd 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -640,6 +640,18 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, convergence_check_frequency = 25; while (true) { + /* check for maxsteps */ + if (sim_steps >= solver.getMaxSteps()) { + throw NewtonFailure(AMICI_TOO_MUCH_WORK, + "exceeded maximum number of steps"); + } + if (state_.t >= 1e200) { + throw NewtonFailure(AMICI_NO_STEADY_STATE, + "simulated to late time" + " point without convergence of RHS"); + } + /* increase counter */ + sim_steps++; /* One step of ODE integration reason for tout specification: max with 1 ensures correct direction (any positive value would do) @@ -648,6 +660,7 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, only direction w.r.t. current t */ solver.step(std::max(state_.t, 1.0) * 10); + if (backward) { solver.writeSolution(&state_.t, xB_, state_.dx, state_.sx, xQ_); } else { @@ -670,17 +683,6 @@ void SteadystateProblem::runSteadystateSimulation(const Solver &solver, if (wrms_ < conv_thresh) break; // converged - /* increase counter, check for maxsteps */ - sim_steps++; - if (sim_steps >= solver.getMaxSteps()) { - throw NewtonFailure(AMICI_TOO_MUCH_WORK, - "exceeded maximum number of steps"); - } - if (state_.t >= 1e200) { - throw NewtonFailure(AMICI_NO_STEADY_STATE, - "simulated to late time" - " point without convergence of RHS"); - } } // if check_sensi_conv_ is deactivated, we still have to update sensis diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 14f319d77b..6f3f1b3f09 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -65,9 +65,12 @@ def _test_case(case, model_type): # compile amici model if case.startswith('0006') and model_type != "pysb": petab.flatten_timepoint_specific_output_overrides(problem) - model_output_dir = f'amici_models/model_{case}' + model_name = f"petab_{model_type}_test_case_{case}" + model_output_dir = f'amici_models/{model_name}' model = import_petab_problem( - problem, model_output_dir=model_output_dir, + petab_problem=problem, + model_output_dir=model_output_dir, + model_name=model_name, force_compile=True) solver = model.getSolver() solver.setSteadyStateToleranceFactor(1.0) diff --git a/version.txt b/version.txt index 36e40e13e0..f41aa7fad3 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.11.29 +0.11.30