diff --git a/.github/workflows/test_python_cplusplus.yml b/.github/workflows/test_python_cplusplus.yml index e651ff01e8..49e300c03e 100644 --- a/.github/workflows/test_python_cplusplus.yml +++ b/.github/workflows/test_python_cplusplus.yml @@ -288,7 +288,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 @@ -350,7 +351,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/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/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/solver_cvodes.cpp b/src/solver_cvodes.cpp index dfccbb2389..0b60304c55 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -34,8 +34,10 @@ static_assert((int)InterpolationType::polynomial == CV_POLYNOMIAL, ""); static_assert((int)LinearMultistepMethod::adams == CV_ADAMS, ""); static_assert((int)LinearMultistepMethod::BDF == CV_BDF, ""); -#define STATIC_ASSERT_EQUAL(amici_constant, cv_constant) \ -static_assert(amici_constant == cv_constant, #amici_constant " != " #cv_constant) +#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); 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/version.txt b/version.txt index 1b58cc1018..697f087f37 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -0.27.0 +0.28.0