diff --git a/.github/workflows/deploy_dockerhub.yml b/.github/workflows/deploy_dockerhub.yml new file mode 100644 index 000000000..a8900d166 --- /dev/null +++ b/.github/workflows/deploy_dockerhub.yml @@ -0,0 +1,18 @@ +# https://github.com/marketplace/actions/publish-docker +name: Deploy to dockerhub +on: [push] +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@master + - run: git archive -v -o container/charliecloud/parpe_base/parpe.tar.gz --format=tar.gz HEAD + - name: Publish to Registry + uses: elgohr/Publish-Docker-Github-Action@2.8 + with: + name: dweindl/parpe + username: ${{ secrets.DOCKER_USERNAME }} + password: ${{ secrets.DOCKER_PASSWORD }} + workdir: container/charliecloud/parpe_base/ + dockerfile: Dockerfile + tag_names: true diff --git a/CMakeLists.txt b/CMakeLists.txt index 70fd2da42..594a01e54 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ -cmake_minimum_required(VERSION 3.6) -cmake_policy(VERSION 3.6) +cmake_minimum_required(VERSION 3.7) +cmake_policy(VERSION 3.7) project(parpe) diff --git a/CMakeModules/BuildType.cmake b/CMakeModules/BuildType.cmake index e5ff5f9e4..3a4918500 100644 --- a/CMakeModules/BuildType.cmake +++ b/CMakeModules/BuildType.cmake @@ -1,6 +1,6 @@ # Set a default build type if none was specified. # If inside git repository build debug, otherwise release. -set(default_build_type "Release") +set(default_build_type "RelWithDebInfo") if(EXISTS "${CMAKE_SOURCE_DIR}/.git") set(default_build_type "Debug") endif() diff --git a/README.md b/README.md index 3a06e178a..4edc27281 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,8 @@ [![Run Status](https://api.shippable.com/projects/59463d3e8993d7070010407b/badge?branch=master)](https://app.shippable.com/github/dweindl/parPE) [![Coverage Badge](https://api.shippable.com/projects/59463d3e8993d7070010407b/coverageBadge?branch=master)](https://app.shippable.com/github/dweindl/parPE) [![Codacy Badge](https://api.codacy.com/project/badge/Grade/1f1ee5a0d90d431499f200a148fb7fdc)](https://www.codacy.com?utm_source=github.com&utm_medium=referral&utm_content=ICB-DCM/parPE&utm_campaign=Badge_Grade) +[![DOI](https://zenodo.org/badge/92953596.svg)](https://zenodo.org/badge/latestdoi/92953596) + # parPE The *parPE* library provides functionality for solving large-scale parameter @@ -110,6 +112,12 @@ parPE is being used or has been used in the following projects: Bioinformatics, btz581, [doi:10.1093/bioinformatics/btz581](https://doi.org/10.1093/bioinformatics/btz581) (preprint: [doi:10.1101/579045](https://www.biorxiv.org/content/10.1101/579045v1)). +- Paul Stapor, 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. + preprint: [doi:10.1101/859884](https://doi.org/10.1101/859884). + - [CanPathPro](http://canpathpro.eu/) diff --git a/benchmark_collection/import_and_run.sh b/benchmark_collection/import_and_run.sh index 0d0e57cf7..477873872 100755 --- a/benchmark_collection/import_and_run.sh +++ b/benchmark_collection/import_and_run.sh @@ -54,7 +54,7 @@ CMD="${PARPE_DIR}/misc/generateHDF5DataFileFromText.py \ -d ${AMICI_MODEL_DIR} \ -m ${PETAB_MODEL_DIR}/measurementData_${MODEL_NAME}.tsv \ -c ${PETAB_MODEL_DIR}/experimentalCondition_${MODEL_NAME}.tsv \ - -n model_${MODEL_NAME} \ + -n ${MODEL_NAME} \ -p ${PETAB_MODEL_DIR}/parameters_${MODEL_NAME}.tsv" echo $CMD $CMD diff --git a/container/charliecloud/parpe_base/install.sh b/container/charliecloud/parpe_base/install.sh index 53fe47f40..ae4cf2b1e 100755 --- a/container/charliecloud/parpe_base/install.sh +++ b/container/charliecloud/parpe_base/install.sh @@ -96,19 +96,20 @@ apt-get install -q -y git #echo "================= Adding awscli 1.14.64 ============" #sudo pip install -q 'awscli==1.14.64' -#echo "================= parPE requirements ============" -apt-get install gfortran libmpich-dev libatlas-base-dev libboost-all-dev libhdf5-dev cmake libceres-dev coinor-libipopt-dev swig3.0 python3-venv hdf5-tools libpython-dev +echo "================= parPE requirements ============" +# using openmpi coming with libboost-all-dev instead of libmpich-dev +apt-get install gfortran libatlas-base-dev libboost-all-dev libhdf5-dev cmake libceres-dev coinor-libipopt-dev swig3.0 python3-venv hdf5-tools libpython-dev # for setuptools to find: ln -s /usr/bin/swig3.0 /usr/bin/swig python3 -m pip install --upgrade pip pip3 install -U setuptools pkgconfig wheel # echo "================= Intalling Shippable CLIs =================" -# +# # git clone https://github.com/Shippable/node.git nodeRepo # ./nodeRepo/shipctl/x86_64/Ubuntu_16.04/install.sh # rm -rf nodeRepo -# +# # echo "Installed Shippable CLIs successfully" # echo "-------------------------------------" diff --git a/container/charliecloud/parpe_base/install_parpe.sh b/container/charliecloud/parpe_base/install_parpe.sh index 1907777b1..74aef1775 100755 --- a/container/charliecloud/parpe_base/install_parpe.sh +++ b/container/charliecloud/parpe_base/install_parpe.sh @@ -3,10 +3,10 @@ set -e cd -# Clone clean repository -git clone --single-branch --branch feature_charlie --depth=1 https://github.com/ICB-DCM/parPE.git +# unpack git archive +mkdir parPE && cd parPE +tar -xzf /u18/parpe.tar.gz -cd parPE export PARPE_BASE=$(pwd) # Build dependencies @@ -36,6 +36,7 @@ CC=mpicc CXX=mpiCC cmake \ -DCERES_INCLUDE_DIRS="/usr/include/;/usr/include/eigen3" \ -DMPI_INCLUDE_DIRS=/usr/include/openmpi-x86_64/ \ -DBUILD_TESTS=ON \ + "-DTESTS_MPIEXEC_COMMAND=mpiexec;--allow-run-as-root;-n;4;--oversubscribe;--mca;btl_vader_single_copy_mechanism;none;--mca;btl;^openib;--mca;oob_tcp_if_include;lo;--mca;btl_tcp_if_include;lo;--mca;orte_base_help_aggregate;0" \ .. make -j12 VERBOSE=1 diff --git a/deps/AMICI/.github/workflows/sbml-semantic-test-suite.yml b/deps/AMICI/.github/workflows/sbml-semantic-test-suite.yml index ad7d27e48..1cbdf6c34 100644 --- a/deps/AMICI/.github/workflows/sbml-semantic-test-suite.yml +++ b/deps/AMICI/.github/workflows/sbml-semantic-test-suite.yml @@ -4,7 +4,6 @@ on: branches: - develop - master - - test_action pull_request: branches: - master diff --git a/deps/AMICI/.github/workflows/test-large-model.yml b/deps/AMICI/.github/workflows/test-large-model.yml new file mode 100644 index 000000000..f77c10ed8 --- /dev/null +++ b/deps/AMICI/.github/workflows/test-large-model.yml @@ -0,0 +1,82 @@ +name: Performance test +on: + push: + branches: + - develop + - master + - feature_sparse_quadratures + + 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 + - name: forward_simulation_non_optimal_parameters + run: | + check_time.sh forward_simulation_non_optimal_parameters tests/performance/test.py forward_simulation_non_optimal_parameters + - name: adjoint_sensitivities_non_optimal_parameters + run: | + check_time.sh adjoint_sensitivities_non_optimal_parameters tests/performance/test.py adjoint_sensitivities_non_optimal_parameters \ No newline at end of file diff --git a/deps/AMICI/.gitrepo b/deps/AMICI/.gitrepo index 0dc0875d7..fa1627219 100644 --- a/deps/AMICI/.gitrepo +++ b/deps/AMICI/.gitrepo @@ -5,8 +5,8 @@ ; [subrepo] remote = git@github.com:ICB-DCM/AMICI.git - branch = v0.10.13 - commit = 76162dfc09d74d439939ab92eac7a2f7b80f0873 - parent = 3fe186d2f265eb58466b7b0e25c56bc9c3a71abf + branch = v0.10.16 + commit = a4e077f317d2d3eb6f93fcb1a849ddb4ffa7a175 + parent = 6f77a9f618dc859380d4dcb6ee73676d2b388b62 cmdver = 0.4.0 method = merge diff --git a/deps/AMICI/.travis.yml b/deps/AMICI/.travis.yml index a11d417e5..36af7b7c0 100644 --- a/deps/AMICI/.travis.yml +++ b/deps/AMICI/.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/deps/AMICI/INSTALL.md b/deps/AMICI/INSTALL.md index e8ddcee09..dd802ccf6 100755 --- a/deps/AMICI/INSTALL.md +++ b/deps/AMICI/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/deps/AMICI/README.md b/deps/AMICI/README.md index efac331da..3270e6608 100644 --- a/deps/AMICI/README.md +++ b/deps/AMICI/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/deps/AMICI/documentation/FAQ.md b/deps/AMICI/documentation/FAQ.md index cbdaf3035..ea79fe25c 100644 --- a/deps/AMICI/documentation/FAQ.md +++ b/deps/AMICI/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/deps/AMICI/documentation/amici_refs.bib b/deps/AMICI/documentation/amici_refs.bib index 8366aecbe..0f3c404f3 100644 --- a/deps/AMICI/documentation/amici_refs.bib +++ b/deps/AMICI/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/deps/AMICI/documentation/development.md b/deps/AMICI/documentation/development.md index f8b951bc9..b4463bdc8 100644 --- a/deps/AMICI/documentation/development.md +++ b/deps/AMICI/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/deps/AMICI/documentation/references.md b/deps/AMICI/documentation/references.md index 759461b43..d4739c704 100644 --- a/deps/AMICI/documentation/references.md +++ b/deps/AMICI/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/deps/AMICI/include/amici/abstract_model.h b/deps/AMICI/include/amici/abstract_model.h index c6cadc3f9..3f6457382 100644 --- a/deps/AMICI/include/amici/abstract_model.h +++ b/deps/AMICI/include/amici/abstract_model.h @@ -108,7 +108,8 @@ class AbstractModel { const AmiVector &dx) = 0; /** - * @brief Parameter derivative of residual function + * @brief Model specific sparse implementation of + explicit parameter derivative of right hand side * @param t time * @param x state * @param dx time derivative of state (DAE only) @@ -668,6 +669,18 @@ class AbstractModel { const realtype *p, const realtype *k, const realtype *h, const realtype *w, const realtype *tcl, const realtype *stcl); + + /** + * @brief Model specific implementation for dwdp, column pointers + * @param indexptrs column pointers + **/ + virtual void fdwdp_colptrs(sunindextype *indexptrs); + + /** + * @brief Model specific implementation for dwdp, row values + * @param indexvals row values + **/ + virtual void fdwdp_rowvals(sunindextype *indexvals); /** * @brief Model specific sensitivity implementation of dwdp diff --git a/deps/AMICI/include/amici/amici.h b/deps/AMICI/include/amici/amici.h index a0848e692..89b4bbca9 100644 --- a/deps/AMICI/include/amici/amici.h +++ b/deps/AMICI/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/deps/AMICI/include/amici/defines.h b/deps/AMICI/include/amici/defines.h index 748257864..1c90baa92 100644 --- a/deps/AMICI/include/amici/defines.h +++ b/deps/AMICI/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/deps/AMICI/include/amici/misc.h b/deps/AMICI/include/amici/misc.h index 627799669..6284a1d8f 100644 --- a/deps/AMICI/include/amici/misc.h +++ b/deps/AMICI/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/deps/AMICI/include/amici/model.h b/deps/AMICI/include/amici/model.h index 351c46f9b..21733ea9b 100644 --- a/deps/AMICI/include/amici/model.h +++ b/deps/AMICI/include/amici/model.h @@ -15,6 +15,9 @@ namespace amici { class ExpData; class Model; class Solver; +class AmiciApplication; + +extern AmiciApplication defaultContext; } // namespace amici @@ -68,14 +71,18 @@ class Model : public AbstractModel { * @param plist indexes wrt to which sensitivities are to be computed * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events + * @param pythonGenerated flag indicating matlab or python wrapping + * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit + * @param ndxdotdp_implicit number of nonzero elements dxdotdp_implicit */ Model(int nx_rdata, int nxtrue_rdata, int nx_solver, int nxtrue_solver, int ny, int nytrue, int nz, int nztrue, int ne, int nJ, int nw, - int ndwdx, int ndwdp, int ndxdotdw, std::vector ndJydy, int nnz, - int ubw, int lbw, amici::SecondOrderMode o2mode, + int ndwdx, int ndwdp, int ndxdotdw, std::vector ndJydy, + int nnz, int ubw, int lbw, amici::SecondOrderMode o2mode, const std::vector &p, std::vector k, const std::vector &plist, std::vector idlist, - std::vector z2event); + std::vector z2event, bool pythonGenerated=false, + int ndxdotdp_explicit=0, int ndxdotdp_implicit=0); /** destructor */ ~Model() override = default; @@ -127,6 +134,8 @@ class Model : public AbstractModel { using AbstractModel::fdsigmaydp; using AbstractModel::fdsigmazdp; using AbstractModel::fdwdp; + using AbstractModel::fdwdp_colptrs; + using AbstractModel::fdwdp_rowvals; using AbstractModel::fdwdx; using AbstractModel::fdwdx_colptrs; using AbstractModel::fdwdx_rowvals; @@ -1030,12 +1039,6 @@ class Model : public AbstractModel { */ bool getAlwaysCheckFinite() const; - /** - * @brief check whether the model was generated from python - * @return that - */ - virtual bool wasPythonGenerated() const { return false; } - /** * @brief Initial states * @param x pointer to state variables @@ -1127,7 +1130,6 @@ class Model : public AbstractModel { /** number of nonzero entries in dxdotdw */ int ndxdotdw{0}; - /** number of nonzero entries in dJydy */ std::vector ndJydy; @@ -1142,7 +1144,16 @@ class Model : public AbstractModel { /** lower bandwith of the jacobian */ int lbw{0}; - + + /** flag indicating Matlab or python based model generation */ + bool pythonGenerated; + + /** number of nonzero entries in ndxdotdp_explicit */ + int ndxdotdp_explicit{0}; + + /** number of nonzero entries in ndxdotdp_implicit */ + int ndxdotdp_implicit{0}; + /** flag indicating whether for sensi == AMICI_SENSI_ORDER_SECOND * directional or full second order derivative will be computed */ SecondOrderMode o2mode{SecondOrderMode::none}; @@ -1150,9 +1161,19 @@ class Model : public AbstractModel { /** flag array for DAE equations */ std::vector idlist; - /** temporary storage of dxdotdp data across functions (dimension: nplist x - * nx_solver, row-major) */ + /** temporary storage of dxdotdp data across functions, Python only + (dimension: nplist x nx_solver, nnz: ndxdotdp_explicit, type CSC_MAT) */ + mutable SUNMatrixWrapper dxdotdp_explicit; + + /** temporary storage of dxdotdp_implicit data across functions, Python only + (dimension: nplist x * nx_solver, nnz: ndxdotdp_implicit, type CSC_MAT) */ + mutable SUNMatrixWrapper dxdotdp_implicit; + + /** temporary storage of dxdotdp data across functions, Matlab only + (dimension: nplist x nx_solver, row-major) */ AmiVectorArray dxdotdp; + /** AMICI context */ + AmiciApplication *app = &defaultContext; protected: /** @@ -1252,7 +1273,7 @@ class Model : public AbstractModel { virtual void fdJydy_colptrs(sunindextype *indexptrs, int index); /** - * @brief Model specific implementation of fdxdotdw row vals + * @brief Model specific implementation of fdJydy row vals * @param indexptrs row val pointers * @param index ytrue index */ @@ -1554,6 +1575,9 @@ class Model : public AbstractModel { /** Sparse dxdotdw temporary storage (dimension: ndxdotdw) */ mutable SUNMatrixWrapper dxdotdw; + /** Sparse dwdp temporary storage (dimension: ndwdp) */ + mutable SUNMatrixWrapper dwdp; + /** Sparse dwdx temporary storage (dimension: ndwdx) */ mutable SUNMatrixWrapper dwdx; @@ -1567,12 +1591,12 @@ class Model : public AbstractModel { mutable std::vector mz; /** Sparse observable derivative of data likelihood, - * only used if wasPythonGenerated()==true + * only used if pythonGenerated==true * (dimension nytrue, nJ x ny, row-major) */ mutable std::vector dJydy; /** observable derivative of data likelihood, - * only used if wasPythonGenerated()==false + * only used if pythonGenerated==false * (dimension nJ x ny x nytrue, row-major) */ mutable std::vector dJydy_matlab; @@ -1655,11 +1679,6 @@ class Model : public AbstractModel { /** tempory storage of w data across functions (dimension: nw) */ mutable std::vector w; - /** tempory storage of sparse/dense dwdp data across functions - * (dimension: ndwdp) - */ - mutable std::vector dwdp; - /** tempory storage for flattened sx, * (dimension: nx_solver x nplist, row-major) */ diff --git a/deps/AMICI/include/amici/model_dae.h b/deps/AMICI/include/amici/model_dae.h index 321851990..c0be1e2ea 100644 --- a/deps/AMICI/include/amici/model_dae.h +++ b/deps/AMICI/include/amici/model_dae.h @@ -13,7 +13,6 @@ #include namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; class ExpData; class IDASolver; @@ -59,21 +58,24 @@ class Model_DAE : public Model { * @param plist indexes wrt to which sensitivities are to be computed * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events + * @param pythonGenerated flag indicating matlab or python wrapping + * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit + * @param ndxdotdp_implicit number of nonzero elements dxdotdp_implicit */ Model_DAE(const int nx_rdata, const int nxtrue_rdata, const int nx_solver, const int nxtrue_solver, const int ny, const int nytrue, const int nz, const int nztrue, const int ne, const int nJ, - const int nw, const int ndwdx, const int ndwdp, - const int ndxdotdw, std::vector ndJydy, const int nnz, - const int ubw, const int lbw, const SecondOrderMode o2mode, + const int nw, const int ndwdx, const int ndwdp, const int ndxdotdw, + std::vector ndJydy, const int nnz, const int ubw, + const int lbw, const SecondOrderMode o2mode, std::vector const &p, std::vector const &k, - std::vector const &plist, - std::vector const &idlist, - std::vector const &z2event) + std::vector const &plist, std::vector const &idlist, + std::vector const &z2event, const bool pythonGenerated=false, + const int ndxdotdp_explicit=0, const int ndxdotdp_implicit=0) : Model(nx_rdata, nxtrue_rdata, nx_solver, nxtrue_solver, ny, nytrue, - nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, - std::move(ndJydy), nnz, ubw, lbw, o2mode, p, k, plist, idlist, - z2event) {} + nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, std::move(ndJydy), + nnz, ubw, lbw, o2mode, p, k, plist, idlist, z2event, + pythonGenerated, ndxdotdp_explicit, ndxdotdp_implicit) {} void fJ(realtype t, realtype cj, const AmiVector &x, const AmiVector &dx, const AmiVector &xdot, SUNMatrix J) override; diff --git a/deps/AMICI/include/amici/model_ode.h b/deps/AMICI/include/amici/model_ode.h index 41aa2be37..c18db306f 100644 --- a/deps/AMICI/include/amici/model_ode.h +++ b/deps/AMICI/include/amici/model_ode.h @@ -14,7 +14,6 @@ #include namespace amici { -extern msgIdAndTxtFp warnMsgIdAndTxt; class CVodeSolver; @@ -59,21 +58,25 @@ class Model_ODE : public Model { * @param plist indexes wrt to which sensitivities are to be computed * @param idlist indexes indicating algebraic components (DAE only) * @param z2event mapping of event outputs to events + * @param pythonGenerated flag indicating matlab or python wrapping + * @param ndxdotdp_explicit number of nonzero elements dxdotdp_explicit + * @param ndxdotdp_implicit number of nonzero elements dxdotdp_implicit */ Model_ODE(const int nx_rdata, const int nxtrue_rdata, const int nx_solver, const int nxtrue_solver, const int ny, const int nytrue, const int nz, const int nztrue, const int ne, const int nJ, const int nw, const int ndwdx, const int ndwdp, - const int ndxdotdw, std::vector ndJydy, const int nnz, - const int ubw, const int lbw, const SecondOrderMode o2mode, - std::vector const &p, std::vector const &k, - std::vector const &plist, + const int ndxdotdw, std::vector ndJydy, + const int nnz, const int ubw, const int lbw, + const SecondOrderMode o2mode, std::vector const &p, + std::vector const &k, std::vector const &plist, std::vector const &idlist, - std::vector const &z2event) + std::vector const &z2event, const bool pythonGenerated=false, + const int ndxdotdp_explicit=0, const int ndxdotdp_implicit=0) : Model(nx_rdata, nxtrue_rdata, nx_solver, nxtrue_solver, ny, nytrue, - nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, - std::move(ndJydy), nnz, ubw, lbw, o2mode, p, k, plist, idlist, - z2event) {} + nz, nztrue, ne, nJ, nw, ndwdx, ndwdp, ndxdotdw, std::move(ndJydy), + nnz, ubw, lbw, o2mode, p, k, plist, idlist, z2event, + pythonGenerated, ndxdotdp_explicit, ndxdotdp_implicit) {} void fJ(realtype t, realtype cj, const AmiVector &x, const AmiVector &dx, const AmiVector &xdot, SUNMatrix J) override; @@ -219,7 +222,7 @@ class Model_ODE : public Model { */ void fdxdotdw(realtype t, const N_Vector x); - /** Sensitivity of dx/dt wrt model parameters p + /** Explicit sensitivity of dx/dt wrt model parameters p * @param t timepoint * @param x Vector with the states * @return status flag indicating successful execution @@ -403,7 +406,7 @@ class Model_ODE : public Model { const realtype *p, const realtype *k, const realtype *h, const realtype *w) = 0; - /** model specific implementation of fdxdotdp, with w chainrule + /** model specific implementation of fdxdotdp, with w chainrule (Matlab) * @param dxdotdp partial derivative xdot wrt p * @param t timepoint * @param x Vector with the states @@ -419,19 +422,39 @@ class Model_ODE : public Model { const realtype *h, int ip, const realtype *w, const realtype *dwdp); - /** model specific implementation of fdxdotdp, without w chainrule - * @param dxdotdp partial derivative xdot wrt p + /** model specific implementation of fdxdotdp_explicit, no w chainrule (Py) + * @param dxdotdp_explicit partial derivative xdot wrt p * @param t timepoint * @param x Vector with the states * @param p parameter vector * @param k constants vector * @param h heavyside vector - * @param ip parameter index * @param w vector with helper variables */ - virtual void fdxdotdp(realtype *dxdotdp, realtype t, const realtype *x, - const realtype *p, const realtype *k, - const realtype *h, int ip, const realtype *w); + virtual void fdxdotdp_explicit(realtype *dxdotdp_explicit, realtype t, + const realtype *x, const realtype *p, + const realtype *k, const realtype *h, + const realtype *w); + + /** model specific implementation of fdxdotdp_explicit, colptrs part + * @param indexptrs column pointers + */ + virtual void fdxdotdp_explicit_colptrs(sunindextype *indexptrs); + + /** model specific implementation of fdxdotdp_explicit, rowvals part + * @param indexvals row values + */ + virtual void fdxdotdp_explicit_rowvals(sunindextype *indexvals); + + /** model specific implementation of fdxdotdp_implicit, colptrs part + * @param indexptrs column pointers + */ + virtual void fdxdotdp_implicit_colptrs(sunindextype *indexptrs); + + /** model specific implementation of fdxdotdp_implicit, rowvals part + * @param indexvals row values + */ + virtual void fdxdotdp_implicit_rowvals(sunindextype *indexvals); /** model specific implementation of fdxdotdw, data part * @param dxdotdw partial derivative xdot wrt w @@ -451,12 +474,11 @@ class Model_ODE : public Model { */ virtual void fdxdotdw_colptrs(sunindextype *indexptrs); - /** model specific implementation of fdxdotdw, colptrs part + /** model specific implementation of fdxdotdw, rowvals part * @param indexvals row values */ virtual void fdxdotdw_rowvals(sunindextype *indexvals); }; - } // namespace amici #endif // MODEL_H diff --git a/deps/AMICI/include/amici/newton_solver.h b/deps/AMICI/include/amici/newton_solver.h index c0c1232b3..5482e0a20 100644 --- a/deps/AMICI/include/amici/newton_solver.h +++ b/deps/AMICI/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/deps/AMICI/include/amici/serialization.h b/deps/AMICI/include/amici/serialization.h index 15dee2c22..b83e2118a 100644 --- a/deps/AMICI/include/amici/serialization.h +++ b/deps/AMICI/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; @@ -109,6 +111,9 @@ void serialize(Archive &ar, amici::Model &u, const unsigned int version) { ar &u.pscale; ar &u.tstart; ar &u.stateIsNonNegative; + ar &u.pythonGenerated; + ar &u.ndxdotdp_explicit; + ar &u.ndxdotdp_implicit; } diff --git a/deps/AMICI/include/amici/solver.h b/deps/AMICI/include/amici/solver.h index 1e9c6dfb4..8ab0b75c1 100644 --- a/deps/AMICI/include/amici/solver.h +++ b/deps/AMICI/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/deps/AMICI/include/amici/solver_cvodes.h b/deps/AMICI/include/amici/solver_cvodes.h index ce2278c3d..72e012aca 100644 --- a/deps/AMICI/include/amici/solver_cvodes.h +++ b/deps/AMICI/include/amici/solver_cvodes.h @@ -26,8 +26,6 @@ namespace amici { class CVodeSolver : public Solver { public: - CVodeSolver() = default; - ~CVodeSolver() override = default; /** diff --git a/deps/AMICI/include/amici/solver_idas.h b/deps/AMICI/include/amici/solver_idas.h index 21648ddd8..5a9bf100d 100644 --- a/deps/AMICI/include/amici/solver_idas.h +++ b/deps/AMICI/include/amici/solver_idas.h @@ -24,7 +24,6 @@ namespace amici { class IDASolver : public Solver { public: - IDASolver() = default; ~IDASolver() override = default; /** diff --git a/deps/AMICI/include/amici/sundials_matrix_wrapper.h b/deps/AMICI/include/amici/sundials_matrix_wrapper.h index 59751425c..35fd16f63 100644 --- a/deps/AMICI/include/amici/sundials_matrix_wrapper.h +++ b/deps/AMICI/include/amici/sundials_matrix_wrapper.h @@ -116,6 +116,12 @@ class SUNMatrixWrapper { */ sunindextype columns() const; + /** + * @brief Get the number of non-zero elements (sparse matrices only) + * @return number + */ + sunindextype nonzeros() const; + /** * @brief Get the index values of a sparse matrix * @return index array @@ -153,15 +159,53 @@ class SUNMatrixWrapper { */ void multiply(gsl::span c, gsl::span b) const; + /** + * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b + * @param c output vector, may already contain values + * @param b multiplication vector + * @param cols int vector for column reordering + * @param transpose bool transpose A before multiplication + */ + void multiply(N_Vector c, + const N_Vector b, + gsl::span cols, + bool transpose) const; + + /** + * @brief Perform reordered matrix vector multiplication c += A[:,cols]*b + * @param c output vector, may already contain values + * @param b multiplication vector + * @param cols int vector for column reordering + * @param transpose bool transpose A before multiplication + */ + void multiply(gsl::span c, + gsl::span b, + gsl::span cols, + bool transpose) const; + + /** + * @brief Perform matrix matrix multiplication + C[:, :] += A * B + for sparse A, B, C + * @param C output matrix, may already contain values + * @param B multiplication matrix + */ + void sparse_multiply(SUNMatrixWrapper *C, + SUNMatrixWrapper *B) const; + /** * @brief Set to 0.0 */ void zero(); + /** + * @brief CSC matrix to which all methods are applied + */ + SUNMatrix matrix = nullptr; + private: void update_ptrs(); - SUNMatrix matrix = nullptr; realtype *data_ptr = nullptr; sunindextype *indexptrs_ptr = nullptr; sunindextype *indexvals_ptr = nullptr; diff --git a/deps/AMICI/matlab/examples/amiExamples.m b/deps/AMICI/matlab/examples/amiExamples.m index 9dd4b5bb3..82ee79c76 100644 --- a/deps/AMICI/matlab/examples/amiExamples.m +++ b/deps/AMICI/matlab/examples/amiExamples.m @@ -13,5 +13,5 @@ example_jakstat_adjoint_hvp example_robertson example_neuron -example_dae_events -path(oldpath); \ No newline at end of file +example_calvetti +path(oldpath); diff --git a/deps/AMICI/models/model_dirac/CMakeLists.txt b/deps/AMICI/models/model_dirac/CMakeLists.txt index 2fbbab79a..f73a17155 100644 --- a/deps/AMICI/models/model_dirac/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_dirac/main.cpp b/deps/AMICI/models/model_dirac/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_dirac/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_dirac/model_dirac.h b/deps/AMICI/models/model_dirac/model_dirac.h index 839b28bb7..389dce41c 100644 --- a/deps/AMICI/models/model_dirac/model_dirac.h +++ b/deps/AMICI/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/deps/AMICI/models/model_events/CMakeLists.txt b/deps/AMICI/models/model_events/CMakeLists.txt index c60998f9e..150dd43b0 100644 --- a/deps/AMICI/models/model_events/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_events/main.cpp b/deps/AMICI/models/model_events/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_events/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_events/model_events.h b/deps/AMICI/models/model_events/model_events.h index 1ea53269c..f3d592eb5 100644 --- a/deps/AMICI/models/model_events/model_events.h +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint/CMakeLists.txt b/deps/AMICI/models/model_jakstat_adjoint/CMakeLists.txt index 27f5fd7c6..76484be70 100644 --- a/deps/AMICI/models/model_jakstat_adjoint/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint/main.cpp b/deps/AMICI/models/model_jakstat_adjoint/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_jakstat_adjoint/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint/model_jakstat_adjoint.h b/deps/AMICI/models/model_jakstat_adjoint/model_jakstat_adjoint.h index e512ee939..3cf5fdec2 100644 --- a/deps/AMICI/models/model_jakstat_adjoint/model_jakstat_adjoint.h +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint_o2/CMakeLists.txt b/deps/AMICI/models/model_jakstat_adjoint_o2/CMakeLists.txt index 3fad28004..92516553d 100644 --- a/deps/AMICI/models/model_jakstat_adjoint_o2/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint_o2/main.cpp b/deps/AMICI/models/model_jakstat_adjoint_o2/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_jakstat_adjoint_o2/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h b/deps/AMICI/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h index 52f25fccd..ef9f6a163 100644 --- a/deps/AMICI/models/model_jakstat_adjoint_o2/model_jakstat_adjoint_o2.h +++ b/deps/AMICI/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/deps/AMICI/models/model_nested_events/CMakeLists.txt b/deps/AMICI/models/model_nested_events/CMakeLists.txt index b5fa8de31..120e5a08e 100644 --- a/deps/AMICI/models/model_nested_events/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_nested_events/main.cpp b/deps/AMICI/models/model_nested_events/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_nested_events/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_nested_events/model_nested_events.h b/deps/AMICI/models/model_nested_events/model_nested_events.h index 76a3c9b80..5ca7091d6 100644 --- a/deps/AMICI/models/model_nested_events/model_nested_events.h +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron/CMakeLists.txt b/deps/AMICI/models/model_neuron/CMakeLists.txt index b7af8a3d7..0c63349ca 100644 --- a/deps/AMICI/models/model_neuron/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron/main.cpp b/deps/AMICI/models/model_neuron/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_neuron/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron/model_neuron.h b/deps/AMICI/models/model_neuron/model_neuron.h index ddf2f6e01..fa33f518a 100644 --- a/deps/AMICI/models/model_neuron/model_neuron.h +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron_o2/CMakeLists.txt b/deps/AMICI/models/model_neuron_o2/CMakeLists.txt index 7ad5944d5..7c2121cb0 100644 --- a/deps/AMICI/models/model_neuron_o2/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron_o2/main.cpp b/deps/AMICI/models/model_neuron_o2/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_neuron_o2/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_neuron_o2/model_neuron_o2.h b/deps/AMICI/models/model_neuron_o2/model_neuron_o2.h index ef16e543d..a4fb5d03d 100644 --- a/deps/AMICI/models/model_neuron_o2/model_neuron_o2.h +++ b/deps/AMICI/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/deps/AMICI/models/model_robertson/CMakeLists.txt b/deps/AMICI/models/model_robertson/CMakeLists.txt index 072011db8..7af14aba5 100644 --- a/deps/AMICI/models/model_robertson/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_robertson/main.cpp b/deps/AMICI/models/model_robertson/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_robertson/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_robertson/model_robertson.h b/deps/AMICI/models/model_robertson/model_robertson.h index 556cea6e7..c526ad14a 100644 --- a/deps/AMICI/models/model_robertson/model_robertson.h +++ b/deps/AMICI/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/deps/AMICI/models/model_steadystate/CMakeLists.txt b/deps/AMICI/models/model_steadystate/CMakeLists.txt index 0537b4b5b..0b1458426 100644 --- a/deps/AMICI/models/model_steadystate/CMakeLists.txt +++ b/deps/AMICI/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/deps/AMICI/models/model_steadystate/main.cpp b/deps/AMICI/models/model_steadystate/main.cpp index 6d6caf54a..07ff62a32 100644 --- a/deps/AMICI/models/model_steadystate/main.cpp +++ b/deps/AMICI/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/deps/AMICI/models/model_steadystate/model_steadystate.h b/deps/AMICI/models/model_steadystate/model_steadystate.h index 06a1c2773..4489e3a66 100644 --- a/deps/AMICI/models/model_steadystate/model_steadystate.h +++ b/deps/AMICI/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/deps/AMICI/python/amici/__init__.py b/deps/AMICI/python/amici/__init__.py index c26c9ff28..4ae68672e 100644 --- a/deps/AMICI/python/amici/__init__.py +++ b/deps/AMICI/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/deps/AMICI/python/amici/__init__.template.py b/deps/AMICI/python/amici/__init__.template.py index b68a5bae9..5c42765ec 100644 --- a/deps/AMICI/python/amici/__init__.template.py +++ b/deps/AMICI/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/deps/AMICI/python/amici/gradient_check.py b/deps/AMICI/python/amici/gradient_check.py index 1a6f73483..b56fca96b 100644 --- a/deps/AMICI/python/amici/gradient_check.py +++ b/deps/AMICI/python/amici/gradient_check.py @@ -64,8 +64,6 @@ def check_derivatives(model, solver, edata, assert_fun, rtol: relative tolerance epsilon: finite difference step-size """ - from scipy.optimize import check_grad - p = np.array(model.getParameters()) rdata = runAmiciSimulation(model, solver, edata) diff --git a/deps/AMICI/python/amici/ode_export.py b/deps/AMICI/python/amici/ode_export.py index cfcb1fd97..de7833321 100644 --- a/deps/AMICI/python/amici/ode_export.py +++ b/deps/AMICI/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__ @@ -55,44 +58,35 @@ '(realtype *J, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w, const realtype *dwdx)', - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity'] }, 'JB': { 'signature': '(realtype *JB, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *xB, const realtype *w, const realtype *dwdx)', - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity'] }, 'JDiag': { 'signature': '(realtype *JDiag, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w, const realtype *dwdx)', - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity'] }, 'JSparse': { 'signature': '(realtype *JSparse, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w, const realtype *dwdx)', - 'sparse': - True, - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity', 'sparse'] }, 'JSparseB': { 'signature': '(realtype *JSparseB, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *xB, const realtype *w, const realtype *dwdx)', - 'sparse': - True, - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity', 'sparse'] }, 'Jy': { 'signature': @@ -111,45 +105,43 @@ '(realtype *dJydy, const int iy, const realtype *p, ' 'const realtype *k, const realtype *y, ' 'const realtype *sigmay, const realtype *my)', - 'sparse': - True, + 'flags': ['sparse'] }, 'dwdp': { 'signature': '(realtype *dwdp, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' - 'const realtype *w, const realtype *tcl, const realtype *dtcldp, ' - 'const int ip)', - 'assume_pow_positivity': - True, + 'const realtype *w, const realtype *tcl, const realtype *dtcldp)', + 'flags': ['assume_pow_positivity', 'sparse'] }, 'dwdx': { 'signature': '(realtype *dwdx, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w, const realtype *tcl)', - 'sparse': - True, - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity', 'sparse'] }, 'dxdotdw': { 'signature': '(realtype *dxdotdw, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w)', - 'sparse': - True, - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity', 'sparse'] }, - 'dxdotdp': { + 'dxdotdp_explicit': { 'signature': - '(realtype *dxdotdp, const realtype t, const realtype *x, ' - 'const realtype *p, const realtype *k, const realtype *h, ' - 'const int ip, const realtype *w)', - 'assume_pow_positivity': - True, + '(realtype *dxdotdp_explicit, const realtype t, ' + 'const realtype *x, const realtype *p, ' + 'const realtype *k, const realtype *h, ' + 'const realtype *w)', + 'flags': ['assume_pow_positivity', 'sparse'] + }, + 'dxdotdp_implicit': { + 'signature': + '(realtype *dxdotdp_implicit, const realtype t, ' + 'const realtype *x, const realtype *p, const realtype *k, ' + 'const realtype *h, const int ip, const realtype *w)', + 'flags': ['assume_pow_positivity', 'sparse', 'dont_generate_body'] }, 'dydx': { 'signature': @@ -178,8 +170,7 @@ '(realtype *w, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, ' 'const realtype *h, const realtype *tcl)', - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity'] }, 'x0': { 'signature': @@ -207,8 +198,7 @@ '(realtype *xdot, const realtype t, const realtype *x, ' 'const realtype *p, const realtype *k, const realtype *h, ' 'const realtype *w)', - 'assume_pow_positivity': - True, + 'flags': ['assume_pow_positivity'] }, 'y': { 'signature': @@ -233,8 +223,12 @@ ## list of sparse functions sparse_functions = [ function for function in functions - if 'sparse' in functions[function] - and functions[function]['sparse'] + if 'sparse' in functions[function].get('flags', []) +] +## list of nobody functions +nobody_functions = [ + function for function in functions + if 'dont_generate_body' in functions[function].get('flags', []) ] ## list of sensitivity functions sensi_functions = [ @@ -712,12 +706,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: @@ -801,9 +800,14 @@ def __init__(self): 'x': 'JB', 'y': 'xB', }, + 'dxdotdp_implicit': { + 'x': 'dxdotdw', + 'y': 'dwdp', + }, } 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 +829,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 +1165,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 +1183,12 @@ 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())}__' + f'{strip_pysb(par.get_id())}', + real=True) + for par in self._parameters + ] 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): @@ -1273,8 +1283,8 @@ def _generateSparseSymbol(self, name): base_index = 0 for iy in range(self.ny()): symbolColPtrs, symbolRowVals, sparseList, symbolList, \ - sparseMatrix = csc_matrix(matrix[iy, :], name, - base_index=base_index) + sparseMatrix = csc_matrix(matrix[iy, :], name, + base_index=base_index) base_index += len(symbolList) self._colptrs[name].append(symbolColPtrs) self._rowvals[name].append(symbolRowVals) @@ -1283,7 +1293,9 @@ def _generateSparseSymbol(self, name): self._syms[name].append(sparseMatrix) else: symbolColPtrs, symbolRowVals, sparseList, symbolList, \ - sparseMatrix = csc_matrix(matrix, name) + sparseMatrix = csc_matrix( + matrix, name, pattern_only=name in nobody_functions + ) self._colptrs[name] = symbolColPtrs self._rowvals[name] = symbolRowVals @@ -1382,6 +1394,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]) @@ -1417,6 +1433,10 @@ def _compute_equation(self, name): # force symbols self._eqs[name] = self.sym(name) + elif name == 'dxdotdp_explicit': + # force symbols + self._derivative('xdot', 'p', name=name) + elif match_deriv: self._derivative(match_deriv.group(1), match_deriv.group(2)) @@ -1429,6 +1449,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 +1579,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) @@ -1631,11 +1654,11 @@ def _multiplication(self, name, x, y, else: xx = variables[x] - if xx.is_zero is not True and variables[y].is_zero is not True \ - and len(xx) and len(variables[y]): - self._eqs[name] = sign * xx * variables[y] - else: + if not len(xx) or not len(variables[y]) or xx.is_zero is True or \ + variables[y].is_zero is True: self._eqs[name] = sp.zeros(len(xx), len(variables[y])) + else: + self._eqs[name] = sign * xx * variables[y] def _equationFromComponent(self, name, component): """Generates the formulas of a symbolic variable from the attributes @@ -1799,7 +1822,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 """ @@ -1908,7 +1931,9 @@ def _generateCCode(self): """ for function in self.functions.keys(): - self._writeFunctionFile(function) + if 'dont_generate_body' not in \ + self.functions[function].get('flags', []): + self._writeFunctionFile(function) if function in sparse_functions: self._write_function_index(function, 'colptrs') self._write_function_index(function, 'rowvals') @@ -2065,8 +2090,7 @@ def _writeFunctionFile(self, function): # first generate the equations to make sure we have everything we # need in subsequent steps - if 'sparse' in self.functions[function] and \ - self.functions[function]['sparse']: + if function in sparse_functions: symbol = self.model.sparseeq(function) elif not self.allow_reinit_fixpar_initcond \ and function == 'sx0_fixedParameters': @@ -2103,9 +2127,8 @@ def _writeFunctionFile(self, function): # function body body = self._getFunctionBody(function, symbol) - if self.assume_pow_positivity \ - and 'assume_pow_positivity' in self.functions[function].keys()\ - and self.functions[function]['assume_pow_positivity']: + if self.assume_pow_positivity and 'assume_pow_positivity' \ + in self.functions[function].get('flags', []): body = [re.sub(r'(^|\W)pow\(', r'\1amici::pos_pow(', line) for line in body] # execute this twice to catch cases where the ending ( would be the @@ -2293,9 +2316,13 @@ def _writeModelHeader(self): 'NEVENT': '0', 'NOBJECTIVE': '1', 'NW': str(len(self.model.sym('w'))), - 'NDWDP': str(len(self.model.eq('dwdp'))), + 'NDWDP': str(len(self.model.sparsesym('dwdp'))), 'NDWDX': str(len(self.model.sparsesym('dwdx'))), 'NDXDOTDW': str(len(self.model.sparsesym('dxdotdw'))), + 'NDXDOTDP_EXPLICIT': str(len(self.model.sparsesym( + 'dxdotdp_explicit'))), + 'NDXDOTDP_IMPLICIT': str(len(self.model.sparsesym( + 'dxdotdp_implicit'))), 'NDJYDY': 'std::vector{%s}' % ','.join(str(len(x)) for x in self.model.sparsesym('dJydy')), @@ -2332,7 +2359,8 @@ def _writeModelHeader(self): for fun in [ 'w', 'dwdp', 'dwdx', 'x_rdata', 'x_solver', 'total_cl', 'dxdotdw', - 'dxdotdp', 'JSparse', 'JSparseB', 'dJydy' + 'dxdotdp_explicit', 'dxdotdp_implicit', 'JSparse', 'JSparseB', + 'dJydy' ]: tplData[f'{fun.upper()}_DEF'] = \ get_function_extern_declaration(fun, self.modelName) @@ -2408,9 +2436,10 @@ def _writeCMakeFile(self): Raises: """ + sources = [self.modelName + '_' + function + '.cpp ' for function in self.functions.keys() - if self.functions[function]['body'] is not None] + if self.functions[function].get('body', None) is not None] # add extra source files for sparse matrices for function in sparse_functions: @@ -2628,7 +2657,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 @@ -2807,7 +2836,7 @@ def getSwitchStatement(condition, cases, return lines -def csc_matrix(matrix, name, base_index=0): +def csc_matrix(matrix, name, base_index=0, pattern_only=False): """Generates the sparse symbolic identifiers, symbolic identifiers, sparse matrix, column pointers and row values for a symbolic variable @@ -2816,6 +2845,7 @@ def csc_matrix(matrix, name, base_index=0): matrix: dense matrix to be sparsified @type sp.Matrix name: name of the symbolic variable @type str base_index: index for first symbol name, defaults to 0 + pattern_only: flag for computing sparsity pattern without whole matrix Returns: symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix @@ -2824,24 +2854,37 @@ def csc_matrix(matrix, name, base_index=0): """ idx = 0 symbol_name_idx = base_index - sparseMatrix = sp.zeros(matrix.rows, matrix.cols) + + nrows, ncols = matrix.shape + + if not pattern_only: + sparseMatrix = sp.zeros(nrows, ncols) symbolList = [] sparseList = [] symbolColPtrs = [] symbolRowVals = [] - for col in range(0, matrix.cols): + + for col in range(0, ncols): symbolColPtrs.append(idx) - for row in range(0, matrix.rows): + for row in range(0, nrows): if not (matrix[row, col] == 0): - symbolName = f'{name}{symbol_name_idx}' - sparseMatrix[row, col] = sp.Symbol(symbolName) - symbolList.append(symbolName) - sparseList.append(matrix[row, col]) symbolRowVals.append(row) idx += 1 + symbolName = f'{name}{symbol_name_idx}' + symbolList.append(symbolName) symbol_name_idx += 1 + if not pattern_only: + sparseMatrix[row, col] = sp.Symbol(symbolName, real=True) + sparseList.append(matrix[row, col]) - symbolColPtrs.append(idx) - sparseList = sp.Matrix(sparseList) + if idx == 0: + symbolColPtrs = [] # avoid bad memory access for empty matrices + else: + symbolColPtrs.append(idx) + + if not pattern_only: + sparseList = sp.Matrix(sparseList) + else: + sparseMatrix = None return symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix diff --git a/deps/AMICI/python/amici/petab_import.py b/deps/AMICI/python/amici/petab_import.py index 8a0bb8f5b..a5c4c51a2 100644 --- a/deps/AMICI/python/amici/petab_import.py +++ b/deps/AMICI/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/deps/AMICI/python/amici/pysb_import.py b/deps/AMICI/python/amici/pysb_import.py index f7a0e6bd5..9cb578184 100644 --- a/deps/AMICI/python/amici/pysb_import.py +++ b/deps/AMICI/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( @@ -235,6 +235,12 @@ def _process_pysb_expressions(model, ODE, observables, sigmas): """ for exp in model.expressions: + ODE.add_component( + Expression( + exp, + f'{exp.name}', + exp.expand_expr(expand_observables=True)) + ) if exp.name in observables: # here we do not define a new Expression from the # pysb.Expression but define an observable, so we do not need @@ -244,7 +250,8 @@ def _process_pysb_expressions(model, ODE, observables, sigmas): Observable( y, f'{exp.name}', - exp.expand_expr(expand_observables=False)) + exp + ) ) sigma_name, sigma_value = _get_sigma_name_and_value( @@ -273,16 +280,6 @@ def _process_pysb_expressions(model, ODE, observables, sigmas): elif exp.name in sigmas.values(): # do nothing pass - else: - # here we do define a new Expression from the pysb.Expression - # so we do need to expand observables as these would otherwise - # lead to dependencies between different Expressions - ODE.add_component( - Expression( - exp, - f'{exp.name}', - exp.expand_expr(expand_observables=True)) - ) def _get_sigma_name_and_value(model, name, sigmas): diff --git a/deps/AMICI/python/amici/sbml_import.py b/deps/AMICI/python/amici/sbml_import.py index 2e2a74f92..a4c2d460a 100644 --- a/deps/AMICI/python/amici/sbml_import.py +++ b/deps/AMICI/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) @@ -916,7 +934,7 @@ def processObservables(self, observables: Dict[str, Dict[str, str]], # set user-provided sigmas for iy, obsName in enumerate(observables): if obsName in sigmas: - sigmaYValues[iy] = sigmas[obsName] + sigmaYValues[iy] = sp.Symbol(sigmas[obsName], real=True) measurementYSyms = sp.Matrix( [sp.symbols(f'm{symbol}', real=True) for symbol in 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/deps/AMICI/python/amici/setuptools.py b/deps/AMICI/python/amici/setuptools.py index cb13caff3..0d66522c0 100644 --- a/deps/AMICI/python/amici/setuptools.py +++ b/deps/AMICI/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/deps/AMICI/python/bin/amici_import_petab.py b/deps/AMICI/python/bin/amici_import_petab.py deleted file mode 100755 index 2aa7174cb..000000000 --- a/deps/AMICI/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/deps/AMICI/python/sdist/setup.py b/deps/AMICI/python/sdist/setup.py index 8373809c3..20462a32d 100755 --- a/deps/AMICI/python/sdist/setup.py +++ b/deps/AMICI/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/deps/AMICI/src/abstract_model.cpp b/deps/AMICI/src/abstract_model.cpp index ba89b0f53..6b1076e20 100644 --- a/deps/AMICI/src/abstract_model.cpp +++ b/deps/AMICI/src/abstract_model.cpp @@ -291,6 +291,10 @@ void AbstractModel::fdwdp(realtype *dwdp, const realtype t, const realtype *x, const realtype *tcl, const realtype *stcl) { } +void AbstractModel::fdwdp_colptrs(sunindextype *indexptrs) {} + +void AbstractModel::fdwdp_rowvals(sunindextype *indexvals) {} + void AbstractModel::fdwdp(realtype *dwdp, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, const realtype *w, diff --git a/deps/AMICI/src/amici.cpp b/deps/AMICI/src/amici.cpp index d71a03c9e..15e8a9885 100644 --- a/deps/AMICI/src/amici.cpp +++ b/deps/AMICI/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/deps/AMICI/src/forwardproblem.cpp b/deps/AMICI/src/forwardproblem.cpp index 52112e567..2b1681c38 100644 --- a/deps/AMICI/src/forwardproblem.cpp +++ b/deps/AMICI/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)!"); @@ -381,7 +379,7 @@ void ForwardProblem::storeJacobianAndDerivativeInReturnData() { } void ForwardProblem::getEventOutput() { - if (t == model->getTimepoint(edata->nt() - 1)) { + if (t == model->getTimepoint(model->nt() - 1)) { // call from fillEvent at last timepoint model->froot(t, x, dx, rootvals); } diff --git a/deps/AMICI/src/hdf5.cpp b/deps/AMICI/src/hdf5.cpp index 510703a0f..003b1460d 100644 --- a/deps/AMICI/src/hdf5.cpp +++ b/deps/AMICI/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/deps/AMICI/src/interface_matlab.cpp b/deps/AMICI/src/interface_matlab.cpp index 5f785d37e..800bbc061 100644 --- a/deps/AMICI/src/interface_matlab.cpp +++ b/deps/AMICI/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/deps/AMICI/src/misc.cpp b/deps/AMICI/src/misc.cpp index 3876c7cbf..731559c96 100644 --- a/deps/AMICI/src/misc.cpp +++ b/deps/AMICI/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/deps/AMICI/src/model.cpp b/deps/AMICI/src/model.cpp index 9a35c36f8..977c56937 100644 --- a/deps/AMICI/src/model.cpp +++ b/deps/AMICI/src/model.cpp @@ -138,16 +138,19 @@ Model::Model(const int nx_rdata, const int nxtrue_rdata, const int nx_solver, const int lbw, SecondOrderMode o2mode, const std::vector &p, std::vector k, const std::vector &plist, std::vector idlist, - std::vector z2event) + std::vector z2event, const bool pythonGenerated, + const int ndxdotdp_explicit, const int ndxdotdp_implicit) : nx_rdata(nx_rdata), nxtrue_rdata(nxtrue_rdata), nx_solver(nx_solver), nxtrue_solver(nxtrue_solver), ny(ny), nytrue(nytrue), nz(nz), nztrue(nztrue), ne(ne), nw(nw), ndwdx(ndwdx), ndwdp(ndwdp), - ndxdotdw(ndxdotdw), ndJydy(std::move(ndJydy)), nnz(nnz), nJ(nJ), ubw(ubw), - lbw(lbw), o2mode(o2mode), idlist(std::move(idlist)), + ndxdotdw(ndxdotdw), ndJydy(std::move(ndJydy)), + nnz(nnz), nJ(nJ), ubw(ubw), lbw(lbw), pythonGenerated(pythonGenerated), + ndxdotdp_explicit(ndxdotdp_explicit), ndxdotdp_implicit(ndxdotdp_implicit), + o2mode(o2mode), idlist(std::move(idlist)), J(nx_solver, nx_solver, nnz, CSC_MAT), - dxdotdw(nx_solver, nw, ndxdotdw, CSC_MAT), + dxdotdw(nx_solver, nw, ndxdotdw, CSC_MAT), dwdp(nw, p.size(), ndwdp, CSC_MAT), dwdx(nw, nx_solver, ndwdx, CSC_MAT), M(nx_solver, nx_solver), w(nw), - dwdp(ndwdp), x_rdata(nx_rdata, 0.0), sx_rdata(nx_rdata, 0.0), h(ne, 0.0), + x_rdata(nx_rdata, 0.0), sx_rdata(nx_rdata, 0.0), h(ne, 0.0), total_cl(nx_rdata - nx_solver), stotal_cl((nx_rdata - nx_solver) * p.size()), x_pos_tmp(nx_solver), unscaledParameters(p), originalParameters(p), @@ -155,21 +158,22 @@ Model::Model(const int nx_rdata, const int nxtrue_rdata, const int nx_solver, stateIsNonNegative(nx_solver, false), pscale(std::vector(p.size(), ParameterScaling::none)) { - // Can't use derivedClass::wasPythonGenerated() in ctor. - // Guess we are using Python if ndJydy is not empty - if (!this->ndJydy.empty()) { - if (static_cast(nytrue) != this->ndJydy.size()) - throw std::runtime_error( - "Number of elements in ndJydy is not equal " - " nytrue."); + /* If Matlab wrapped: dxdotdp is a full AmiVector, + if Python wrapped: dxdotdp_explicit and dxdotdp_implicit are CSC matrices */ + if (pythonGenerated) { + dxdotdp_explicit = SUNMatrixWrapper(nx_solver, p.size(), ndxdotdp_explicit, CSC_MAT); + dxdotdp_implicit = SUNMatrixWrapper(nx_solver, p.size(), ndxdotdp_implicit, CSC_MAT); + // also dJydy depends on the way of wrapping + if (static_cast(nytrue) != this->ndJydy.size()) + throw std::runtime_error("Number of elements in ndJydy is not equal " + " nytrue."); + for (int iytrue = 0; iytrue < nytrue; ++iytrue) - dJydy.emplace_back( - SUNMatrixWrapper(nJ, ny, this->ndJydy[iytrue], CSC_MAT)); + dJydy.emplace_back(SUNMatrixWrapper(nJ, ny, this->ndJydy[iytrue], CSC_MAT)); } else { dJydy_matlab = std::vector(nJ * nytrue * ny, 0.0); } - requireSensitivitiesForAllParameters(); } @@ -177,6 +181,13 @@ bool operator==(const Model &a, const Model &b) { if (typeid(a) != typeid(b)) return false; + bool bool_dxdotdp = true; + if (a.pythonGenerated && b.pythonGenerated) + bool_dxdotdp = (a.ndxdotdp_explicit == b.ndxdotdp_explicit) && + (a.ndxdotdp_implicit == b.ndxdotdp_implicit); + if (a.pythonGenerated != b.pythonGenerated) + bool_dxdotdp = false; + return (a.nx_rdata == b.nx_rdata) && (a.nxtrue_rdata == b.nxtrue_rdata) && (a.nx_solver == b.nx_solver) && (a.nxtrue_solver == b.nxtrue_solver) && (a.ny == b.ny) && @@ -193,7 +204,7 @@ bool operator==(const Model &a, const Model &b) { (a.ts == b.ts) && (a.nmaxevent == b.nmaxevent) && (a.pscale == b.pscale) && (a.stateIsNonNegative == b.stateIsNonNegative) && - (a.tstart == b.tstart); + (a.tstart == b.tstart) && bool_dxdotdp; } void Model::initialize(AmiVector &x, AmiVector &dx, AmiVectorArray &sx, @@ -1002,7 +1013,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 +1038,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); @@ -1039,7 +1050,7 @@ void Model::addAdjointStateEventUpdate(AmiVector &xB, const int ie, const AmiVector &xdot, const AmiVector &xdot_old) { - deltasx.assign(nx_solver, 0.0); + deltaxB.assign(nx_solver, 0.0); // compute update fdeltaxB(deltaxB.data(), t, x.data(), unscaledParameters.data(), @@ -1047,7 +1058,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 +1083,7 @@ void Model::addAdjointQuadratureEventUpdate( } if (alwaysCheckFinite) { - amici::checkFinite(deltaqB, "deltaqB"); + app->checkFinite(deltaqB, "deltaqB"); } } @@ -1089,12 +1100,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; @@ -1225,8 +1236,9 @@ void Model::checkLLHBufferSize(std::vector &sllh, } void Model::initializeVectors() { - dxdotdp = AmiVectorArray(nx_solver, nplist()); sx0data.clear(); + if (!pythonGenerated) + dxdotdp = AmiVectorArray(nx_solver, nplist()); } void Model::fy(const realtype t, const AmiVector &x) { @@ -1240,7 +1252,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"); } } @@ -1249,23 +1261,17 @@ void Model::fdydp(const realtype t, const AmiVector &x) { return; dydp.assign(ny * nplist(), 0.0); - fw(t, x.data()); fdwdp(t, x.data()); - - // if dwdp is not dense, fdydp will expect the full sparse array - realtype *dwdp_tmp = dwdp.data(); - for (int ip = 0; ip < nplist(); ip++) { - // get dydp slice (ny) for current time and parameter - if (wasPythonGenerated() && nw) - dwdp_tmp = &dwdp.at(nw * ip); - + + /* get dydp slice (ny) for current time and parameter */ + for (int ip = 0; ip < nplist(); ip++) fdydp(&dydp.at(ip * ny), t, x.data(), unscaledParameters.data(), - fixedParameters.data(), h.data(), plist(ip), w.data(), dwdp_tmp); - } + fixedParameters.data(), h.data(), plist(ip), w.data(), + dwdp.data()); if (alwaysCheckFinite) { - amici::checkFinite(dydp, "dydp"); + app->checkFinite(dydp, "dydp"); } } @@ -1281,7 +1287,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 +1345,7 @@ void Model::fdsigmaydp(const int it, const ExpData *edata) { } if (alwaysCheckFinite) { - amici::checkFinite(dsigmaydp, "dsigmaydp"); + app->checkFinite(dsigmaydp, "dsigmaydp"); } } @@ -1360,7 +1366,7 @@ void Model::fdJydy(const int it, const AmiVector &x, const ExpData &edata) { fy(edata.getTimepoint(it), x); fsigmay(it, &edata); - if (wasPythonGenerated()) { + if (pythonGenerated) { for (int iyt = 0; iyt < nytrue; iyt++) { dJydy[iyt].zero(); fdJydy_colptrs(dJydy[iyt].indexptrs(), iyt); @@ -1375,7 +1381,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 +1395,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 +1416,7 @@ void Model::fdJydsigma(const int it, const AmiVector &x, const ExpData &edata) { } if (alwaysCheckFinite) { - amici::checkFinite(dJydsigma, "dJydsigma"); + app->checkFinite(dJydsigma, "dJydsigma"); } } @@ -1432,7 +1438,7 @@ void Model::fdJydp(const int it, const AmiVector &x, const ExpData &edata) { if (!edata.isSetObservedData(it, iyt)) continue; - if (wasPythonGenerated()) { + if (pythonGenerated) { // dJydp = 1.0 * dJydp + 1.0 * dJydy * dydp for (int iplist = 0; iplist < nplist(); ++iplist) { dJydy[iyt].multiply( @@ -1471,7 +1477,7 @@ void Model::fdJydx(const int it, const AmiVector &x, const ExpData &edata) { // M K K N M N // lda ldb ldc - if (wasPythonGenerated()) { + if (pythonGenerated) { for (int ix = 0; ix < nx_solver; ++ix) { dJydy[iyt].multiply( gsl::span(&dJydx.at(ix * nJ), nJ), @@ -1486,7 +1492,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 +1514,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 +1526,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 +1548,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 +1560,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 +1619,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 +1640,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 +1662,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 +1765,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 +1786,7 @@ void Model::fdJrzdsigma(const int ie, const int nroots, const realtype t, } if (alwaysCheckFinite) { - amici::checkFinite(dJrzdsigma, "dJrzdsigma"); + app->checkFinite(dJrzdsigma, "dJrzdsigma"); } } @@ -1790,27 +1796,24 @@ void Model::fw(const realtype t, const realtype *x) { h.data(), total_cl.data()); if (alwaysCheckFinite) { - amici::checkFinite(w, "w"); + app->checkFinite(w, "w"); } } void Model::fdwdp(const realtype t, const realtype *x) { fw(t, x); - std::fill(dwdp.begin(), dwdp.end(), 0.0); - if (wasPythonGenerated()) { - realtype *stcl = nullptr; + if (pythonGenerated) { + dwdp.reset(); // avoid bad memory access when slicing if (!nw) return; - - for (int ip = 0; ip < nplist(); ++ip) { - if (ncl()) - stcl = &stotal_cl.at(plist(ip) * ncl()); - fdwdp(&dwdp.at(nw * ip), t, x, unscaledParameters.data(), - fixedParameters.data(), h.data(), w.data(), total_cl.data(), - stcl, plist_[ip]); - } + + fdwdp_colptrs(dwdp.indexptrs()); + fdwdp_rowvals(dwdp.indexvals()); + fdwdp(dwdp.data(), t, x, unscaledParameters.data(), fixedParameters.data(), + h.data(), w.data(), total_cl.data(), stotal_cl.data()); + } else { // matlab generated fdwdp(dwdp.data(), t, x, unscaledParameters.data(), @@ -1819,20 +1822,21 @@ void Model::fdwdp(const realtype t, const realtype *x) { } if (alwaysCheckFinite) { - amici::checkFinite(dwdp, "dwdp"); + app->checkFinite(gsl::make_span(dwdp.get()), "dwdp"); } } void Model::fdwdx(const realtype t, const realtype *x) { fw(t, x); dwdx.reset(); + + fdwdx_colptrs(dwdx.indexptrs()); + fdwdx_rowvals(dwdx.indexvals()); fdwdx(dwdx.data(), t, x, unscaledParameters.data(), fixedParameters.data(), h.data(), w.data(), total_cl.data()); - fdwdx_colptrs(dwdx.indexptrs()); - 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/deps/AMICI/src/model_dae.cpp b/deps/AMICI/src/model_dae.cpp index 339a802b7..d0384ae0b 100644 --- a/deps/AMICI/src/model_dae.cpp +++ b/deps/AMICI/src/model_dae.cpp @@ -34,14 +34,14 @@ void Model_DAE::fJSparse(realtype t, realtype cj, N_Vector x, N_Vector dx, fixedParameters.data(), h.data(), cj, N_VGetArrayPointer(dx), w.data(), dwdx.data()); } - + void Model_DAE::fJSparseB(SUNMatrixContent_Sparse /*JSparseB*/, - const realtype /*t*/, const realtype * /*x*/, - const double * /*p*/, const double * /*k*/, - const realtype * /*h*/, const realtype /*cj*/, - const realtype * /*xB*/, const realtype * /*dx*/, - const realtype * /*dxB*/, const realtype * /*w*/, - const realtype * /*dwdx*/) { + const realtype /*t*/, const realtype * /*x*/, + const double * /*p*/, const double * /*k*/, + const realtype * /*h*/, const realtype /*cj*/, + const realtype * /*xB*/, const realtype * /*dx*/, + const realtype * /*dxB*/, const realtype * /*w*/, + const realtype * /*dwdx*/) { throw AmiException("Requested functionality is not supported as %s " "is not implemented for this model!", __func__); @@ -104,12 +104,20 @@ void Model_DAE::fJDiag(const realtype t, AmiVector &JDiag, void Model_DAE::fdxdotdp(const realtype t, const N_Vector x, const N_Vector dx) { auto x_pos = computeX_pos(x); - fdwdp(t, N_VGetArrayPointer(x_pos)); - for (int ip = 0; ip < nplist(); ip++) { - N_VConst(0.0, dxdotdp.getNVector(ip)); - fdxdotdp(dxdotdp.data(ip), t, N_VGetArrayPointer(x_pos), - unscaledParameters.data(), fixedParameters.data(), h.data(), - plist_[ip], N_VGetArrayPointer(dx), w.data(), dwdp.data()); + + if (pythonGenerated) { + // python generated, not yet implemented for DAEs + throw AmiException("Wrapping of DAEs is not yet implemented from Python"); + } else { + // matlab generated + fdwdp(t, N_VGetArrayPointer(x_pos)); // Why is it x_pos here and x ind model_ode.cpp? + + for (int ip = 0; ip < nplist(); ip++) { + N_VConst(0.0, dxdotdp.getNVector(ip)); + fdxdotdp(dxdotdp.data(ip), t, N_VGetArrayPointer(x_pos), + unscaledParameters.data(), fixedParameters.data(), h.data(), + plist_[ip], N_VGetArrayPointer(dx), w.data(), dwdp.data()); + } } } @@ -243,7 +251,15 @@ void Model_DAE::fsxdot(realtype t, N_Vector x, N_Vector dx, int ip, N_Vector sx, fdxdotdp(t, x, dx); fJSparse(t, 0.0, x, dx, J.get()); } - N_VScale(1.0, dxdotdp.getNVector(ip), sxdot); + + if (pythonGenerated) { + // python generated, not yet implemented for DAEs + throw AmiException("Wrapping of DAEs is not yet implemented from Python"); + } else { + /* copy dxdotdp over */ + N_VScale(1.0, dxdotdp.getNVector(ip), sxdot); + } + J.multiply(sxdot, sx); N_VScale(-1.0, sdx, sdx); M.multiply(sxdot, sdx); diff --git a/deps/AMICI/src/model_header.ODE_template.h b/deps/AMICI/src/model_header.ODE_template.h index bde8831f8..df0523811 100644 --- a/deps/AMICI/src/model_header.ODE_template.h +++ b/deps/AMICI/src/model_header.ODE_template.h @@ -46,13 +46,20 @@ TPL_DJYDY_DEF TPL_DJYDY_COLPTRS_DEF TPL_DJYDY_ROWVALS_DEF TPL_DWDP_DEF +TPL_DWDP_COLPTRS_DEF +TPL_DWDP_ROWVALS_DEF TPL_DWDX_DEF TPL_DWDX_COLPTRS_DEF TPL_DWDX_ROWVALS_DEF TPL_DXDOTDW_DEF TPL_DXDOTDW_COLPTRS_DEF TPL_DXDOTDW_ROWVALS_DEF -TPL_DXDOTDP_DEF +TPL_DXDOTDP_EXPLICIT_DEF +TPL_DXDOTDP_EXPLICIT_COLPTRS_DEF +TPL_DXDOTDP_EXPLICIT_ROWVALS_DEF +TPL_DXDOTDP_IMPLICIT_COLPTRS_DEF +TPL_DXDOTDP_IMPLICIT_ROWVALS_DEF + extern void dydx_TPL_MODELNAME(realtype *dydx, const realtype t, const realtype *x, const realtype *p, const realtype *k, const realtype *h, @@ -124,7 +131,10 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { std::vector{TPL_FIXED_PARAMETERS}, // fixedParameters std::vector{}, // plist std::vector(TPL_NX_SOLVER, 0.0), // idlist - std::vector{} // z2event + std::vector{}, // z2event + true, // pythonGenerated + TPL_NDXDOTDP_EXPLICIT, // ndxdotdp_explicit + TPL_NDXDOTDP_IMPLICIT // ndxdotdp_implicit ) {} /** @@ -444,17 +454,31 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { TPL_DJYDY_ROWVALS_IMPL TPL_DWDP_IMPL + + TPL_DWDP_COLPTRS_IMPL + + TPL_DWDP_ROWVALS_IMPL TPL_DWDX_IMPL + TPL_DWDX_COLPTRS_IMPL + + TPL_DWDX_ROWVALS_IMPL + TPL_DXDOTDW_IMPL - TPL_DXDOTDW_COLPTRS_IMPL - TPL_DXDOTDW_ROWVALS_IMPL - TPL_DXDOTDP_IMPL + TPL_DXDOTDP_EXPLICIT_IMPL + + TPL_DXDOTDP_EXPLICIT_COLPTRS_IMPL + + TPL_DXDOTDP_EXPLICIT_ROWVALS_IMPL + TPL_DXDOTDP_IMPLICIT_COLPTRS_IMPL + + TPL_DXDOTDP_IMPLICIT_ROWVALS_IMPL + /** model specific implementation of fdydx * @param dydx partial derivative of observables y w.r.t. model states x * @param t current time @@ -801,10 +825,6 @@ class Model_TPL_MODELNAME : public amici::Model_ODE { virtual const std::string getAmiciCommit() const override { return "TPL_AMICI_COMMIT_STRING"; } - - virtual bool wasPythonGenerated() const override { - return true; - } }; #endif /* _amici_TPL_MODELNAME_h */ diff --git a/deps/AMICI/src/model_ode.cpp b/deps/AMICI/src/model_ode.cpp index 4284c801d..02389c534 100644 --- a/deps/AMICI/src/model_ode.cpp +++ b/deps/AMICI/src/model_ode.cpp @@ -1,3 +1,4 @@ +#include #include "amici/model_ode.h" #include "amici/solver_cvodes.h" @@ -27,7 +28,7 @@ void Model_ODE::fJSparse(realtype t, N_Vector x, SUNMatrix J) { auto x_pos = computeX_pos(x); fdwdx(t, N_VGetArrayPointer(x_pos)); SUNMatZero(J); - if (wasPythonGenerated()) { + if (pythonGenerated) { fJSparse(SM_DATA_S(J), t, N_VGetArrayPointer(x_pos), unscaledParameters.data(), fixedParameters.data(), h.data(), w.data(), dwdx.data()); @@ -87,30 +88,41 @@ void Model_ODE::fJDiag(const realtype t, AmiVector &JDiag, } void Model_ODE::fdxdotdw(const realtype t, const N_Vector x) { - dxdotdw.reset(); - auto x_pos = computeX_pos(x); - fdxdotdw(dxdotdw.data(), t, N_VGetArrayPointer(x_pos), - unscaledParameters.data(), fixedParameters.data(), h.data(), - w.data()); - fdxdotdw_colptrs(dxdotdw.indexptrs()); - fdxdotdw_rowvals(dxdotdw.indexvals()); + if (nw > 0 && ndxdotdw > 0) { + auto x_pos = computeX_pos(x); + dxdotdw.reset(); + fdxdotdw_colptrs(dxdotdw.indexptrs()); + fdxdotdw_rowvals(dxdotdw.indexvals()); + fdxdotdw(dxdotdw.data(), t, N_VGetArrayPointer(x_pos), + unscaledParameters.data(), fixedParameters.data(), h.data(), + w.data()); + } } void Model_ODE::fdxdotdp(const realtype t, const N_Vector x) { auto x_pos = computeX_pos(x); - fdwdp(t, N_VGetArrayPointer(x)); - if (wasPythonGenerated()) { + fdwdp(t, N_VGetArrayPointer(x_pos)); + fdxdotdw(t, x_pos); + + if (pythonGenerated) { // python generated - fdxdotdw(t, x); - for (int ip = 0; ip < nplist(); ip++) { - N_VConst(0.0, dxdotdp.getNVector(ip)); - fdxdotdp(dxdotdp.data(ip), t, N_VGetArrayPointer(x_pos), - unscaledParameters.data(), fixedParameters.data(), - h.data(), plist_[ip], w.data()); - if (nw > 0) - dxdotdw.multiply( - gsl::span(dxdotdp.data(ip), nx_solver), - gsl::span(&dwdp.at(nw * ip), nw)); + if (ndxdotdp_explicit > 0) { + dxdotdp_explicit.reset(); + fdxdotdp_explicit_colptrs(dxdotdp_explicit.indexptrs()); + fdxdotdp_explicit_rowvals(dxdotdp_explicit.indexvals()); + fdxdotdp_explicit(dxdotdp_explicit.data(), + t, N_VGetArrayPointer(x_pos), + unscaledParameters.data(), + fixedParameters.data(), + h.data(), w.data()); + } + if (nw > 0 && ndxdotdp_implicit > 0) { + /* Sparse matrix multiplication + dxdotdp_implicit += dxdotdw * dwdp */ + dxdotdp_implicit.reset(); + fdxdotdp_implicit_colptrs(dxdotdp_implicit.indexptrs()); + fdxdotdp_implicit_rowvals(dxdotdp_implicit.indexvals()); + dxdotdw.sparse_multiply(&dxdotdp_implicit, &dwdp); } } else { // matlab generated @@ -232,10 +244,34 @@ void Model_ODE::fdxdotdp(realtype * /*dxdotdp*/, const realtype /*t*/, __func__); // not implemented } -void Model_ODE::fdxdotdp(realtype * /*dxdotdp*/, const realtype /*t*/, - const realtype * /*x*/, const realtype * /*p*/, - const realtype * /*k*/, const realtype * /*h*/, - const int /*ip*/, const realtype * /*w*/) { +void Model_ODE::fdxdotdp_explicit(realtype * /*dxdotdp_explicit*/, const realtype /*t*/, + const realtype * /*x*/, const realtype * /*p*/, + const realtype * /*k*/, const realtype * /*h*/, + const realtype * /*w*/) { + throw AmiException("Requested functionality is not supported as %s " + "is not implemented for this model!", + __func__); // not implemented +} + +void Model_ODE::fdxdotdp_explicit_colptrs(sunindextype * /*indexptrs*/) { + throw AmiException("Requested functionality is not supported as %s " + "is not implemented for this model!", + __func__); // not implemented +} + +void Model_ODE::fdxdotdp_explicit_rowvals(sunindextype * /*indexvals*/) { + throw AmiException("Requested functionality is not supported as %s " + "is not implemented for this model!", + __func__); // not implemented +} + +void Model_ODE::fdxdotdp_implicit_colptrs(sunindextype * /*indexptrs*/) { + throw AmiException("Requested functionality is not supported as %s " + "is not implemented for this model!", + __func__); // not implemented +} + +void Model_ODE::fdxdotdp_implicit_rowvals(sunindextype * /*indexvals*/) { throw AmiException("Requested functionality is not supported as %s " "is not implemented for this model!", __func__); // not implemented @@ -277,7 +313,7 @@ void Model_ODE::fJSparseB(realtype t, N_Vector x, N_Vector xB, auto x_pos = computeX_pos(x); fdwdx(t, N_VGetArrayPointer(x_pos)); SUNMatZero(JB); - if (wasPythonGenerated()) { + if (pythonGenerated) { fJSparseB(SM_DATA_S(JB), t, N_VGetArrayPointer(x_pos), unscaledParameters.data(), fixedParameters.data(), h.data(), N_VGetArrayPointer(xB), w.data(), dwdx.data()); @@ -314,17 +350,29 @@ void Model_ODE::fxBdot(realtype t, N_Vector x, N_Vector xB, N_Vector xBdot) { } void Model_ODE::fqBdot(realtype t, N_Vector x, N_Vector xB, N_Vector qBdot) { + /* initialize with zeros */ N_VConst(0.0, qBdot); fdxdotdp(t, x); - for (int ip = 0; ip < nplist(); ip++) { - for (int ix = 0; ix < nxtrue_solver; ix++) - NV_Ith_S(qBdot, ip * nJ) -= NV_Ith_S(xB, ix) * dxdotdp.at(ix, ip); - // second order part - for (int iJ = 1; iJ < nJ; iJ++) + + if (pythonGenerated) { + /* call multiplication */ + if (ndxdotdp_explicit > 0) + dxdotdp_explicit.multiply(qBdot, xB, plist_, true); + if (ndxdotdp_implicit > 0) + dxdotdp_implicit.multiply(qBdot, xB, plist_, true); + N_VScale(-1.0, qBdot, qBdot); + } else { + /* was matlab generated */ + for (int ip = 0; ip < nplist(); ip++) { for (int ix = 0; ix < nxtrue_solver; ix++) - NV_Ith_S(qBdot, ip * nJ + iJ) -= + NV_Ith_S(qBdot, ip * nJ) -= NV_Ith_S(xB, ix) * dxdotdp.at(ix, ip); + // second order part + for (int iJ = 1; iJ < nJ; iJ++) + for (int ix = 0; ix < nxtrue_solver; ix++) + NV_Ith_S(qBdot, ip * nJ + iJ) -= NV_Ith_S(xB, ix) * dxdotdp.at(ix + iJ * nxtrue_solver, ip) + NV_Ith_S(xB, ix + iJ * nxtrue_solver) * dxdotdp.at(ix, ip); + } } } @@ -337,13 +385,43 @@ void Model_ODE::fsxdot(const realtype t, const AmiVector &x, void Model_ODE::fsxdot(realtype t, N_Vector x, int ip, N_Vector sx, N_Vector sxdot) { + + /* sxdot is just the total derivative d(xdot)dp, + so we just call dxdotdp and copy the stuff over */ if (ip == 0) { // we only need to call this for the first parameter index will be // the same for all remaining fdxdotdp(t, x); fJSparse(t, x, J.get()); } - N_VScale(1.0, dxdotdp.getNVector(ip), sxdot); + if (pythonGenerated) { + /* copy dxdotdp and the implicit version over */ + // initialize + N_VConst(0.0, sxdot); + realtype *sxdot_tmp = N_VGetArrayPointer(sxdot); + + // copy explicit version + if (ndxdotdp_explicit > 0) { + auto col_exp = dxdotdp_explicit.indexptrs(); + auto row_exp = dxdotdp_explicit.indexvals(); + auto data_exp_ptr = dxdotdp_explicit.data(); + for (sunindextype i = col_exp[plist(ip)]; i < col_exp[plist(ip) + 1]; ++i) + sxdot_tmp[row_exp[i]] += data_exp_ptr[i]; + } + + // copy implicit version + if (ndxdotdp_implicit > 0) { + auto col_imp = dxdotdp_implicit.indexptrs(); + auto row_imp = dxdotdp_implicit.indexvals(); + auto data_imp_ptr = dxdotdp_implicit.data(); + for (sunindextype i = col_imp[plist(ip)]; i < col_imp[plist(ip) + 1]; ++i) + sxdot_tmp[row_imp[i]] += data_imp_ptr[i]; + } + + } else { + /* copy dxdotdp over */ + N_VScale(1.0, dxdotdp.getNVector(ip), sxdot); + } J.multiply(sxdot, sx); } diff --git a/deps/AMICI/src/newton_solver.cpp b/deps/AMICI/src/newton_solver.cpp index 26d1247ba..198a4a1c9 100644 --- a/deps/AMICI/src/newton_solver.cpp +++ b/deps/AMICI/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; } @@ -93,14 +96,41 @@ void NewtonSolver::getStep(int ntry, int nnewt, AmiVector &delta) { void NewtonSolver::computeNewtonSensis(AmiVectorArray &sx) { prepareLinearSystem(0, -1); - model->fdxdotdp(*t, *x, dx); - for (int ip = 0; ip < model->nplist(); ip++) { - for (int ix = 0; ix < model->nx_solver; ix++) { - sx.at(ix,ip) = -model->dxdotdp.at(ix, ip); + if (model->pythonGenerated) { + for (int ip = 0; ip < model->nplist(); ip++) { + N_VConst(0.0, sx.getNVector(ip)); + + // copy explicit version + if (model->ndxdotdp_explicit > 0) { + auto col = model->dxdotdp_explicit.indexptrs(); + auto row = model->dxdotdp_explicit.indexvals(); + auto data_ptr = model->dxdotdp_explicit.data(); + for (sunindextype iCol = col[model->plist(ip)]; + iCol < col[model->plist(ip) + 1]; ++iCol) + sx.at(row[iCol], ip) -= data_ptr[iCol]; + } + + // copy implicit version + if (model->ndxdotdp_implicit > 0) { + auto col = model->dxdotdp_implicit.indexptrs(); + auto row = model->dxdotdp_implicit.indexvals(); + auto data_ptr = model->dxdotdp_implicit.data(); + for (sunindextype iCol = col[model->plist(ip)]; + iCol < col[model->plist(ip) + 1]; ++iCol) + sx.at(row[iCol], ip) -= data_ptr[iCol]; + } + + solveLinearSystem(sx[ip]); + } + } else { + for (int ip = 0; ip < model->nplist(); ip++) { + for (int ix = 0; ix < model->nx_solver; ix++) + sx.at(ix,ip) = -model->dxdotdp.at(ix, ip); + + solveLinearSystem(sx[ip]); } - solveLinearSystem(sx[ip]); } } /* ------------------------------------------------------------------------- */ @@ -212,7 +242,8 @@ void NewtonSolverIterative::prepareLinearSystem(int ntry, int nnewt) { newton_try = ntry; i_newton = nnewt; if (nnewt == -1) { - throw AmiException("Linear solver SPBCG does not support sensitivity computation for steady state problems."); + throw AmiException("Linear solver SPBCG does not support sensitivity " + "computation for steady state problems."); } } diff --git a/deps/AMICI/src/solver.cpp b/deps/AMICI/src/solver.cpp index cad80533a..f86a896fa 100644 --- a/deps/AMICI/src/solver.cpp +++ b/deps/AMICI/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/deps/AMICI/src/solver_cvodes.cpp b/deps/AMICI/src/solver_cvodes.cpp index 56441f782..34fbfe9ef 100644 --- a/deps/AMICI/src/solver_cvodes.cpp +++ b/deps/AMICI/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/deps/AMICI/src/solver_idas.cpp b/deps/AMICI/src/solver_idas.cpp index 670b22805..a8b9b51e4 100644 --- a/deps/AMICI/src/solver_idas.cpp +++ b/deps/AMICI/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/deps/AMICI/src/steadystateproblem.cpp b/deps/AMICI/src/steadystateproblem.cpp index 8ea73c352..3b3890d77 100644 --- a/deps/AMICI/src/steadystateproblem.cpp +++ b/deps/AMICI/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/deps/AMICI/src/sundials_matrix_wrapper.cpp b/deps/AMICI/src/sundials_matrix_wrapper.cpp index dd107c5ef..496972f08 100644 --- a/deps/AMICI/src/sundials_matrix_wrapper.cpp +++ b/deps/AMICI/src/sundials_matrix_wrapper.cpp @@ -142,6 +142,19 @@ sunindextype SUNMatrixWrapper::columns() const { } } +sunindextype SUNMatrixWrapper::nonzeros() const { + if (!matrix) + return 0; + + switch (SUNMatGetID(matrix)) { + case SUNMATRIX_SPARSE: + return SM_NNZ_S(matrix); + default: + throw std::domain_error("Non-zeros property only available for " + "sparse matrices"); + } +} + sunindextype *SUNMatrixWrapper::indexvals() const { return indexvals_ptr; } @@ -166,32 +179,36 @@ void SUNMatrixWrapper::multiply(N_Vector c, const_N_Vector b) const { gsl::make_span(NV_DATA_S(b), NV_LENGTH_S(b))); } -void SUNMatrixWrapper::multiply(gsl::span c, gsl::span b) const { +void SUNMatrixWrapper::multiply(gsl::span c, + gsl::span b) const { 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 +216,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]]; @@ -216,6 +233,135 @@ void SUNMatrixWrapper::multiply(gsl::span c, gsl::span } } +void SUNMatrixWrapper::multiply(N_Vector c, + const N_Vector b, + gsl::span cols, + bool transpose) const { + multiply(gsl::make_span(NV_DATA_S(c), NV_LENGTH_S(c)), + gsl::make_span(NV_DATA_S(b), NV_LENGTH_S(b)), + cols, transpose); +} + +void SUNMatrixWrapper::multiply(gsl::span c, + gsl::span b, + gsl::span cols, + bool transpose) const { + if (!matrix) + return; + + sunindextype nrows = rows(); + sunindextype ncols = columns(); + + if (transpose) { + if (c.size() != cols.size()) + throw std::invalid_argument("Dimension mismatch between number of cols " + "in index vector cols (" + std::to_string(ncols) + + ") and elements in c (" + std::to_string(c.size()) + + " ), when using transposed A"); + + if (static_cast(b.size()) != nrows) + throw std::invalid_argument("Dimension mismatch between number of rows " + "in A (" + std::to_string(nrows) + ") and " + "elements in b (" + std::to_string(b.size()) + + "), when using transposed A"); + } else { + if (static_cast(c.size()) != nrows) + throw std::invalid_argument("Dimension mismatch between number of rows " + "in A (" + std::to_string(nrows) + ") and " + "elements in c (" + std::to_string(c.size()) + + ")"); + + if (static_cast(b.size()) != ncols) + throw std::invalid_argument("Dimension mismatch between number of cols " + "in A (" + std::to_string(ncols) + + ") and elements in b (" + + std::to_string(b.size()) + ")"); + } + + if (SUNMatGetID(matrix) != SUNMATRIX_SPARSE) + throw std::invalid_argument("Reordered multiply only implemented for " + "sparse matrices, but A is not sparse"); + + if (sparsetype() != CSC_MAT) + throw std::invalid_argument("Reordered multiply only implemented for " + "matrix type CSC, but A is not of type CSC"); + + /* Carry out actual multiplication */ + if (transpose) { + for (int i = 0; i < (int)cols.size(); ++i) + for (sunindextype k = indexptrs_ptr[cols[i]]; + k < indexptrs_ptr[cols[i] + 1]; ++k) + c[i] += data_ptr[k] * b[indexvals_ptr[k]]; + } else { + for (sunindextype i = 0; i < ncols; ++i) + for (sunindextype k = indexptrs_ptr[cols[i]]; + k < indexptrs_ptr[cols[i] + 1]; ++k) + c[indexvals_ptr[k]] += data_ptr[k] * b[i]; + } +} + + +void SUNMatrixWrapper::sparse_multiply(SUNMatrixWrapper *C, + SUNMatrixWrapper *B) const { + if (!matrix) + return; + + sunindextype nrows = rows(); + sunindextype ncols = columns(); + + if (SUNMatGetID(matrix) != SUNMATRIX_SPARSE) + throw std::invalid_argument("Matrix A not sparse in sparse_multiply"); + + if (sparsetype() != CSC_MAT) + throw std::invalid_argument("Matrix A not of type CSC_MAT"); + + if (SUNMatGetID(B->matrix) != SUNMATRIX_SPARSE) + throw std::invalid_argument("Matrix B not sparse in sparse_multiply"); + + if (B->sparsetype() != CSC_MAT) + throw std::invalid_argument("Matrix B not of type CSC_MAT"); + + if (SUNMatGetID(C->matrix) != SUNMATRIX_SPARSE) + throw std::invalid_argument("Matrix C not sparse in sparse_multiply"); + + if (C->sparsetype() != CSC_MAT) + throw std::invalid_argument("Matrix C not of type CSC_MAT"); + + if (C->rows() != nrows) + throw std::invalid_argument("Dimension mismatch between number of rows " + "in A (" + std::to_string(nrows) + ") and " + "number of rows in C (" + + std::to_string((int)C->rows()) + ")"); + + if (B->rows() != ncols) + throw std::invalid_argument("Dimension mismatch between number of rows " + "in A (" + std::to_string(ncols) + + ") and number of cols in B (" + + std::to_string((int)B->rows()) + ")"); + + /* Carry out actual multiplication */ + sunindextype iC_data = 0; + for (sunindextype iC_col = 0; iC_col < C->columns(); ++iC_col) { + for(sunindextype iC_row = C->indexptrs_ptr[iC_col]; + iC_row < C->indexptrs_ptr[iC_col + 1]; ++iC_row) { + // Current entry in C: (C.indexvals[iC_row], iC_col) + for(sunindextype iB_row = B->indexptrs_ptr[iC_col]; + iB_row < B->indexptrs_ptr[iC_col + 1]; ++iB_row) { + // Loop over column iC_col in B + sunindextype iA_col = B->indexvals_ptr[iB_row]; + for (sunindextype iA_row = indexptrs_ptr[iA_col]; + iA_row < indexptrs_ptr[iA_col + 1]; ++iA_row) + // loop over entries in column col_in_A + // if two entries match together, carry out multiplication + if (indexvals_ptr[iA_row] == C->indexvals_ptr[iC_row]) + C->data_ptr[iC_data] += + data_ptr[iA_row] * B->data_ptr[iB_row]; + } + iC_data++; + } + } +} + void SUNMatrixWrapper::zero() { if(int res = SUNMatZero(matrix)) @@ -244,7 +390,7 @@ void SUNMatrixWrapper::update_ptrs() { data_ptr = SM_DATA_B(matrix); break; case SUNMATRIX_CUSTOM: - throw std::domain_error("Amici currently does not support" + throw std::domain_error("Amici currently does not support " "custom matrix types."); } } diff --git a/deps/AMICI/swig/amici.i b/deps/AMICI/swig/amici.i index f57eeb8e5..359520677 100644 --- a/deps/AMICI/swig/amici.i +++ b/deps/AMICI/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/deps/AMICI/tests/performance/check_time.sh b/deps/AMICI/tests/performance/check_time.sh new file mode 100755 index 000000000..d8f440e1f --- /dev/null +++ b/deps/AMICI/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/deps/AMICI/tests/performance/reference.yml b/deps/AMICI/tests/performance/reference.yml new file mode 100644 index 000000000..f6ab156b9 --- /dev/null +++ b/deps/AMICI/tests/performance/reference.yml @@ -0,0 +1,10 @@ +# Reference wall times (seconds) with some buffer +create_sdist: 25 +install_sdist: 140 +petab_import: 1920 # ^= 32' +install_model: 400 +forward_simulation: 2 +forward_sensitivities: 6 +adjoint_sensitivities: 4 +forward_simulation_non_optimal_parameters: 2 +adjoint_sensitivities_non_optimal_parameters: 5 diff --git a/deps/AMICI/tests/performance/test.py b/deps/AMICI/tests/performance/test.py new file mode 100755 index 000000000..535a18f72 --- /dev/null +++ b/deps/AMICI/tests/performance/test.py @@ -0,0 +1,52 @@ +#!/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) + elif arg == 'forward_simulation_non_optimal_parameters': + tmpPar = model.getParameters() + model.setParameters([0.1 for iPar in tmpPar]) + solver.setSensitivityMethod(amici.SensitivityMethod_none) + solver.setSensitivityOrder(amici.SensitivityOrder_none) + elif arg == 'adjoint_sensitivities_non_optimal_parameters': + tmpPar = model.getParameters() + model.setParameters([0.1 for iPar in tmpPar]) + 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/deps/AMICI/tests/testMisc.py b/deps/AMICI/tests/testMisc.py index 877afceec..7aa21e756 100755 --- a/deps/AMICI/tests/testMisc.py +++ b/deps/AMICI/tests/testMisc.py @@ -56,7 +56,7 @@ def test_csc_matrix_empty(self): symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix = \ amici.ode_export.csc_matrix(matrix, 'a') print(symbolColPtrs, symbolRowVals, sparseList, symbolList, sparseMatrix) - assert symbolColPtrs == [0] + assert symbolColPtrs == [] assert symbolRowVals == [] assert sparseList == sp.Matrix(0, 0, []) assert symbolList == [] @@ -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/deps/AMICI/tests/testPYSB.py b/deps/AMICI/tests/testPYSB.py index 1092fb201..a60670b24 100755 --- a/deps/AMICI/tests/testPYSB.py +++ b/deps/AMICI/tests/testPYSB.py @@ -93,10 +93,12 @@ def test_compare_to_sbml_import(self): rdata_sbml = get_results(model_sbml, edata) # check if preequilibration fixed parameters are correctly applied: - for rdata, model in zip([rdata_sbml, rdata_pysb], - [model_sbml, model_pysb]): + for rdata, model, importer in zip([rdata_sbml, rdata_pysb], + [model_sbml, model_pysb], + ['sbml', 'pysb']): # check equilibrium fixed parameters - with self.subTest(fixed_pars='preequilibration'): + with self.subTest(fixed_pars='preequilibration', + importer=importer): self.assertTrue(np.isclose( [ sum(rdata["x_ss"][[1, 3]]), @@ -111,7 +113,8 @@ def test_compare_to_sbml_import(self): model.getParameterByName('PROT_0'), atol=1e-6, rtol=1e-6 )) - with self.subTest(fixed_pars='simulation'): + with self.subTest(fixed_pars='presimulation', + importer=importer): # check reinitialization with fixed parameter after # presimulation self.assertTrue(np.isclose( diff --git a/deps/AMICI/tests/testSBMLSuite.py b/deps/AMICI/tests/testSBMLSuite.py index 0b8c07c52..0a1816c8e 100755 --- a/deps/AMICI/tests/testSBMLSuite.py +++ b/deps/AMICI/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/deps/AMICI/version.txt b/deps/AMICI/version.txt index f25c43cdc..55a6d615b 100644 --- a/deps/AMICI/version.txt +++ b/deps/AMICI/version.txt @@ -1 +1 @@ -0.10.13 +0.10.16 diff --git a/doc/parpe_with_charliecloud.md b/doc/parpe_with_charliecloud.md index a2c556e99..7da6ee4ad 100644 --- a/doc/parpe_with_charliecloud.md +++ b/doc/parpe_with_charliecloud.md @@ -21,6 +21,7 @@ This will create the parPE base image *from parPE from github* (takes about 10'): cd container/charliecloud/parpe_base + git archive -v -o parpe.tar.gz --format=tar.gz HEAD ch-build -t parpe . Export image to charliecloud archive in the current directory: diff --git a/examples/parpeamici/steadystate/run-examples.sh b/examples/parpeamici/steadystate/run-examples.sh index 844ce247c..b0e92ac1e 100755 --- a/examples/parpeamici/steadystate/run-examples.sh +++ b/examples/parpeamici/steadystate/run-examples.sh @@ -10,6 +10,10 @@ HDF5_FILE_TEST=$2 MPIEXEC="mpiexec --oversubscribe -n 5" # Allow running as root in docker grep docker /proc/1/cgroup -qa && MPIEXEC="${MPIEXEC} --allow-run-as-root" +# If we are running in docker, we generally don't have SYS_PTRACE permissions +# and thus, cannot use vader. Also disable Infiniband. +grep docker /proc/1/cgroup -qa && mpiexec --version | grep open-mpi && MPIEXEC="${MPIEXEC} --oversubscribe --mca btl_vader_single_copy_mechanism none --mca btl ^openib -mca oob_tcp_if_include lo --mca btl_tcp_if_include lo --mca orte_base_help_aggregate 0" + rm -f test.log diff --git a/include/parpeamici/amiciMisc.h b/include/parpeamici/amiciMisc.h index f18abadfd..cedb4b903 100644 --- a/include/parpeamici/amiciMisc.h +++ b/include/parpeamici/amiciMisc.h @@ -3,11 +3,6 @@ namespace parpe { - -void printAmiciErrMsgIdAndTxt(const char *identifier, const char *format, ...); - -void printAmiciWarnMsgIdAndTxt(const char *identifier, const char *format, ...); - using amici::getUnscaledParameter; using amici::getScaledParameter; diff --git a/include/parpeamici/amiciSimulationRunner.h b/include/parpeamici/amiciSimulationRunner.h index eab00f36e..1b3afb1df 100644 --- a/include/parpeamici/amiciSimulationRunner.h +++ b/include/parpeamici/amiciSimulationRunner.h @@ -73,10 +73,13 @@ class AmiciSimulationRunner { /** * @brief SimulationRunner + * @param optimizationParameters + * @param sensitivityOrder + * @param conditionIndices * @param callbackJobFinished Function which is called after any finished simulation. May be nullptr. * @param aggregate Function which is called after all simulations are completed. May be nullptr. + * @param logPrefix */ - AmiciSimulationRunner(const std::vector &optimizationParameters, amici::SensitivityOrder sensitivityOrder, const std::vector &conditionIndices, @@ -90,6 +93,7 @@ class AmiciSimulationRunner { /** * @brief Dispatch simulation jobs using LoadBalancerMaster * @param loadBalancer + * @param maxSimulationsPerPackage * @return */ int runDistributedMemory(LoadBalancerMaster *loadBalancer, const int maxSimulationsPerPackage = 1); @@ -99,6 +103,7 @@ class AmiciSimulationRunner { * @brief Runs simulations within the same thread. Mostly intended for * debugging. * @param messageHandler + * @param sequential Run sequential (not in parallel) * @return */ int runSharedMemory(const messageHandlerFunc& messageHandler, bool sequential = false); diff --git a/include/parpeamici/hierarchicalOptimization.h b/include/parpeamici/hierarchicalOptimization.h index 1213adec0..8216a4f53 100644 --- a/include/parpeamici/hierarchicalOptimization.h +++ b/include/parpeamici/hierarchicalOptimization.h @@ -457,10 +457,14 @@ double getDefaultOffsetParameter(amici::ParameterScaling scaling); * See Supplement 1.1 of [1]. * * [1] Loos, Krause, Hasenauer. Hierarchical optimization for the efficient parametrization of ODE models. - * @param + * + * @param scalingIdx + * @param modelOutputsUnscaled + * @param measurements + * @param scalingReader + * @param numObservables * @return */ - double computeAnalyticalScalings(int scalingIdx, const std::vector > &modelOutputsUnscaled, const std::vector > &measurements, @@ -506,7 +510,9 @@ std::vector spliceParameters(const gsl::span reducedParame /** * @brief Compute negative log-likelihood for normal distribution based on the model outputs and measurements for multiple conditions. + * @param measurements * @param modelOutputsScaled + * @param sigmas * @return */ double computeNegLogLikelihood(std::vector > const& measurements, @@ -515,7 +521,9 @@ double computeNegLogLikelihood(std::vector > const& measurem /** * @brief Compute negative log-likelihood for normal distribution based on the model outputs and measurements for a single condition. + * @param measurements * @param modelOutputsScaled + * @param sigmas * @return */ double computeNegLogLikelihood(std::vector const& measurements, diff --git a/include/parpeamici/multiConditionProblem.h b/include/parpeamici/multiConditionProblem.h index e2f366aed..7c9c1c516 100644 --- a/include/parpeamici/multiConditionProblem.h +++ b/include/parpeamici/multiConditionProblem.h @@ -50,10 +50,17 @@ AmiciSimulationRunner::AmiciResultPackageSimple runAndLogSimulation(amici::Solve /** * @brief Run simulations (no gradient) with given parameters and collect * model outputs + * @param dataProvider + * @param loadBalancer + * @param maxSimulationsPerPackage + * @param resultWriter + * @param logLineSearch * @param parameters Model parameters for simulation * @param modelOutput in: some vector reference, will be resized. * output: Vector of double vectors containing AMICI ReturnData::y * (nt x ny, column-major) + * @param logger + * @param cpuTime * @return Simulation status */ FunctionEvaluationStatus getModelOutputs( @@ -68,6 +75,9 @@ FunctionEvaluationStatus getModelOutputs( /** * @brief Callback function for LoadBalancer + * @param dataProvider + * @param resultWriter + * @param logLineSearch * @param buffer In/out: message buffer * @param jobId: In: Identifier of the job (unique up to INT_MAX) */ diff --git a/include/parpeamici/optimizationApplication.h b/include/parpeamici/optimizationApplication.h index 7b35cc5ef..6048a7ef9 100644 --- a/include/parpeamici/optimizationApplication.h +++ b/include/parpeamici/optimizationApplication.h @@ -166,6 +166,7 @@ class OptimizationApplication { /** * @brief CPU time for whole application run + * @param file_id * @param timeInSeconds */ void saveTotalCpuTime(hid_t file_id, const double timeInSeconds); diff --git a/include/parpecommon/misc.h b/include/parpecommon/misc.h index d32940c83..0998dae94 100644 --- a/include/parpecommon/misc.h +++ b/include/parpecommon/misc.h @@ -97,11 +97,10 @@ std::string getBacktrace(int nMaxFrames = 20); double randDouble(double min, double max); /** - * @brief fillArrayRandomDoubleIndividualInterval Fill "buffer" with "length" + * @brief fillArrayRandomDoubleIndividualInterval Fill "buffer" with * random double values, drawn from an interval [min, max] given for each value. * @param min * @param max - * @param length * @param buffer */ void fillArrayRandomDoubleIndividualInterval(gsl::span min, diff --git a/include/parpeloadbalancer/loadBalancerMaster.h b/include/parpeloadbalancer/loadBalancerMaster.h index 4f80c4099..f13ecedd1 100644 --- a/include/parpeloadbalancer/loadBalancerMaster.h +++ b/include/parpeloadbalancer/loadBalancerMaster.h @@ -27,7 +27,7 @@ struct JobData { } /** auto-assigned (unique number up to MAX_INT) */ - int jobId; + int jobId = -1; /** data to send */ std::vector sendBuffer; diff --git a/include/parpeoptimization/minibatchOptimization.h b/include/parpeoptimization/minibatchOptimization.h index d60f8edf9..1cdaf0ee5 100755 --- a/include/parpeoptimization/minibatchOptimization.h +++ b/include/parpeoptimization/minibatchOptimization.h @@ -221,6 +221,39 @@ class ParameterUpdaterRmsProp: public ParameterUpdater { std::vector oldGradientNormCache; }; +/** + * @brief Minibatch optimizer: Momentum Updater + * A classical gradient based optimizer using a vanilla momentum formula + */ +class ParameterUpdaterMomentum: public ParameterUpdater { +public: + ParameterUpdaterMomentum() = default; + + void updateParameters(double learningRate, + int iteration, + gsl::span gradient, + gsl::span parameters, + gsl::span lowerBounds = gsl::span(), + gsl::span upperBounds = gsl::span()) override; + + void undoLastStep() override; + + void clearCache() override; + + void initialize(unsigned int numParameters) override; + +private: + + /** Rate for memorizing gradient norms (between 0 and 1, high rates mean long memory) */ + double decayRate = 0.8; + + /** Accumulated momentum (decaying average) from last steps */ + std::vector momentum; + + /** Accumulated momentum (decaying average), one step back (if one step must be undone) */ + std::vector oldMomentum; +}; + /** * @brief Minibatch optimizer: Adam Updater * A momentum-based and so-called adaptive mini batching algorithm @@ -266,6 +299,52 @@ class ParameterUpdaterAdam: public ParameterUpdater { std::vector oldGradientCache; }; +/** + * @brief Minibatch optimizer: Adam Classic Updater + * A momentum-based and so-called adaptive mini batching algorithm, + * using the original settings from the literature + */ +class ParameterUpdaterAdamClassic: public ParameterUpdater { +public: + ParameterUpdaterAdamClassic() = default; + + void updateParameters(double learningRate, + int iteration, + gsl::span gradient, + gsl::span parameters, + gsl::span lowerBounds = gsl::span(), + gsl::span upperBounds = gsl::span()) override; + + void undoLastStep() override; + + void clearCache() override; + + void initialize(unsigned int numParameters) override; + +private: + + /** Rate for memorizing gradients (between 0 and 1, high rates mean long memory) */ + double decayRateGradient = 0.9; + + /** Rate for memorizing gradient norms (between 0 and 1, high rates mean long memory) */ + double decayRateGradientNorm = 0.999; + + /** Stabilization factor for gradient normalization (avoid dividing by 0) */ + double delta = 1e-7; + + /** Memorized gradient norms (decaying average) from last steps */ + std::vector gradientNormCache; + + /** Memorized gradient norms (decaying average), one step back (if one step must be undone) */ + std::vector oldGradientNormCache; + + /** Memorized gradients (decaying average) from last steps */ + std::vector gradientCache; + + /** Memorized gradients (decaying average), one step back (if one step must be undone) */ + std::vector oldGradientCache; +}; + /** * @brief Split a vector into batches of the given size. * @@ -515,7 +594,7 @@ class MinibatchOptimizer { Logger *logger, OptimizationReporter *reporter) { - // initialize disgnostic variables + // initialize diagnostic variables int maxSubsequentFails = 10; bool finalFail = false; bool initialFail = false; diff --git a/include/parpeoptimization/optimizationProblem.h b/include/parpeoptimization/optimizationProblem.h index d7c5b76f5..69eb86b7a 100644 --- a/include/parpeoptimization/optimizationProblem.h +++ b/include/parpeoptimization/optimizationProblem.h @@ -219,7 +219,7 @@ int getLocalOptimum(OptimizationProblem *problem); /** * @brief getLocalOptimumThreadWrapper wrapper for using getLocalOptimum with * pThreads. - * @param problem + * @param optimizationProblemVp * @return Pointer to int indicating status. 0: success, != 0: failure */ diff --git a/misc/generateHDF5DataFileFromText.py b/misc/generateHDF5DataFileFromText.py index 87208dabf..45278ec7b 100755 --- a/misc/generateHDF5DataFileFromText.py +++ b/misc/generateHDF5DataFileFromText.py @@ -1135,22 +1135,22 @@ def write_bounds(self): def write_starting_points(self): """ - Write a list of random starting points uniformly sampled from the - parameter bounds. - Parameter bounds need to be written beforehand. + Write a list of random starting points sampled as specified in PEtab + file. """ num_params = self.f['/parameters/parameterNames'].shape[0] num_starting_points = 100 np.random.seed(0) + starting_points = self.f.require_dataset( '/optimizationOptions/randomStarts', [num_params, num_starting_points], 'f8') - lower = self.f['/parameters/lowerBound'][:] - upper = self.f['/parameters/upperBound'][:] - starting_points[:] = np.transpose( - np.random.rand(num_starting_points, num_params) * ( - upper - lower) + lower) + starting_points[:] = \ + self.petab_problem.sample_parameter_startpoints( + num_starting_points) + + # Write nominal values for testing purposes if 'nominalValue' in self.parameter_df: self.f['/parameters/nominalValues'] = \ self.parameter_df.nominalValue[ diff --git a/misc/venv.sh b/misc/venv.sh index 8231ba5e2..fca421674 100755 --- a/misc/venv.sh +++ b/misc/venv.sh @@ -2,7 +2,6 @@ # # Setup virtual environment for building/testing parPE # -set -e SCRIPT_PATH=$(dirname $BASH_SOURCE) PARPE_ROOT=$(cd ${SCRIPT_PATH}/.. && pwd) @@ -15,8 +14,20 @@ if [[ -z "${BUILD_DIR}" ]]; then BUILD_DIR="${PARPE_ROOT}/build"; fi # NOTE: Must remove folder if AMICI is updated if [[ ! -d ${BUILD_DIR}/venv ]]; then # create venv - python3 -m venv ${BUILD_DIR}/venv - source ${BUILD_DIR}/venv/bin/activate + python3 -m venv ${BUILD_DIR}/venv 2>/dev/null + # in case this fails (usually due to missing ensurepip), try getting pip + # manually + if [[ $? ]]; then + set -e + python3 -m venv ${BUILD_DIR}/venv --clear --without-pip + source ${BUILD_DIR}/venv/bin/activate + curl https://bootstrap.pypa.io/get-pip.py -o ${BUILD_DIR}/get-pip.py + python3 ${BUILD_DIR}/get-pip.py + else + set -e + source ${BUILD_DIR}/venv/bin/activate + fi + pip3 install wheel # install AMICI @@ -26,5 +37,6 @@ if [[ ! -d ${BUILD_DIR}/venv ]]; then # install parPE cd ${PARPE_ROOT}/python pip3 install -e . - pip3 install -U git+https://github.com/ICB-DCM/PEtab.git@develop + #pip3 install https://github.com/ICB-DCM/PEtab/archive/develop.zip + pip3 install -U petab==0.0.0a17 fi diff --git a/python/setup.py b/python/setup.py index 3659b87c1..315852137 100644 --- a/python/setup.py +++ b/python/setup.py @@ -11,8 +11,8 @@ install_requires=['numpy', 'termcolor', 'colorama', - 'petab>=0.0.0a14', - 'amici>=0.10.13', + 'petab>=0.0.0a17', + 'amici>=0.10.16', 'h5py', 'python-libsbml>=5.17.0', 'jinja2', diff --git a/shippable.yml b/shippable.yml index af566019f..e2fb04ce6 100644 --- a/shippable.yml +++ b/shippable.yml @@ -44,7 +44,7 @@ build: - CTEST_OUTPUT_ON_FAILURE=1 make ExperimentalMemCheck; cat Testing/Temporary/MemoryChecker.*.log # coverage report - - cd $PARPE_BASE/build && make parpe_coverage_cobertura + - cd $PARPE_BASE/build && CTEST_OUTPUT_ON_FAILURE=1 make parpe_coverage_cobertura - mkdir -p $PARPE_BASE/shippable/codecoverage && cp $PARPE_BASE/build/parpe_coverage_cobertura.xml $PARPE_BASE/shippable/codecoverage # on_failure: diff --git a/src/parpeamici/amiciMisc.cpp b/src/parpeamici/amiciMisc.cpp index a4aa1bf7d..d970a5ff3 100644 --- a/src/parpeamici/amiciMisc.cpp +++ b/src/parpeamici/amiciMisc.cpp @@ -1,33 +1,6 @@ #include -#include - -#include - namespace parpe { -void printAmiciErrMsgIdAndTxt(const char *identifier, const char *format, ...) { - std::stringstream ss; - if(identifier) { - ss <<"["< -#include -#include #include #include + +#include +#include + +#include #include #include -#include +#include #include @@ -274,18 +277,21 @@ void saveSimulation(const H5::H5File &file, std::string const& pathStr, file.getId(), fullGroupPath, "simulationWallTimeInSec", gsl::make_span(&timeElapsedInSeconds, 1)); - if (!states.empty()) - hdf5CreateOrExtendAndWriteToDouble2DArray( - file.getId(), fullGroupPath, "simulationStates", states); + // TODO: This was broken by allowing different numbers of timepoints + // for different simulation conditions. Vector lengths now may differ and + // lead to crashes. + // if (!states.empty()) + // hdf5CreateOrExtendAndWriteToDouble2DArray( + // file.getId(), fullGroupPath, "simulationStates", states); - if (!outputs.empty()) - hdf5CreateOrExtendAndWriteToDouble2DArray( - file.getId(), fullGroupPath, "simulationObservables", outputs); + // if (!outputs.empty()) + // hdf5CreateOrExtendAndWriteToDouble2DArray( + // file.getId(), fullGroupPath, "simulationObservables", outputs); - if (!stateSensi.empty()) - hdf5CreateOrExtendAndWriteToDouble3DArray( - file.getId(), fullGroupPath, "simulationStateSensitivities", stateSensi, - stateSensi.size() / parameters.size(), parameters.size()); + // if (!stateSensi.empty()) + // hdf5CreateOrExtendAndWriteToDouble3DArray( + // file.getId(), fullGroupPath, "simulationStateSensitivities", stateSensi, + // stateSensi.size() / parameters.size(), parameters.size()); hdf5CreateOrExtendAndWriteToInt2DArray( file.getId(), fullGroupPath, "simulationStatus", @@ -322,10 +328,34 @@ AmiciSimulationRunner::AmiciResultPackageSimple runAndLogSimulation( constexpr double errorRelaxation = 1e2; std::unique_ptr rdata; + // redirect AMICI output to parPE logging + amici::AmiciApplication amiciApp; + amiciApp.error = [logger]( + std::string const& identifier, + std::string const& message){ + if(!identifier.empty()) { + logger->logmessage(LOGLVL_ERROR, "[" + identifier + "] " + message); + } else { + logger->logmessage(LOGLVL_ERROR, message); + } + }; + amiciApp.warning = [logger]( + std::string const& identifier, + std::string const& message){ + if(!identifier.empty()) { + logger->logmessage(LOGLVL_WARNING, + "[" + identifier + "] " + message); + } else { + logger->logmessage(LOGLVL_WARNING, message); + } + }; + model.app = &amiciApp; // TODO: may dangle need to unset on exit + for(int trial = 1; trial <= maxNumTrials; ++trial) { /* It is currently not safe to reuse solver if an exception has * occurred,so clone every time */ auto solver = std::unique_ptr(solverTemplate.clone()); + solver->app = &amiciApp; if(trial - 1 == maxNumTrials) { logger->logmessage(LOGLVL_ERROR, @@ -383,7 +413,7 @@ AmiciSimulationRunner::AmiciResultPackageSimple runAndLogSimulation( } try { - rdata = amici::runAmiciSimulation(*solver, edata.get(), model); + rdata = amiciApp.runAmiciSimulation(*solver, edata.get(), model); } catch (std::exception const& e) { logger->logmessage( LOGLVL_WARNING, "Error during simulation: %s (%d)", @@ -439,7 +469,7 @@ FunctionEvaluationStatus getModelOutputs( modelOutput.resize(dataIndices.size()); auto parameterVector = std::vector(parameters.begin(), parameters.end()); - auto jobFinished = [&](JobData *job, int dataIdx) { // jobFinished + auto jobFinished = [&](JobData *job, int /*dataIdx*/) { // jobFinished // deserialize auto results = amici::deserializeFromChar ( diff --git a/src/parpeamici/optimizationApplication.cpp b/src/parpeamici/optimizationApplication.cpp index 85a2ebd78..2fb72998c 100644 --- a/src/parpeamici/optimizationApplication.cpp +++ b/src/parpeamici/optimizationApplication.cpp @@ -59,10 +59,6 @@ int OptimizationApplication::init(int argc, char **argv) { logmessage(LOGLVL_DEBUG, "Seeding RNG with %u", seed); srand(seed); // TODO to CLI - // redirect AMICI output to parPE logging - amici::errMsgIdAndTxt = printAmiciErrMsgIdAndTxt; - amici::warnMsgIdAndTxt = printAmiciWarnMsgIdAndTxt; - return status; } diff --git a/src/parpeamici/standaloneSimulator.cpp b/src/parpeamici/standaloneSimulator.cpp index b0dd3a70b..d278edb9c 100644 --- a/src/parpeamici/standaloneSimulator.cpp +++ b/src/parpeamici/standaloneSimulator.cpp @@ -272,7 +272,33 @@ AmiciSimulationRunner::AmiciResultPackageSimple StandaloneSimulator::runSimulati // currently requires edata, since all condition specific parameters are set via edata auto edata = dataProvider->getExperimentalDataForCondition(conditionIdx); - auto rdata = amici::runAmiciSimulation(solver, edata.get(), model); + // redirect AMICI output to parPE logging + Logger logger("c" + std::to_string(conditionIdx)); + amici::AmiciApplication amiciApp; + amiciApp.error = [&logger]( + std::string const& identifier, + std::string const& message){ + if(!identifier.empty()) { + logger.logmessage(LOGLVL_ERROR, "[" + identifier + "] " + message); + } else { + logger.logmessage(LOGLVL_ERROR, message); + } + }; + amiciApp.warning = [&logger]( + std::string const& identifier, + std::string const& message){ + if(!identifier.empty()) { + logger.logmessage(LOGLVL_WARNING, + "[" + identifier + "] " + message); + } else { + logger.logmessage(LOGLVL_WARNING, message); + } + }; + model.app = &amiciApp; // TODO: may dangle need to unset on exit + solver.app = &amiciApp; // TODO: may dangle need to unset on exit + + + auto rdata = amiciApp.runAmiciSimulation(solver, edata.get(), model); RELEASE_ASSERT(rdata != nullptr, ""); diff --git a/src/parpeloadbalancer/loadBalancerMaster.cpp b/src/parpeloadbalancer/loadBalancerMaster.cpp index 282789188..5602ca27d 100644 --- a/src/parpeloadbalancer/loadBalancerMaster.cpp +++ b/src/parpeloadbalancer/loadBalancerMaster.cpp @@ -140,7 +140,7 @@ int LoadBalancerMaster::handleFinishedJobs() { // loadBalancer.recvRequests pthread_testcancel(); - int messageWaiting; + int messageWaiting = 0; MPI_Iprobe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &messageWaiting, &status); diff --git a/src/parpeoptimization/minibatchOptimization.cpp b/src/parpeoptimization/minibatchOptimization.cpp index acbe76715..ee9101b00 100755 --- a/src/parpeoptimization/minibatchOptimization.cpp +++ b/src/parpeoptimization/minibatchOptimization.cpp @@ -22,7 +22,7 @@ std::vector getVectorDifference(gsl::span v, std::vector difference(v.size(), 0.0); for (unsigned int i = 0; i < v.size(); ++i) difference[i] = v[i] - w[i]; - + return difference; } @@ -57,6 +57,9 @@ void setMinibatchOption(const std::pair &p } else if (val == "Adam" && !dynamic_cast(optimizer->parameterUpdater.get())) { // this might have been set previously if there was an updater-specific option before optimizer->parameterUpdater = std::make_unique(); + } else if (val == "AdamClassic" && !dynamic_cast(optimizer->parameterUpdater.get())) { + // this might have been set previously if there was an updater-specific option before + optimizer->parameterUpdater = std::make_unique(); } else { logmessage(LOGLVL_WARNING, "Ignoring unknown Minibatch parameterUpdater %s.", val.c_str()); } @@ -162,6 +165,8 @@ void LearningRateUpdater::setEndLearningRate(double learningRate) void ParameterUpdaterRmsProp::initialize(unsigned int numParameters) { gradientNormCache.resize(numParameters); oldGradientNormCache.resize(numParameters); + std::fill(gradientNormCache.begin(), gradientNormCache.end(), 0.0); + std::fill(oldGradientNormCache.begin(), oldGradientNormCache.end(), 0.0); } void ParameterUpdaterRmsProp::updateParameters(double learningRate, @@ -195,11 +200,54 @@ void ParameterUpdaterRmsProp::clearCache() { std::fill(oldGradientNormCache.begin(), oldGradientNormCache.end(), 0.0); } +void ParameterUpdaterMomentum::initialize(unsigned int numParameters) { + momentum.resize(numParameters); + oldMomentum.resize(numParameters); + std::fill(momentum.begin(), momentum.end(), 0.0); + std::fill(oldMomentum.begin(), oldMomentum.end(), 0.0); +} + +void ParameterUpdaterMomentum::updateParameters(double learningRate, + int /*iteration*/, + gsl::span gradient, + gsl::span parameters, + gsl::span lowerBounds, + gsl::span upperBounds) { + + int numParameters = gradient.size(); + oldMomentum = momentum; + + for (int i = 0; i < numParameters; ++i) { + momentum[i] = decayRate * momentum[i] + (1 - decayRate) * gradient[i]; + + double momentumNormalize = std::max(getVectorNorm(momentum), 1.0); + + parameters[i] += -learningRate * momentum[i] / momentumNormalize; + } + + clipToBounds(lowerBounds, upperBounds, parameters); +} + +void ParameterUpdaterMomentum::undoLastStep() { + // The cached gradient norm needs to be restored, since the new one is probably NaN + momentum = oldMomentum; +} + +void ParameterUpdaterMomentum::clearCache() { + // Reset all cached information + std::fill(momentum.begin(), momentum.end(), 0.0); +} + + void ParameterUpdaterAdam::initialize(unsigned int numParameters) { gradientCache.resize(numParameters); gradientNormCache.resize(numParameters); oldGradientCache.resize(numParameters); oldGradientNormCache.resize(numParameters); + std::fill(gradientNormCache.begin(), gradientNormCache.end(), 0.0); + std::fill(oldGradientNormCache.begin(), oldGradientNormCache.end(), 0.0); + std::fill(gradientCache.begin(), gradientCache.end(), 0.0); + std::fill(oldGradientCache.begin(), oldGradientCache.end(), 0.0); } void ParameterUpdaterAdam::updateParameters(double learningRate, @@ -215,8 +263,8 @@ void ParameterUpdaterAdam::updateParameters(double learningRate, oldGradientNormCache = gradientNormCache; oldGradientCache = gradientCache; - - for (int i = 0; i < numParameters; ++i) { + + for (int i = 0; i < numParameters; ++i) { // compute new steps from last gradient information gradientCache[i] = decayRateGradient * gradientCache[i] + (1 - decayRateGradient) * gradient[i]; gradientNormCache[i] = decayRateGradientNorm * gradientNormCache[i] @@ -246,6 +294,62 @@ void ParameterUpdaterAdam::clearCache() { std::fill(oldGradientCache.begin(), oldGradientCache.end(), 0.0); } +void ParameterUpdaterAdamClassic::initialize(unsigned int numParameters) { + gradientCache.resize(numParameters); + gradientNormCache.resize(numParameters); + oldGradientCache.resize(numParameters); + oldGradientNormCache.resize(numParameters); + std::fill(gradientNormCache.begin(), gradientNormCache.end(), 0.0); + std::fill(oldGradientNormCache.begin(), oldGradientNormCache.end(), 0.0); + std::fill(gradientCache.begin(), gradientCache.end(), 0.0); + std::fill(oldGradientCache.begin(), oldGradientCache.end(), 0.0); +} + +void ParameterUpdaterAdamClassic::updateParameters(double learningRate, + int iteration, + gsl::span gradient, + gsl::span parameters, + gsl::span lowerBounds, + gsl::span upperBounds) { + + int numParameters = gradient.size(); + double tmpNumerator; + double tmpDenominator; + + oldGradientNormCache = gradientNormCache; + oldGradientCache = gradientCache; + + for (int i = 0; i < numParameters; ++i) { + // compute new steps from last gradient information + gradientCache[i] = decayRateGradient * gradientCache[i] + (1 - decayRateGradient) * gradient[i]; + gradientNormCache[i] = decayRateGradientNorm * gradientNormCache[i] + + (1 - decayRateGradientNorm) * gradient[i] * gradient[i]; + + tmpNumerator = gradientCache[i] / (1 - std::pow(decayRateGradient, (double) iteration)); + tmpDenominator = std::sqrt(gradientNormCache[i] / (1 - std::pow(decayRateGradientNorm, (double) iteration))) + + delta; + + parameters[i] += -learningRate * tmpNumerator / tmpDenominator; + } + + clipToBounds(lowerBounds, upperBounds, parameters); +} + +void ParameterUpdaterAdamClassic::undoLastStep() { + // The cached gradient norm needs to be restored, since the new one is probably NaN + gradientNormCache = oldGradientNormCache; + gradientCache = oldGradientCache; +} + +void ParameterUpdaterAdamClassic::clearCache() { + // Reset all cached information + std::fill(gradientNormCache.begin(), gradientNormCache.end(), 0.0); + std::fill(oldGradientNormCache.begin(), oldGradientNormCache.end(), 0.0); + std::fill(gradientCache.begin(), gradientCache.end(), 0.0); + std::fill(oldGradientCache.begin(), oldGradientCache.end(), 0.0); +} + + void ParameterUpdaterVanilla::updateParameters(double learningRate, int /*iteration*/, gsl::span gradient, diff --git a/src/parpeoptimization/optimizationProblem.cpp b/src/parpeoptimization/optimizationProblem.cpp index 3724c9a24..9d6dc1406 100644 --- a/src/parpeoptimization/optimizationProblem.cpp +++ b/src/parpeoptimization/optimizationProblem.cpp @@ -1,7 +1,8 @@ +#include + #include #include #include -#include #include #include #include