Skip to content

Commit

Permalink
Remove broken error example
Browse files Browse the repository at this point in the history
  • Loading branch information
@dweindl
dweindl committed Dec 12, 2023
1 parent bbfe5f7 commit fd3f076
Showing 1 changed file with 0 additions and 245 deletions.
245 changes: 0 additions & 245 deletions python/examples/example_errors.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -840,251 +840,6 @@
")\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": {
Expand Down

0 comments on commit fd3f076

Please sign in to comment.