diff --git a/.github/workflows/sbml-semantic-test-suite.yml b/.github/workflows/sbml-semantic-test-suite.yml index ad7d27e48d..1cbdf6c34c 100644 --- a/.github/workflows/sbml-semantic-test-suite.yml +++ b/.github/workflows/sbml-semantic-test-suite.yml @@ -4,7 +4,6 @@ on: branches: - develop - master - - test_action pull_request: branches: - master diff --git a/.github/workflows/test-large-model.yml b/.github/workflows/test-large-model.yml new file mode 100644 index 0000000000..07445a43df --- /dev/null +++ b/.github/workflows/test-large-model.yml @@ -0,0 +1,75 @@ +name: Performance test +on: + push: + branches: + - develop + - master + - fix_862_powsimp + pull_request: + branches: + - master + +jobs: + build: + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v1 + with: + fetch-depth: 20 + + # install dependencies + - name: apt + run: | + sudo apt-get install -y swig3.0 libatlas-base-dev \ + && sudo ln -s /usr/bin/swig3.0 /usr/bin/swig + - name: pip + run: | + pip3 install --upgrade --user wheel \ + && pip3 install --upgrade --user setuptools + - run: pip3 install petab shyaml + - run: | + echo ::add-path::${HOME}/.local/bin/ + echo ::add-path::${GITHUB_WORKSPACE}/tests/performance/ + # install AMICI + - name: Create AMICI sdist + run: | + cd python/sdist \ + && check_time.sh create_sdist /usr/bin/python3 setup.py sdist + - name: Install AMICI sdist + run: | + AMICI_PARALLEL_COMPILE=2 check_time.sh \ + install_sdist pip3 install -v --user \ + $(ls -t python/sdist/dist/amici-*.tar.gz | head -1) + + # retrieve test model + - run: git clone --depth 1 https://github.com/ICB-DCM/CS_Signalling_ERBB_RAS_AKT + + # import test model + - name: Import test model + run: | + cd CS_Signalling_ERBB_RAS_AKT \ + && check_time.sh \ + petab_import amici_import_petab.py -v \ + -s 'PEtab/CS_Signalling_ERBB_RAS_AKT_petab.xml' \ + -c 'PEtab/conditions_petab.tsv' \ + -m 'PEtab/measurements_petab.tsv' \ + -p 'PEtab/parameters_petab.tsv' --no-compile + + # install model package + - name: Install test model + run: | + check_time.sh install_model \ + python3 CS_Signalling_ERBB_RAS_AKT/CS_Signalling_ERBB_RAS_AKT_petab/setup.py install --user + + # run simulations + - name: forward_simulation + run: | + check_time.sh forward_simulation tests/performance/test.py forward_simulation + - name: forward_sensitivities + run: | + check_time.sh forward_sensitivities tests/performance/test.py forward_sensitivities + - name: adjoint_sensitivities + run: | + check_time.sh adjoint_sensitivities tests/performance/test.py adjoint_sensitivities diff --git a/.travis.yml b/.travis.yml index a11d417e50..36af7b7c0e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,8 @@ branches: only: - master - develop + # Release branches, otherwise they cannot be merged into protected master + - /^release.*/ # Version tags - /^v\d+\.\d+(\.\d+)?(-\S*)?$/ @@ -21,7 +23,6 @@ matrix: - libatlas-base-dev - lcov - libboost-serialization-dev - - swig3.0 - g++-5 - libc6-dbg env: @@ -95,7 +96,8 @@ matrix: install: - export BASE_DIR=`pwd` - - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then mkdir -p ~/bin/ && ln -s /usr/bin/swig3.0 ~/bin/swig && export PATH=~/bin/:$PATH; fi # Python distutils only looks for `swig` and does not find `swig3.0` + # Build swig4.0 (not yet available with apt) to include pydoc in source distribution for pypi + - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then scripts/downloadAndBuildSwig.sh && export PATH=${BASE_DIR}/ThirdParty/swig-4.0.1/install/bin:${PATH}; fi - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then export PYTHON_EXECUTABLE=$(which python3); fi # cmake wont be able to find python3 on its own ... - if [[ "$TRAVIS_OS_NAME" == "linux" ]] && [[ "$CI_MODE" == "test" ]]; then pip3 install --upgrade pip==9.0.3 setuptools wheel pkgconfig scipy; fi - if [[ "$TRAVIS_OS_NAME" != "linux" ]] && [[ "$CI_MODE" == "test" ]]; then pip3 install --user --upgrade pip==9.0.3 setuptools wheel pkgconfig scipy; fi diff --git a/INSTALL.md b/INSTALL.md index e8ddcee092..dd802ccf61 100755 --- a/INSTALL.md +++ b/INSTALL.md @@ -52,6 +52,20 @@ and custom installations. For Python-AMICI usage see [https://github.com/ICB-DCM/AMICI/blob/master/documentation/PYTHON.md](https://github.com/ICB-DCM/AMICI/blob/master/documentation/PYTHON.md). + +### Installation of development versions + +To install development versions which have not been released to pypi yet, +you can install AMICI with pip directly from GitHub using: + + pip3 install -e https://github.com/icb-dcm/amici/archive/develop.zip#egg=amici\&subdirectory=python/sdist + +Replace `develop` by the branch or commit you want to install. + +Note that this will probably not work on Windows which does not support +symlinks by default +(https://stackoverflow.com/questions/5917249/git-symlinks-in-windows/49913019#49913019). + ### Light installation In case you only want to use the AMICI Python package for generating model code @@ -308,6 +322,10 @@ implementations such as Accelerate, Intel MKL, cblas, openblas, atlas. The matlab interface uses the MATLAB MKL, which requires no separate installation. +On Ubuntu, this requirement can be satisfied with + + apt install libatlas-base-dev + #### C++ compiler All AMICI installations require a C++11-compatible C++ compiler. diff --git a/README.md b/README.md index efac331da8..3270e6608d 100644 --- a/README.md +++ b/README.md @@ -128,6 +128,8 @@ However, the following features are unlikely to be supported: - initial assignments for parameters - `delay()` due to missing SUNDIALS solver support +In addition to SBML, we also plan to implement support for the [Simulation Experiment Description Markup Language (SED-ML)](https://sed-ml.org/). + ## Current build status diff --git a/documentation/FAQ.md b/documentation/FAQ.md index cbdaf30353..ea79fe25cf 100644 --- a/documentation/FAQ.md +++ b/documentation/FAQ.md @@ -42,3 +42,17 @@ __Q__: The simulation/sensitivities I get are incorrect. __A__: There are some known issues, especially with adjoint sensitivities, events and DAEs. If your particular problem is not featured in the [issues](https://github.com/ICB-DCM/AMICI/issues) list, please add it! +--- + +__Q__: I am trying to install the AMICI Python package, but installation fails +with something like + + amici/src/cblas.cpp:16:13: fatal error: cblas.h: No such file or directory + #include + ^~~~~~~~~ + compilation terminated. + error: command 'x86_64-linux-gnu-gcc' failed with exit status 1 + +__A__: You will have to install a CBLAS-compatible BLAS library and/or set +`BLAS_CFLAGS` as described in the [installation guide](INSTALL.md). + diff --git a/documentation/amici_refs.bib b/documentation/amici_refs.bib index 8366aecbe1..0f3c404f39 100644 --- a/documentation/amici_refs.bib +++ b/documentation/amici_refs.bib @@ -145,21 +145,6 @@ @article{FroehlichKal2017 Year = {2017}, Bdsk-Url-1 = {https://doi.org/10.1371/journal.pcbi.1005331}} -@article{FroehlichKes2017, - Author = {Fr\"ohlich, Fabian and Kessler, Thomas and Weindl, Daniel and Shadrin, Alexey and Schmiester, Leonard and Hache, Hendrik and Muradyan, Artur and Schuette, Moritz and Lim, Ji-Hyun and Heinig, Matthias and Theis, Fabian and Lehrach, Hans and Wierling, Christoph and Lange, Bodo and Hasenauer, Jan}, - Date-Added = {2019-03-19 23:04:09 +0100}, - Date-Modified = {2019-03-19 23:04:09 +0100}, - Doi = {10.1101/174094}, - Eprint = {http://www.biorxiv.org/content/early/2017/08/09/174094.full.pdf}, - Journal = {bioRxiv}, - Keywords = {drug response; cancer; adjoint sensitivity; ccle}, - Publisher = {Cold Spring Harbor Labs Journals}, - Title = {Efficient parameterization of large-scale mechanistic models enables drug response prediction for cancer cell lines}, - Url = {http://www.biorxiv.org/content/early/2017/08/09/174094}, - Year = {2017}, - Bdsk-Url-1 = {http://www.biorxiv.org/content/early/2017/08/09/174094}, - Bdsk-Url-2 = {https://doi.org/10.1101/174094}} - @article{FroehlichRei2018, Abstract = {Single-cell time-lapse studies have advanced the quantitative understanding of cellular pathways and their inherent cell-to-cell variability. However, parameters retrieved from individual experiments are model dependent and their estimation is limited, if based on solely one kind of experiment. Hence, methods to integrate data collected under different conditions are expected to improve model validation and information content. Here we present a multi-experiment nonlinear mixed effect modeling approach for mechanistic pathway models, which allows the integration of multiple single-cell perturbation experiments. We apply this approach to the translation of green fluorescent protein after transfection using a massively parallel read-out of micropatterned single-cell arrays. We demonstrate that the integration of data from perturbation experiments allows the robust reconstruction of cell-to-cell variability, i.e., parameter densities, while each individual experiment provides insufficient information. Indeed, we show that the integration of the datasets on the population level also improves the estimates for individual cells by breaking symmetries, although each of them is only measured in one experiment. Moreover, we confirmed that the suggested approach is robust with respect to batch effects across experimental replicates and can provide mechanistic insights into the nature of batch effects. We anticipate that the proposed multi-experiment nonlinear mixed effect modeling approach will serve as a basis for the analysis of cellular heterogeneity in single-cell dynamics.}, Author = {Fr{\"o}hlich, Fabian and Reiser, Anita and Fink, Laura and Wosch{\'e}e, Daniel and Ligon, Thomas and Theis, Fabian Joachim and R{\"a}dler, Joachim Oskar and Hasenauer, Jan}, @@ -478,4 +463,61 @@ @Article{SchmiesterSch2019 url = {https://doi.org/10.1093/bioinformatics/btz581}, } +@Article{SchmiesterWei2019, + author = {Schmiester, Leonard and Weindl, Daniel and Hasenauer, Jan}, + title = {Statistical inference of mechanistic models from qualitative data using an efficient optimal scaling approach}, + journal = {bioRxiv}, + year = {2019}, + abstract = {Quantitative dynamical models facilitate the understanding of biological processes and the prediction of their dynamics. These models usually comprise unknown parameters, which have to be inferred from experimental data. For quantitative experimental data, there are several methods and software tools available. However, for qualitative data the available approaches are limited and computationally demanding.Here, we consider the optimal scaling method which has been developed in statistics for categorical data and has been applied to dynamical systems. This approach turns qualitative variables into quantitative ones, accounting for constraints on their relation. We derive a reduced formulation for the optimization problem defining the optimal scaling. The reduced formulation possesses the same optimal points as the established formulation but requires less degrees of freedom. Parameter estimation for dynamical models of cellular pathways revealed that the reduced formulation improves the robustness and convergence of optimizers. This resulted in substantially reduced computation times.We implemented the proposed approach in the open-source Python Parameter EStimation TOolbox (pyPESTO) to facilitate reuse and extension. The proposed approach enables efficient parameterization of quantitative dynamical models using qualitative data.}, + doi = {10.1101/848648}, + elocation-id = {848648}, + eprint = {https://www.biorxiv.org/content/early/2019/11/20/848648.full.pdf}, + publisher = {Cold Spring Harbor Laboratory}, + url = {https://www.biorxiv.org/content/early/2019/11/20/848648}, +} + +@Article{PedretscherKal2019, + author = {B. Pedretscher and B. Kaltenbacher and O. Pfeiler}, + title = {Parameter identification and uncertainty quantification in stochastic state space models and its application to texture analysis}, + journal = {Applied Numerical Mathematics}, + year = {2019}, + volume = {146}, + pages = {38 - 54}, + issn = {0168-9274}, + abstract = {In this paper, a computational framework, which enables efficient and robust parameter identification, as well as uncertainty quantification in state space models based on Itô stochastic processes, is presented. For optimization, a Maximum Likelihood approach based on the system's corresponding Fokker-Planck equation is followed. Gradient information is included by means of an adjoint approach, which is based on the Lagrangian of the optimization problem. To quantify the uncertainty of the Maximum-A-Posteriori estimates of the model parameters, a Bayesian inference approach based on Markov Chain Monte Carlo simulations, as well as profile likelihoods are implemented and compared in terms of runtime and accuracy. The framework is applied to experimental electron backscatter diffraction data of a fatigued metal film, where the aim is to develop a model, which consistently and physically meaningfully captures the metal's microstructural changes that are caused by external loading.}, + doi = {https://doi.org/10.1016/j.apnum.2019.06.020}, + keywords = {Stochastic state space model, Parameter identification, Uncertainty quantification, Profile likelihood, Adjoint approach, Fokker-Planck equation, Thermo-mechanical fatigue, Texture analysis}, + url = {http://www.sciencedirect.com/science/article/pii/S0168927419301722}, +} + +@Article{FroehlichKes2018, + author = {Fröhlich, Fabian and Kessler, Thomas and Weindl, Daniel and Shadrin, Alexey and Schmiester, Leonard and Hache, Hendrik and Muradyan, Artur and Schütte, Moritz and Lim, Ji-Hyun and Heinig, Matthias and Theis, Fabian J. and Lehrach, Hans and Wierling, Christoph and Lange, Bodo and Hasenauer, Jan}, + title = {Efficient Parameter Estimation Enables the Prediction of Drug Response Using a Mechanistic Pan-Cancer Pathway Model}, + journal = {Cell Systems}, + year = {2018}, + volume = {7}, + number = {6}, + pages = {567--579.e6}, + issn = {2405-4712}, + abstract = {Mechanistic models are essential to deepen the understanding of complex diseases at the molecular level. Nowadays, high-throughput molecular and phenotypic characterizations are possible, but the integration of such data with prior knowledge on signaling pathways is limited by the availability of scalable computational methods. Here, we present a computational framework for the parameterization of large-scale mechanistic models and its application to the prediction of drug response of cancer cell lines from exome and transcriptome sequencing data. This framework is over 104 times faster than state-of-the-art methods, which enables modeling at previously infeasible scales. By applying the framework to a model describing major cancer-associated pathways (>1,200 species and >2,600 reactions), we could predict the effect of drug combinations from single drug data. This is the first integration of high-throughput datasets using large-scale mechanistic models. We anticipate this to be the starting point for development of more comprehensive models allowing a deeper mechanistic insight. +Mechanistic models are essential to deepen the understanding of complex diseases at the molecular level. Nowadays, high-throughput molecular and phenotypic characterizations are possible, but the integration of such data with prior knowledge on signaling pathways is limited by the availability of scalable computational methods. Here, we present a computational framework for the parameterization of large-scale mechanistic models and its application to the prediction of drug response of cancer cell lines from exome and transcriptome sequencing data. This framework is over 104 times faster than state-of-the-art methods, which enables modeling at previously infeasible scales. By applying the framework to a model describing major cancer-associated pathways (>1,200 species and >2,600 reactions), we could predict the effect of drug combinations from single drug data. This is the first integration of high-throughput datasets using large-scale mechanistic models. We anticipate this to be the starting point for development of more comprehensive models allowing a deeper mechanistic insight.}, + comment = {doi: 10.1016/j.cels.2018.10.013}, + doi = {10.1016/j.cels.2018.10.013}, + publisher = {Elsevier}, + url = {https://doi.org/10.1016/j.cels.2018.10.013}, +} + +@Article{StaporSch2019, + author = {Stapor, Paul and Schmiester, Leonard and Wierling, Christoph and Lange, Bodo and Weindl, Daniel and Hasenauer, Jan}, + title = {Mini-batch optimization enables training of ODE models on large-scale datasets}, + journal = {bioRxiv}, + year = {2019}, + abstract = {Quantitative dynamical models are widely used to study cellular signal processing. A critical step in modeling is the estimation of unknown model parameters from experimental data. As model sizes and datasets are steadily growing, established parameter optimization approaches for mechanistic models become computationally extremely challenging. However, mini-batch optimization methods, as employed in deep learning, have better scaling properties. In this work, we adapt, apply, and benchmark mini-batch optimization for ordinary differential equation (ODE) models thereby establishing a direct link between dynamic modeling and machine learning. On our main application example, a large-scale model of cancer signaling, we benchmark mini-batch optimization against established methods, achieving better optimization results and reducing computation by more than an order of magnitude. We expect that our work will serve as a first step towards mini-batch optimization tailored to ODE models and enable modeling of even larger and more complex systems than what is currently possible.}, + doi = {10.1101/859884}, + elocation-id = {859884}, + eprint = {https://www.biorxiv.org/content/early/2019/11/30/859884.full.pdf}, + publisher = {Cold Spring Harbor Laboratory}, + url = {https://www.biorxiv.org/content/early/2019/11/30/859884}, +} + @Comment{jabref-meta: databaseType:bibtex;} diff --git a/documentation/development.md b/documentation/development.md index f8b951bc91..b4463bdc85 100644 --- a/documentation/development.md +++ b/documentation/development.md @@ -98,8 +98,10 @@ described below: case for all existing code, any new contributions should do so. * We use Python [type hints](https://docs.python.org/3/library/typing.html) - for all functions. In Python code type hints should be used instead of - doxygen `@type`. (All legacy `@type` attributes are to be removed.) + for all functions (but not for class attributes, since they are not supported + by the current Python doxygen filter). In Python code type hints should be + used instead of doxygen `@type`. (All legacy `@type` attributes are to be + removed.) For function docstrings, follow this format: diff --git a/documentation/references.md b/documentation/references.md index 759461b437..d4739c7040 100644 --- a/documentation/references.md +++ b/documentation/references.md @@ -1,6 +1,6 @@ # References -List of publications using AMICI. Total number is 32. +List of publications using AMICI. Total number is 35. If you applied AMICI in your work and your publication is missing, please let us know via a new Github issue. @@ -24,12 +24,21 @@ If you applied AMICI in your work and your publication is missing, please let us

Nousiainen, Kari, Jukka Intosalmi, and Harri Lähdesmäki. 2019. “A Mathematical Model for Enhancer Activation Kinetics During Cell Differentiation.” In Algorithms for Computational Biology, edited by Ian Holmes, Carlos Martín-Vide, and Miguel A. Vega-Rodríguez, 191–202. Cham: Springer International Publishing.

+
+

Pedretscher, B., B. Kaltenbacher, and O. Pfeiler. 2019. “Parameter Identification and Uncertainty Quantification in Stochastic State Space Models and Its Application to Texture Analysis.” Applied Numerical Mathematics 146: 38–54. https://doi.org/https://doi.org/10.1016/j.apnum.2019.06.020.

+

Pitt, Jake Alan, and Julio R Banga. 2019. “Parameter Estimation in Models of Biological Oscillators: An Automated Regularised Estimation Approach.” BMC Bioinformatics 20 (1): 82. https://doi.org/10.1186/s12859-019-2630-y.

Schmiester, Leonard, Yannik Schälte, Fabian Fröhlich, Jan Hasenauer, and Daniel Weindl. 2019. “Efficient parameterization of large-scale dynamic models based on relative measurements.” Bioinformatics, July. https://doi.org/10.1093/bioinformatics/btz581.

+
+

Schmiester, Leonard, Daniel Weindl, and Jan Hasenauer. 2019. “Statistical Inference of Mechanistic Models from Qualitative Data Using an Efficient Optimal Scaling Approach.” bioRxiv. Cold Spring Harbor Laboratory. https://doi.org/10.1101/848648.

+
+
+

Stapor, Paul, Leonard Schmiester, Christoph Wierling, Bodo Lange, Daniel Weindl, and Jan Hasenauer. 2019. “Mini-Batch Optimization Enables Training of Ode Models on Large-Scale Datasets.” bioRxiv. Cold Spring Harbor Laboratory. https://doi.org/10.1101/859884.

+

Villaverde, Alejandro F., Elba Raimúndez, Jan Hasenauer, and Julio R. Banga. 2019. “A Comparison of Methods for Quantifying Prediction Uncertainty in Systems Biology.” Accepted for Publication in Proc. Of the Foundations of Syst. Biol. In Engin. (FOSBE).

@@ -45,6 +54,9 @@ If you applied AMICI in your work and your publication is missing, please let us

Bast, Lisa, Filippo Calzolari, Michael Strasser, Jan Hasenauer, Fabian J. Theis, Jovica Ninkovic, and Carsten Marr. 2018. “Subtle Changes in Clonal Dynamics Underlie the Age-Related Decline in Neurogenesis.” Cell Reports.

+
+

Fröhlich, Fabian, Thomas Kessler, Daniel Weindl, Alexey Shadrin, Leonard Schmiester, Hendrik Hache, Artur Muradyan, et al. 2018. “Efficient Parameter Estimation Enables the Prediction of Drug Response Using a Mechanistic Pan-Cancer Pathway Model.” Cell Systems 7 (6). Elsevier: 567–579.e6. https://doi.org/10.1016/j.cels.2018.10.013.

+

Fröhlich, Fabian, Anita Reiser, Laura Fink, Daniel Woschée, Thomas Ligon, Fabian Joachim Theis, Joachim Oskar Rädler, and Jan Hasenauer. 2018. “Multi-Experiment Nonlinear Mixed Effect Modeling of Single-Cell Translation Kinetics After Transfection.” Npj Systems Biology and Applications 5 (1): 1. https://doi.org/10.1038/s41540-018-0079-7.

@@ -72,9 +84,6 @@ If you applied AMICI in your work and your publication is missing, please let us

Ballnus, B., S. Hug, K. Hatz, L. Görlitz, J. Hasenauer, and F. J. Theis. 2017. “Comprehensive Benchmarking of Markov Chain Monte Carlo Methods for Dynamical Systems.” BMC Syst. Biol. 11 (63). https://doi.org/10.1186/s12918-017-0433-1.

-
-

Fröhlich, Fabian, Thomas Kessler, Daniel Weindl, Alexey Shadrin, Leonard Schmiester, Hendrik Hache, Artur Muradyan, et al. 2017. “Efficient Parameterization of Large-Scale Mechanistic Models Enables Drug Response Prediction for Cancer Cell Lines.” bioRxiv. Cold Spring Harbor Labs Journals. https://doi.org/10.1101/174094.

-

Fröhlich, F., B. Kaltenbacher, F. J. Theis, and J. Hasenauer. 2017. “Scalable Parameter Estimation for Genome-Scale Biochemical Reaction Networks.” PLoS Comput. Biol. 13 (1): e1005331. https://doi.org/10.1371/journal.pcbi.1005331.

diff --git a/include/amici/amici.h b/include/amici/amici.h index a0848e692a..89b4bbca90 100644 --- a/include/amici/amici.h +++ b/include/amici/amici.h @@ -13,32 +13,105 @@ namespace amici { /*! - * @brief Prints a specified error message associated to the specified + * @brief Prints a specified error message associated with the specified * identifier * - * @param identifier error identifier - * @param format string with error message printf-style format - * @param ... arguments to be formatted + * @param id error identifier + * @param message error message */ void -printErrMsgIdAndTxt(const char* identifier, const char* format, ...); +printErrMsgIdAndTxt(std::string const& id, std::string const& message); /*! - * @brief Prints a specified warning message associated to the specified + * @brief Prints a specified warning message associated with the specified * identifier * - * @param identifier warning identifier - * @param format string with error message printf-style format - * @param ... arguments to be formatted + * @param id warning identifier + * @param message warning message */ void -printWarnMsgIdAndTxt(const char* identifier, const char* format, ...); +printWarnMsgIdAndTxt(std::string const& id, std::string const& message); -/** errMsgIdAndTxt is a function pointer for printErrMsgIdAndTxt */ -extern msgIdAndTxtFp errMsgIdAndTxt; +/** + * @brief Main class for making calls to AMICI. + * + * This class is used to provide separate AMICI contexts, for example, for use + * in multi-threaded applications where different threads want to use AMICI with + * different settings, such custom logging functions. + * + * NOTE: For this moment, the context object needs to be set manually to any + * Model and Solver object. If not set, they will use the default output + * channel. + */ +class AmiciApplication +{ + public: + AmiciApplication() = default; + + /** + * Core integration routine. Initializes the solver and runs the forward + * and backward problem. + * + * @param solver Solver instance + * @param edata pointer to experimental data object + * @param model model specification object + * @param rethrow rethrow integration exceptions? + * @return rdata pointer to return data object + */ + std::unique_ptr runAmiciSimulation(Solver& solver, + const ExpData* edata, + Model& model, + bool rethrow = false); + + /** + * Same as runAmiciSimulation, but for multiple ExpData instances. + * + * @param solver Solver instance + * @param edatas experimental data objects + * @param model model specification object + * @param failfast flag to allow early termination + * @param num_threads number of threads for parallel execution + * @return vector of pointers to return data objects + */ + std::vector> runAmiciSimulations( + Solver const& solver, + const std::vector& edatas, + Model const& model, + bool failfast, + int num_threads); + + /** Function to process warnings */ + outputFunctionType warning = printWarnMsgIdAndTxt; + + /** Function to process errors */ + outputFunctionType error = printErrMsgIdAndTxt; + + /** + * @brief printf interface to warning() + * @param identifier warning identifier + * @param format string with warning message printf-style format + * @param ... arguments to be formatted + */ + void warningF(const char* identifier, const char* format, ...); + + /** + * @brief printf interface to error() + * @param identifier warning identifier + * @param format string with error message printf-style format + * @param ... arguments to be formatted + */ + void errorF(const char* identifier, const char* format, ...); -/** warnMsgIdAndTxt is a function pointer for printWarnMsgIdAndTxt */ -extern msgIdAndTxtFp warnMsgIdAndTxt; + /** + * @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); +}; /*! * runAmiciSimulation is the core integration routine. It initializes the solver diff --git a/include/amici/defines.h b/include/amici/defines.h index 7482578645..1c90baa92c 100644 --- a/include/amici/defines.h +++ b/include/amici/defines.h @@ -2,6 +2,7 @@ #define AMICI_DEFINES_H #include +#include namespace amici { @@ -148,13 +149,17 @@ enum class NewtonStatus { newt_sim_newt=3, }; +/** Damping factor flag for the Newton method */ +enum class NewtonDampingFactorMode { + off = 0, + on = 1 +}; + /** - * @brief msgIdAndTxtFp - * @param identifier string with error message identifier - * @param format string with error message printf-style format - * @param ... arguments to be formatted + * Type for function to process warnings or error messages. */ -using msgIdAndTxtFp = void (*)(const char *, const char *, ...); +using outputFunctionType = std::function; // clang-format on diff --git a/include/amici/misc.h b/include/amici/misc.h index 6277996699..6284a1d8f9 100644 --- a/include/amici/misc.h +++ b/include/amici/misc.h @@ -25,14 +25,6 @@ namespace amici { gsl::span slice(std::vector &data, int index, unsigned size); -/** - * @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); /** @@ -92,6 +84,15 @@ std::string backtraceString(int maxFrames); */ std::string regexErrorToString(std::regex_constants::error_type err_type); +/** + * @brief Format printf-style arguments to std::string + * @param fmt Format string + * @param ap Argument list pointer + * @return Formatted String + */ +std::string printfToString(const char *fmt, va_list ap); + + } // namespace amici #ifndef __cpp_lib_make_unique diff --git a/include/amici/model.h b/include/amici/model.h index 351c46f9b0..b08e5492cd 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -15,6 +15,9 @@ namespace amici { class ExpData; class Model; class Solver; +class AmiciApplication; + +extern AmiciApplication defaultContext; } // namespace amici @@ -1154,6 +1157,9 @@ class Model : public AbstractModel { * nx_solver, row-major) */ AmiVectorArray dxdotdp; + /** AMICI context */ + AmiciApplication *app = &defaultContext; + protected: /** * @brief Writes part of a slice to a buffer according to indices specified diff --git a/include/amici/model_dae.h b/include/amici/model_dae.h index 3218519900..bb8ed38cd2 100644 --- a/include/amici/model_dae.h +++ b/include/amici/model_dae.h @@ -13,7 +13,6 @@ #include namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; class ExpData; class IDASolver; diff --git a/include/amici/model_ode.h b/include/amici/model_ode.h index 41aa2be377..5e223af243 100644 --- a/include/amici/model_ode.h +++ b/include/amici/model_ode.h @@ -14,7 +14,6 @@ #include namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; class CVodeSolver; diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h index c0c1232b31..5482e0a208 100644 --- a/include/amici/newton_solver.h +++ b/include/amici/newton_solver.h @@ -46,6 +46,8 @@ class NewtonSolver { * @param maxsteps maximum number of allowed Newton steps for steady state computation * @param atol absolute tolerance * @param rtol relative tolerance + * @param dampingFactorMode to switch on/off a use of a damping factor + * @param dampingFactorLowerBound a lower bound for the damping factor * @return solver NewtonSolver according to the specified linsolType */ static std::unique_ptr getSolver(realtype *t, AmiVector *x, @@ -54,7 +56,9 @@ class NewtonSolver { ReturnData *rdata, int maxlinsteps, int maxsteps, - double atol, double rtol); + double atol, double rtol, + NewtonDampingFactorMode dampingFactorMode, + double dampingFactorLowerBound); /** * Computes the solution of one Newton iteration @@ -103,6 +107,10 @@ class NewtonSolver { double atol = 1e-16; /** relative tolerance */ double rtol = 1e-8; + /** damping factor flag */ + NewtonDampingFactorMode dampingFactorMode = NewtonDampingFactorMode::on; + /** damping factor lower bound */ + double dampingFactorLowerBound = 1e-8; protected: /** time variable */ diff --git a/include/amici/serialization.h b/include/amici/serialization.h index 15dee2c22d..6fb4779e8d 100644 --- a/include/amici/serialization.h +++ b/include/amici/serialization.h @@ -57,6 +57,8 @@ void serialize(Archive &ar, amici::Solver &u, const unsigned int version) { ar &u.requires_preequilibration; ar &u.newton_maxsteps; ar &u.newton_maxlinsteps; + ar &u.newton_damping_factor_mode; + ar &u.newton_damping_factor_lower_bound; ar &u.ism; ar &u.sensi_meth; ar &u.linsol; diff --git a/include/amici/solver.h b/include/amici/solver.h index 1e9c6dfb41..8ab0b75c17 100644 --- a/include/amici/solver.h +++ b/include/amici/solver.h @@ -5,6 +5,7 @@ #include "amici/sundials_linsol_wrapper.h" #include "amici/symbolic_functions.h" #include "amici/vector.h" +#include "amici/amici.h" #include #include @@ -16,6 +17,9 @@ class ForwardProblem; class BackwardProblem; class Model; class Solver; +class AmiciApplication; + +extern AmiciApplication defaultContext; } // namespace amici // for serialization friend in Solver @@ -44,6 +48,12 @@ class Solver { public: Solver() = default; + /** + * @brief Constructor + * @param app AMICI application context + */ + Solver(AmiciApplication *app); + /** * @brief Solver copy constructor * @param other @@ -219,6 +229,30 @@ class Solver { */ void setNewtonMaxLinearSteps(int newton_maxlinsteps); + /** + * @brief Get a state of the damping factor used in the Newton solver + * @return + */ + NewtonDampingFactorMode getNewtonDampingFactorMode() const; + + /** + * @brief Turn on/off a damping factor in the Newton method + * @param dampingFactorMode + */ + void setNewtonDampingFactorMode(NewtonDampingFactorMode dampingFactorMode); + + /** + * @brief Get a lower bound of the damping factor used in the Newton solver + * @return + */ + double getNewtonDampingFactorLowerBound() const; + + /** + * @brief Set a lower bound of the damping factor in the Newton solver + * @param dampingFactorLowerBound + */ + void setNewtonDampingFactorLowerBound(double dampingFactorLowerBound); + /** * @brief Get sensitvity order * @return sensitivity order @@ -678,6 +712,9 @@ class Solver { */ friend bool operator==(const Solver &a, const Solver &b); + /** AMICI context */ + AmiciApplication *app = &defaultContext; + protected: /** * @brief Sets a timepoint at which the simulation will be stopped @@ -1381,6 +1418,12 @@ class Solver { * computation */ long int newton_maxlinsteps = 0; + /** Damping factor state used int the Newton method */ + NewtonDampingFactorMode newton_damping_factor_mode = NewtonDampingFactorMode::on; + + /** Lower bound of the damping factor. */ + double newton_damping_factor_lower_bound = 1e-8; + /** Enable model preequilibration */ bool requires_preequilibration = false; @@ -1460,13 +1503,13 @@ bool operator==(const Solver &a, const Solver &b); /** * @brief Extracts diagnosis information from solver memory block and - * writes them into the return data object for the backward problem + * passes them to the specified output function * * @param error_code error identifier * @param module name of the module in which the error occured * @param function name of the function in which the error occured * @param msg error message - * @param eh_data unused input + * @param eh_data amici::Solver as void* */ void wrapErrHandlerFn(int error_code, const char *module, const char *function, char *msg, void *eh_data); diff --git a/include/amici/solver_cvodes.h b/include/amici/solver_cvodes.h index ce2278c3d8..72e012aca1 100644 --- a/include/amici/solver_cvodes.h +++ b/include/amici/solver_cvodes.h @@ -26,8 +26,6 @@ namespace amici { class CVodeSolver : public Solver { public: - CVodeSolver() = default; - ~CVodeSolver() override = default; /** diff --git a/include/amici/solver_idas.h b/include/amici/solver_idas.h index 21648ddd88..5a9bf100d8 100644 --- a/include/amici/solver_idas.h +++ b/include/amici/solver_idas.h @@ -24,7 +24,6 @@ namespace amici { class IDASolver : public Solver { public: - IDASolver() = default; ~IDASolver() override = default; /** diff --git a/models/model_dirac/CMakeLists.txt b/models/model_dirac/CMakeLists.txt index 2fbbab79a1..f73a171551 100644 --- a/models/model_dirac/CMakeLists.txt +++ b/models/model_dirac/CMakeLists.txt @@ -52,6 +52,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_dirac/main.cpp b/models/model_dirac/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_dirac/main.cpp +++ b/models/model_dirac/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_dirac/model_dirac.h b/models/model_dirac/model_dirac.h index 839b28bb78..389dce41c7 100644 --- a/models/model_dirac/model_dirac.h +++ b/models/model_dirac/model_dirac.h @@ -1,6 +1,6 @@ #ifndef _amici_model_dirac_h #define _amici_model_dirac_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -61,7 +61,7 @@ class Model_model_dirac : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_dirac(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_dirac(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_events/CMakeLists.txt b/models/model_events/CMakeLists.txt index c60998f9e6..150dd43b02 100644 --- a/models/model_events/CMakeLists.txt +++ b/models/model_events/CMakeLists.txt @@ -66,6 +66,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_events/main.cpp b/models/model_events/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_events/main.cpp +++ b/models/model_events/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_events/model_events.h b/models/model_events/model_events.h index 1ea53269cf..f3d592eb5e 100644 --- a/models/model_events/model_events.h +++ b/models/model_events/model_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_events_h #define _amici_model_events_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -75,7 +75,7 @@ class Model_model_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_events(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint/CMakeLists.txt b/models/model_jakstat_adjoint/CMakeLists.txt index 27f5fd7c62..76484be708 100644 --- a/models/model_jakstat_adjoint/CMakeLists.txt +++ b/models/model_jakstat_adjoint/CMakeLists.txt @@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_jakstat_adjoint/main.cpp b/models/model_jakstat_adjoint/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_jakstat_adjoint/main.cpp +++ b/models/model_jakstat_adjoint/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/models/model_jakstat_adjoint/model_jakstat_adjoint.h index e512ee9394..3cf5fdec25 100644 --- a/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/models/model_jakstat_adjoint/model_jakstat_adjoint.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_h #define _amici_model_jakstat_adjoint_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -64,7 +64,7 @@ class Model_model_jakstat_adjoint : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_jakstat_adjoint_o2/CMakeLists.txt b/models/model_jakstat_adjoint_o2/CMakeLists.txt index 3fad280043..92516553d0 100644 --- a/models/model_jakstat_adjoint_o2/CMakeLists.txt +++ b/models/model_jakstat_adjoint_o2/CMakeLists.txt @@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_jakstat_adjoint_o2/main.cpp b/models/model_jakstat_adjoint_o2/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_jakstat_adjoint_o2/main.cpp +++ b/models/model_jakstat_adjoint_o2/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index 52f25fccd9..ef9f6a1638 100644 --- a/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_jakstat_adjoint_o2_h #define _amici_model_jakstat_adjoint_o2_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -64,7 +64,7 @@ class Model_model_jakstat_adjoint_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_jakstat_adjoint_o2(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_jakstat_adjoint_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_nested_events/CMakeLists.txt b/models/model_nested_events/CMakeLists.txt index b5fa8de311..120e5a08ed 100644 --- a/models/model_nested_events/CMakeLists.txt +++ b/models/model_nested_events/CMakeLists.txt @@ -55,6 +55,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_nested_events/main.cpp b/models/model_nested_events/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_nested_events/main.cpp +++ b/models/model_nested_events/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_nested_events/model_nested_events.h b/models/model_nested_events/model_nested_events.h index 76a3c9b80f..5ca7091d6f 100644 --- a/models/model_nested_events/model_nested_events.h +++ b/models/model_nested_events/model_nested_events.h @@ -1,6 +1,6 @@ #ifndef _amici_model_nested_events_h #define _amici_model_nested_events_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -64,7 +64,7 @@ class Model_model_nested_events : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_nested_events(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_nested_events(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron/CMakeLists.txt b/models/model_neuron/CMakeLists.txt index b7af8a3d7d..0c63349ca8 100644 --- a/models/model_neuron/CMakeLists.txt +++ b/models/model_neuron/CMakeLists.txt @@ -69,6 +69,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_neuron/main.cpp b/models/model_neuron/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_neuron/main.cpp +++ b/models/model_neuron/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_neuron/model_neuron.h b/models/model_neuron/model_neuron.h index ddf2f6e014..fa33f518ab 100644 --- a/models/model_neuron/model_neuron.h +++ b/models/model_neuron/model_neuron.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_h #define _amici_model_neuron_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */ #include #include #include "amici/defines.h" @@ -78,7 +78,7 @@ class Model_model_neuron : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_neuron_o2/CMakeLists.txt b/models/model_neuron_o2/CMakeLists.txt index 7ad5944d56..7c2121cb0f 100644 --- a/models/model_neuron_o2/CMakeLists.txt +++ b/models/model_neuron_o2/CMakeLists.txt @@ -71,6 +71,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_neuron_o2/main.cpp b/models/model_neuron_o2/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_neuron_o2/main.cpp +++ b/models/model_neuron_o2/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_neuron_o2/model_neuron_o2.h b/models/model_neuron_o2/model_neuron_o2.h index ef16e543d1..a4fb5d03d9 100644 --- a/models/model_neuron_o2/model_neuron_o2.h +++ b/models/model_neuron_o2/model_neuron_o2.h @@ -1,6 +1,6 @@ #ifndef _amici_model_neuron_o2_h #define _amici_model_neuron_o2_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */ #include #include #include "amici/defines.h" @@ -80,7 +80,7 @@ class Model_model_neuron_o2 : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_neuron_o2(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_neuron_o2(J, t, x, p, k, h, w, dwdx); diff --git a/models/model_robertson/CMakeLists.txt b/models/model_robertson/CMakeLists.txt index 072011db85..7af14aba5b 100644 --- a/models/model_robertson/CMakeLists.txt +++ b/models/model_robertson/CMakeLists.txt @@ -53,6 +53,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_robertson/main.cpp b/models/model_robertson/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_robertson/main.cpp +++ b/models/model_robertson/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_robertson/model_robertson.h b/models/model_robertson/model_robertson.h index 556cea6e75..c526ad14ad 100644 --- a/models/model_robertson/model_robertson.h +++ b/models/model_robertson/model_robertson.h @@ -1,6 +1,6 @@ #ifndef _amici_model_robertson_h #define _amici_model_robertson_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5193a397de39f5bd4b7d98718067dc71f02f916e */ #include #include #include "amici/defines.h" @@ -62,7 +62,7 @@ class Model_model_robertson : public amici::Model_DAE { virtual amici::Model* clone() const override { return new Model_model_robertson(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5193a397de39f5bd4b7d98718067dc71f02f916e"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype cj, const realtype *dx, const realtype *w, const realtype *dwdx) override { J_model_robertson(J, t, x, p, k, h, cj, dx, w, dwdx); diff --git a/models/model_steadystate/CMakeLists.txt b/models/model_steadystate/CMakeLists.txt index 0537b4b5b1..0b14584260 100644 --- a/models/model_steadystate/CMakeLists.txt +++ b/models/model_steadystate/CMakeLists.txt @@ -52,6 +52,13 @@ ${MODEL_DIR}/wrapfunctions.cpp add_library(${PROJECT_NAME} ${SRC_LIST_LIB}) add_library(model ALIAS ${PROJECT_NAME}) +# This option can be helpful when using the Intel compiler and compilation of +# wrapfunctions.cpp fails due to insufficient memory. +option(ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS "Enable compiler optimizations for wrapfunctions.cpp?" ON) +if(NOT ENABLE_WRAPFUNCTIONS_OPTIMIZATIONS) + set_source_files_properties(wrapfunctions.cpp PROPERTIES COMPILE_FLAGS -O0) +endif() + target_include_directories(${PROJECT_NAME} PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}") target_link_libraries(${PROJECT_NAME} diff --git a/models/model_steadystate/main.cpp b/models/model_steadystate/main.cpp index 6d6caf54a8..07ff62a322 100644 --- a/models/model_steadystate/main.cpp +++ b/models/model_steadystate/main.cpp @@ -10,8 +10,7 @@ #include "wrapfunctions.h" /* model-provided functions */ /* This is a scaffold for a stand-alone AMICI simulation executable - * demonstrating - * use of the AMICI C++ API. + * demonstrating the use of the AMICI C++ API. * * This program reads AMICI options from an HDF5 file, prints some results * and writes additional results to an HDF5 file. The name of the HDF5 file diff --git a/models/model_steadystate/model_steadystate.h b/models/model_steadystate/model_steadystate.h index 06a1c2773b..4489e3a664 100644 --- a/models/model_steadystate/model_steadystate.h +++ b/models/model_steadystate/model_steadystate.h @@ -1,6 +1,6 @@ #ifndef _amici_model_steadystate_h #define _amici_model_steadystate_h -/* Generated by amiwrap (R2017b) 5afd3fccf92ca20486e1bbaa4a50343b107e239c */ +/* Generated by amiwrap (R2017b) 5dfc4abddf5cb50611e9881fcab074d8859ccc6b */ #include #include #include "amici/defines.h" @@ -61,7 +61,7 @@ class Model_model_steadystate : public amici::Model_ODE { virtual amici::Model* clone() const override { return new Model_model_steadystate(*this); }; - const std::string getAmiciCommit() const override { return "5afd3fccf92ca20486e1bbaa4a50343b107e239c"; }; + const std::string getAmiciCommit() const override { return "5dfc4abddf5cb50611e9881fcab074d8859ccc6b"; }; virtual void fJ(realtype *J, const realtype t, const realtype *x, const double *p, const double *k, const realtype *h, const realtype *w, const realtype *dwdx) override { J_model_steadystate(J, t, x, p, k, h, w, dwdx); diff --git a/python/amici/__init__.py b/python/amici/__init__.py index c26c9ff281..4ae68672ef 100644 --- a/python/amici/__init__.py +++ b/python/amici/__init__.py @@ -157,6 +157,8 @@ def ExpData(*args): # the constructor signature may have changed and you are glad this # wrapper did not break. return amici.ExpData(args[0].get(), *args[:1]) + elif isinstance(args[0], ModelPtr): + return amici.ExpData(args[0].get()) else: return amici.ExpData(*args) diff --git a/python/amici/__init__.template.py b/python/amici/__init__.template.py index b68a5bae90..5c42765ec7 100644 --- a/python/amici/__init__.template.py +++ b/python/amici/__init__.template.py @@ -11,6 +11,6 @@ 'model version or regenerate the model with the AMICI ' 'currently in your path.') -from TPL_MODELNAME.TPL_MODELNAME import * +from TPL_MODELNAME._TPL_MODELNAME import * __version__ = 'TPL_PACKAGE_VERSION' diff --git a/python/amici/ode_export.py b/python/amici/ode_export.py index cfcb1fd97d..552a43466f 100644 --- a/python/amici/ode_export.py +++ b/python/amici/ode_export.py @@ -18,8 +18,11 @@ pysb = None +from typing import Callable, Optional from string import Template import sympy.printing.ccode as ccode +from sympy.matrices.immutable import ImmutableDenseMatrix +from sympy.matrices.dense import MutableDenseMatrix from . import ( amiciSwigPath, amiciSrcPath, amiciModulePath, __version__, __commit__ @@ -712,12 +715,17 @@ class ODEModel: derivative from a partial derivative call to enforce a partial derivative in the next recursion. prevents infinite recursion + _simplify: If not None, this function will be used to simplify symbolic + derivative expressions. Receives sympy expressions as only argument. + To apply multiple simplifications, wrap them in a lambda expression. + NOTE: This does currently not work with PySB symbols. """ - def __init__(self): + def __init__(self, simplify: Optional[Callable] = sp.powsimp): """Create a new ODEModel instance. Arguments: + simplify: see ODEModel._simplify Raises: @@ -804,6 +812,7 @@ def __init__(self): } self._lock_total_derivative = False + self._simplify = simplify def import_from_sbml_importer(self, si): """Imports a model specification from a amici.SBMLImporter instance. @@ -825,7 +834,8 @@ def import_from_sbml_importer(self, si): self._eqs['dxdotdw'] = si.stoichiometricMatrix self._eqs['w'] = si.fluxVector self._syms['w'] = sp.Matrix( - [sp.Symbol(f'flux_r{idx}') for idx in range(len(si.fluxVector))] + [sp.Symbol(f'flux_r{idx}', real=True) + for idx in range(len(si.fluxVector))] ) self._eqs['dxdotdx'] = sp.zeros(si.stoichiometricMatrix.shape[0]) if len(si.stoichiometricMatrix): @@ -1160,12 +1170,12 @@ def _generateSymbol(self, name): for comp in getattr(self, component) ]) self._strippedsyms[name] = sp.Matrix([ - sp.Symbol(comp.get_name()) + sp.Symbol(comp.get_name(), real=True) for comp in getattr(self, component) ]) if name == 'y': self._syms['my'] = sp.Matrix([ - sp.Symbol(f'm{strip_pysb(comp.get_id())}') + sp.Symbol(f'm{strip_pysb(comp.get_id())}', real=True) for comp in getattr(self, component) ]) return @@ -1178,7 +1188,7 @@ def _generateSymbol(self, name): return elif name == 'dtcldp': self._syms[name] = sp.Matrix([ - sp.Symbol(f's{strip_pysb(tcl.get_id())}') + sp.Symbol(f's{strip_pysb(tcl.get_id())}', real=True) for tcl in self._conservationlaws ]) return @@ -1193,7 +1203,7 @@ def _generateSymbol(self, name): length = len(self.eq(name)) self._syms[name] = sp.Matrix([ - sp.Symbol(f'{name}{i}') for i in range(length) + sp.Symbol(f'{name}{i}', real=True) for i in range(length) ]) def generateBasicVariables(self): @@ -1382,6 +1392,10 @@ def _compute_equation(self, name): for index, formula in enumerate(self.eq('x0_fixedParameters')): if formula == 0 or formula == 0.0: + # sp.simplify returns ImmutableDenseMatrix, if we need to + # change them, they need to be made mutable + if isinstance(self._eqs[name], ImmutableDenseMatrix): + self._eqs[name] = MutableDenseMatrix(self._eqs[name]) self._eqs[name][index, :] = \ sp.zeros(1, self._eqs[name].shape[1]) @@ -1429,6 +1443,9 @@ def _compute_equation(self, name): if not self._lock_total_derivative: self._eqs[name] = self._eqs[name].transpose() + if self._simplify: + self._eqs[name] = self._simplify(self._eqs[name]) + def symNames(self): """Returns a list of names of generated symbolic variables @@ -1556,7 +1573,7 @@ def _total_derivative(self, name, eq, chainvars, var, if dydx_name is None: dydx_name = f'd{eq}d{chainvar}' if dxdz_name is None: - dxdz_name = f'd{chainvar}d{var}' + dxdz_name = f'd{chainvar}d{var}' dydx = self.sym_or_eq(name, dydx_name) dxdz = self.sym_or_eq(name, dxdz_name) @@ -1799,7 +1816,7 @@ class ODEExporter: modelSwigPath: path to the generated swig files @type str allow_reinit_fixpar_initcond: indicates whether reinitialization of - initial states depending on fixedParmeters is allowed for this model + initial states depending on fixedParameters is allowed for this model @type bool """ @@ -2628,7 +2645,7 @@ def strip_pysb(symbol): # this ensures that the pysb type specific __repr__ is used when converting # to string if pysb and isinstance(symbol, pysb.Component): - return sp.Symbol(symbol.name) + return sp.Symbol(symbol.name, real=True) else: # in this case we will use sympy specific transform anyways return symbol @@ -2834,7 +2851,7 @@ def csc_matrix(matrix, name, base_index=0): for row in range(0, matrix.rows): if not (matrix[row, col] == 0): symbolName = f'{name}{symbol_name_idx}' - sparseMatrix[row, col] = sp.Symbol(symbolName) + sparseMatrix[row, col] = sp.Symbol(symbolName, real=True) symbolList.append(symbolName) sparseList.append(matrix[row, col]) symbolRowVals.append(row) diff --git a/python/amici/petab_import.py b/python/amici/petab_import.py index 8a0bb8f5ba..a5c4c51a22 100644 --- a/python/amici/petab_import.py +++ b/python/amici/petab_import.py @@ -10,6 +10,7 @@ import logging from typing import List, Dict import pandas as pd +import argparse import petab from colorama import Fore @@ -81,7 +82,7 @@ def get_fixed_parameters(condition_df: pd.DataFrame, # check global parameters if not sbml_model.getParameter(fixed_parameter) \ and not sbml_model.getSpecies(fixed_parameter): - logger.log(logging.warning, + logger.log(logging.WARN, f"{Fore.YELLOW}Parameter or species '{fixed_parameter}'" " provided in condition table but not present in" " model.") @@ -286,3 +287,88 @@ def show_model_info(sbml_model: 'libsbml.Model'): + str(len(sbml_model.getListOfParameters()))) logger.log(logging.INFO, f'Reactions: {len(sbml_model.getListOfReactions())}') + + +def parse_cli_args(): + """Parse command line arguments""" + + parser = argparse.ArgumentParser( + description='Import PEtab-format model into AMICI.') + + # General options: + parser.add_argument('-v', '--verbose', dest='verbose', action='store_true', + help='More verbose output') + parser.add_argument('-o', '--output-dir', dest='model_output_dir', + help='Name of the model directory to create') + parser.add_argument('--no-compile', action='store_false', + dest='compile', + help='Only generate model code, do not compile') + + # Call with set of files + parser.add_argument('-s', '--sbml', dest='sbml_file_name', + help='SBML model filename') + parser.add_argument('-m', '--measurements', dest='measurement_file_name', + help='Measurement table') + parser.add_argument('-c', '--conditions', dest='condition_file_name', + help='Conditions table') + parser.add_argument('-p', '--parameters', dest='parameter_file_name', + help='Parameter table') + + # or with model name, following default naming + parser.add_argument('-n', '--model-name', dest='model_name', + help='Model name where all files are in the working ' + 'directory and follow PEtab naming convention. ' + 'Specifying -[smcp] will override defaults') + + args = parser.parse_args() + + if args.model_name: + if not args.sbml_file_name: + args.sbml_file_name = petab.get_default_sbml_file_name( + args.model_name) + if not args.measurement_file_name: + args.measurement_file_name = \ + petab.get_default_measurement_file_name(args.model_name) + if not args.condition_file_name: + args.condition_file_name = petab.get_default_condition_file_name( + args.model_name) + if not args.parameter_file_name: + args.parameter_file_name = petab.get_default_parameter_file_name( + args.model_name) + + if not args.model_name and \ + (not args.sbml_file_name + or not args.condition_file_name + or not args.measurement_file_name): + parser.error('When not specifying a model name, sbml, ' + 'condition and measurement file must be specified') + + return args + + +def main(): + """ + Command line interface to import a model in the PEtab + (https://github.com/ICB-DCM/PEtab/) format into sAMICI + """ + args = parse_cli_args() + + # First check for valid PEtab + pp = petab.Problem.from_files( + sbml_file=args.sbml_file_name, + condition_file=args.condition_file_name, + measurement_file=args.measurement_file_name, + parameter_file=args.parameter_file_name) + petab.lint_problem(pp) + + import_model(model_name=args.model_name, + sbml_file=args.sbml_file_name, + condition_file=args.condition_file_name, + measurement_file=args.measurement_file_name, + model_output_dir=args.model_output_dir, + compile=args.compile, + verbose=True) + + +if __name__ == '__main__': + main() diff --git a/python/amici/pysb_import.py b/python/amici/pysb_import.py index f7a0e6bd59..579f5c76d0 100644 --- a/python/amici/pysb_import.py +++ b/python/amici/pysb_import.py @@ -115,7 +115,7 @@ def ODEModel_from_pysb_importer(model, constants=None, """ - ODE = ODEModel() + ODE = ODEModel(simplify=None) if not pysb_available: raise ImportError( diff --git a/python/amici/sbml_import.py b/python/amici/sbml_import.py index 2e2a74f926..ab7b58caef 100644 --- a/python/amici/sbml_import.py +++ b/python/amici/sbml_import.py @@ -244,7 +244,7 @@ def sbml2amici(self, self.reset_symbols() self.processSBML(constantParameters) self.processObservables(observables, sigmas, noise_distributions) - ode_model = ODEModel() + ode_model = ODEModel(simplify=sp.powsimp) ode_model.import_from_sbml_importer(self) exporter = ODEExporter( ode_model, @@ -342,10 +342,21 @@ def _gather_locals(self): """ for s in self.sbml.getListOfSpecies(): - self.local_symbols[s.getId()] = sp.Symbol(s.getId()) + self.local_symbols[s.getId()] = sp.Symbol(s.getId(), real=True) for p in self.sbml.getListOfParameters(): - self.local_symbols[p.getId()] = sp.Symbol(p.getId()) + self.local_symbols[p.getId()] = sp.Symbol(p.getId(), real=True) + + for c in self.sbml.getListOfCompartments(): + self.local_symbols[c.getId()] = sp.Symbol(c.getId(), real=True) + + for r in self.sbml.getListOfRules(): + self.local_symbols[r.getVariable()] = sp.Symbol(r.getVariable(), + real=True) + + # SBML time symbol + constants + self.local_symbols['time'] = sp.Symbol('time', real=True) + self.local_symbols['avogadro'] = sp.Symbol('avogadro', real=True) def processSpecies(self): """Get species information from SBML model. @@ -365,12 +376,12 @@ def processSpecies(self): } self.symbols['species']['identifier'] = sp.Matrix( - [sp.Symbol(spec.getId()) for spec in species] + [sp.Symbol(spec.getId(), real=True) for spec in species] ) self.symbols['species']['name'] = [spec.getName() for spec in species] self.speciesCompartment = sp.Matrix( - [sp.Symbol(spec.getCompartment()) for spec in species] + [sp.Symbol(spec.getCompartment(), real=True) for spec in species] ) self.constantSpecies = [species_element.getId() @@ -448,7 +459,8 @@ def get_species_initial(index, conc): self.symbols['species']['value'] = species_initial if self.sbml.isSetConversionFactor(): - conversion_factor = sp.Symbol(self.sbml.getConversionFactor()) + conversion_factor = sp.Symbol(self.sbml.getConversionFactor(), + real=True) else: conversion_factor = 1.0 @@ -516,7 +528,7 @@ def processParameters(self, constantParameters: List[str] = None): settings = loop_settings[partype] self.symbols[partype]['identifier'] = sp.Matrix( - [sp.Symbol(par.getId()) for par in settings['var']] + [sp.Symbol(par.getId(), real=True) for par in settings['var']] ) self.symbols[partype]['name'] = [ par.getName() for par in settings['var'] @@ -547,7 +559,7 @@ def processCompartments(self): """ compartments = self.sbml.getListOfCompartments() self.compartmentSymbols = sp.Matrix( - [sp.Symbol(comp.getId()) for comp in compartments] + [sp.Symbol(comp.getId(), real=True) for comp in compartments] ) self.compartmentVolume = sp.Matrix( [sp.sympify(comp.getVolume()) if comp.isSetVolume() @@ -565,10 +577,7 @@ def processCompartments(self): locals=self.local_symbols ) - for comp, vol in zip(self.compartmentSymbols, self.compartmentVolume): - self.replaceInAllExpressions( - comp, vol - ) + def processReactions(self): """Get reactions from SBML model. @@ -620,7 +629,7 @@ def getElementStoichiometry(element): if symMath is None: symMath = sp.sympify(element.getStoichiometry()) elif element.getId() in rulevars: - return sp.Symbol(element.getId()) + return sp.Symbol(element.getId(), real=True) else: # dont put the symbol if it wont get replaced by a # rule @@ -712,7 +721,8 @@ def processRules(self): volumevars = self.compartmentVolume.free_symbols compartmentvars = self.compartmentSymbols.free_symbols parametervars = sp.Matrix([ - sp.Symbol(par.getId()) for par in self.sbml.getListOfParameters() + sp.Symbol(par.getId(), real=True) + for par in self.sbml.getListOfParameters() ]) stoichvars = self.stoichiometricMatrix.free_symbols @@ -777,6 +787,10 @@ def processRules(self): sp.sympify(variable, locals=self.local_symbols), assignments[variable] ) + for comp, vol in zip(self.compartmentSymbols, self.compartmentVolume): + self.replaceInAllExpressions( + comp, vol + ) def processVolumeConversion(self): """Convert equations from amount to volume. @@ -810,8 +824,8 @@ def processTime(self): Raises: """ - sbmlTimeSymbol = sp.Symbol('time') - amiciTimeSymbol = sp.Symbol('t') + sbmlTimeSymbol = sp.Symbol('time', real=True) + amiciTimeSymbol = sp.Symbol('t', real=True) self.replaceInAllExpressions(sbmlTimeSymbol, amiciTimeSymbol) @@ -897,17 +911,21 @@ def processObservables(self, observables: Dict[str, Dict[str, str]], observableSyms = sp.Matrix([ sp.symbols(obs, real=True) for obs in observables.keys() ]) + observable_ids = observables.keys() else: observableValues = speciesSyms - observableNames = [ + observable_ids = [ f'x{index}' for index in range(len(speciesSyms)) ] + observableNames = observable_ids[:] observableSyms = sp.Matrix( - [sp.symbols(f'y{index}', real=True) for index in range(len(speciesSyms))] + [sp.symbols(f'y{index}', real=True) + for index in range(len(speciesSyms))] ) sigmaYSyms = sp.Matrix( - [sp.symbols(f'sigma{symbol}', real=True) for symbol in observableSyms] + [sp.symbols(f'sigma{symbol}', real=True) + for symbol in observableSyms] ) sigmaYValues = sp.Matrix( [1.0] * len(observableSyms) @@ -927,7 +945,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]], # set cost functions llhYStrings = [] - for y_name in observableNames: + for y_name in observable_ids: llhYStrings.append(noise_distribution_to_cost_function( noise_distributions.get(y_name, 'normal'))) @@ -942,7 +960,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]], llhYValues = sp.Matrix(llhYValues) llhYSyms = sp.Matrix( - [sp.Symbol(f'J{symbol}') for symbol in observableSyms] + [sp.Symbol(f'J{symbol}', real=True) for symbol in observableSyms] ) # set symbols @@ -1002,8 +1020,8 @@ def cleanReservedSymbols(self): """ reservedSymbols = ['k','p','y','w'] for str in reservedSymbols: - old_symbol = sp.Symbol(str) - new_symbol = sp.Symbol('amici_' + str) + old_symbol = sp.Symbol(str, real=True) + new_symbol = sp.Symbol('amici_' + str, real=True) self.replaceInAllExpressions(old_symbol, new_symbol) for symbol in self.symbols.keys(): if 'identifier' in self.symbols[symbol].keys(): @@ -1022,7 +1040,7 @@ def replaceSpecialConstants(self): """ constants = [ - (sp.Symbol('avogadro'), sp.Symbol('6.02214179e23')), + (sp.Symbol('avogadro', real=True), sp.Symbol('6.02214179e23')), ] for constant, value in constants: # do not replace if any symbol is shadowing default definition diff --git a/python/amici/setuptools.py b/python/amici/setuptools.py index cb13caff39..0d66522c01 100644 --- a/python/amici/setuptools.py +++ b/python/amici/setuptools.py @@ -207,6 +207,8 @@ def generateSwigInterfaceFiles(): swig_exe = find_swig() swig_version = get_swig_version(swig_exe) + print(f"Found SWIG version {swig_version}") + # Swig AMICI interface without HDF5 dependency swig_cmd = [swig_exe, '-c++', diff --git a/python/bin/amici_import_petab.py b/python/bin/amici_import_petab.py deleted file mode 100755 index 2aa7174cb5..0000000000 --- a/python/bin/amici_import_petab.py +++ /dev/null @@ -1,92 +0,0 @@ -#!/usr/bin/env python3 - -""" -Command line interface to import a model in the PEtab -(https://github.com/ICB-DCM/PEtab/) format into sAMICI -""" - -import petab -import argparse - -from amici.petab_import import import_model - - -def parse_cli_args(): - """Parse command line arguments""" - - parser = argparse.ArgumentParser( - description='Import PEtab-format model into AMICI.') - - # General options: - parser.add_argument('-v', '--verbose', dest='verbose', action='store_true', - help='More verbose output') - parser.add_argument('-o', '--output-dir', dest='model_output_dir', - help='Name of the model directory to create') - parser.add_argument('--no-compile', action='store_false', - dest='compile', - help='Only generate model code, do not compile') - - # Call with set of files - parser.add_argument('-s', '--sbml', dest='sbml_file_name', - help='SBML model filename') - parser.add_argument('-m', '--measurements', dest='measurement_file_name', - help='Measurement table') - parser.add_argument('-c', '--conditions', dest='condition_file_name', - help='Conditions table') - parser.add_argument('-p', '--parameters', dest='parameter_file_name', - help='Parameter table') - - # or with model name, following default naming - parser.add_argument('-n', '--model-name', dest='model_name', - help='Model name where all files are in the working ' - 'directory and follow PEtab naming convention. ' - 'Specifying -[smcp] will override defaults') - - args = parser.parse_args() - - if args.model_name: - if not args.sbml_file_name: - args.sbml_file_name = petab.get_default_sbml_file_name( - args.model_name) - if not args.measurement_file_name: - args.measurement_file_name = \ - petab.get_default_measurement_file_name(args.model_name) - if not args.condition_file_name: - args.condition_file_name = petab.get_default_condition_file_name( - args.model_name) - if not args.parameter_file_name: - args.parameter_file_name = petab.get_default_parameter_file_name( - args.model_name) - - if not args.model_name and \ - (not args.sbml_file_name - or not args.condition_file_name - or not args.measurement_file_name): - parser.error('When not specifying a model name, sbml, ' - 'condition and measurement file must be specified') - - return args - - -def main(): - args = parse_cli_args() - - # First check for valid PEtab - pp = petab.Problem.from_files( - sbml_file=args.sbml_file_name, - condition_file=args.condition_file_name, - measurement_file=args.measurement_file_name, - parameter_file=args.parameter_file_name) - petab.lint_problem(pp) - - import_model(model_name=args.model_name, - sbml_file=args.sbml_file_name, - condition_file=args.condition_file_name, - measurement_file=args.measurement_file_name, - model_output_dir=args.model_output_dir, - compile=args.compile, - verbose=True) - - -if __name__ == '__main__': - main() diff --git a/python/sdist/setup.py b/python/sdist/setup.py index 8373809c3f..20462a32dd 100755 --- a/python/sdist/setup.py +++ b/python/sdist/setup.py @@ -177,7 +177,13 @@ def main(): ], packages=find_packages(), package_dir={'amici': 'amici'}, - scripts=['bin/amici_import_petab.py'], + entry_points={ + 'console_scripts': [ + 'amici_import_petab = amici.petab_import:main', + # for backwards compatibility + 'amici_import_petab.py = amici.petab_import:main' + ] + }, install_requires=['sympy', 'python-libsbml', 'h5py', diff --git a/src/amici.cpp b/src/amici.cpp index d71a03c9eb..15e8a98852 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -9,15 +9,16 @@ #include "amici/forwardproblem.h" #include "amici/misc.h" +#include //return codes #include //realtype -#include //return codes -#include #include +#include #include #include -#include +#include #include +#include // ensure definitions are in sync static_assert(AMICI_SUCCESS == CV_SUCCESS, "AMICI_SUCCESS != CV_SUCCESS"); @@ -32,53 +33,120 @@ static_assert(AMICI_ONE_STEP == CV_ONE_STEP, "AMICI_ONE_STEP != CV_ONE_STEP"); static_assert(std::is_same::value, "Definition of realtype does not match"); - namespace amici { -msgIdAndTxtFp errMsgIdAndTxt = &printErrMsgIdAndTxt; -msgIdAndTxtFp warnMsgIdAndTxt = &printWarnMsgIdAndTxt; +/** AMICI default application context, kept around for convenience for using + * amici::runAmiciSimulation or instantiating Solver and Model without special + * needs. + */ +AmiciApplication defaultContext = AmiciApplication(); + +std::unique_ptr +runAmiciSimulation(Solver& solver, + const ExpData* edata, + Model& model, + bool rethrow) +{ + return defaultContext.runAmiciSimulation(solver, edata, model, rethrow); +} + +void +printErrMsgIdAndTxt(std::string const& id, std::string const& message) +{ + std::cerr << "[Error] "; + if (!id.empty()) { + std::cerr << id << ": "; + } + std::cerr << message << std::endl; +} + +void +printWarnMsgIdAndTxt(std::string const& id, std::string const& message) +{ + std::cerr << "[Warning] "; + if (!id.empty()) { + std::cerr << id << ": "; + } + std::cerr << message << std::endl; +} +std::vector> +runAmiciSimulations(const Solver& solver, + const std::vector& edatas, + const Model& model, + const bool failfast, +#if defined(_OPENMP) + int num_threads +#else + int /* num_threads */ +#endif +) +{ +#if defined(_OPENMP) + return defaultContext.runAmiciSimulations( + solver, edatas, model, failfast, num_threads); +#else + return defaultContext.runAmiciSimulations(solver, edatas, model, failfast, 1); +#endif +} -std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *edata, Model &model, bool rethrow) { +std::unique_ptr +AmiciApplication::runAmiciSimulation(Solver& solver, + const ExpData* edata, + Model& model, + bool rethrow) +{ std::unique_ptr rdata; /* Applies condition-specific model settings and restores them when going * out of scope */ ConditionContext conditionContext(&model, edata); - try{ + try { rdata = std::unique_ptr(new ReturnData(solver, model)); if (model.nx_solver <= 0) { return rdata; } - auto fwd = std::unique_ptr(new ForwardProblem(rdata.get(),edata,&model,&solver)); + auto fwd = std::unique_ptr( + new ForwardProblem(rdata.get(), edata, &model, &solver)); fwd->workForwardProblem(); - auto bwd = std::unique_ptr(new BackwardProblem(fwd.get())); + auto bwd = + std::unique_ptr(new BackwardProblem(fwd.get())); bwd->workBackwardProblem(); rdata->status = AMICI_SUCCESS; } catch (amici::IntegrationFailure const& ex) { rdata->invalidate(ex.time); rdata->status = ex.error_code; - if(rethrow) throw; - amici::warnMsgIdAndTxt("AMICI:mex:simulation","AMICI forward simulation failed at t = %f:\n%s\n",ex.time,ex.what()); + if (rethrow) + throw; + warningF("AMICI:mex:simulation", + "AMICI forward simulation failed at t = %f:\n%s\n", + ex.time, + ex.what()); } catch (amici::IntegrationFailureB const& ex) { rdata->invalidateSLLH(); rdata->status = ex.error_code; - if(rethrow) throw; - amici::warnMsgIdAndTxt( - "AMICI:mex:simulation", - "AMICI backward simulation failed when trying to solve until t = %f" - " (see message above):\n%s\n", - ex.time, ex.what()); + if (rethrow) + throw; + warningF( + "AMICI:mex:simulation", + "AMICI backward simulation failed when trying to solve until t = %f" + " (see message above):\n%s\n", + ex.time, + ex.what()); } catch (amici::AmiException const& ex) { rdata->invalidate(model.t0()); rdata->status = AMICI_ERROR; - if(rethrow) throw; - amici::warnMsgIdAndTxt("AMICI:mex:simulation","AMICI simulation failed:\n%s\nError occured in:\n%s",ex.what(),ex.getBacktrace()); + if (rethrow) + throw; + warningF("AMICI:mex:simulation", + "AMICI simulation failed:\n%s\nError occured in:\n%s", + ex.what(), + ex.getBacktrace()); } rdata->applyChainRuleFactorToSimulationResults(&model); @@ -86,50 +154,28 @@ std::unique_ptr runAmiciSimulation(Solver &solver, const ExpData *ed return rdata; } -void printErrMsgIdAndTxt(const char *identifier, const char *format, ...) { - if(identifier && *identifier != '\0') - fprintf(stderr, "[Error] %s: ", identifier); - else - fprintf(stderr, "[Error] "); - va_list argptr; - va_start(argptr,format); - vfprintf(stderr, format, argptr); - va_end(argptr); - fprintf(stderr, "\n"); -} - -void printWarnMsgIdAndTxt(const char *identifier, const char *format, ...) { - if(identifier && *identifier != '\0') - printf("[Warning] %s: ", identifier); - else - printf("[Warning] "); - va_list argptr; - va_start(argptr,format); - vprintf(format, argptr); - va_end(argptr); - printf("\n"); -} - -std::vector > runAmiciSimulations(const Solver &solver, - const std::vector &edatas, - const Model &model, - const bool failfast, +std::vector> +AmiciApplication::runAmiciSimulations(const Solver& solver, + const std::vector& edatas, + const Model& model, + bool failfast, #if defined(_OPENMP) - int num_threads + int num_threads #else - int /* num_threads */ + int /* num_threads */ #endif + ) { - std::vector > results(edatas.size()); + std::vector> results(edatas.size()); // is set to true if one simulation fails and we should skip the rest. // shared across threads. bool skipThrough = false; #if defined(_OPENMP) - #pragma omp parallel for num_threads(num_threads) +#pragma omp parallel for num_threads(num_threads) #endif - for(int i = 0; i < (int)edatas.size(); ++i) { + for (int i = 0; i < (int)edatas.size(); ++i) { auto mySolver = std::unique_ptr(solver.clone()); auto myModel = std::unique_ptr(model.clone()); @@ -138,7 +184,7 @@ std::vector > runAmiciSimulations(const Solver &solv if (skipThrough) { ConditionContext conditionContext(myModel.get(), edatas[i]); results[i] = - std::unique_ptr(new ReturnData(solver, model)); + std::unique_ptr(new ReturnData(solver, model)); } else { results[i] = runAmiciSimulation(*mySolver, edatas[i], *myModel); } @@ -149,4 +195,49 @@ std::vector > runAmiciSimulations(const Solver &solv return results; } +void +AmiciApplication::warningF(const char* identifier, const char* format, ...) +{ + va_list argptr; + va_start(argptr, format); + auto str = printfToString(format, argptr); + va_end(argptr); + warning(identifier, str); +} + +void +AmiciApplication::errorF(const char* identifier, const char* format, ...) +{ + va_list argptr; + va_start(argptr, format); + auto str = printfToString(format, argptr); + va_end(argptr); + error(identifier, str); +} + +int +AmiciApplication::checkFinite(gsl::span array, const char* fun) +{ + + for (int idx = 0; idx < (int)array.size(); idx++) { + if (isNaN(array[idx])) { + warningF("AMICI:NaN", + "AMICI encountered a NaN value at index %i of %i in %s!", + idx, + (int)array.size(), + fun); + return AMICI_RECOVERABLE_ERROR; + } + if (isInf(array[idx])) { + warningF("AMICI:Inf", + "AMICI encountered an Inf value at index %i of %i in %s!", + idx, + (int)array.size(), + fun); + return AMICI_RECOVERABLE_ERROR; + } + } + return AMICI_SUCCESS; +} + } // namespace amici diff --git a/src/forwardproblem.cpp b/src/forwardproblem.cpp index 52112e5670..5b8ff882f3 100644 --- a/src/forwardproblem.cpp +++ b/src/forwardproblem.cpp @@ -15,9 +15,6 @@ namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; - - ForwardProblem::ForwardProblem(ReturnData *rdata, const ExpData *edata, Model *model, Solver *solver) : model(model), @@ -310,7 +307,8 @@ void ForwardProblem::handleEvent(realtype *tlastroot, const bool seflag) { discs[iroot] = t; ++iroot; } else { - warnMsgIdAndTxt("AMICI:mex:TOO_MUCH_EVENT", + solver->app->warning( + "AMICI:mex:TOO_MUCH_EVENT", "Event was recorded but not reported as the number of " "occured events exceeded (nmaxevents)*(number of " "events in model definition)!"); diff --git a/src/hdf5.cpp b/src/hdf5.cpp index 510703a0fe..003b1460d0 100644 --- a/src/hdf5.cpp +++ b/src/hdf5.cpp @@ -628,6 +628,17 @@ void readSolverSettingsFromHDF5(H5::H5File const& file, Solver &solver, getIntScalarAttribute(file, datasetPath, "newton_preeq")); } + if(attributeExists(file, datasetPath, "newton_damping_factor_mode")) { + solver.setNewtonDampingFactorMode( + static_cast( + getIntScalarAttribute(file, datasetPath, "newton_damping_factor_mode"))); + } + + if(attributeExists(file, datasetPath, "newton_damping_factor_lower_bound")) { + solver.setNewtonDampingFactorLowerBound( + getDoubleScalarAttribute(file, datasetPath, "newton_damping_factor_lower_bound")); + } + if(attributeExists(file, datasetPath, "newton_maxlinsteps")) { solver.setNewtonMaxLinearSteps( getIntScalarAttribute(file, datasetPath, diff --git a/src/interface_matlab.cpp b/src/interface_matlab.cpp index 5f785d37ee..800bbc061c 100644 --- a/src/interface_matlab.cpp +++ b/src/interface_matlab.cpp @@ -462,17 +462,31 @@ void setModelData(const mxArray *prhs[], int nrhs, Model &model) */ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { // use matlab error reporting - amici::warnMsgIdAndTxt = &mexWarnMsgIdAndTxt; - amici::errMsgIdAndTxt = &mexErrMsgIdAndTxt; + amici::AmiciApplication amiciApp; + amiciApp.warning = []( + std::string const& identifier, + std::string const& message){ + mexWarnMsgIdAndTxt(identifier.c_str(), message.c_str()); + }; + amiciApp.error = []( + std::string const& identifier, + std::string const& message){ + mexErrMsgIdAndTxt(identifier.c_str(), message.c_str()); + }; if (nlhs != 1) { - amici::errMsgIdAndTxt("AMICI:mex:setup","Incorrect number of output arguments (must be 1)!"); + amiciApp.errorF("AMICI:mex:setup", + "Incorrect number of output arguments (must be 1)!"); } else if(nrhs < amici::RHS_NUMARGS_REQUIRED) { - amici::errMsgIdAndTxt("AMICI:mex:setup", "Incorrect number of input arguments (must be at least 7)!"); + amiciApp.errorF("AMICI:mex:setup", + "Incorrect number of input arguments (must be at least 7)!"); }; auto model = getModel(); + model->app = &amiciApp; + auto solver = model->getSolver(); + solver->app = &amiciApp; setModelData(prhs, nrhs, *model); setSolverOptions(prhs, nrhs, *solver); @@ -481,11 +495,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { try { edata = amici::expDataFromMatlabCall(prhs, *model); } catch (amici::AmiException const& ex) { - amici::errMsgIdAndTxt("AMICI:mex:setup","Failed to read experimental data:\n%s",ex.what()); + amiciApp.errorF("AMICI:mex:setup","Failed to read experimental data:\n%s",ex.what()); } } else if (solver->getSensitivityOrder() >= amici::SensitivityOrder::first && solver->getSensitivityMethod() == amici::SensitivityMethod::adjoint) { - amici::errMsgIdAndTxt("AMICI:mex:setup","No data provided!"); + amiciApp.errorF("AMICI:mex:setup","No data provided!"); } /* ensures that plhs[0] is available */ diff --git a/src/misc.cpp b/src/misc.cpp index 3876c7cbfd..731559c96b 100644 --- a/src/misc.cpp +++ b/src/misc.cpp @@ -5,6 +5,8 @@ #include #include #include +#include + #if defined(_WIN32) #define PLATFORM_WINDOWS // Windows #elif defined(_WIN64) @@ -86,27 +88,6 @@ void scaleParameters(gsl::span bufferUnscaled, } -int checkFinite(gsl::span array, const char *fun) -{ - for (int idx = 0; idx < (int) array.size(); idx++) { - if (isNaN(array[idx])) { - warnMsgIdAndTxt( - "AMICI:NaN", - "AMICI encountered a NaN value at index %i of %i in %s!", idx, - (int) array.size(), fun); - return AMICI_RECOVERABLE_ERROR; - } - if (isInf(array[idx])) { - warnMsgIdAndTxt( - "AMICI:Inf", - "AMICI encountered an Inf value at index %i of %i in %s!", idx, - (int) array.size(), fun); - return AMICI_RECOVERABLE_ERROR; - } - } - return AMICI_SUCCESS; -} - std::string backtraceString(const int maxFrames) { std::ostringstream trace_buf; @@ -186,4 +167,21 @@ std::string regexErrorToString(std::regex_constants::error_type err_type) } } +std::string printfToString(const char *fmt, va_list ap) { + // Get size of string + va_list ap_count; + va_copy(ap_count, ap); + auto size = vsnprintf(nullptr, 0, fmt, ap_count); + va_end(ap_count); + ++size; + + // actual formatting + auto buf = new char[size]; + size = vsnprintf(buf, size, fmt, ap); + std::string str(buf, size); + delete[] buf; + + return str; +} + } // namespace amici diff --git a/src/model.cpp b/src/model.cpp index 9a35c36f83..24c14ee728 100644 --- a/src/model.cpp +++ b/src/model.cpp @@ -1002,7 +1002,7 @@ void Model::addStateEventUpdate(AmiVector &x, const int ie, const realtype t, fixedParameters.data(), h.data(), ie, xdot.data(), xdot_old.data()); if (alwaysCheckFinite) { - amici::checkFinite(deltax, "deltax"); + app->checkFinite(deltax, "deltax"); } // update @@ -1027,7 +1027,7 @@ void Model::addStateSensitivityEventUpdate(AmiVectorArray &sx, const int ie, xdot.data(), xdot_old.data(), sx.data(ip), &stau.at(ip)); if (alwaysCheckFinite) { - amici::checkFinite(deltasx, "deltasx"); + app->checkFinite(deltasx, "deltasx"); } amici_daxpy(nx_solver, 1.0, deltasx.data(), 1, sx.data(ip), 1); @@ -1047,7 +1047,7 @@ void Model::addAdjointStateEventUpdate(AmiVector &xB, const int ie, xB.data()); if (alwaysCheckFinite) { - amici::checkFinite(deltaxB, "deltaxB"); + app->checkFinite(deltaxB, "deltaxB"); } // apply update @@ -1072,7 +1072,7 @@ void Model::addAdjointQuadratureEventUpdate( } if (alwaysCheckFinite) { - amici::checkFinite(deltaqB, "deltaqB"); + app->checkFinite(deltaqB, "deltaqB"); } } @@ -1089,12 +1089,12 @@ void Model::updateHeavisideB(const int *rootsfound) { } int Model::checkFinite(gsl::span array, const char *fun) const { - auto result = amici::checkFinite(array, fun); + auto result = app->checkFinite(array, fun); if (result != AMICI_SUCCESS) { - amici::checkFinite(fixedParameters, "k"); - amici::checkFinite(unscaledParameters, "p"); - amici::checkFinite(w, "w"); + app->checkFinite(fixedParameters, "k"); + app->checkFinite(unscaledParameters, "p"); + app->checkFinite(w, "w"); } return result; @@ -1240,7 +1240,7 @@ void Model::fy(const realtype t, const AmiVector &x) { h.data(), w.data()); if (alwaysCheckFinite) { - amici::checkFinite(gsl::make_span(y.data(), ny), "y"); + app->checkFinite(gsl::make_span(y.data(), ny), "y"); } } @@ -1265,7 +1265,7 @@ void Model::fdydp(const realtype t, const AmiVector &x) { } if (alwaysCheckFinite) { - amici::checkFinite(dydp, "dydp"); + app->checkFinite(dydp, "dydp"); } } @@ -1281,7 +1281,7 @@ void Model::fdydx(const realtype t, const AmiVector &x) { fixedParameters.data(), h.data(), w.data(), dwdx.data()); if (alwaysCheckFinite) { - amici::checkFinite(dydx, "dydx"); + app->checkFinite(dydx, "dydx"); } } @@ -1339,7 +1339,7 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) { } if (alwaysCheckFinite) { - amici::checkFinite(dsigmaydp, "dsigmaydp"); + app->checkFinite(dsigmaydp, "dsigmaydp"); } } @@ -1375,7 +1375,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { edata.getObservedDataPtr(it)); if (alwaysCheckFinite) { - amici::checkFinite(gsl::make_span(dJydy[iyt].get()), "dJydy"); + app->checkFinite(gsl::make_span(dJydy[iyt].get()), "dJydy"); } } } else { @@ -1389,7 +1389,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { } if (alwaysCheckFinite) { // get dJydy slice (ny) for current timepoint and observable - amici::checkFinite(dJydy_matlab, "dJydy"); + app->checkFinite(dJydy_matlab, "dJydy"); } } } @@ -1410,7 +1410,7 @@ void Model::fdJydsigma(const int it, const AmiVector &x, const ExpData &edata) { } if (alwaysCheckFinite) { - amici::checkFinite(dJydsigma, "dJydsigma"); + app->checkFinite(dJydsigma, "dJydsigma"); } } @@ -1486,7 +1486,7 @@ void Model::fdJydx(const int it, const AmiVector &x, const ExpData &edata) { } if (alwaysCheckFinite) { - amici::checkFinite(dJydx, "dJydx"); + app->checkFinite(dJydx, "dJydx"); } } @@ -1508,7 +1508,7 @@ void Model::fdzdp(const int ie, const realtype t, const AmiVector &x) { } if (alwaysCheckFinite) { - amici::checkFinite(dzdp, "dzdp"); + app->checkFinite(dzdp, "dzdp"); } } @@ -1520,7 +1520,7 @@ void Model::fdzdx(const int ie, const realtype t, const AmiVector &x) { fixedParameters.data(), h.data()); if (alwaysCheckFinite) { - amici::checkFinite(dzdx, "dzdx"); + app->checkFinite(dzdx, "dzdx"); } } @@ -1542,7 +1542,7 @@ void Model::fdrzdp(const int ie, const realtype t, const AmiVector &x) { } if (alwaysCheckFinite) { - amici::checkFinite(drzdp, "drzdp"); + app->checkFinite(drzdp, "drzdp"); } } @@ -1554,7 +1554,7 @@ void Model::fdrzdx(const int ie, const realtype t, const AmiVector &x) { fixedParameters.data(), h.data()); if (alwaysCheckFinite) { - amici::checkFinite(drzdx, "drzdx"); + app->checkFinite(drzdx, "drzdx"); } } @@ -1613,7 +1613,7 @@ void Model::fdsigmazdp(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dsigmazdp, "dsigmazdp"); + app->checkFinite(dsigmazdp, "dsigmazdp"); } } @@ -1634,7 +1634,7 @@ void Model::fdJzdz(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dJzdz, "dJzdz"); + app->checkFinite(dJzdz, "dJzdz"); } } @@ -1656,7 +1656,7 @@ void Model::fdJzdsigma(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dJzdsigma, "dJzdsigma"); + app->checkFinite(dJzdsigma, "dJzdsigma"); } } @@ -1759,7 +1759,7 @@ void Model::fdJrzdz(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dJrzdz, "dJrzdz"); + app->checkFinite(dJrzdz, "dJrzdz"); } } @@ -1780,7 +1780,7 @@ void Model::fdJrzdsigma(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dJrzdsigma, "dJrzdsigma"); + app->checkFinite(dJrzdsigma, "dJrzdsigma"); } } @@ -1790,7 +1790,7 @@ void Model::fw(const realtype t, const realtype *x) { h.data(), total_cl.data()); if (alwaysCheckFinite) { - amici::checkFinite(w, "w"); + app->checkFinite(w, "w"); } } @@ -1819,7 +1819,7 @@ void Model::fdwdp(const realtype t, const realtype *x) { } if (alwaysCheckFinite) { - amici::checkFinite(dwdp, "dwdp"); + app->checkFinite(dwdp, "dwdp"); } } @@ -1832,7 +1832,7 @@ void Model::fdwdx(const realtype t, const realtype *x) { fdwdx_rowvals(dwdx.indexptrs()); if (alwaysCheckFinite) { - amici::checkFinite(gsl::make_span(dwdx.get()), "dwdx"); + app->checkFinite(gsl::make_span(dwdx.get()), "dwdx"); } } diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index 26d1247baa..ca86694b44 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -28,7 +28,8 @@ NewtonSolver::NewtonSolver(realtype *t, AmiVector *x, Model *model, std::unique_ptr NewtonSolver::getSolver( realtype *t, AmiVector *x, LinearSolver linsolType, Model *model, - ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol) { + ReturnData *rdata, int maxlinsteps, int maxsteps, double atol, double rtol, + NewtonDampingFactorMode dampingFactorMode, double dampingFactorLowerBound) { std::unique_ptr solver; @@ -76,6 +77,8 @@ std::unique_ptr NewtonSolver::getSolver( solver->rtol = rtol; solver->maxlinsteps = maxlinsteps; solver->maxsteps = maxsteps; + solver->dampingFactorMode = dampingFactorMode; + solver->dampingFactorLowerBound = dampingFactorLowerBound; return solver; } diff --git a/src/solver.cpp b/src/solver.cpp index cad80533a6..f86a896fa9 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -12,7 +12,10 @@ namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; +Solver::Solver(AmiciApplication *app) : app(app) +{ + +} Solver::Solver(const Solver &other) : ism(other.ism), lmm(other.lmm), iter(other.iter), @@ -20,6 +23,8 @@ Solver::Solver(const Solver &other) sensi_meth(other.sensi_meth), stldet(other.stldet), ordering(other.ordering), newton_maxsteps(other.newton_maxsteps), newton_maxlinsteps(other.newton_maxlinsteps), + newton_damping_factor_mode(other.newton_damping_factor_mode), + newton_damping_factor_lower_bound(other.newton_damping_factor_lower_bound), requires_preequilibration(other.requires_preequilibration), linsol(other.linsol), atol(other.atol), rtol(other.rtol), atol_fsa(other.atol_fsa), rtol_fsa(other.rtol_fsa), atolB(other.atolB), rtolB(other.rtolB), @@ -378,6 +383,8 @@ bool operator==(const Solver &a, const Solver &b) { (a.ordering == b.ordering) && (a.newton_maxsteps == b.newton_maxsteps) && (a.newton_maxlinsteps == b.newton_maxlinsteps) && + (a.newton_damping_factor_mode == b.newton_damping_factor_mode) && + (a.newton_damping_factor_lower_bound == b.newton_damping_factor_lower_bound) && (a.requires_preequilibration == b.requires_preequilibration) && (a.ism == b.ism) && (a.linsol == b.linsol) && (a.atol == b.atol) && (a.rtol == b.rtol) && (a.maxsteps == b.maxsteps) && (a.maxstepsB == b.maxstepsB) && @@ -493,6 +500,18 @@ void Solver::setNewtonMaxLinearSteps(const int newton_maxlinsteps) { this->newton_maxlinsteps = newton_maxlinsteps; } +NewtonDampingFactorMode Solver::getNewtonDampingFactorMode() const { return newton_damping_factor_mode; } + +void Solver::setNewtonDampingFactorMode(NewtonDampingFactorMode dampingFactorMode) { + this->newton_damping_factor_mode = dampingFactorMode; +} + +double Solver::getNewtonDampingFactorLowerBound() const { return newton_damping_factor_lower_bound; } + +void Solver::setNewtonDampingFactorLowerBound(double dampingFactorLowerBound) { + this->newton_damping_factor_lower_bound = dampingFactorLowerBound; +} + SensitivityOrder Solver::getSensitivityOrder() const { return sensi; } void Solver::setSensitivityOrder(const SensitivityOrder sensi) { @@ -1000,38 +1019,44 @@ const AmiVector &Solver::getAdjointQuadrature(const int which, realtype Solver::gett() const { return t; } void wrapErrHandlerFn(int error_code, const char *module, - const char *function, char *msg, void * /*eh_data*/) { - char buffer[250]; - char buffid[250]; - sprintf(buffer, "AMICI ERROR: in module %s in function %s : %s ", module, + const char *function, char *msg, void * eh_data) { +#define BUF_SIZE 250 + char buffer[BUF_SIZE]; + char buffid[BUF_SIZE]; + snprintf(buffer, BUF_SIZE, "AMICI ERROR: in module %s in function %s : %s ", module, function, msg); switch (error_code) { case 99: - sprintf(buffid, "AMICI:mex:%s:%s:WARNING", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:WARNING", module, function); break; case -1: - sprintf(buffid, "AMICI:mex:%s:%s:TOO_MUCH_WORK", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:TOO_MUCH_WORK", module, function); break; case -2: - sprintf(buffid, "AMICI:mex:%s:%s:TOO_MUCH_ACC", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:TOO_MUCH_ACC", module, function); break; case -3: - sprintf(buffid, "AMICI:mex:%s:%s:ERR_FAILURE", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:ERR_FAILURE", module, function); break; case -4: - sprintf(buffid, "AMICI:mex:%s:%s:CONV_FAILURE", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:CONV_FAILURE", module, function); break; default: - sprintf(buffid, "AMICI:mex:%s:%s:OTHER", module, function); + snprintf(buffid, BUF_SIZE, "AMICI:mex:%s:%s:OTHER", module, function); break; } - warnMsgIdAndTxt(buffid, buffer); + + if(!eh_data) { + throw std::runtime_error("eh_data unset"); + } + auto solver = static_cast(eh_data); + solver->app->warning(buffid, buffer); } } // namespace amici diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index 56441f7828..34fbfe9ef1 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -307,7 +307,9 @@ void CVodeSolver::setNonLinearSolverB(const int which) const { void CVodeSolver::setErrHandlerFn() const { int status = - CVodeSetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, nullptr); + CVodeSetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, + reinterpret_cast( + const_cast(this))); if (status != CV_SUCCESS) throw CvodeException(status, "CVodeSetErrHandlerFn"); } @@ -342,9 +344,9 @@ void CVodeSolver::setStabLimDetB(const int which, const int stldet) const { throw CvodeException(status, "CVodeSetStabLimDetB"); } -void CVodeSolver::setId(const Model *model) const {} +void CVodeSolver::setId(const Model */*model*/) const {} -void CVodeSolver::setSuppressAlg(const bool flag) const {} +void CVodeSolver::setSuppressAlg(const bool /*flag*/) const {} void CVodeSolver::resetState(void *ami_mem, const_N_Vector y0) const { @@ -570,7 +572,7 @@ void CVodeSolver::allocateSolverB(int *which) const { solverMemoryB.resize(*which + 1); solverMemoryB.at(*which) = std::unique_ptr>( - getAdjBmem(solverMemory.get(), *which), [](void *ptr) {}); + getAdjBmem(solverMemory.get(), *which), [](void */*ptr*/) {}); if (status != CV_SUCCESS) throw CvodeException(status, "CVodeCreateB"); } @@ -686,9 +688,9 @@ void *CVodeSolver::getAdjBmem(void *ami_mem, int which) const { return CVodeGetAdjCVodeBmem(ami_mem, which); } -void CVodeSolver::calcIC(const realtype tout1) const {}; +void CVodeSolver::calcIC(const realtype /*tout1*/) const {}; -void CVodeSolver::calcICB(const int which, const realtype tout1) const {}; +void CVodeSolver::calcICB(const int /*which*/, const realtype /*tout1*/) const {}; void CVodeSolver::setStopTime(const realtype tstop) const { int status = CVodeSetStopTime(solverMemory.get(), tstop); @@ -903,7 +905,7 @@ int froot(realtype t, N_Vector x, realtype *root, int fxdot(realtype t, N_Vector x, N_Vector xdot, void *user_data) { auto model = static_cast(user_data); - if (t > 1e200 && !amici::checkFinite(gsl::make_span(x), "fxdot")) { + if (t > 1e200 && !model->checkFinite(gsl::make_span(x), "fxdot")) { /* 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/solver_idas.cpp b/src/solver_idas.cpp index 670b22805b..a8b9b51e4f 100644 --- a/src/solver_idas.cpp +++ b/src/solver_idas.cpp @@ -244,7 +244,9 @@ void IDASolver::getRootInfo(int *rootsfound) const { void IDASolver::setErrHandlerFn() const { int status = - IDASetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, nullptr); + IDASetErrHandlerFn(solverMemory.get(), wrapErrHandlerFn, + reinterpret_cast( + const_cast(this))); if (status != IDA_SUCCESS) throw IDAException(status, "IDASetErrHandlerFn"); } @@ -267,9 +269,10 @@ void IDASolver::setMaxNumSteps(const long int mxsteps) const { throw IDAException(status, "IDASetMaxNumSteps"); } -void IDASolver::setStabLimDet(const int stldet) const {} +void IDASolver::setStabLimDet(const int /*stldet*/) const {} -void IDASolver::setStabLimDetB(const int which, const int stldet) const {} +void IDASolver::setStabLimDetB(const int /*which*/, + const int /*stldet*/) const {} void IDASolver::setId(const Model *model) const { @@ -511,7 +514,7 @@ void IDASolver::allocateSolverB(int *which) const { solverMemoryB.resize(*which + 1); solverMemoryB.at(*which) = std::unique_ptr>( - getAdjBmem(solverMemory.get(), *which), [](void *ptr) {}); + getAdjBmem(solverMemory.get(), *which), [](void */*ptr*/) {}); if (status != IDA_SUCCESS) throw IDAException(status, "IDACreateB"); } @@ -938,7 +941,7 @@ int fxdot(realtype t, N_Vector x, N_Vector dx, N_Vector xdot, void *user_data) { auto model = static_cast(user_data); - if (t > 1e200 && !amici::checkFinite(gsl::make_span(x), "fxdot")) { + if (t > 1e200 && !model->app->checkFinite(gsl::make_span(x), "fxdot")) { /* 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 8ea73c352a..3b3890d77d 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -16,7 +16,7 @@ #include namespace amici { - + SteadystateProblem::SteadystateProblem(const Solver *solver, const AmiVector &x0): t(solver->gett()), delta(solver->nx()), ewt(solver->nx()), @@ -48,7 +48,9 @@ void SteadystateProblem::workSteadyStateProblem(ReturnData *rdata, auto newtonSolver = NewtonSolver::getSolver( &t, &x, solver->getLinearSolver(), model, rdata, solver->getNewtonMaxLinearSteps(), solver->getNewtonMaxSteps(), - solver->getAbsoluteTolerance(), solver->getRelativeTolerance()); + solver->getAbsoluteTolerance(), solver->getRelativeTolerance(), + solver->getNewtonDampingFactorMode(), + solver->getNewtonDampingFactorLowerBound()); auto newton_status = NewtonStatus::failed; try { @@ -207,20 +209,30 @@ void SteadystateProblem::applyNewtonsMethod(ReturnData *rdata, Model *model, converged = wrms < RCONST(1.0); if (converged) { - /* Ensure positivity of the found state */ + /* Ensure positivity of the found state and recheck if + the convergence still holds */ + bool recheck_convergence = false; for (ix = 0; ix < model->nx_solver; ix++) { if (x[ix] < 0.0) { x[ix] = 0.0; - converged = FALSE; + recheck_convergence = true; } } - } else { + if (recheck_convergence) { + model->fxdot(t, x, dx, xdot); + wrms = getWrmsNorm(x_newton, xdot, newtonSolver->atol, newtonSolver->rtol); + converged = wrms < RCONST(1.0); + } + } else if (newtonSolver->dampingFactorMode==NewtonDampingFactorMode::on) { /* increase dampening factor (superfluous, if converged) */ gamma = fmin(1.0, 2.0 * gamma); } - } else { - /* Reduce dampening factor */ + } else if (newtonSolver->dampingFactorMode==NewtonDampingFactorMode::on) { + /* Reduce dampening factor and raise an error when becomes too small */ gamma = gamma / 4.0; + if (gamma < newtonSolver->dampingFactorLowerBound) + throw AmiException("Newton solver failed: a damping factor reached its lower bound"); + /* No new linear solve, only try new dampening */ compNewStep = FALSE; } @@ -325,7 +337,7 @@ std::unique_ptr SteadystateProblem::createSteadystateSimSolver( return newton_solver; } - + void SteadystateProblem::writeSolution(realtype *t, AmiVector &x, AmiVectorArray &sx) const { *t = this->t; diff --git a/src/sundials_matrix_wrapper.cpp b/src/sundials_matrix_wrapper.cpp index dd107c5efa..efd17ba9da 100644 --- a/src/sundials_matrix_wrapper.cpp +++ b/src/sundials_matrix_wrapper.cpp @@ -170,28 +170,31 @@ void SUNMatrixWrapper::multiply(gsl::span c, gsl::span if (!matrix) return; - if (static_cast(c.size()) != rows()) + sunindextype nrows = rows(); + sunindextype ncols = columns(); + + if (static_cast(c.size()) != nrows) throw std::invalid_argument("Dimension mismatch between number of rows " - "in A (" + std::to_string(rows()) + ") and " + "in A (" + std::to_string(nrows) + ") and " "elements in c (" + std::to_string(c.size()) + ")"); - if (static_cast(b.size()) != columns()) + if (static_cast(b.size()) != ncols) throw std::invalid_argument("Dimension mismatch between number of cols " - "in A (" + std::to_string(columns()) + "in A (" + std::to_string(ncols) + ") and elements in b (" + std::to_string(b.size()) + ")"); switch (SUNMatGetID(matrix)) { case SUNMATRIX_DENSE: - amici_dgemv(BLASLayout::colMajor, BLASTranspose::noTrans, rows(), - columns(), 1.0, data(), rows(), b.data(), 1, 1.0, c.data(), 1); + amici_dgemv(BLASLayout::colMajor, BLASTranspose::noTrans, nrows, + ncols, 1.0, data(), nrows, b.data(), 1, 1.0, c.data(), 1); break; case SUNMATRIX_SPARSE: switch (sparsetype()) { case CSC_MAT: - for (sunindextype i = 0; i < columns(); ++i) { + for (sunindextype i = 0; i < ncols; ++i) { for (sunindextype k = indexptrs_ptr[i]; k < indexptrs_ptr[i + 1]; ++k) { c[indexvals_ptr[k]] += data_ptr[k] * b[i]; @@ -199,7 +202,7 @@ void SUNMatrixWrapper::multiply(gsl::span c, gsl::span } break; case CSR_MAT: - for (sunindextype i = 0; i < rows(); ++i) { + for (sunindextype i = 0; i < nrows; ++i) { for (sunindextype k = indexptrs_ptr[i]; k < indexptrs_ptr[i + 1]; ++k) { c[i] += data_ptr[k] * b[indexvals_ptr[k]]; diff --git a/swig/amici.i b/swig/amici.i index f57eeb8e50..359520677c 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -87,8 +87,8 @@ wrap_unique_ptr(ExpDataPtr, amici::ExpData) // Add necessary symbols to generated header // Ignore due to https://github.com/swig/swig/issues/1643 -%ignore printErrMsgIdAndTxt; -%ignore printWarnMsgIdAndTxt; +%ignore amici::AmiciApplication::warningF; +%ignore amici::AmiciApplication::errorF; %{ #include "amici/amici.h" using namespace amici; diff --git a/tests/performance/check_time.sh b/tests/performance/check_time.sh new file mode 100755 index 0000000000..d8f440e1f2 --- /dev/null +++ b/tests/performance/check_time.sh @@ -0,0 +1,32 @@ +#!/bin/bash +# Execute command and compare wall time to a reference +set -e + +SCRIPT_PATH=$(dirname "${BASH_SOURCE}") +SCRIPT_PATH=$(cd "${SCRIPT_PATH}" && pwd) + +# YAML file with reference wall times +REFERENCE_FILE="${SCRIPT_PATH}"/reference.yml + +# Reference key +REF="$1" +# Command to time +CMD="${@:2}" +# Logfile +LOG=$(tempfile) + +# Run and time +/usr/bin/time -f %e ${CMD} 2>&1 | tee "$LOG" +RET=${PIPESTATUS[0]} +test "$RET" != 0 && exit "$RET" +TIME=$(tail -n1 "$LOG") + +# Read reference +REF_TIME=$(shyaml get-value $REF < $REFERENCE_FILE) + +if (( $(echo "$TIME > $REF_TIME" | bc) )); then + echo "TOO LONG: $TIME > $REF_TIME" + exit 1 +else + echo "OKAY: $TIME < $REF_TIME" +fi diff --git a/tests/performance/reference.yml b/tests/performance/reference.yml new file mode 100644 index 0000000000..b392b8bb6f --- /dev/null +++ b/tests/performance/reference.yml @@ -0,0 +1,8 @@ +# Reference wall times (seconds) with some buffer +create_sdist: 25 +install_sdist: 120 +petab_import: 1680 # ^= 28' +install_model: 400 +forward_simulation: 2 +forward_sensitivities: 6 +adjoint_sensitivities: 120 # TODO reduce when #849 fixed diff --git a/tests/performance/test.py b/tests/performance/test.py new file mode 100755 index 0000000000..8ca43c4733 --- /dev/null +++ b/tests/performance/test.py @@ -0,0 +1,42 @@ +#!/usr/bin/env python3 + +import amici +import sys + +import CS_Signalling_ERBB_RAS_AKT_petab as model_module + + +def main(): + arg = sys.argv[1] + + model = model_module.getModel() + solver = model.getSolver() + # TODO + edata = amici.ExpData(model) + edata.setTimepoints([1e8]) + edata.setObservedData([1.0]) + + if arg == 'forward_simulation': + solver.setSensitivityMethod(amici.SensitivityMethod_none) + solver.setSensitivityOrder(amici.SensitivityOrder_none) + elif arg == 'forward_sensitivities': + solver.setSensitivityMethod(amici.SensitivityMethod_forward) + solver.setSensitivityOrder(amici.SensitivityOrder_first) + elif arg == 'adjoint_sensitivities': + solver.setSensitivityMethod(amici.SensitivityMethod_adjoint) + solver.setSensitivityOrder(amici.SensitivityOrder_first) + else: + print("Unknown argument:", arg) + sys.exit(1) + rdata = amici.runAmiciSimulation(model, solver, edata) + + diagnostics = ['numsteps', 'numstepsB', 'numrhsevals', 'numrhsevalsB', + 'numerrtestfails', 'numerrtestfailsB', + 'numnonlinsolvconvfails', 'numnonlinsolvconvfailsB'] + for d in diagnostics: + print(d, rdata[d]) + assert rdata['status'] == amici.AMICI_SUCCESS + + +if __name__ == '__main__': + main() diff --git a/tests/testMisc.py b/tests/testMisc.py index 877afceec4..76c7eecc53 100755 --- a/tests/testMisc.py +++ b/tests/testMisc.py @@ -153,6 +153,26 @@ def test_constant_species_to_parameters(self): assert len(list(r.getListOfProducts())) == 0 assert len(list(r.getListOfModifiers())) == 0 + def test_hill_function_dwdx(self): + """Kinetic laws with Hill functions, may lead to NaNs in the Jacobian + if involved states are zero if not properly arranged symbolically. + Test that what we are applying the right sympy simplification.""" + + w = sp.sympify('Pow(x1, p1) / (Pow(x1, p1) + a)') + dwdx = w.diff(sp.Symbol('x1')) + + # Verify that without simplification we fail + with self.assertRaises(ZeroDivisionError): + with sp.evaluate(False): + res = dwdx.subs({'x1': 0.0}) + _ = str(res) + + # Test that powsimp does the job + dwdx = sp.powsimp(dwdx) + with sp.evaluate(False): + res = dwdx.subs({'x1': 0.0}) + _ = str(res) + if __name__ == '__main__': suite = unittest.TestSuite() diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index 0b8c07c52b..0a1816c8e9 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -118,7 +118,7 @@ def concentrations_to_amounts(amount_species, wrapper, model, simulated_x): ) }) volume = volume.subs({ - sp.Symbol(name): value + sp.Symbol(name, real=True): value for name, value in zip( model.getParameterIds(), model.getParameters() @@ -127,7 +127,7 @@ def concentrations_to_amounts(amount_species, wrapper, model, simulated_x): # required for 525-527, 530 as k is renamed to amici_k volume = volume.subs({ - sp.Symbol(name): value + sp.Symbol(name, real=True): value for name, value in zip( model.getParameterNames(), model.getParameters() diff --git a/version.txt b/version.txt index f25c43cdc2..c70613aa09 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.10.13 +0.10.14