diff --git a/doc/examples/example_errors.ipynb b/doc/examples/example_errors.ipynb index a34a95c93b..87da94c912 100644 --- a/doc/examples/example_errors.ipynb +++ b/doc/examples/example_errors.ipynb @@ -26,9 +26,19 @@ "import numpy as np\n", "\n", "import amici\n", + "from amici import (\n", + " simulation_status_to_str,\n", + " SteadyStateSensitivityMode,\n", + " SensitivityMethod,\n", + " SensitivityOrder,\n", + " runAmiciSimulation,\n", + " SteadyStateStatus,\n", + ")\n", "from amici.petab.petab_import import import_petab_problem\n", "from amici.petab.simulations import simulate_petab, RDATAS, EDATAS\n", "from amici.plotting import plot_state_trajectories, plot_jacobian\n", + "from petab.v1.sbml import get_sbml_model\n", + "from amici import SbmlImporter, import_model_module\n", "\n", "try:\n", " import benchmark_models_petab\n", @@ -48,11 +58,11 @@ "source": [ "## Overview\n", "\n", - "In the following, we will simulate models contained in the [PEtab Benchmark Collection](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/) to demonstrate a number of simulation failures to analyze and fix them. We use the PEtab format, as it makes model import and simulation much easier, but everything illustrated here, also applies to plain SBML or PySB import.\n", + "In the following, we will simulate models contained in the [PEtab Benchmark Collection](https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/) and the [SBML test suite](https://github.com/sbmlteam/sbml-test-suite/) to demonstrate a number of simulation failures to analyze and fix them. We use the PEtab format, as it makes model import and simulation much easier, but everything illustrated here, also applies to plain SBML or PySB import.\n", "\n", "Note that, due to numerical issues, the examples below may not be fully reproducible on every system.\n", "\n", - "If any simulation failures occur, they will be printed via Python logging.\n", + "If any simulation failures occur, they will be printed via Python logging. Additionally, they will be stored in `ReturnData.messages`.\n", "\n", "Programmatically, simulation success can be checked via `ReturnDataView.status`. In case of a successful simulation, and only then, this value corresponds to `amici.AMICI_SUCCESS`.\n", "In case of a simulation error, all quantities in `ReturnData`/`ReturnDataView` will be reported up to the time of failure, the rest will be `NaN`. The likelihood and it's gradient will always be `NaN` in case of failure." @@ -93,11 +103,9 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", - "] == [\n", + "assert [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", " \"AMICI_SUCCESS\",\n", " \"AMICI_SUCCESS\",\n", " \"AMICI_SUCCESS\",\n", @@ -149,7 +157,7 @@ "\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])\n", "print(\"Simulations finished successfully.\")\n", @@ -170,7 +178,7 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])\n", "print(\"Simulations finished successfully.\")" @@ -211,11 +219,11 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", - "] == [\"AMICI_TOO_MUCH_WORK\"]" + "assert [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", + " \"AMICI_TOO_MUCH_WORK\"\n", + "]" ] }, { @@ -245,7 +253,7 @@ "edata = amici.ExpData(res[EDATAS][0])\n", "edata.setTimepoints(np.linspace(0, 0.33011, 5000))\n", "amici_solver = amici_model.getSolver()\n", - "rdata = amici.runAmiciSimulation(amici_model, amici_solver, edata)\n", + "rdata = runAmiciSimulation(amici_model, amici_solver, edata)\n", "\n", "# Visualize state trajectories\n", "plot_state_trajectories(rdata, model=amici_model)\n", @@ -315,12 +323,10 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", - "assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", - "] == [\n", + "assert [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", " \"AMICI_SUCCESS\",\n", " \"AMICI_ERR_FAILURE\",\n", " \"AMICI_NOT_RUN\",\n", @@ -396,7 +402,7 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])\n", "print(\"Simulations finished successfully.\")" @@ -412,6 +418,14 @@ "Let's run a simulation:" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "d91e78142032ba01", + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, @@ -419,32 +433,28 @@ "metadata": {}, "outputs": [], "source": [ - "petab_problem = benchmark_models_petab.get_problem(\"Weber_BMC2015\")\n", - "amici_model = import_petab_problem(petab_problem, verbose=False, compile_=None)\n", - "\n", - "np.random.seed(4)\n", - "problem_parameters = dict(\n", - " zip(\n", - " petab_problem.x_free_ids,\n", - " petab_problem.sample_parameter_startpoints(n_starts=1)[0],\n", + "# retrieve and import the SBML model\n", + "if \"imported_00885\" not in globals():\n", + " sbml_reader, sbml_document, sbml_model = get_sbml_model(\n", + " \"https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/7ab011efe471877987b5d6dcba8cc19a6ff71254/cases/semantic/00885/00885-sbml-l3v2.xml\"\n", " )\n", - ")\n", - "res = simulate_petab(\n", - " petab_problem=petab_problem,\n", - " amici_model=amici_model,\n", - " problem_parameters=problem_parameters,\n", - " scaled_parameters=True,\n", - ")\n", - "print(\n", - " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", - ")\n", - "assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", - "] == [\n", - " \"AMICI_ERROR\",\n", - " \"AMICI_NOT_RUN\",\n", - "]" + " SbmlImporter(sbml_model).sbml2amici(\n", + " model_name=\"sbml_00885\", output_dir=\"amici_models/sbml_00885\"\n", + " )\n", + " imported_00885 = True\n", + "\n", + "amici_model = import_model_module(\n", + " \"sbml_00885\", \"amici_models/sbml_00885\"\n", + ").get_model()\n", + "\n", + "# run the simulation\n", + "amici_solver = amici_model.getSolver()\n", + "amici_solver.setSensitivityMethod(SensitivityMethod.forward)\n", + "amici_solver.setSensitivityOrder(SensitivityOrder.first)\n", + "amici_model.setTimepoints(np.linspace(0, 6, 51))\n", + "rdata = runAmiciSimulation(amici_model, amici_solver)\n", + "print(\"Status:\", simulation_status_to_str(rdata.status))\n", + "assert simulation_status_to_str(rdata.status) == \"AMICI_ERROR\"" ] }, { @@ -454,7 +464,7 @@ "source": [ "**What happened?**\n", "\n", - "The simulation failed because the initial step-size after an event or heaviside function was too small. The error occurred during simulation of condition `model1_data1` after successful preequilibration (`model1_data2`)." + "The simulation failed because the initial step-size after an event or Heaviside function was too small. The error occurred around `t=2.1`." ] }, { @@ -464,7 +474,7 @@ "source": [ "**How to address?**\n", "\n", - "The error message already suggests a fix for this situation, so let's try increasing the relative tolerance:" + "The error message already suggests a fix for this situation, so let's try decreasing the relative tolerance:" ] }, { @@ -475,29 +485,21 @@ "outputs": [], "source": [ "amici_solver = amici_model.getSolver()\n", - "amici_solver.setRelativeTolerance(200 * amici_solver.getRelativeTolerance())\n", - "\n", - "np.random.seed(4)\n", - "problem_parameters = dict(\n", - " zip(\n", - " petab_problem.x_free_ids,\n", - " petab_problem.sample_parameter_startpoints(n_starts=1)[0],\n", - " )\n", - ")\n", - "res = simulate_petab(\n", - " petab_problem=petab_problem,\n", - " amici_model=amici_model,\n", - " problem_parameters=problem_parameters,\n", - " scaled_parameters=True,\n", - " solver=amici_solver,\n", - ")\n", - "print(\n", - " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", - ")\n", - "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])" + "amici_solver.setSensitivityMethod(SensitivityMethod.forward)\n", + "amici_solver.setSensitivityOrder(SensitivityOrder.first)\n", + "amici_solver.setRelativeTolerance(1e-15)\n", + "amici_model.setTimepoints(np.linspace(0, 6, 51))\n", + "rdata = runAmiciSimulation(amici_model, amici_solver)\n", + "print(\"Status:\", simulation_status_to_str(rdata.status))\n", + "assert simulation_status_to_str(rdata.status) == \"AMICI_SUCCESS\"" ] }, + { + "cell_type": "markdown", + "id": "4f8557aa8e05f6cc", + "metadata": {}, + "source": "Sometimes this problem can also be solved be *in*creasing tolerances or by slightly changing the output timepoints." + }, { "cell_type": "markdown", "id": "86f4db29", @@ -533,11 +535,11 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", - "] == [\"AMICI_FIRST_RHSFUNC_ERR\"]" + "assert [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", + " \"AMICI_FIRST_RHSFUNC_ERR\"\n", + "]" ] }, { @@ -700,10 +702,10 @@ ")\n", "\n", "amici_solver = amici_model.getSolver()\n", - "amici_solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", - "amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)\n", + "amici_solver.setSensitivityMethod(SensitivityMethod.forward)\n", + "amici_solver.setSensitivityOrder(SensitivityOrder.first)\n", "amici_model.setSteadyStateSensitivityMode(\n", - " amici.SteadyStateSensitivityMode.newtonOnly\n", + " SteadyStateSensitivityMode.newtonOnly\n", ")\n", "\n", "np.random.seed(2020)\n", @@ -722,13 +724,13 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "# hard to reproduce on GHA\n", "if os.getenv(\"GITHUB_ACTIONS\") is None:\n", " assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + " simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", " ] == [\"AMICI_ERROR\"]" ] }, @@ -781,7 +783,7 @@ "source": [ "# use numerical integration\n", "amici_model.setSteadyStateSensitivityMode(\n", - " amici.SteadyStateSensitivityMode.integrationOnly\n", + " SteadyStateSensitivityMode.integrationOnly\n", ")\n", "\n", "res = simulate_petab(\n", @@ -793,7 +795,7 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])" ] @@ -820,8 +822,8 @@ "del os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"]\n", "\n", "amici_solver = amici_model.getSolver()\n", - "amici_solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", - "amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)\n", + "amici_solver.setSensitivityMethod(SensitivityMethod.forward)\n", + "amici_solver.setSensitivityOrder(SensitivityOrder.first)\n", "\n", "res = simulate_petab(\n", " petab_problem=petab_problem,\n", @@ -832,7 +834,7 @@ ")\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])" ] @@ -880,13 +882,13 @@ "\n", "print(\n", " \"Status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "# hard to reproduce on GHA\n", "if os.getenv(\"GITHUB_ACTIONS\") is None:\n", " assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + " simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", " ] == [\n", " \"AMICI_ERROR\",\n", " \"AMICI_NOT_RUN\",\n", @@ -919,7 +921,7 @@ "outputs": [], "source": [ "rdata = res[RDATAS][0]\n", - "list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))" + "list(map(SteadyStateStatus, rdata.preeq_status.flatten()))" ] }, { @@ -962,12 +964,12 @@ ")\n", "print(\n", " \"status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", "print(\n", - " f\"preeq_status={list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", + " f\"preeq_status={list(map(SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", ")\n", "print(f\"{rdata.preeq_numsteps=}\")\n", "\n", @@ -1009,12 +1011,12 @@ "\n", "print(\n", " \"status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", "print(\n", - " f\"preeq_status={list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", + " f\"preeq_status={list(map(SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", ")\n", "print(f\"{rdata.preeq_numsteps=}\")\n", "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])" @@ -1050,18 +1052,18 @@ ")\n", "print(\n", " \"status:\",\n", - " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", + " [simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", "print(\n", - " f\"preeq_status={list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", + " f\"preeq_status={list(map(SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n", ")\n", "print(f\"{rdata.preeq_numsteps=}\")\n", "# hard to reproduce on GHA\n", "if os.getenv(\"GITHUB_ACTIONS\") is None:\n", " assert [\n", - " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + " simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", " ] == [\n", " \"AMICI_ERROR\",\n", " \"AMICI_NOT_RUN\",\n", diff --git a/src/solver_cvodes.cpp b/src/solver_cvodes.cpp index 4a4f0dc87f..89eb45b48b 100644 --- a/src/solver_cvodes.cpp +++ b/src/solver_cvodes.cpp @@ -541,8 +541,8 @@ void CVodeSolver::reInitPostProcess( = std::string("CVode returned a root after reinitialization at t=") + std::to_string(*t) + ". The initial step-size after the event or " - "Heaviside function is too small. To fix this, increase " - "absolute and relative tolerances!"; + "Heaviside function is too small. To fix this, adjust " + "absolute or relative tolerances!"; throw CvodeException(status, message.c_str()); } if (status != CV_SUCCESS) {