From e0fbeec982c8a3389a3c7c300f4524f38e09f5d2 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 21:45:33 +0100 Subject: [PATCH 01/12] Fix piecewise/Heavisides handling * Substitute non-time-dependent Heavisides * Substitute inside the expanded expression, otherwise not all substitution targets are found Closes #2231 --- python/sdist/amici/de_export.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index ac857b366c..2943b000d3 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2757,19 +2757,22 @@ def _process_heavisides( for tmp_old in tmp_roots_old: # we want unique identifiers for the roots tmp_new = self._get_unique_root(tmp_old, roots) + heavisides.append((sp.Heaviside(tmp_old), tmp_new)) + # `tmp_new` is None if the root is not time-dependent. if tmp_new is None: continue # For Heavisides, we need to add the negative function as well self._get_unique_root(sp.sympify(-tmp_old), roots) - heavisides.append((sp.Heaviside(tmp_old), tmp_new)) if heavisides: # only apply subs if necessary for heaviside_sympy, heaviside_amici in heavisides: - dxdt = dxdt.subs(heaviside_sympy, heaviside_amici) + dt_expanded = dt_expanded.subs( + heaviside_sympy, heaviside_amici + ) - return dxdt + return dt_expanded class DEExporter: From fd227ecbf728626857f588cf6edfd45194980524 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 21:50:59 +0100 Subject: [PATCH 02/12] perf: Handle unique roots only --- python/sdist/amici/de_export.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 2943b000d3..a4328a610e 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2754,7 +2754,7 @@ def _process_heavisides( heavisides = [] # run through the expression tree and get the roots tmp_roots_old = self._collect_heaviside_roots(dt_expanded.args) - for tmp_old in tmp_roots_old: + for tmp_old in set(tmp_roots_old): # we want unique identifiers for the roots tmp_new = self._get_unique_root(tmp_old, roots) heavisides.append((sp.Heaviside(tmp_old), tmp_new)) From fcc36ccd359be554ec574b9c2ffdd2ced7b5d25f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 21:56:52 +0100 Subject: [PATCH 03/12] .. --- python/sdist/amici/de_export.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index a4328a610e..f9d553f06c 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2757,11 +2757,12 @@ def _process_heavisides( for tmp_old in set(tmp_roots_old): # we want unique identifiers for the roots tmp_new = self._get_unique_root(tmp_old, roots) - heavisides.append((sp.Heaviside(tmp_old), tmp_new)) - # `tmp_new` is None if the root is not time-dependent. if tmp_new is None: continue + + heavisides.append((sp.Heaviside(tmp_old), tmp_new)) + # For Heavisides, we need to add the negative function as well self._get_unique_root(sp.sympify(-tmp_old), roots) From bbfe5f79fdd63a30baf2dbf888ac4a313fdd47bb Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 22:06:42 +0100 Subject: [PATCH 04/12] .. --- python/sdist/amici/de_export.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index f9d553f06c..ffcd0f3993 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2760,11 +2760,9 @@ def _process_heavisides( # `tmp_new` is None if the root is not time-dependent. if tmp_new is None: continue - - heavisides.append((sp.Heaviside(tmp_old), tmp_new)) - # For Heavisides, we need to add the negative function as well self._get_unique_root(sp.sympify(-tmp_old), roots) + heavisides.append((sp.Heaviside(tmp_old), tmp_new)) if heavisides: # only apply subs if necessary From fd3f076b94d4f772f5f52704aea2e59d757107a9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 22:40:11 +0100 Subject: [PATCH 05/12] Remove broken error example --- python/examples/example_errors.ipynb | 245 --------------------------- 1 file changed, 245 deletions(-) diff --git a/python/examples/example_errors.ipynb b/python/examples/example_errors.ipynb index 2b35964d8b..01e477b67d 100644 --- a/python/examples/example_errors.ipynb +++ b/python/examples/example_errors.ipynb @@ -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": { From 31518e2b0f224a856358dae58fc63817d1807d52 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 12 Dec 2023 23:07:40 +0100 Subject: [PATCH 06/12] deterministic unique --- python/sdist/amici/de_export.py | 4 ++-- python/sdist/amici/import_utils.py | 15 +++++++++++++++ python/sdist/amici/sbml_import.py | 2 +- 3 files changed, 18 insertions(+), 3 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index ffcd0f3993..3770068a5e 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -58,8 +58,8 @@ generate_flux_symbol, smart_subs_dict, strip_pysb, - symbol_with_assumptions, toposort_symbols, + unique_preserve_order, ) from .logging import get_logger, log_execution_time, set_log_level @@ -2754,7 +2754,7 @@ def _process_heavisides( heavisides = [] # run through the expression tree and get the roots tmp_roots_old = self._collect_heaviside_roots(dt_expanded.args) - for tmp_old in set(tmp_roots_old): + for tmp_old in unique_preserve_order(tmp_roots_old): # we want unique identifiers for the roots tmp_new = self._get_unique_root(tmp_old, roots) # `tmp_new` is None if the root is not time-dependent. diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 77a2add60b..2dfc60d0c5 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -734,5 +734,20 @@ def strip_pysb(symbol: sp.Basic) -> sp.Basic: return symbol +def unique_preserve_order(seq: Sequence) -> list: + """Return a list of unique elements in Sequence, keeping only the first + occurrence of each element + + Parameters: + seq: Sequence to prune + + Returns: + List of unique elements in ``seq`` + """ + seen = set() + seen_add = seen.add + return [x for x in seq if not (x in seen or seen_add(x))] + + sbml_time_symbol = symbol_with_assumptions("time") amici_time_symbol = symbol_with_assumptions("t") diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index dd24b98cf8..b6a7b37ed4 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -38,7 +38,6 @@ DEModel, _default_simplify, smart_is_zero_matrix, - symbol_with_assumptions, ) from .import_utils import ( RESERVED_SYMBOLS, @@ -54,6 +53,7 @@ sbml_time_symbol, smart_subs, smart_subs_dict, + symbol_with_assumptions, toposort_symbols, ) from .logging import get_logger, log_execution_time, set_log_level From 1f7dc9a18a68a836dd32e45fc931b930583c5fe2 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 09:25:12 +0100 Subject: [PATCH 07/12] seed++? --- tests/benchmark-models/test_petab_benchmark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) mode change 100755 => 100644 tests/benchmark-models/test_petab_benchmark.py diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py old mode 100755 new mode 100644 index d0a783e2c4..08008d0660 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -120,7 +120,7 @@ def test_benchmark_gradient(model, scale): noise_level = 0.1 - np.random.seed(0) + np.random.seed(1) if scale: point = np.asarray( list( From f76678e0824d66285e635aa277e7af1114194c1c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 10:11:39 +0100 Subject: [PATCH 08/12] seed++ --- tests/benchmark-models/test_petab_benchmark.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 08008d0660..b18de3830e 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -120,7 +120,7 @@ def test_benchmark_gradient(model, scale): noise_level = 0.1 - np.random.seed(1) + np.random.seed(2) if scale: point = np.asarray( list( From c854930de560fb6b77a878d4375ac88793402f79 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 10:52:00 +0100 Subject: [PATCH 09/12] skip --- tests/benchmark-models/test_petab_benchmark.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index b18de3830e..92492738c9 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -38,6 +38,8 @@ # excluded due to excessive numerical failures "Crauste_CellSystems2017", "Fujita_SciSignal2010", + # fails for unclear reasons on GHA + "Giordano_Nature2020", ) ] @@ -120,7 +122,7 @@ def test_benchmark_gradient(model, scale): noise_level = 0.1 - np.random.seed(2) + np.random.seed(0) if scale: point = np.asarray( list( From 2771ba0b1adbac377cb1a8c1a98dfae7f9ef8f7e Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 15:59:58 +0100 Subject: [PATCH 10/12] don't expand --- python/sdist/amici/de_export.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index 3770068a5e..2694f753ad 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -2745,15 +2745,11 @@ def _process_heavisides( :returns: dxdt with Heaviside functions replaced by amici helper variables """ - - # expanding the rhs will in general help to collect the same - # heaviside function - dt_expanded = dxdt.expand() # track all the old Heaviside expressions in tmp_roots_old # replace them later by the new expressions heavisides = [] # run through the expression tree and get the roots - tmp_roots_old = self._collect_heaviside_roots(dt_expanded.args) + tmp_roots_old = self._collect_heaviside_roots(dxdt.args) for tmp_old in unique_preserve_order(tmp_roots_old): # we want unique identifiers for the roots tmp_new = self._get_unique_root(tmp_old, roots) @@ -2767,11 +2763,9 @@ def _process_heavisides( if heavisides: # only apply subs if necessary for heaviside_sympy, heaviside_amici in heavisides: - dt_expanded = dt_expanded.subs( - heaviside_sympy, heaviside_amici - ) + dxdt = dxdt.subs(heaviside_sympy, heaviside_amici) - return dt_expanded + return dxdt class DEExporter: From bcbea7e6be5d1655529029d1995ff5ccf41666e9 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 16:00:57 +0100 Subject: [PATCH 11/12] try unskip Giordano_Nature2020 --- tests/benchmark-models/test_petab_benchmark.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index d3ea3e62a1..753c88e500 100644 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -38,8 +38,6 @@ # excluded due to excessive numerical failures "Crauste_CellSystems2017", "Fujita_SciSignal2010", - # fails for unclear reasons on GHA - "Giordano_Nature2020", ) ] From 4e727f6f4cd020b20073d3ad468011de1f3c75de Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 13 Dec 2023 16:01:24 +0100 Subject: [PATCH 12/12] Revert "Remove broken error example" This reverts commit fd3f076b94d4f772f5f52704aea2e59d757107a9. --- python/examples/example_errors.ipynb | 245 +++++++++++++++++++++++++++ 1 file changed, 245 insertions(+) diff --git a/python/examples/example_errors.ipynb b/python/examples/example_errors.ipynb index 01e477b67d..2b35964d8b 100644 --- a/python/examples/example_errors.ipynb +++ b/python/examples/example_errors.ipynb @@ -840,6 +840,251 @@ ")\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": {