python/examples/example_errors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d8e8890",
   "metadata": {},
   "source": [
    "# Debugging simulation failures\n",
    "\n",
    "**Objective:** Demonstrate common simulation failures and give some hints for interpreting, debugging, and fixing them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87383e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "from contextlib import suppress\n",
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import amici\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",
    "\n",
    "try:\n",
    "    import benchmark_models_petab\n",
    "except ModuleNotFoundError:\n",
    "    # install `benchmark_models_petab` if necessary\n",
    "    %pip install -q -e \"git+https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab.git@master#subdirectory=src/python&egg=benchmark_models_petab\"\n",
    "    try:\n",
    "        import benchmark_models_petab\n",
    "    except ModuleNotFoundError:\n",
    "        print(\"** Please restart the kernel. **\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46500fb0",
   "metadata": {},
   "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",
    "\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",
    "\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."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ad882ac",
   "metadata": {},
   "source": [
    "## `AMICI_TOO_MUCH_WORK` - `mxstep steps taken before reaching tout`\n",
    "\n",
    "Let's run a simulation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae3cb45e",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Fujita_SciSignal2010\")\n",
    "amici_model = import_petab_problem(petab_problem, verbose=False, compile_=None)\n",
    "\n",
    "np.random.seed(2991)\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",
    ")\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_SUCCESS\",\n",
    "    \"AMICI_SUCCESS\",\n",
    "    \"AMICI_SUCCESS\",\n",
    "    \"AMICI_TOO_MUCH_WORK\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e75a6dcc",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "AMICI failed to integrate the forward problem. The problem occurred for only one simulation condition, `condition_step_03_0`. The issue occurred at $t = 3031.8$, where the CVODES reached the maximum number of steps.\n",
    "\n",
    "**How to address?**\n",
    "\n",
    "The number of steps the solver has to take is closely related to the chosen error tolerance. More accurate results, more steps. Therefore, this problem can be solved in two ways:\n",
    "\n",
    "1. Increasing the maximum number of steps via [amici.Solver.setMaxSteps](https://amici.readthedocs.io/en/latest/generated/amici.amici.Solver.html#amici.amici.Solver.setMaxSteps). Note that this will increase the time required for simulation, and that simulation may still fail eventually. Sometimes it may be preferable to not increase this limit but rather fail fast. Also note that increasing the number of allowed steps increase RAM requirements (even if fewer steps are actually taken), so don't set this to ridiculously large values in order to avoid this error.\n",
    "\n",
    "2. Reducing the number of steps CVODES has to take. This is determined by the required error tolerance. There are various solver error tolerances than can be adjusted. The most relevant ones are those controlled via [amici.Solver.setRelativeTolerance()](https://amici.readthedocs.io/en/latest/generated/amici.amici.Solver.html#amici.amici.Solver.setRelativeTolerance) and [amici.Solver.setAbsoluteTolerance()](https://amici.readthedocs.io/en/latest/generated/amici.amici.Solver.html#amici.amici.Solver.setAbsoluteTolerance).\n",
    "\n",
    "So, let's fix that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9667f55",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's increase the allowed number of steps by 10x:\n",
    "print(\"Increasing allowed number of steps ...\")\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setMaxSteps(10 * amici_solver.getMaxSteps())\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",
    "\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])\n",
    "print(\"Simulations finished successfully.\")\n",
    "print()\n",
    "\n",
    "\n",
    "# let's relax the relative error tolerance by a factor of 50\n",
    "print(\"Relaxing relative error tolerance ...\")\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setRelativeTolerance(50 * amici_solver.getRelativeTolerance())\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])\n",
    "print(\"Simulations finished successfully.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18fd3fa6",
   "metadata": {},
   "source": [
    "## `Internal t = [...] and h = [...] are such that t + h = t on the next step`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f78179dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Crauste_CellSystems2017\")\n",
    "amici_model = import_petab_problem(petab_problem, verbose=False)\n",
    "\n",
    "np.random.seed(1)\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",
    ")\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",
    "] == [\"AMICI_TOO_MUCH_WORK\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53e4b822",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "The forward simulation failed because the AMICI solver exceeded the maximum number of steps. Unlike in the previous case of  `mxstep steps taken before reaching tout` (see above), here we got several additional warnings that the current step size $h$ is numerically zero.\n",
    "\n",
    "**How to address?**\n",
    "\n",
    "The warning `Internal t = [...] and h = [...] are such that t + h = t on the next step` tells us that the solver is not able to move forward. The solver may be able to recover from that, but not always.\n",
    "\n",
    "Let's look at the state trajectories to see what's going on. Such a tiny step size is usually related to very fast dynamics. We repeat the simulation with additional timepoints before the point of failure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a6794d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a copy of this simulation condition\n",
    "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",
    "\n",
    "# Visualize state trajectories\n",
    "plot_state_trajectories(rdata, model=amici_model)\n",
    "plt.yscale(\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7fe6d08",
   "metadata": {},
   "source": [
    "We can see a steep increase for `Pathogen` just before the error occurs. Let's zoom in:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a645efb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(rdata.t, rdata.by_id(\"Pathogen\"))\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"Pathogen\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b85cd813",
   "metadata": {},
   "source": [
    "The solver is unable to handle such a steep increase. There is not much we can do. Increasing the tolerances will let the solver proceed a bit further, but this is usually not enough. Most likely there is a problem in the model or in the choice of parameter values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5afbd242",
   "metadata": {},
   "source": [
    "## `the error test failed repeatedly or with |h| = hmin`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb8910f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Fujita_SciSignal2010\")\n",
    "amici_model = import_petab_problem(petab_problem, verbose=False)\n",
    "\n",
    "np.random.seed(4920)\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",
    ")\n",
    "print(\n",
    "    \"Status:\",\n",
    "    [amici.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",
    "    \"AMICI_SUCCESS\",\n",
    "    \"AMICI_ERR_FAILURE\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "    \"AMICI_NOT_RUN\",\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b3a9904",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "AMICI failed to integrate the forward problem. The problem occurred for only one simulation condition, `condition_step_00_3`. The issue occurred at $t = 429.232$, where the error test failed.\n",
    "This means, the solver is unable to take a step of non-zero size without violating the chosen error tolerances."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c16ac6c8",
   "metadata": {},
   "source": [
    "**How to address?**\n",
    "\n",
    "The step size is computed based on the Jacobian. Inspecting `ReturnData.J` shows us that we have rather large values in the Jacobian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467c3d36",
   "metadata": {},
   "outputs": [],
   "source": [
    "rdata = res[RDATAS][1]\n",
    "\n",
    "# Show Jacobian as heatmap\n",
    "plot_jacobian(rdata)\n",
    "\n",
    "print(f\"largest absolute Jacobian value: {np.max(np.abs(rdata.J)):.3g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f6ec0e3",
   "metadata": {},
   "source": [
    "In this case, the default relative error tolerance may be too high and lead too large absolute errors. \n",
    "\n",
    "Let's retry simulation using stricter tolerances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "600ae826",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set stricter relative error tolerance\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setRelativeTolerance(amici_solver.getRelativeTolerance() / 10)\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])\n",
    "print(\"Simulations finished successfully.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "616710b6",
   "metadata": {},
   "source": [
    "## `Cvode routine CVode returned a root after reinitialization`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cec31332",
   "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",
    "    )\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",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4386603",
   "metadata": {},
   "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`)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "02f2dd5d",
   "metadata": {},
   "source": [
    "**How to address?**\n",
    "\n",
    "The error message already suggests a fix for this situation, so let's try increasing the relative tolerance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d1552e3",
   "metadata": {},
   "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])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f4db29",
   "metadata": {},
   "source": [
    "## `AMICI encountered a NaN / Inf value for [...]`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d97349ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Borghans_BiophysChem1997\")\n",
    "amici_model = import_petab_problem(petab_problem, verbose=False)\n",
    "\n",
    "np.random.seed(18)\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",
    ")\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",
    "] == [\"AMICI_FIRST_RHSFUNC_ERR\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63641cff",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "The forward simulation failed because AMICI encountered a `NaN` value when simulating condition `model1_data1`.\n",
    "Then `NaN`s occurred in $\\dot x$ and $w$ (model expressions, such as reaction fluxes or assignment rules). Furthermore, the failure occurred at the first call, so at $t = t_0$ (here: $t = 0$).\n",
    "\n",
    "**How to address?**\n",
    "\n",
    "The `NaN` in $\\dot x$ is most likely a consequence of the one in $w$. (A subset of) the dependency tree looks something like:\n",
    "\n",
    "[![](https://mermaid.ink/img/pako:eNpdkrFuwyAQhl8FIWVL1MwMndKFtd1wBmJIg2IDwucCivLuxQp2zvaA_P1398Od7kFbpzRl9DdIfyM_p8aS8t3FnZHW2QGkheH8Er3wjHgZZK9Bh1kFAYyA6XXldBTpyIixBozsSHGAJSQSWwlccEa4bN3FSFu1KCIjOvmgh8GUF8y1yoGYDkZU-lBQ5SwyI-4y6PAnL52esl-B3a68pDZDDofPhfyKYEVLacSVERdGXFchzYAu59iBYweOHTh2qBAxvLupI3s9eJ41pndqGdOqvYXThv2G7xuOiBf7jHMzNsr41oyvzNgv0z3tdeilUWXzHlOooXDTvW4oK79KX-XYQUMb-yypcgT3nW1LGYRR7-noVVmhk5FlZ3vKrrIbFvVLGXChis9_j9jNUw?type=png)](https://mermaid.live/edit#pako:eNpdkrFuwyAQhl8FIWVL1MwMndKFtd1wBmJIg2IDwucCivLuxQp2zvaA_P1398Od7kFbpzRl9DdIfyM_p8aS8t3FnZHW2QGkheH8Er3wjHgZZK9Bh1kFAYyA6XXldBTpyIixBozsSHGAJSQSWwlccEa4bN3FSFu1KCIjOvmgh8GUF8y1yoGYDkZU-lBQ5SwyI-4y6PAnL52esl-B3a68pDZDDofPhfyKYEVLacSVERdGXFchzYAu59iBYweOHTh2qBAxvLupI3s9eJ41pndqGdOqvYXThv2G7xuOiBf7jHMzNsr41oyvzNgv0z3tdeilUWXzHlOooXDTvW4oK79KX-XYQUMb-yypcgT3nW1LGYRR7-noVVmhk5FlZ3vKrrIbFvVLGXChis9_j9jNUw)\n",
    "\n",
    "Always look for the most basic (furthest up) model quantities.\n",
    "In cases where non-finite values occur in expressions further down, rerunning the simulation after calling `Model.setAlwaysCheckFinite(True)` may give some further hints on where the issue originates.\n",
    "\n",
    "The `NaN` in $w$ occurred for `flux_v7_v_6` (see error log), i.e., when computing the reaction flux for reaction `v7_v_6`. As $w$ only depends on $(t, p, k, x)$ and no non-finite values have been reported for those, the issue has to be in the respective flux equation.\n",
    "\n",
    "Let's look at that expression. This can either be done by inspecting the underlying SBML model (e.g., using COPASI), or by checking the generated model code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fca95143",
   "metadata": {},
   "outputs": [],
   "source": [
    "# model source code location\n",
    "model_src_dir = Path(amici_model.module.__file__).parents[1]\n",
    "\n",
    "# find the problematic expression in the model source code\n",
    "!grep flux_v7_v_6 {model_src_dir}/w.cpp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f49a00a",
   "metadata": {},
   "source": [
    "What could go wrong? We can obtain `NaN` from any of these symbols being `NaN`, or through division by zero.\n",
    "\n",
    "Let's let's check the denominator first: $$(A\\_state^2 + Kp^2)*(Kd^{n\\_par} + Z\\_state^{n\\_par})$$\n",
    "\n",
    "\n",
    "`A_state` and `Z_state` are state variables, `Kd`, `K_p`, and `n_par` are parameters.\n",
    "\n",
    "As the error occurred at $t = t_0$, let's ensure the initial state is non-zero and finite:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eec6fe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "rdata = res[RDATAS][0]\n",
    "edata = res[EDATAS][0]\n",
    "# check initial states\n",
    "x0 = dict(zip(amici_model.getStateIds(), rdata.x0))\n",
    "print(f\"{x0=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "050c4c38",
   "metadata": {},
   "source": [
    "The initial states are fine - the first multiplicand is non-zero, as $x_0$ was non-zero. \n",
    "\n",
    "So let's check the parameter values occurring in the second multiplicand:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec724ce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# we have to account for the chosen parameter scale\n",
    "from itertools import starmap\n",
    "\n",
    "unscaled_parameter = dict(\n",
    "    zip(\n",
    "        amici_model.getParameterIds(),\n",
    "        starmap(\n",
    "            amici.getUnscaledParameter, zip(edata.parameters, edata.pscale)\n",
    "        ),\n",
    "    )\n",
    ")\n",
    "print(dict((p, unscaled_parameter[p]) for p in (\"Kd\", \"Kp\", \"n_par\")))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d82971",
   "metadata": {},
   "source": [
    "Considering that `n_par` occurs as exponent, it's magnitude looks pretty high.\n",
    "This term is very likely causing the problem - let's check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f3f8bdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"{x0['Z_state']**unscaled_parameter['n_par'] + unscaled_parameter['Kd']**unscaled_parameter['n_par']=}\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6616c2ff",
   "metadata": {},
   "source": [
    "Indeed, no way we can fix this for the given model.\n",
    "This was most likely an unrealistic parameter value, originating from a too high upper parameter bound for `n_par`.\n",
    "Therefore, if this error occurs during optimization, a first step could be adapting the respective parameter bounds.\n",
    "In other cases, this may be a result of unfortunate arrangement of model expressions, which can sometimes be solved by passing a suitable simplification function to the model import."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22cfbbbf",
   "metadata": {},
   "source": [
    "<a id='unsuccessful_factorization'></a>\n",
    "\n",
    "## `Steady state sensitivity computation failed due to unsuccessful factorization of RHS Jacobian`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b41a5017",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Blasi_CellSystems2016\")\n",
    "with suppress(KeyError):\n",
    "    del os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"]\n",
    "amici_model = import_petab_problem(\n",
    "    petab_problem,\n",
    "    verbose=False,\n",
    "    compile_=True,\n",
    "    model_name=\"Blasi_CellSystems2016_1\",\n",
    ")\n",
    "\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n",
    "amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)\n",
    "amici_model.setSteadyStateSensitivityMode(\n",
    "    amici.SteadyStateSensitivityMode.newtonOnly\n",
    ")\n",
    "\n",
    "np.random.seed(2020)\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",
    "\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",
    "    ] == [\"AMICI_ERROR\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82267e47",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "AMICI failed to compute steady-state sensitivities, because it was not able to factorize the Jacobian.\n",
    "\n",
    "**How to address?**\n",
    "\n",
    "This is most likely a result of a singular Jacobian. Let's check the condition number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cd349b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "rdata = res[RDATAS][0]\n",
    "np.linalg.cond(rdata.J)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da131d3e",
   "metadata": {},
   "source": [
    "Indeed, the condition number shows that the Jacobian is numerically singular. If this happens consistently, it is usually due to conserved quantities in the model.\n",
    "\n",
    "There are two ways we can address that:\n",
    "\n",
    "1. Use numerical integration to compute sensitivities, for which a singular Jacobian is not an issue. This is, usually, slower, though.\n",
    "2. Remove any conserved quantities.\n",
    "\n",
    "Let's try both approaches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f82078e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# use numerical integration\n",
    "amici_model.setSteadyStateSensitivityMode(\n",
    "    amici.SteadyStateSensitivityMode.integrationOnly\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])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d7be541",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove conserved quantities - this requires re-importing the model\n",
    "\n",
    "# this is enabled by the `AMICI_EXPERIMENTAL_SBML_NONCONST_CLS` environment variable\n",
    "os.environ[\"AMICI_EXPERIMENTAL_SBML_NONCONST_CLS\"] = \"1\"\n",
    "amici_model = import_petab_problem(\n",
    "    petab_problem,\n",
    "    verbose=False,\n",
    "    # we need a different model name if we import the model again\n",
    "    # we cannot load a model with the same name as an already loaded model\n",
    "    model_name=\"Blasi_CellSystems2016_2\",\n",
    "    compile_=True,\n",
    ")\n",
    "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",
    "\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])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b977c0b",
   "metadata": {},
   "source": [
    "## `Steady state computation failed`\n",
    "\n",
    "Let's run a simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97f797dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_problem = benchmark_models_petab.get_problem(\"Brannmark_JBC2010\")\n",
    "amici_model = import_petab_problem(\n",
    "    petab_problem,\n",
    "    verbose=False,\n",
    ")\n",
    "\n",
    "amici_solver = amici_model.getSolver()\n",
    "\n",
    "\n",
    "np.random.seed(1851)\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",
    "\n",
    "print(\n",
    "    \"Status:\",\n",
    "    [amici.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",
    "    ] == [\n",
    "        \"AMICI_ERROR\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0713b830",
   "metadata": {},
   "source": [
    "**What happened?**\n",
    "\n",
    "All given experimental conditions require pre-equilibration, i.e., finding a steady state. AMICI first tries to find a steady state using the Newton solver, if that fails, it tries simulating until steady state, if that also fails, it tries the Newton solver from the end of the simulation. In this case, all three failed. Neither Newton's method nor simulation yielded a steady state satisfying the required tolerances.\n",
    "\n",
    "This can also be seen in `ReturnDataView.preeq_status` (the three statuses corresponds to Newton \\#1, Simulation, Newton \\#2):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffdc8e82",
   "metadata": {},
   "outputs": [],
   "source": [
    "rdata = res[RDATAS][0]\n",
    "list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "189dd964",
   "metadata": {},
   "source": [
    "**How to address?**\n",
    "\n",
    "There are several ways to address that:\n",
    "\n",
    "1. Stricter integration tolerances (preferred if affordable - higher accuracy, but generally slower)\n",
    "\n",
    "2. Looser steady-state tolerances (lower accuracy, generally faster)\n",
    "\n",
    "3. Increase the number of allowed steps for Newton's method\n",
    "\n",
    "Let's try all of them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28fada9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reduce relative tolerance for integration\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setRelativeTolerance(\n",
    "    1 / 100 * amici_solver.getRelativeTolerance()\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",
    "\n",
    "rdata = res[RDATAS][0]\n",
    "print(\n",
    "    f\"preeq_status={list(map(amici.SteadyStateStatus, rdata.preeq_status.flatten()))}\"\n",
    ")\n",
    "print(f\"{rdata.preeq_numsteps=}\")\n",
    "\n",
    "# hard to reproduce on GHA\n",
    "if os.getenv(\"GITHUB_ACTIONS\") is None:\n",
    "    assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6d34467",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase relative steady state tolerance\n",
    "for log10_relaxation_factor in range(1, 10):\n",
    "    print(f\"Relaxing tolerances by factor {10 ** log10_relaxation_factor}\")\n",
    "    amici_solver = amici_model.getSolver()\n",
    "    amici_solver.setRelativeToleranceSteadyState(\n",
    "        amici_solver.getRelativeToleranceSteadyState()\n",
    "        * 10**log10_relaxation_factor\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",
    "    if all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS]):\n",
    "        print(\n",
    "            f\"-> Succeeded with relative steady state tolerance {amici_solver.getRelativeToleranceSteadyState()}\\n\"\n",
    "        )\n",
    "        break\n",
    "    else:\n",
    "        print(\"-> Failed.\\n\")\n",
    "\n",
    "print(\n",
    "    \"status:\",\n",
    "    [amici.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",
    ")\n",
    "print(f\"{rdata.preeq_numsteps=}\")\n",
    "assert all(rdata.status == amici.AMICI_SUCCESS for rdata in res[RDATAS])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d1b1835",
   "metadata": {},
   "source": [
    "That fixed the error, and took only a quarter of the number steps as the previous run, but at the cost of much lower accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df1ee3fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's try increasing the number of Newton steps\n",
    "# (this is 0 by default, so the Newton solver wasn't used before,\n",
    "#  as can be seen from the 0 in `rdata.preeq_numsteps[0]`)\n",
    "amici_solver = amici_model.getSolver()\n",
    "amici_solver.setNewtonMaxSteps(10**4)\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",
    "\n",
    "rdata = res[RDATAS][0]\n",
    "print(\n",
    "    f\"preeq_status={list(map(amici.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",
    "    ] == [\n",
    "        \"AMICI_ERROR\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "        \"AMICI_NOT_RUN\",\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c746982",
   "metadata": {},
   "source": [
    "Increasing the maximum number of Newton steps doesn't seem to help here. The Jacobian was numerically singular and its factorization failed. This can be a result of conserved quantities in the model. Section [Steady state sensitivity computation failed due to unsuccessful factorization of RHS Jacobian](#unsuccessful_factorization) shows how to address that."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "335.6px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}