diff --git a/.github/workflows/test_benchmark_collection_models.yml b/.github/workflows/test_benchmark_collection_models.yml index 0b14d73716..3e46af8992 100644 --- a/.github/workflows/test_benchmark_collection_models.yml +++ b/.github/workflows/test_benchmark_collection_models.yml @@ -52,25 +52,28 @@ jobs: AMICI_PARALLEL_COMPILE="" pip3 install -v --user \ $(ls -t python/sdist/dist/amici-*.tar.gz | head -1)[petab,test,vis] - - run: | + - name: Install test dependencies + run: | python3 -m pip uninstall -y petab && python3 -m pip install git+https://github.com/petab-dev/libpetab-python.git@develop \ - && python3 -m pip install -U sympy + && python3 -m pip install -U sympy \ + && python3 -m pip install git+https://github.com/ICB-DCM/fiddy.git - # retrieve test models - - name: Download and test benchmark collection + - name: Download benchmark collection run: | git clone --depth 1 https://github.com/benchmarking-initiative/Benchmark-Models-PEtab.git \ - && export BENCHMARK_COLLECTION="$(pwd)/Benchmark-Models-PEtab/Benchmark-Models/" \ - && pip3 install -e $BENCHMARK_COLLECTION/../src/python \ - && AMICI_PARALLEL_COMPILE="" tests/benchmark-models/test_benchmark_collection.sh + && python3 -m pip install -e Benchmark-Models-PEtab/src/python + + - name: Run tests + env: + AMICI_PARALLEL_COMPILE: "" + run: | + cd tests/benchmark-models && pytest --durations=10 - # run gradient checks - - name: Run Gradient Checks + # collect & upload results + - name: Aggregate results run: | - pip install git+https://github.com/ICB-DCM/fiddy.git \ - && cd tests/benchmark-models && pytest --durations=10 ./test_petab_benchmark.py + cd tests/benchmark-models && python3 evaluate_benchmark.py - # upload results - uses: actions/upload-artifact@v4 with: name: computation-times-${{ matrix.python-version }}-${{ matrix.extract_subexpressions }} diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index 23337986db..fcb1067b85 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -247,6 +247,8 @@ jobs: macos_cpp_py: name: Tests MacOS C++/Python runs-on: macos-latest + env: + PYTHONFAULTHANDLER: "1" steps: - name: Set up Python @@ -291,7 +293,8 @@ jobs: - name: Python tests run: | - scripts/run-python-tests.sh \ + # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved + scripts/run-python-tests.sh -W ignore:: \ test_pregenerated_models.py \ test_splines_short.py \ test_misc.py @@ -300,6 +303,8 @@ jobs: macos_python: name: Tests MacOS Python runs-on: macos-latest + env: + PYTHONFAULTHANDLER: "1" steps: - name: Cache @@ -351,7 +356,8 @@ jobs: - name: Python tests run: | - scripts/run-python-tests.sh \ + # ignore warnings until https://github.com/swig/swig/issues/3061 is resolved + scripts/run-python-tests.sh -W ignore:: \ --ignore=test_pregenerated_models.py \ --ignore=test_splines_short.py \ --ignore=test_misc.py diff --git a/.github/workflows/test_python_ver_matrix.yml b/.github/workflows/test_python_ver_matrix.yml index 8ccf2c1558..a7d2d0e514 100644 --- a/.github/workflows/test_python_ver_matrix.yml +++ b/.github/workflows/test_python_ver_matrix.yml @@ -33,7 +33,7 @@ jobs: - python-version: '3.12' experimental: false - python-version: '3.13' - experimental: true + experimental: false steps: - run: echo "AMICI_DIR=$(pwd)" >> $GITHUB_ENV diff --git a/CHANGELOG.md b/CHANGELOG.md index 8eb0abd8ae..d4c667b7c7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,63 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni ## v0.X Series +### v0.28.0 (2024-11-11) + +**Breaking changes** + +* Changed the default steady-state method to `integrationOnly` + (by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2574) + + The default mode for computing steady states and sensitivities at steady + state was changed to `integrationOnly` + (from previously `integrateIfNewtonFails`). + + This was done for a more robust default behavior. For example, the evaluation + in https://doi.org/10.1371/journal.pone.0312148 shows that - at least for + some models - Newton's method may easily lead to physically impossible + solutions. + + To keep the previous behavior, use: + ```python + amici_model.setSteadyStateComputationMode(amici.SteadyStateComputationMode.integrateIfNewtonFails) + amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.integrateIfNewtonFails) + ``` + +**Fixes** + +* PEtab import: **Fixed potentially incorrect sensitivities** with + observable/state-dependent sigmas. + This was fixed for all cases amici can handle, others cases will now result + in `ValueError`s (https://github.com/AMICI-dev/AMICI/pull/2563). + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2562 + +* Fixed potentially incorrect disabling of Newton's method + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2576 + +* Fixed `ModelStateDerived` copy ctor, where previously dangling pointers + could lead to crashes in some situations + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2583 + +* Added missing simulation status codes + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2560 + +* Check for unsupported observable IDs in sigma expressions + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2563 + + +**Features** + +* Optional warning in `fill_in_parameters` + + by @dweindl in https://github.com/AMICI-dev/AMICI/pull/2578 + +**Full Changelog**: https://github.com/AMICI-dev/AMICI/compare/v0.27.0...v0.28.0 + ### v0.27.0 (2024-10-21) This release comes with an **updated version of the SUNDIALS package (7.1.1)** (https://github.com/AMICI-dev/AMICI/pull/2513). diff --git a/codecov.yml b/codecov.yml index 15beea9dfc..2779b486d1 100644 --- a/codecov.yml +++ b/codecov.yml @@ -1,8 +1,10 @@ -fixes: - - "build/venv/lib/python3.9/site-packages/::python/" - - "build/venv/lib/python3.10/site-packages/::python/" - - "build/venv/lib/python3.11/site-packages/::python/" +# see https://docs.codecov.com/docs/codecovyml-reference +fixes: + # https://docs.codecov.com/docs/fixing-paths + - "build/venv/lib/python[0-9.]+/site-packages/::python/" + - "python/sdist/build/temp_amici/CMakeFiles/amici.dir/src/::src/" + - "build/CMakeFiles/amici.dir/src/::src/" codecov: require_ci_to_pass: yes @@ -27,6 +29,8 @@ comment: ignore: - "tests/*" - "tests/**/*" + - "build/tests/**" + - "amici_models/**" flags: python: diff --git a/include/amici/defines.h b/include/amici/defines.h index 1c6741f8d2..b49ebc6a83 100644 --- a/include/amici/defines.h +++ b/include/amici/defines.h @@ -74,7 +74,13 @@ constexpr int AMICI_CONSTR_FAIL= -15; constexpr int AMICI_CVODES_CONSTR_FAIL= -15; constexpr int AMICI_IDAS_CONSTR_FAIL= -11; constexpr int AMICI_ILL_INPUT= -22; +constexpr int AMICI_BAD_T= -25; +constexpr int AMICI_BAD_DKY= -26; constexpr int AMICI_FIRST_QRHSFUNC_ERR= -32; +constexpr int AMICI_SRHSFUNC_FAIL= -41; +constexpr int AMICI_FIRST_SRHSFUNC_ERR= -42; +constexpr int AMICI_REPTD_SRHSFUNC_ERR= -43; +constexpr int AMICI_UNREC_SRHSFUNC_ERR= -44; constexpr int AMICI_ERROR= -99; constexpr int AMICI_NO_STEADY_STATE= -81; constexpr int AMICI_DAMPING_FACTOR_ERROR= -86; diff --git a/include/amici/model.h b/include/amici/model.h index 533f0139dc..57591dd2f7 100644 --- a/include/amici/model.h +++ b/include/amici/model.h @@ -2063,12 +2063,12 @@ class Model : public AbstractModel, public ModelDimensions { /** method for steady-state computation */ SteadyStateComputationMode steadystate_computation_mode_{ - SteadyStateComputationMode::integrateIfNewtonFails + SteadyStateComputationMode::integrationOnly }; /** method for steadystate sensitivities computation */ SteadyStateSensitivityMode steadystate_sensitivity_mode_{ - SteadyStateSensitivityMode::integrateIfNewtonFails + SteadyStateSensitivityMode::integrationOnly }; /** diff --git a/include/amici/model_state.h b/include/amici/model_state.h index 149c5757b8..defb12d4c0 100644 --- a/include/amici/model_state.h +++ b/include/amici/model_state.h @@ -146,71 +146,37 @@ struct ModelStateDerived { , dwdw_(other.dwdw_) , dwdx_hierarchical_(other.dwdx_hierarchical_) , dJydy_dense_(other.dJydy_dense_) { - // Update the SUNContext of all matrices - if (J_.data()) { - J_.get()->sunctx = sunctx_; + // Update the SUNContext of all SUNDIALS objects + J_.set_ctx(sunctx_); + JB_.set_ctx(sunctx_); + dxdotdw_.set_ctx(sunctx_); + dwdx_.set_ctx(sunctx_); + dwdp_.set_ctx(sunctx_); + M_.set_ctx(sunctx_); + MSparse_.set_ctx(sunctx_); + dfdx_.set_ctx(sunctx_); + dxdotdp_full.set_ctx(sunctx_); + dxdotdp_explicit.set_ctx(sunctx_); + dxdotdp_implicit.set_ctx(sunctx_); + dxdotdx_explicit.set_ctx(sunctx_); + dxdotdx_implicit.set_ctx(sunctx_); + dx_rdatadx_solver.set_ctx(sunctx_); + dx_rdatadtcl.set_ctx(sunctx_); + dtotal_cldx_rdata.set_ctx(sunctx_); + dxdotdp.set_ctx(sunctx_); + + for (auto& dJydy : dJydy_) { + dJydy.set_ctx(sunctx_); } - if (JB_.data()) { - JB_.get()->sunctx = sunctx_; + for (auto& dwdp : dwdp_hierarchical_) { + dwdp.set_ctx(sunctx_); } - if (dxdotdw_.data()) { - dxdotdw_.get()->sunctx = sunctx_; - } - if (dwdx_.data()) { - dwdx_.get()->sunctx = sunctx_; - } - if (dwdp_.data()) { - dwdp_.get()->sunctx = sunctx_; - } - if (M_.data()) { - M_.get()->sunctx = sunctx_; - } - if (MSparse_.data()) { - MSparse_.get()->sunctx = sunctx_; - } - if (dfdx_.data()) { - dfdx_.get()->sunctx = sunctx_; - } - if (dxdotdp_full.data()) { - dxdotdp_full.get()->sunctx = sunctx_; - } - if (dxdotdp_explicit.data()) { - dxdotdp_explicit.get()->sunctx = sunctx_; - } - if (dxdotdp_implicit.data()) { - dxdotdp_implicit.get()->sunctx = sunctx_; - } - if (dxdotdx_explicit.data()) { - dxdotdx_explicit.get()->sunctx = sunctx_; - } - if (dxdotdx_implicit.data()) { - dxdotdx_implicit.get()->sunctx = sunctx_; - } - if (dx_rdatadx_solver.data()) { - dx_rdatadx_solver.get()->sunctx = sunctx_; - } - if (dx_rdatadtcl.data()) { - dx_rdatadtcl.get()->sunctx = sunctx_; - } - if (dtotal_cldx_rdata.data()) { - dtotal_cldx_rdata.get()->sunctx = sunctx_; - } - for (auto const& dwdp : dwdp_hierarchical_) { - if (dwdp.data()) { - dwdp.get()->sunctx = sunctx_; - } - } - for (auto const& dwdx : dwdx_hierarchical_) { - if (dwdx.data()) { - dwdx.get()->sunctx = sunctx_; - } - } - if (dwdw_.data()) { - dwdw_.get()->sunctx = sunctx_; - } - if (dJydy_dense_.data()) { - dJydy_dense_.get()->sunctx = sunctx_; + for (auto& dwdx : dwdx_hierarchical_) { + dwdx.set_ctx(sunctx_); } + sspl_.set_ctx(sunctx_); + dwdw_.set_ctx(sunctx_); + dJydy_dense_.set_ctx(sunctx_); } /** diff --git a/include/amici/sundials_matrix_wrapper.h b/include/amici/sundials_matrix_wrapper.h index 942069a803..a7ad361acb 100644 --- a/include/amici/sundials_matrix_wrapper.h +++ b/include/amici/sundials_matrix_wrapper.h @@ -506,6 +506,19 @@ class SUNMatrixWrapper { */ SUNContext get_ctx() const; + /** + * @brief Set SUNContext + * + * Update the SUNContext of the wrapped SUNMatrix. + * + * @param ctx SUNDIALS context to set + */ + void set_ctx(SUNContext ctx) { + if (matrix_) { + matrix_->sunctx = ctx; + } + } + private: /** * @brief SUNMatrix to which all methods are applied diff --git a/include/amici/vector.h b/include/amici/vector.h index 32c436fbda..0a7648e460 100644 --- a/include/amici/vector.h +++ b/include/amici/vector.h @@ -413,6 +413,20 @@ class AmiVectorArray { */ void copy(AmiVectorArray const& other); + /** + * @brief Set SUNContext + * + * If any AmiVector is non-empty, this changes the current SUNContext of the + * associated N_Vector. If empty, do nothing. + * + * @param ctx SUNDIALS context to set + */ + void set_ctx(SUNContext ctx) { + for (auto& vec : vec_array_) { + vec.set_ctx(ctx); + } + } + private: /** main data storage */ std::vector vec_array_; diff --git a/python/sdist/amici/petab/conditions.py b/python/sdist/amici/petab/conditions.py index 08c2f90302..ab06e8850d 100644 --- a/python/sdist/amici/petab/conditions.py +++ b/python/sdist/amici/petab/conditions.py @@ -41,6 +41,7 @@ def fill_in_parameters( scaled_parameters: bool, parameter_mapping: ParameterMapping, amici_model: AmiciModel, + warn_unused: bool = True, ) -> None: """Fill fixed and dynamic parameters into the edatas (in-place). @@ -59,9 +60,15 @@ def fill_in_parameters( Parameter mapping for all conditions. :param amici_model: AMICI model. + :param warn_unused: + Whether a warning should be emitted if not all problem parameters + were used. I.e., if there are parameters in `problem_parameters` + that are not in `parameter_mapping`. """ - if unused_parameters := ( - set(problem_parameters.keys()) - parameter_mapping.free_symbols + if warn_unused and ( + unused_parameters := ( + set(problem_parameters.keys()) - parameter_mapping.free_symbols + ) ): warnings.warn( "The following problem parameters were not used: " @@ -223,6 +230,7 @@ def create_parameterized_edatas( scaled_parameters: bool = False, parameter_mapping: ParameterMapping = None, simulation_conditions: pd.DataFrame | dict = None, + warn_unused: bool = True, ) -> list[amici.ExpData]: """Create list of :class:amici.ExpData objects with parameters filled in. @@ -244,6 +252,11 @@ def create_parameterized_edatas( :param simulation_conditions: Result of :func:`petab.get_simulation_conditions`. Can be provided to save time if this has been obtained before. + :param warn_unused: + Whether a warning should be emitted if not all problem parameters + were used. I.e., if there are parameters in `problem_parameters` + that are not in `parameter_mapping` or in the generated parameter + mapping. :return: List with one :class:`amici.amici.ExpData` per simulation condition, @@ -282,6 +295,7 @@ def create_parameterized_edatas( scaled_parameters=scaled_parameters, parameter_mapping=parameter_mapping, amici_model=amici_model, + warn_unused=warn_unused, ) return edatas @@ -387,7 +401,9 @@ def create_edatas( :return: List with one :class:`amici.amici.ExpData` per simulation condition, - with filled in timepoints and data. + with filled in timepoints and data, but without parameter values + (see :func:`create_parameterized_edatas` or + :func:`fill_in_parameters` for that). """ if simulation_conditions is None: simulation_conditions = ( diff --git a/python/sdist/amici/petab/import_helpers.py b/python/sdist/amici/petab/import_helpers.py index b52d74004d..70af87c3b3 100644 --- a/python/sdist/amici/petab/import_helpers.py +++ b/python/sdist/amici/petab/import_helpers.py @@ -47,8 +47,8 @@ def get_observation_model( observables = {} sigmas = {} - nan_pat = r"^[nN]a[nN]$" + for _, observable in observable_df.iterrows(): oid = str(observable.name) # need to sanitize due to https://github.com/PEtab-dev/PEtab/issues/447 @@ -56,20 +56,18 @@ def get_observation_model( formula_obs = re.sub(nan_pat, "", str(observable[OBSERVABLE_FORMULA])) formula_noise = re.sub(nan_pat, "", str(observable[NOISE_FORMULA])) observables[oid] = {"name": name, "formula": formula_obs} - sigmas[oid] = formula_noise - - # PEtab does currently not allow observables in noiseFormula and AMICI - # cannot handle states in sigma expressions. Therefore, where possible, - # replace species occurring in error model definition by observableIds. - replacements = { - sp.sympify(observable["formula"], locals=_clash): sp.Symbol( - observable_id + formula_noise_sym = sp.sympify(formula_noise, locals=_clash) + formula_obs_sym = sp.sympify(formula_obs, locals=_clash) + # PEtab does currently not allow observables in noiseFormula, and + # AMICI cannot handle state variables in sigma expressions. + # Therefore, where possible, replace state variables occurring in + # the error model definition by equivalent observableIds. + # Do not use observableIds from other rows + # (https://github.com/AMICI-dev/AMICI/issues/2561). + formula_noise_sym = formula_noise_sym.subs( + formula_obs_sym, sp.Symbol(oid) ) - for observable_id, observable in observables.items() - } - for observable_id, formula in sigmas.items(): - repl = sp.sympify(formula, locals=_clash).subs(replacements) - sigmas[observable_id] = str(repl) + sigmas[oid] = str(formula_noise_sym) noise_distrs = petab_noise_distributions_to_amici(observable_df) diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 61ce9a0ee1..fcaa1ed752 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -2052,12 +2052,26 @@ def _process_log_likelihood( obs["reg_symbol"] = generate_regularization_symbol(obs_id) if not event_reg: + sigmas = { + obs_id: self._sympy_from_sbml_math(sigma) + for obs_id, sigma in sigmas.items() + } + obs_syms = set(self.symbols[obs_symbol].keys()) + for obs_id in self.symbols[obs_symbol]: + if (sigma := sigmas.get(str(obs_id))) and sigma.has( + *(obs_syms - {obs_id}) + ): + raise ValueError( + f"Sigma expression for {obs_id} ({sigma}) must not " + f"contain any observable symbols other than {obs_id}." + ) + # check that only the corresponding observable ID is used in the + # sigma formula, but not any other observable ID + # https://github.com/AMICI-dev/AMICI/issues/2561 self.symbols[sigma_symbol] = { symbol_with_assumptions(f"sigma_{obs_id}"): { "name": f'sigma_{obs["name"]}', - "value": self._sympy_from_sbml_math( - sigmas.get(str(obs_id), "1.0") - ), + "value": sigmas.get(str(obs_id), sp.Float(1.0)), } for obs_id, obs in self.symbols[obs_symbol].items() } diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index 640d2dd988..a19a7565d2 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -102,6 +102,7 @@ def get_results( sensi_order=0, sensi_meth=amici.SensitivityMethod.forward, sensi_meth_preeq=amici.SensitivityMethod.forward, + stst_mode=amici.SteadyStateComputationMode.integrateIfNewtonFails, stst_sensi_mode=amici.SteadyStateSensitivityMode.newtonOnly, reinitialize_states=False, ): @@ -115,6 +116,7 @@ def get_results( solver.setSensitivityMethodPreequilibration(sensi_meth_preeq) solver.setSensitivityMethod(sensi_meth) model.setSteadyStateSensitivityMode(stst_sensi_mode) + model.setSteadyStateComputationMode(stst_mode) if edata is None: model.setTimepoints(np.linspace(0, 5, 101)) else: diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index bf81e0256b..c7e254114e 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -18,7 +18,12 @@ def preeq_fixture(pysb_example_presimulation_module): model = pysb_example_presimulation_module.getModel() model.setReinitializeFixedParameterInitialStates(True) - + model.setSteadyStateComputationMode( + amici.SteadyStateComputationMode.integrateIfNewtonFails + ) + model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) solver = model.getSolver() solver.setSensitivityOrder(amici.SensitivityOrder.first) solver.setSensitivityMethod(amici.SensitivityMethod.forward) diff --git a/python/tests/test_pregenerated_models.py b/python/tests/test_pregenerated_models.py index 5a110cdfc2..e9daf13ba8 100755 --- a/python/tests/test_pregenerated_models.py +++ b/python/tests/test_pregenerated_models.py @@ -63,6 +63,13 @@ def test_pregenerated_model(sub_test, case): amici.readModelDataFromHDF5( options_file, model.get(), f"/{sub_test}/{case}/options" ) + if model_name == "model_steadystate": + model.setSteadyStateComputationMode( + amici.SteadyStateComputationMode.integrateIfNewtonFails + ) + model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode.integrateIfNewtonFails + ) amici.readSolverSettingsFromHDF5( options_file, solver.get(), f"/{sub_test}/{case}/options" ) diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py index f61b28ea4d..9686a25d94 100644 --- a/python/tests/test_swig_interface.py +++ b/python/tests/test_swig_interface.py @@ -105,12 +105,12 @@ def test_copy_constructors(pysb_example_presimulation_module): # `pysb_example_presimulation_module.getModel()`. "StateIsNonNegative": None, "SteadyStateComputationMode": [ - 2, - 1, + amici.SteadyStateComputationMode.integrationOnly, + amici.SteadyStateComputationMode.integrateIfNewtonFails, ], "SteadyStateSensitivityMode": [ - 2, - 1, + amici.SteadyStateSensitivityMode.integrationOnly, + amici.SteadyStateSensitivityMode.integrateIfNewtonFails, ], ("t0", "setT0"): [ 0.0, diff --git a/src/amici.cpp b/src/amici.cpp index bba344358f..8312200e0b 100644 --- a/src/amici.cpp +++ b/src/amici.cpp @@ -51,6 +51,12 @@ std::map simulation_status_to_str_map = { {AMICI_LSETUP_FAIL, "AMICI_LSETUP_FAIL"}, {AMICI_FIRST_QRHSFUNC_ERR, "AMICI_FIRST_QRHSFUNC_ERR"}, {AMICI_WARNING, "AMICI_WARNING"}, + {AMICI_BAD_T, "AMICI_BAD_T"}, + {AMICI_BAD_DKY, "AMICI_BAD_DKY"}, + {AMICI_FIRST_SRHSFUNC_ERR, "AMICI_FIRST_SRHSFUNC_ERR"}, + {AMICI_SRHSFUNC_FAIL, "AMICI_SRHSFUNC_FAIL"}, + {AMICI_REPTD_SRHSFUNC_ERR, "AMICI_REPTD_SRHSFUNC_ERR"}, + {AMICI_UNREC_SRHSFUNC_ERR, "AMICI_UNREC_SRHSFUNC_ERR"}, }; std::unique_ptr runAmiciSimulation( diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index f777011329..0b60304c55 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -34,65 +34,32 @@ static_assert((int)InterpolationType::polynomial == CV_POLYNOMIAL, ""); static_assert((int)LinearMultistepMethod::adams == CV_ADAMS, ""); static_assert((int)LinearMultistepMethod::BDF == CV_BDF, ""); -static_assert(AMICI_ROOT_RETURN == CV_ROOT_RETURN, ""); - -static_assert( - amici::AMICI_SUCCESS == CV_SUCCESS, "AMICI_SUCCESS != CV_SUCCESS" -); -static_assert( - amici::AMICI_DATA_RETURN == CV_TSTOP_RETURN, - "AMICI_DATA_RETURN != CV_TSTOP_RETURN" -); -static_assert( - amici::AMICI_ROOT_RETURN == CV_ROOT_RETURN, - "AMICI_ROOT_RETURN != CV_ROOT_RETURN" -); -static_assert( - amici::AMICI_ILL_INPUT == CV_ILL_INPUT, "AMICI_ILL_INPUT != CV_ILL_INPUT" -); -static_assert(amici::AMICI_NORMAL == CV_NORMAL, "AMICI_NORMAL != CV_NORMAL"); -static_assert( - amici::AMICI_ONE_STEP == CV_ONE_STEP, "AMICI_ONE_STEP != CV_ONE_STEP" -); -static_assert( - amici::AMICI_TOO_MUCH_ACC == CV_TOO_MUCH_ACC, - "AMICI_TOO_MUCH_ACC != CV_TOO_MUCH_ACC" -); -static_assert( - amici::AMICI_TOO_MUCH_WORK == CV_TOO_MUCH_WORK, - "AMICI_TOO_MUCH_WORK != CV_TOO_MUCH_WORK" -); -static_assert( - amici::AMICI_ERR_FAILURE == CV_ERR_FAILURE, - "AMICI_ERR_FAILURE != CV_ERR_FAILURE" -); -static_assert( - amici::AMICI_CONV_FAILURE == CV_CONV_FAILURE, - "AMICI_CONV_FAILURE != CV_CONV_FAILURE" -); -static_assert( - amici::AMICI_LSETUP_FAIL == CV_LSETUP_FAIL, - "AMICI_LSETUP_FAIL != CV_LSETUP_FAIL" -); -static_assert( - amici::AMICI_CONSTR_FAIL == CV_CONSTR_FAIL, - "AMICI_CONSTR_FAIL != CV_CONSTR_FAIL" -); -static_assert( - amici::AMICI_RHSFUNC_FAIL == CV_RHSFUNC_FAIL, - "AMICI_RHSFUNC_FAIL != CV_RHSFUNC_FAIL" -); -static_assert( - amici::AMICI_FIRST_RHSFUNC_ERR == CV_FIRST_RHSFUNC_ERR, - "AMICI_FIRST_RHSFUNC_ERR != CV_FIRST_RHSFUNC_ERR" -); -static_assert( - amici::AMICI_FIRST_QRHSFUNC_ERR == CV_FIRST_QRHSFUNC_ERR, - "AMICI_FIRST_QRHSFUNC_ERR != CV_FIRST_QRHSFUNC_ERR" -); -static_assert( - amici::AMICI_WARNING == CV_WARNING, "AMICI_WARNING != CV_WARNING" -); +#define STATIC_ASSERT_EQUAL(amici_constant, cv_constant) \ + static_assert( \ + amici_constant == cv_constant, #amici_constant " != " #cv_constant \ + ) + +STATIC_ASSERT_EQUAL(amici::AMICI_SUCCESS, CV_SUCCESS); +STATIC_ASSERT_EQUAL(amici::AMICI_ROOT_RETURN, CV_ROOT_RETURN); +STATIC_ASSERT_EQUAL(amici::AMICI_DATA_RETURN, CV_TSTOP_RETURN); +STATIC_ASSERT_EQUAL(amici::AMICI_ILL_INPUT, CV_ILL_INPUT); +STATIC_ASSERT_EQUAL(amici::AMICI_NORMAL, CV_NORMAL); +STATIC_ASSERT_EQUAL(amici::AMICI_ONE_STEP, CV_ONE_STEP); +STATIC_ASSERT_EQUAL(amici::AMICI_TOO_MUCH_ACC, CV_TOO_MUCH_ACC); +STATIC_ASSERT_EQUAL(amici::AMICI_TOO_MUCH_WORK, CV_TOO_MUCH_WORK); +STATIC_ASSERT_EQUAL(amici::AMICI_ERR_FAILURE, CV_ERR_FAILURE); +STATIC_ASSERT_EQUAL(amici::AMICI_CONV_FAILURE, CV_CONV_FAILURE); +STATIC_ASSERT_EQUAL(amici::AMICI_LSETUP_FAIL, CV_LSETUP_FAIL); +STATIC_ASSERT_EQUAL(amici::AMICI_CONSTR_FAIL, CV_CONSTR_FAIL); +STATIC_ASSERT_EQUAL(amici::AMICI_RHSFUNC_FAIL, CV_RHSFUNC_FAIL); +STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_RHSFUNC_ERR, CV_FIRST_RHSFUNC_ERR); +STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_QRHSFUNC_ERR, CV_FIRST_QRHSFUNC_ERR); +STATIC_ASSERT_EQUAL(amici::AMICI_BAD_T, CV_BAD_T); +STATIC_ASSERT_EQUAL(amici::AMICI_BAD_DKY, CV_BAD_DKY); +STATIC_ASSERT_EQUAL(amici::AMICI_SRHSFUNC_FAIL, CV_SRHSFUNC_FAIL); +STATIC_ASSERT_EQUAL(amici::AMICI_FIRST_SRHSFUNC_ERR, CV_FIRST_SRHSFUNC_ERR); +STATIC_ASSERT_EQUAL(amici::AMICI_REPTD_SRHSFUNC_ERR, CV_REPTD_SRHSFUNC_ERR); +STATIC_ASSERT_EQUAL(amici::AMICI_UNREC_SRHSFUNC_ERR, CV_UNREC_SRHSFUNC_ERR); /* * The following static members are callback function to CVODES. diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d2923d5bcd..7c3ac9fab3 100644 --- a/src/steadystateproblem.cpp +++ b/src/steadystateproblem.cpp @@ -133,8 +133,9 @@ void SteadystateProblem::findSteadyState( = model.getSteadyStateComputationMode() == SteadyStateComputationMode::integrationOnly || solver.getNewtonMaxSteps() == 0 - || (model.getSteadyStateSensitivityMode() - == SteadyStateSensitivityMode::integrationOnly + || (solver.getSensitivityOrder() >= SensitivityOrder::first + && model.getSteadyStateSensitivityMode() + == SteadyStateSensitivityMode::integrationOnly && ((it == -1 && solver.getSensitivityMethodPreequilibration() == SensitivityMethod::forward) diff --git a/tests/benchmark-models/benchmark_models.yaml b/tests/benchmark-models/benchmark_models.yaml index 8550c9ca38..2fed74cf3f 100644 --- a/tests/benchmark-models/benchmark_models.yaml +++ b/tests/benchmark-models/benchmark_models.yaml @@ -57,7 +57,7 @@ Elowitz_Nature2000: t_adj: 0.11 note: benchmark collection reference value matches up to sign when applying log10-correction +sum(log(meas*log(10))) / 2 -Fiedler_BMC2016: +Fiedler_BMCSystBiol2016: llh: 58.58390161681 t_sim: 0.005 t_fwd: 0.05 diff --git a/tests/benchmark-models/test_benchmark_collection.sh b/tests/benchmark-models/test_benchmark_collection.sh deleted file mode 100755 index 581b8db028..0000000000 --- a/tests/benchmark-models/test_benchmark_collection.sh +++ /dev/null @@ -1,136 +0,0 @@ -#!/bin/bash -# Import and run selected benchmark models with nominal parameters and check -# agreement with reference values -# -# Expects environment variable BENCHMARK_COLLECTION to provide path to -# benchmark collection model directory - -# Confirmed to be working -models=" -Bachmann_MSB2011 -Beer_MolBioSystems2014 -Boehm_JProteomeRes2014 -Borghans_BiophysChem1997 -Brannmark_JBC2010 -Bruno_JExpBot2016 -Crauste_CellSystems2017 -Elowitz_Nature2000 -Fiedler_BMCSystBiol2016 -Fujita_SciSignal2010 -Isensee_JCB2018 -Lucarelli_CellSystems2018 -Schwen_PONE2014 -Smith_BMCSystBiol2013 -Sneyd_PNAS2002 -Weber_BMC2015 -Zheng_PNAS2012" - -# -# not merged: -# Becker_Science2010 (multiple models) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Becker_Science2010 -# Hass_PONE2017 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Hass_PONE2017 -# Korkut_eLIFE2015 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Korkut_eLIFE2015 -# Casaletto_PNAS2019 (yaml missing) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Casaletto_PNAS2019 -# Merkle_PCB2016 (model missing) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Merkle_PCB2016 -# Parmar_PCB2019 (SBML extensions) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Parmar_PCB2019 -# Swameye_PNAS2003 (splines) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Swameye_PNAS2003 -# Sobotta_Frontiers2017 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Sobotta_Frontiers2017 -# Raia_CancerResearch2011 (state dependent sigmas) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Raia_CancerResearch2011 -# -# no reference value: -# Alkan_SciSignal2018 (d2d: Alkan_DDRP_SciSignal2018) -# Bertozzi_PNAS2020 (gh: vanako, doi: https://doi.org/10.1073/pnas.2006520117, code/data: https://github.com/gomohler/pnas2020 (R)) -# Blasi_CellSystems2016 (gh: Leonard Schmiester, doi: https://doi.org/10.1016/j.cels.2016.01.002, code/data: not available) -# Giordano_Nature2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1038/s41591-020-0883-7, code/data: http://users.dimi.uniud.it/~giulia.giordano/docs/papers/SIDARTHEcode.zip (MATLAB)) -# Laske_PLOSComputBiol2019 (gh: Clemens Peiter, doi: https://doi.org/10.1128/JVI.00080-12 (?), code/data: ???) -# Okuonghae_ChaosSolitonsFractals2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1016/j.chaos.2020.110032, code/data: ???) -# Oliveira_NatCommun2021 (gh: lorenapohl, doi: https://doi.org/10.1038/s41467-020-19798-3, code: https://github.com/cidacslab/Mathematical-and-Statistical-Modeling-of-COVID19-in-Brazil (python) data: https://infovis.sei.ba.gov.br/covid19/ ) -# Perelson_Science1996 (gh: Philipp Staedter, doi: https://doi.org/10.1126/science.271.5255.1582, code/data: ???) -# Rahman_MBS2016 (gh: Yannik Schaelte, doi: https://doi.org/10.1016/j.mbs.2016.07.009, code: not available, data: table in paper ...) -# Raimundez_PCB2020 (gh: Elba Raimundez, doi: https://doi.org/10.1371/journal.pcbi.1007147, code/data: https://zenodo.org/record/2908234#.Y5hUUS8w3yw (d2d)) -# SalazarCavazos_MBoC2020 (gh: Dilan Pathirana, doi: https://doi.org/10.1091/mbc.E19-09-0548, code/data: supplement (BNGL)) -# Zhao_QuantBiol2020 (gh: Iva Ewert, doi: https://doi.org/10.1007/s40484-020-0199-0, code: not available, data: table in supp) -# -# covered by performance test: -# Froehlich_CellSystems2018 -# -# Unknown reasons: -# Chen_MSB2009 -# - -set -e - -[[ -n "${BENCHMARK_COLLECTION}" ]] && model_dir="${BENCHMARK_COLLECTION}" - -function show_help() { - echo "-h: this help; -n: dry run, print commands; -b path_to_models_dir" -} - -OPTIND=1 -while getopts "h?nb:" opt; do - case "$opt" in - h | \?) - show_help - exit 0 - ;; - n) - dry_run=1 - ;; - b) - model_dir=$OPTARG - ;; - esac -done - -script_path=$(dirname "$BASH_SOURCE") -script_path=$(cd "$script_path" && pwd) - -for model in $models; do - yaml="${model_dir}"/"${model}"/"${model}".yaml - - # different naming scheme - if [[ "$model" == "Bertozzi_PNAS2020" ]]; then - yaml="${model_dir}"/"${model}"/problem.yaml - fi - - amici_model_dir=test_bmc/"${model}" - mkdir -p "$amici_model_dir" - cmd_import="amici_import_petab ${yaml} -o ${amici_model_dir} -n ${model} --flatten" - cmd_run="$script_path/test_petab_model.py -y ${yaml} -d ${amici_model_dir} -m ${model} -c" - - printf '=%.0s' {1..40} - printf " %s " "${model}" - printf '=%.0s' {1..40} - echo - - if [[ -z "$dry_run" ]]; then - $cmd_import - $cmd_run - else - echo "$cmd_import" - echo "$cmd_run" - fi - - printf '=%.0s' {1..100} - echo - echo -done - -cd "$script_path" && python evaluate_benchmark.py - -# Test deprecated import from individual PEtab files -model="Zheng_PNAS2012" -problem_dir="${model_dir}/${model}" -amici_model_dir=test_bmc/"${model}-deprecated" -cmd_import="amici_import_petab -s "${problem_dir}/model_${model}.xml" \ - -m "${problem_dir}/measurementData_${model}.tsv" \ - -c "${problem_dir}/experimentalCondition_${model}.tsv" \ - -p "${problem_dir}/parameters_${model}.tsv" \ - -b "${problem_dir}/observables_${model}.tsv" \ - -o ${amici_model_dir} -n ${model} --no-compile" - -if [[ -z "$dry_run" ]]; then - $cmd_import -else - echo "$cmd_import" -fi diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 4892100877..b055f8961f 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -1,9 +1,12 @@ -"""Tests for simulate_petab on PEtab benchmark problems.""" +"""Tests based on the PEtab benchmark problems. -from pathlib import Path +Tests simulate_petab, correctness of the log-likelihood computation at nominal +parameters, correctness of the gradient computation, and simulation times +for a subset of the benchmark problems. +""" +from pathlib import Path import fiddy - import amici import numpy as np import pandas as pd @@ -12,19 +15,48 @@ from amici.petab.petab_import import import_petab_problem import benchmark_models_petab from collections import defaultdict -from dataclasses import dataclass +from dataclasses import dataclass, field from amici import SensitivityMethod +from petab.v1.lint import measurement_table_has_timepoint_specific_mappings from fiddy import MethodId, get_derivative from fiddy.derivative_check import NumpyIsCloseDerivativeCheck from fiddy.extensions.amici import simulate_petab_to_cached_functions from fiddy.success import Consistency +import contextlib +import logging +import yaml +from amici.logging import get_logger +from amici.petab.simulations import ( + LLH, + RDATAS, + rdatas_to_measurement_df, + simulate_petab, +) +from petab.v1.visualize import plot_problem + +# Enable various debug output +debug = False -repo_root = Path(__file__).parent.parent.parent +logger = get_logger( + f"amici.{__name__}", logging.DEBUG if debug else logging.INFO +) -# reuse compiled models from test_benchmark_collection.sh +script_dir = Path(__file__).parent.absolute() +repo_root = script_dir.parent.parent benchmark_outdir = repo_root / "test_bmc" -models = set(benchmark_models_petab.MODELS) - { +debug_path = script_dir / "debug" +if debug: + debug_path.mkdir(exist_ok=True, parents=True) + +# reference values for simulation times and log-likelihoods +references_yaml = script_dir / "benchmark_models.yaml" +with open(references_yaml) as f: + reference_values = yaml.full_load(f) + +# problem IDs for which to check the gradient +# TODO: extend +problems_for_gradient_check = set(benchmark_models_petab.MODELS) - { # excluded due to excessive runtime "Bachmann_MSB2011", "Chen_MSB2009", @@ -38,15 +70,75 @@ "Smith_BMCSystBiol2013", # excluded due to excessive numerical failures "Crauste_CellSystems2017", - "Fujita_SciSignal2010", - # FIXME: re-enable - "Raia_CancerResearch2011", } -models = list(sorted(models)) +problems_for_gradient_check = list(sorted(problems_for_gradient_check)) + + +# Problems for checking the log-likelihood computation at nominal parameters +# +# not merged: +# Becker_Science2010 (multiple models) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Becker_Science2010 +# Hass_PONE2017 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Hass_PONE2017 +# Korkut_eLIFE2015 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Korkut_eLIFE2015 +# Casaletto_PNAS2019 (yaml missing) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Casaletto_PNAS2019 +# Merkle_PCB2016 (model missing) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Merkle_PCB2016 +# Parmar_PCB2019 (SBML extensions) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Parmar_PCB2019 +# Swameye_PNAS2003 (splines) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Swameye_PNAS2003 +# Sobotta_Frontiers2017 (???) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Sobotta_Frontiers2017 +# Raia_CancerResearch2011 (state dependent sigmas) https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/model_Raia_CancerResearch2011 +# +# no reference value: +# Alkan_SciSignal2018 (d2d: Alkan_DDRP_SciSignal2018) +# Bertozzi_PNAS2020 (gh: vanako, doi: https://doi.org/10.1073/pnas.2006520117, code/data: https://github.com/gomohler/pnas2020 (R)) +# Blasi_CellSystems2016 (gh: Leonard Schmiester, doi: https://doi.org/10.1016/j.cels.2016.01.002, code/data: not available) +# Giordano_Nature2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1038/s41591-020-0883-7, code/data: http://users.dimi.uniud.it/~giulia.giordano/docs/papers/SIDARTHEcode.zip (MATLAB)) +# Laske_PLOSComputBiol2019 (gh: Clemens Peiter, doi: https://doi.org/10.1128/JVI.00080-12 (?), code/data: ???) +# Okuonghae_ChaosSolitonsFractals2020 (gh: Paul Jonas Jost, doi: https://doi.org/10.1016/j.chaos.2020.110032, code/data: ???) +# Oliveira_NatCommun2021 (gh: lorenapohl, doi: https://doi.org/10.1038/s41467-020-19798-3, code: https://github.com/cidacslab/Mathematical-and-Statistical-Modeling-of-COVID19-in-Brazil (python) data: https://infovis.sei.ba.gov.br/covid19/ ) +# Perelson_Science1996 (gh: Philipp Staedter, doi: https://doi.org/10.1126/science.271.5255.1582, code/data: ???) +# Rahman_MBS2016 (gh: Yannik Schaelte, doi: https://doi.org/10.1016/j.mbs.2016.07.009, code: not available, data: table in paper ...) +# Raimundez_PCB2020 (gh: Elba Raimundez, doi: https://doi.org/10.1371/journal.pcbi.1007147, code/data: https://zenodo.org/record/2908234#.Y5hUUS8w3yw (d2d)) +# SalazarCavazos_MBoC2020 (gh: Dilan Pathirana, doi: https://doi.org/10.1091/mbc.E19-09-0548, code/data: supplement (BNGL)) +# Zhao_QuantBiol2020 (gh: Iva Ewert, doi: https://doi.org/10.1007/s40484-020-0199-0, code: not available, data: table in supp) +# +# covered by performance test: +# Froehlich_CellSystems2018 +# +# Unknown reasons: +# Chen_MSB2009 +# +# Confirmed to be working: +# TODO: extend +problems_for_llh_check = [ + "Bachmann_MSB2011", + "Beer_MolBioSystems2014", + "Boehm_JProteomeRes2014", + "Borghans_BiophysChem1997", + "Brannmark_JBC2010", + "Bruno_JExpBot2016", + "Crauste_CellSystems2017", + "Elowitz_Nature2000", + "Fiedler_BMCSystBiol2016", + "Fujita_SciSignal2010", + "Isensee_JCB2018", + "Lucarelli_CellSystems2018", + "Schwen_PONE2014", + "Smith_BMCSystBiol2013", + "Sneyd_PNAS2002", + "Weber_BMC2015", + "Zheng_PNAS2012", +] + +# all PEtab problems we need to import +problems = list( + sorted(set(problems_for_gradient_check + problems_for_llh_check)) +) @dataclass class GradientCheckSettings: + """Problem-specific settings for gradient checks.""" + # Absolute and relative tolerances for simulation atol_sim: float = 1e-16 rtol_sim: float = 1e-12 @@ -58,14 +150,18 @@ class GradientCheckSettings: atol_consistency: float = 1e-5 rtol_consistency: float = 1e-1 # Step sizes for finite difference gradient checks. - step_sizes = [ - 1e-1, - 5e-2, - 1e-2, - 1e-3, - 1e-4, - 1e-5, - ] + step_sizes: list[float] = field( + default_factory=lambda: [ + 2e-1, + 1e-1, + 5e-2, + 1e-2, + 5e-1, + 1e-3, + 1e-4, + 1e-5, + ] + ) rng_seed: int = 0 ss_sensitivity_mode: amici.SteadyStateSensitivityMode = ( amici.SteadyStateSensitivityMode.integrateIfNewtonFails @@ -88,6 +184,10 @@ class GradientCheckSettings: settings["Brannmark_JBC2010"] = GradientCheckSettings( ss_sensitivity_mode=amici.SteadyStateSensitivityMode.integrationOnly, ) +settings["Fujita_SciSignal2010"] = GradientCheckSettings( + atol_check=1e-7, + rtol_check=5e-4, +) settings["Giordano_Nature2020"] = GradientCheckSettings( atol_check=1e-6, rtol_check=1e-3, rng_seed=1 ) @@ -97,12 +197,15 @@ class GradientCheckSettings: noise_level=0.01, atol_consistency=1e-3, ) -settings["Okuonghae_ChaosSolitonsFractals2020"].step_sizes.extend([0.2, 0.005]) settings["Oliveira_NatCommun2021"] = GradientCheckSettings( # Avoid "root after reinitialization" atol_sim=1e-12, rtol_sim=1e-12, ) +settings["Raia_CancerResearch2011"] = GradientCheckSettings( + atol_check=1e-10, + rtol_check=1e-3, +) settings["Smith_BMCSystBiol2013"] = GradientCheckSettings( atol_sim=1e-10, rtol_sim=1e-10, @@ -130,10 +233,173 @@ class GradientCheckSettings: ) -debug = False -if debug: - debug_path = Path(__file__).parent / "debug" - debug_path.mkdir(exist_ok=True, parents=True) +@pytest.fixture(scope="session", params=problems, ids=problems) +def benchmark_problem(request): + """Fixture providing model and PEtab problem for a problem from + the benchmark problem collection.""" + problem_id = request.param + petab_problem = benchmark_models_petab.get_problem(problem_id) + if measurement_table_has_timepoint_specific_mappings( + petab_problem.measurement_df, + ): + petab.flatten_timepoint_specific_output_overrides(petab_problem) + + # Setup AMICI objects. + amici_model = import_petab_problem( + petab_problem, + model_output_dir=benchmark_outdir / problem_id, + ) + return problem_id, petab_problem, amici_model + + +@pytest.mark.filterwarnings( + "ignore:divide by zero encountered in log", + # https://github.com/AMICI-dev/AMICI/issues/18 + "ignore:Adjoint sensitivity analysis for models with discontinuous " + "right hand sides .*:UserWarning", +) +def test_nominal_parameters_llh(benchmark_problem): + """Test the log-likelihood computation at nominal parameters. + + Also check that the simulation time is within the reference range. + """ + problem_id, petab_problem, amici_model = benchmark_problem + if problem_id not in problems_for_llh_check: + pytest.skip("Excluded from log-likelihood check.") + + amici_solver = amici_model.getSolver() + amici_solver.setAbsoluteTolerance(1e-8) + amici_solver.setRelativeTolerance(1e-8) + amici_solver.setMaxSteps(10_000) + if problem_id in ("Brannmark_JBC2010", "Isensee_JCB2018"): + amici_model.setSteadyStateSensitivityMode( + amici.SteadyStateSensitivityMode.integrationOnly + ) + + times = dict() + + for label, sensi_mode in { + "t_sim": amici.SensitivityMethod.none, + "t_fwd": amici.SensitivityMethod.forward, + "t_adj": amici.SensitivityMethod.adjoint, + }.items(): + amici_solver.setSensitivityMethod(sensi_mode) + if sensi_mode == amici.SensitivityMethod.none: + amici_solver.setSensitivityOrder(amici.SensitivityOrder.none) + else: + amici_solver.setSensitivityOrder(amici.SensitivityOrder.first) + + res_repeats = [ + simulate_petab( + petab_problem=petab_problem, + amici_model=amici_model, + solver=amici_solver, + log_level=logging.DEBUG, + ) + for _ in range(3) # repeat to get more stable timings + ] + res = res_repeats[0] + + times[label] = np.min( + [ + sum(r.cpu_time + r.cpu_timeB for r in res[RDATAS]) / 1000 + # only forwards/backwards simulation + for res in res_repeats + ] + ) + + if sensi_mode == amici.SensitivityMethod.none: + rdatas = res[RDATAS] + llh = res[LLH] + # TODO: check that all llh match, check that all sllh match + times["np"] = sum(petab_problem.parameter_df[petab.ESTIMATE]) + + pd.Series(times).to_csv(script_dir / f"{problem_id}_benchmark.csv") + for rdata in rdatas: + assert ( + rdata.status == amici.AMICI_SUCCESS + ), f"Simulation failed for {rdata.id}" + + if debug: + # create simulation PEtab table + sim_df = rdatas_to_measurement_df( + rdatas=rdatas, + model=amici_model, + measurement_df=petab_problem.measurement_df, + ) + sim_df.rename( + columns={petab.MEASUREMENT: petab.SIMULATION}, inplace=True + ) + sim_df.to_csv( + benchmark_outdir / problem_id / f"{problem_id}_simulation.tsv", + index=False, + sep="\t", + ) + + # visualize fit, save to file + with contextlib.suppress(NotImplementedError): + axs = plot_problem( + petab_problem=petab_problem, simulations_df=sim_df + ) + for plot_id, ax in axs.items(): + fig_path = ( + benchmark_outdir + / problem_id + / f"{problem_id}_{plot_id}_vis.png" + ) + logger.info(f"Saving figure to {fig_path}") + ax.get_figure().savefig(fig_path, dpi=150) + + try: + ref_llh = reference_values[problem_id]["llh"] + rdiff = np.abs((llh - ref_llh) / ref_llh) + rtol = 1e-3 + adiff = np.abs(llh - ref_llh) + atol = 1e-3 + tolstr = ( + f" Absolute difference is {adiff:.2e} " + f"(tol {atol:.2e}) and relative difference is " + f"{rdiff:.2e} (tol {rtol:.2e})." + ) + + if np.isclose(llh, ref_llh, rtol=rtol, atol=atol): + logger.info( + f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}." + + tolstr + ) + else: + pytest.fail( + f"Computed llh {llh:.4e} does not match reference " + f"{ref_llh:.4e}." + tolstr + ) + except KeyError: + logger.error( + "No reference likelihood found for " + f"{problem_id} in {references_yaml}" + ) + + for label, key in { + "simulation": "t_sim", + "adjoint sensitivity": "t_adj", + "forward sensitivity": "t_fwd", + }.items(): + try: + ref = reference_values[problem_id][key] + if times[key] > ref: + pytest.fail( + f"Computation time for {label} ({times[key]:.2e}) " + f"exceeds reference ({ref:.2e})." + ) + else: + logger.info( + f"Computation time for {label} ({times[key]:.2e}) " + f"within reference ({ref:.2e})." + ) + except KeyError: + logger.error( + f"No reference time for {label} found for " + f"{problem_id} in {references_yaml}" + ) @pytest.mark.filterwarnings( @@ -148,9 +414,14 @@ class GradientCheckSettings: (SensitivityMethod.forward, SensitivityMethod.adjoint), ids=["forward", "adjoint"], ) -@pytest.mark.parametrize("model", models) -def test_benchmark_gradient(model, scale, sensitivity_method, request): - if not scale and model in ( +def test_benchmark_gradient( + benchmark_problem, scale, sensitivity_method, request +): + problem_id, petab_problem, amici_model = benchmark_problem + if problem_id not in problems_for_gradient_check: + pytest.skip("Excluded from gradient check.") + + if not scale and problem_id in ( "Smith_BMCSystBiol2013", "Brannmark_JBC2010", "Elowitz_Nature2000", @@ -161,10 +432,10 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): ): # not really worth the effort trying to fix these cases if they # only fail on linear scale - pytest.skip() + pytest.skip("scale=False disabled for this problem") if ( - model + problem_id in ( # events with parameter-dependent triggers # https://github.com/AMICI-dev/AMICI/issues/18 @@ -172,20 +443,22 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): ) and sensitivity_method == SensitivityMethod.adjoint ): - # FIXME - pytest.skip() + pytest.skip("Unsupported ASA+events") - petab_problem = benchmark_models_petab.get_problem(model) - petab.flatten_timepoint_specific_output_overrides(petab_problem) + petab_problem = benchmark_models_petab.get_problem(problem_id) + if measurement_table_has_timepoint_specific_mappings( + petab_problem.measurement_df, + ): + petab.flatten_timepoint_specific_output_overrides(petab_problem) # Only compute gradient for estimated parameters. parameter_ids = petab_problem.x_free_ids - cur_settings = settings[model] + cur_settings = settings[problem_id] # Setup AMICI objects. amici_model = import_petab_problem( petab_problem, - model_output_dir=benchmark_outdir / model, + model_output_dir=benchmark_outdir / problem_id, ) amici_solver = amici_model.getSolver() amici_solver.setAbsoluteTolerance(cur_settings.atol_sim) @@ -299,9 +572,18 @@ def assert_gradient_check_success( df["rtol_success"] = df["rel_diff"] <= rtol max_adiff = df["abs_diff"].max() max_rdiff = df["rel_diff"].max() - with pd.option_context("display.max_columns", None, "display.width", None): + + success_fail = "succeeded" if check_result.success else "failed" + with pd.option_context( + "display.max_columns", + None, + "display.width", + None, + "display.max_rows", + None, + ): message = ( - f"Gradient check failed:\n{df}\n\n" + f"Gradient check {success_fail}:\n{df}\n\n" f"Maximum absolute difference: {max_adiff} (tolerance: {atol})\n" f"Maximum relative difference: {max_rdiff} (tolerance: {rtol})" ) diff --git a/tests/benchmark-models/test_petab_model.py b/tests/benchmark-models/test_petab_model.py deleted file mode 100755 index c4ec2f5dd2..0000000000 --- a/tests/benchmark-models/test_petab_model.py +++ /dev/null @@ -1,263 +0,0 @@ -#!/usr/bin/env python3 - -""" -Simulate a PEtab problem and compare results to reference values -""" - -import argparse -import contextlib -import importlib -import logging -import os -import sys -from pathlib import Path - -import amici -import numpy as np -import pandas as pd -import petab.v1 as petab -import yaml -from amici.logging import get_logger -from amici.petab.simulations import ( - LLH, - RDATAS, - rdatas_to_measurement_df, - simulate_petab, -) -from petab.v1.visualize import plot_problem - -logger = get_logger(f"amici.{__name__}", logging.WARNING) - - -def parse_cli_args(): - """Parse command line arguments - - Returns: - Parsed CLI arguments from ``argparse``. - """ - - parser = argparse.ArgumentParser( - description="Simulate PEtab-format model using AMICI." - ) - - # General options: - parser.add_argument( - "-v", - "--verbose", - dest="verbose", - action="store_true", - help="More verbose output", - ) - parser.add_argument( - "-c", - "--check", - dest="check", - action="store_true", - help="Compare to reference value", - ) - parser.add_argument( - "-p", - "--plot", - dest="plot", - action="store_true", - help="Plot measurement and simulation results", - ) - - # PEtab problem - parser.add_argument( - "-y", - "--yaml", - dest="yaml_file_name", - required=True, - help="PEtab YAML problem filename", - ) - - # Corresponding AMICI model - parser.add_argument( - "-m", - "--model-name", - dest="model_name", - help="Name of the AMICI module of the model to " "simulate.", - required=True, - ) - parser.add_argument( - "-d", - "--model-dir", - dest="model_directory", - help="Directory containing the AMICI module of the " - "model to simulate. Required if model is not " - "in python path.", - ) - - parser.add_argument( - "-o", - "--simulation-file", - dest="simulation_file", - help="File to write simulation result to, in PEtab" - "measurement table format.", - ) - - return parser.parse_args() - - -def main(): - """Simulate the model specified on the command line""" - script_dir = Path(__file__).parent.absolute() - args = parse_cli_args() - loglevel = logging.DEBUG if args.verbose else logging.INFO - logger.setLevel(loglevel) - - logger.info( - f"Simulating '{args.model_name}' " - f"({args.model_directory}) using PEtab data from " - f"{args.yaml_file_name}" - ) - - # load PEtab files - problem = petab.Problem.from_yaml(args.yaml_file_name) - petab.flatten_timepoint_specific_output_overrides(problem) - - # load model - if args.model_directory: - sys.path.insert(0, args.model_directory) - model_module = importlib.import_module(args.model_name) - amici_model = model_module.getModel() - amici_solver = amici_model.getSolver() - - amici_solver.setAbsoluteTolerance(1e-8) - amici_solver.setRelativeTolerance(1e-8) - amici_solver.setMaxSteps(int(1e4)) - if args.model_name in ("Brannmark_JBC2010", "Isensee_JCB2018"): - amici_model.setSteadyStateSensitivityMode( - amici.SteadyStateSensitivityMode.integrationOnly - ) - - times = dict() - - for label, sensi_mode in { - "t_sim": amici.SensitivityMethod.none, - "t_fwd": amici.SensitivityMethod.forward, - "t_adj": amici.SensitivityMethod.adjoint, - }.items(): - amici_solver.setSensitivityMethod(sensi_mode) - if sensi_mode == amici.SensitivityMethod.none: - amici_solver.setSensitivityOrder(amici.SensitivityOrder.none) - else: - amici_solver.setSensitivityOrder(amici.SensitivityOrder.first) - - res_repeats = [ - simulate_petab( - petab_problem=problem, - amici_model=amici_model, - solver=amici_solver, - log_level=loglevel, - ) - for _ in range(3) # repeat to get more stable timings - ] - res = res_repeats[0] - - times[label] = np.min( - [ - sum(r.cpu_time + r.cpu_timeB for r in res[RDATAS]) / 1000 - # only forwards/backwards simulation - for res in res_repeats - ] - ) - - if sensi_mode == amici.SensitivityMethod.none: - rdatas = res[RDATAS] - llh = res[LLH] - - times["np"] = sum(problem.parameter_df[petab.ESTIMATE]) - - pd.Series(times).to_csv(script_dir / f"{args.model_name}_benchmark.csv") - for rdata in rdatas: - assert ( - rdata.status == amici.AMICI_SUCCESS - ), f"Simulation failed for {rdata.id}" - - # create simulation PEtab table - sim_df = rdatas_to_measurement_df( - rdatas=rdatas, model=amici_model, measurement_df=problem.measurement_df - ) - sim_df.rename(columns={petab.MEASUREMENT: petab.SIMULATION}, inplace=True) - - if args.simulation_file: - sim_df.to_csv(args.simulation_file, index=False, sep="\t") - - if args.plot: - with contextlib.suppress(NotImplementedError): - # visualize fit - axs = plot_problem(petab_problem=problem, simulations_df=sim_df) - - # save figure - for plot_id, ax in axs.items(): - fig_path = os.path.join( - args.model_directory, - f"{args.model_name}_{plot_id}_vis.png", - ) - logger.info(f"Saving figure to {fig_path}") - ax.get_figure().savefig(fig_path, dpi=150) - - if args.check: - references_yaml = script_dir / "benchmark_models.yaml" - with open(references_yaml) as f: - refs = yaml.full_load(f) - - try: - ref_llh = refs[args.model_name]["llh"] - - rdiff = np.abs((llh - ref_llh) / ref_llh) - rtol = 1e-3 - adiff = np.abs(llh - ref_llh) - atol = 1e-3 - tolstr = ( - f" Absolute difference is {adiff:.2e} " - f"(tol {atol:.2e}) and relative difference is " - f"{rdiff:.2e} (tol {rtol:.2e})." - ) - - if np.isclose(llh, ref_llh, rtol=rtol, atol=atol): - logger.info( - f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}." - + tolstr - ) - else: - logger.error( - f"Computed llh {llh:.4e} does not match reference " - f"{ref_llh:.4e}." + tolstr - ) - sys.exit(1) - except KeyError: - logger.error( - "No reference likelihood found for " - f"{args.model_name} in {references_yaml}" - ) - - for label, key in { - "simulation": "t_sim", - "adjoint sensitivity": "t_adj", - "forward sensitivity": "t_fwd", - }.items(): - try: - ref = refs[args.model_name][key] - if times[key] > ref: - logger.error( - f"Computation time for {label} ({times[key]:.2e}) " - f"exceeds reference ({ref:.2e})." - ) - sys.exit(1) - else: - logger.info( - f"Computation time for {label} ({times[key]:.2e}) " - f"within reference ({ref:.2e})." - ) - except KeyError: - logger.error( - f"No reference time for {label} found for " - f"{args.model_name} in {references_yaml}" - ) - - -if __name__ == "__main__": - main() diff --git a/version.txt b/version.txt index 1b58cc1018..697f087f37 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.27.0 +0.28.0