From a20447a0982f014f171985dbd3b0c58e7632e026 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 23 Oct 2024 19:08:20 +0200 Subject: [PATCH 1/3] Add missing simulation status codes (#2560) Closes #2559. --- include/amici/defines.h | 6 +++ src/amici.cpp | 6 +++ src/solver_cvodes.cpp | 83 ++++++++++++----------------------------- 3 files changed, 36 insertions(+), 59 deletions(-) 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/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..dfccbb2389 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -34,65 +34,30 @@ 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. From 9c603e1015ff0f233078f4e741b564bb9e9dc6c8 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 24 Oct 2024 10:53:08 +0200 Subject: [PATCH 2/3] PEtab import: Fix potentially incorrect sensitivities with state-dependent sigmas (#2562) Due to PEtab not allowing and observables in noiseFormula and AMICI not allowing state variables in sigma expressions, we are trying to replace any state variables by matching observables during PEtab import. However, if there were multiple observables with the same observableFormula, or one observableFormula was contained in another one, the substition of observableFormula for observableId could - depending on the ordering - result in noiseFormulas depending on other observables, which could lead to incorrect sensitivities due to https://github.com/AMICI-dev/AMICI/issues/2561. This change ensures that we are only using the observableId from the same row of the observable table as the noiseFormula. --------- Co-authored-by: Dilan Pathirana <59329744+dilpath@users.noreply.github.com> --- python/sdist/amici/petab/import_helpers.py | 26 ++++++++++------------ 1 file changed, 12 insertions(+), 14 deletions(-) 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) From 56e69568ad4c5de3f074b0d9143eeadf5e82367b Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 24 Oct 2024 15:05:55 +0200 Subject: [PATCH 3/3] Fix benchmark collection gradient check (#2564) * Fix instance vs class attribute for step sizes * Only flatten problems where necessary --- .../test_benchmark_collection.sh | 15 +++++++++- .../benchmark-models/test_petab_benchmark.py | 29 ++++++++++++------- tests/benchmark-models/test_petab_model.py | 7 ++++- 3 files changed, 38 insertions(+), 13 deletions(-) diff --git a/tests/benchmark-models/test_benchmark_collection.sh b/tests/benchmark-models/test_benchmark_collection.sh index 581b8db028..d3d1cad712 100755 --- a/tests/benchmark-models/test_benchmark_collection.sh +++ b/tests/benchmark-models/test_benchmark_collection.sh @@ -93,9 +93,22 @@ for model in $models; do yaml="${model_dir}"/"${model}"/problem.yaml fi + # problems we need to flatten + to_flatten=( + "Bruno_JExpBot2016" "Chen_MSB2009" "Crauste_CellSystems2017" + "Fiedler_BMCSystBiol2016" "Fujita_SciSignal2010" "SalazarCavazos_MBoC2020" + ) + flatten="" + for item in "${to_flatten[@]}"; do + if [[ "$item" == "$model" ]]; then + flatten="--flatten" + break + fi + done + 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_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} diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 4892100877..69df16f181 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -12,8 +12,9 @@ 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 @@ -58,14 +59,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 @@ -97,7 +102,6 @@ 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, @@ -176,7 +180,10 @@ def test_benchmark_gradient(model, scale, sensitivity_method, request): pytest.skip() petab_problem = benchmark_models_petab.get_problem(model) - petab.flatten_timepoint_specific_output_overrides(petab_problem) + 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 diff --git a/tests/benchmark-models/test_petab_model.py b/tests/benchmark-models/test_petab_model.py index c4ec2f5dd2..125a046a5e 100755 --- a/tests/benchmark-models/test_petab_model.py +++ b/tests/benchmark-models/test_petab_model.py @@ -25,6 +25,7 @@ simulate_petab, ) from petab.v1.visualize import plot_problem +from petab.v1.lint import measurement_table_has_timepoint_specific_mappings logger = get_logger(f"amici.{__name__}", logging.WARNING) @@ -115,7 +116,11 @@ def main(): # load PEtab files problem = petab.Problem.from_yaml(args.yaml_file_name) - petab.flatten_timepoint_specific_output_overrides(problem) + + if measurement_table_has_timepoint_specific_mappings( + problem.measurement_df + ): + petab.flatten_timepoint_specific_output_overrides(problem) # load model if args.model_directory: