From f7cbca4784486685b7a3020fb737304118d905d0 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 30 Aug 2023 19:41:40 +0200 Subject: [PATCH] Fix pre-commit configuration Relevant changes are only in .pre-commit-config.yaml * run black on notebooks via pre-commit * set proper line width (pyproject.toml is not used if we are running things from the repo root) * compatible line width for black + isort * fix json error in binder/overview.ipynb * re-blacken everything --- .pre-commit-config.yaml | 8 +- binder/overview.ipynb | 3 +- documentation/ExampleJax.ipynb | 16 +- documentation/conf.py | 16 +- documentation/recreate_reference_list.py | 8 +- python/benchmark/benchmark_pysb.py | 14 +- .../ExampleEquilibrationLogic.ipynb | 68 +++-- python/examples/example_errors.ipynb | 109 +++++--- python/examples/example_jax/ExampleJax.ipynb | 22 +- .../example_performance_optimization.ipynb | 9 +- python/examples/example_petab/petab.ipynb | 5 +- .../createModelPresimulation.py | 3 +- .../example_splines/ExampleSplines.ipynb | 39 ++- .../ExampleSplinesSwameye2003.ipynb | 121 ++++++--- .../ExampleSteadystate.ipynb | 20 +- python/sdist/amici/__init__.py | 10 +- python/sdist/amici/__main__.py | 4 +- .../amici/conserved_quantities_demartino.py | 55 +++- python/sdist/amici/custom_commands.py | 3 +- python/sdist/amici/cxxcodeprinter.py | 23 +- python/sdist/amici/de_export.py | 257 +++++++++++++----- python/sdist/amici/de_model.py | 26 +- python/sdist/amici/gradient_check.py | 7 +- python/sdist/amici/import_utils.py | 42 ++- python/sdist/amici/logging.py | 6 +- python/sdist/amici/numpy.py | 26 +- python/sdist/amici/pandas.py | 47 +++- python/sdist/amici/parameter_mapping.py | 37 ++- python/sdist/amici/petab_import.py | 51 +++- python/sdist/amici/petab_import_pysb.py | 23 +- python/sdist/amici/petab_objective.py | 123 ++++++--- python/sdist/amici/petab_simulate.py | 5 +- python/sdist/amici/petab_util.py | 12 +- python/sdist/amici/pysb_import.py | 140 +++++++--- python/sdist/amici/sbml_import.py | 215 +++++++++++---- python/sdist/amici/sbml_utils.py | 8 +- python/sdist/amici/splines.py | 156 ++++++++--- python/sdist/amici/swig.py | 16 +- python/sdist/amici/swig_wrappers.py | 10 +- python/sdist/setup.py | 4 +- python/tests/conftest.py | 4 +- .../bngwiki_egfr_simple_deletemolecules.py | 4 +- python/tests/splines_utils.py | 38 ++- .../test_compare_conservation_laws_sbml.py | 20 +- .../test_conserved_quantities_demartino.py | 27 +- python/tests/test_edata.py | 4 +- python/tests/test_events.py | 37 ++- python/tests/test_hdf5.py | 4 +- python/tests/test_heavisides.py | 20 +- python/tests/test_misc.py | 4 +- python/tests/test_observable_events.py | 4 +- python/tests/test_pandas.py | 8 +- python/tests/test_parameter_mapping.py | 13 +- python/tests/test_petab_import.py | 8 +- python/tests/test_petab_objective.py | 4 +- python/tests/test_petab_simulate.py | 4 +- python/tests/test_preequilibration.py | 35 ++- python/tests/test_pregenerated_models.py | 21 +- python/tests/test_pysb.py | 16 +- python/tests/test_rdata.py | 4 +- python/tests/test_sbml_import.py | 39 ++- .../test_sbml_import_special_functions.py | 12 +- python/tests/test_splines.py | 4 +- python/tests/test_splines_python.py | 12 +- python/tests/test_splines_short.py | 4 +- python/tests/test_swig_interface.py | 23 +- python/tests/util.py | 20 +- tests/benchmark-models/evaluate_benchmark.py | 7 +- .../benchmark-models/test_petab_benchmark.py | 12 +- tests/benchmark-models/test_petab_model.py | 14 +- tests/conftest.py | 21 +- tests/generateTestConfig/example.py | 4 +- tests/generateTestConfig/example_events.py | 8 +- tests/generateTestConfig/example_jakstat.py | 16 +- .../example_nested_events.py | 8 +- tests/generateTestConfig/example_neuron.py | 4 +- tests/performance/test.py | 8 +- tests/petab_test_suite/conftest.py | 13 +- tests/petab_test_suite/test_petab_suite.py | 26 +- tests/testSBMLSuite.py | 37 ++- 80 files changed, 1712 insertions(+), 626 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f73aa39939..521ff54f85 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,7 +6,7 @@ repos: hooks: - id: isort name: isort (python) - args: ["--profile", "black", "--filter-files"] + args: ["--profile", "black", "--filter-files", "--line-length", "79"] - repo: https://github.com/pre-commit/pre-commit-hooks rev: v4.4.0 hooks: @@ -17,16 +17,14 @@ repos: - id: end-of-file-fixer - id: trailing-whitespace - repo: https://github.com/psf/black - rev: 23.3.0 + rev: 23.7.0 hooks: - - id: black + - id: black-jupyter # It is recommended to specify the latest version of Python # supported by your project here, or alternatively use # pre-commit's default_language_version, see # https://pre-commit.com/#top_level-default_language_version language_version: python3.11 args: ["--line-length", "79"] - - id: black-jupyter - language_version: python3.11 exclude: '^(ThirdParty|models)/' diff --git a/binder/overview.ipynb b/binder/overview.ipynb index 9c4959f372..6d98f0f9fb 100644 --- a/binder/overview.ipynb +++ b/binder/overview.ipynb @@ -34,7 +34,8 @@ "\n", "* [Interfacing JAX](../python/examples/example_jax/ExampleJax.ipynb)\n", "\n", - " Provides guidance on how to combine AMICI with differential programming frameworks such as JAX.\n" + " Provides guidance on how to combine AMICI with differential programming frameworks such as JAX.\n", + "\n", "* [Efficient spline interpolation](../python/examples/example_splines/ExampleSplines.ipynb)\n", "\n", " Shows how to add annotated spline formulas to existing SBML models in order to speed up AMICI's model import.\n", diff --git a/documentation/ExampleJax.ipynb b/documentation/ExampleJax.ipynb index af4267b59d..1899305b67 100644 --- a/documentation/ExampleJax.ipynb +++ b/documentation/ExampleJax.ipynb @@ -645,7 +645,9 @@ "\n", "def amici_hcb_sllh(parameters: jnp.array):\n", " sllh = amici_hcb_base(parameters)[\"sllh\"]\n", - " return jnp.asarray(tuple(sllh[par_id] for par_id in petab_problem.x_free_ids))" + " return jnp.asarray(\n", + " tuple(sllh[par_id] for par_id in petab_problem.x_free_ids)\n", + " )" ] }, { @@ -769,7 +771,9 @@ "metadata": {}, "outputs": [], "source": [ - "llh_jax, sllh_jax = jax_objective_with_parameter_transform(petab_problem.x_nominal_free)" + "llh_jax, sllh_jax = jax_objective_with_parameter_transform(\n", + " petab_problem.x_nominal_free\n", + ")" ] }, { @@ -791,7 +795,9 @@ "# TODO remove me as soon as sllh in simulate_petab is fixed\n", "sllh = {\n", " name: value / (np.log(10) * par_value)\n", - " for (name, value), par_value in zip(r[\"sllh\"].items(), petab_problem.x_nominal_free)\n", + " for (name, value), par_value in zip(\n", + " r[\"sllh\"].items(), petab_problem.x_nominal_free\n", + " )\n", "}" ] }, @@ -942,7 +948,9 @@ "outputs": [], "source": [ "jax.config.update(\"jax_enable_x64\", True)\n", - "llh_jax, sllh_jax = jax_objective_with_parameter_transform(petab_problem.x_nominal_free)" + "llh_jax, sllh_jax = jax_objective_with_parameter_transform(\n", + " petab_problem.x_nominal_free\n", + ")" ] }, { diff --git a/documentation/conf.py b/documentation/conf.py index 5b55dd82b9..ba88b25a8d 100644 --- a/documentation/conf.py +++ b/documentation/conf.py @@ -52,7 +52,9 @@ def my_exhale_generate_doxygen(doxygen_input): DomainDirectiveFactory as breathe_DomainDirectiveFactory, ) -old_breathe_DomainDirectiveFactory_create = breathe_DomainDirectiveFactory.create +old_breathe_DomainDirectiveFactory_create = ( + breathe_DomainDirectiveFactory.create +) def my_breathe_DomainDirectiveFactory_create(domain: str, args): @@ -67,7 +69,9 @@ def my_breathe_DomainDirectiveFactory_create(domain: str, args): return cls(domain + ":" + name, *args[1:]) -breathe_DomainDirectiveFactory.create = my_breathe_DomainDirectiveFactory_create +breathe_DomainDirectiveFactory.create = ( + my_breathe_DomainDirectiveFactory_create +) # END Monkeypatch breathe @@ -102,7 +106,9 @@ def install_doxygen(): subprocess.run(cmd, shell=True, check=True) assert os.path.islink(os.path.join(some_dir_on_path, "doxygen")) # verify it's available - res = subprocess.run(["doxygen", "--version"], check=False, capture_output=True) + res = subprocess.run( + ["doxygen", "--version"], check=False, capture_output=True + ) print(res.stdout.decode(), res.stderr.decode()) assert version in res.stdout.decode() @@ -294,7 +300,9 @@ def install_doxygen(): "verboseBuild": True, } -mtocpp_filter = os.path.join(amici_dir, "matlab", "mtoc", "config", "mtocpp_filter.sh") +mtocpp_filter = os.path.join( + amici_dir, "matlab", "mtoc", "config", "mtocpp_filter.sh" +) exhale_projects_args = { "AMICI_CPP": { "exhaleDoxygenStdin": "\n".join( diff --git a/documentation/recreate_reference_list.py b/documentation/recreate_reference_list.py index 1dd1c13b4b..034c884c4b 100755 --- a/documentation/recreate_reference_list.py +++ b/documentation/recreate_reference_list.py @@ -42,7 +42,10 @@ def get_sub_bibliography(year, by_year, bibfile): entries = ",".join(["@" + x for x in by_year[year]]) stdin_input = ( - "---\n" f"bibliography: {bibfile}\n" f'nocite: "{entries}"\n...\n' f"# {year}" + "---\n" + f"bibliography: {bibfile}\n" + f'nocite: "{entries}"\n...\n' + f"# {year}" ) out = subprocess.run( @@ -67,7 +70,8 @@ def main(): with open(outfile, "w") as f: f.write("# References\n\n") f.write( - "List of publications using AMICI. " f"Total number is {num_total}.\n\n" + "List of publications using AMICI. " + f"Total number is {num_total}.\n\n" ) f.write( "If you applied AMICI in your work and your publication is " diff --git a/python/benchmark/benchmark_pysb.py b/python/benchmark/benchmark_pysb.py index 42c45d6274..079041ed4d 100644 --- a/python/benchmark/benchmark_pysb.py +++ b/python/benchmark/benchmark_pysb.py @@ -33,7 +33,9 @@ with amici.add_path(os.path.dirname(pysb.examples.__file__)): with amici.add_path( - os.path.join(os.path.dirname(__file__), "..", "tests", "pysb_test_models") + os.path.join( + os.path.dirname(__file__), "..", "tests", "pysb_test_models" + ) ): pysb.SelfExporter.cleanup() # reset pysb pysb.SelfExporter.do_export = True @@ -52,9 +54,9 @@ integrator_options={"rtol": rtol, "atol": atol}, ) time_pysb = ( - timeit.Timer("pysb_simres = sim.run()", globals={"sim": sim}).timeit( - number=N_REPEATS - ) + timeit.Timer( + "pysb_simres = sim.run()", globals={"sim": sim} + ).timeit(number=N_REPEATS) / N_REPEATS ) @@ -76,7 +78,9 @@ observables=list(pysb_model.observables.keys()), ) - amici_model_module = amici.import_model_module(pysb_model.name, outdir) + amici_model_module = amici.import_model_module( + pysb_model.name, outdir + ) model_pysb = amici_model_module.getModel() diff --git a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb index dcdb4ab04d..5f66ea4db9 100644 --- a/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb +++ b/python/examples/example_constant_species/ExampleEquilibrationLogic.ipynb @@ -61,7 +61,9 @@ "source": [ "from IPython.display import Image\n", "\n", - "fig = Image(filename=(\"../../../documentation/gfx/steadystate_solver_workflow.png\"))\n", + "fig = Image(\n", + " filename=(\"../../../documentation/gfx/steadystate_solver_workflow.png\")\n", + ")\n", "fig" ] }, @@ -252,7 +254,9 @@ ")\n", "model_reduced = model_reduced_module.getModel()\n", "\n", - "model_module = amici.import_model_module(model_name, os.path.abspath(model_output_dir))\n", + "model_module = amici.import_model_module(\n", + " model_name, os.path.abspath(model_output_dir)\n", + ")\n", "model = model_module.getModel()\n", "\n", "\n", @@ -431,7 +435,9 @@ "solver.setMaxSteps(100)\n", "rdata = amici.runAmiciSimulation(model, solver)\n", "print(\"Status of postequilibration:\", rdata[\"posteq_status\"])\n", - "print(\"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"])" + "print(\n", + " \"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"]\n", + ")" ] }, { @@ -470,7 +476,8 @@ "\n", "print(\"Status of postequilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", - " \"Number of steps employed in postequilibration:\", rdata_reduced[\"posteq_numsteps\"]\n", + " \"Number of steps employed in postequilibration:\",\n", + " rdata_reduced[\"posteq_numsteps\"],\n", ")" ] }, @@ -542,7 +549,9 @@ "rdata = amici.runAmiciSimulation(model, solver)\n", "\n", "print(\"Status of postequilibration:\", rdata[\"posteq_status\"])\n", - "print(\"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"])\n", + "print(\n", + " \"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"]\n", + ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" ] @@ -585,7 +594,9 @@ "source": [ "# Call simulation with singular Jacobian and newtonOnly mode (will fail)\n", "model.setTimepoints(np.full(1, np.inf))\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "model.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "solver = model.getSolver()\n", "solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", "solver.setSensitivityOrder(amici.SensitivityOrder.first)\n", @@ -593,7 +604,9 @@ "rdata = amici.runAmiciSimulation(model, solver)\n", "\n", "print(\"Status of postequilibration:\", rdata[\"posteq_status\"])\n", - "print(\"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"])\n", + "print(\n", + " \"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"]\n", + ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata[\"sx\"][0, :, :])" ] @@ -621,7 +634,9 @@ "source": [ "# Call postequilibration by setting an infinity timepoint\n", "model_reduced.setTimepoints(np.full(1, np.inf))\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "model.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "solver_reduced = model_reduced.getSolver()\n", "solver_reduced.setNewtonMaxSteps(10)\n", "solver_reduced.setSensitivityMethod(amici.SensitivityMethod.forward)\n", @@ -631,7 +646,8 @@ "\n", "print(\"Status of postequilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", - " \"Number of steps employed in postequilibration:\", rdata_reduced[\"posteq_numsteps\"]\n", + " \"Number of steps employed in postequilibration:\",\n", + " rdata_reduced[\"posteq_numsteps\"],\n", ")\n", "print(\"Computed state sensitivities:\")\n", "print(rdata_reduced[\"sx\"][0, :, :])" @@ -687,7 +703,9 @@ "edata.setObservedData([1.8] * 2)\n", "edata.fixedParameters = np.array([3.0, 5.0])\n", "\n", - "model_reduced.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "model_reduced.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "solver_reduced = model_reduced.getSolver()\n", "solver_reduced.setNewtonMaxSteps(10)\n", "solver_reduced.setSensitivityMethod(amici.SensitivityMethod.adjoint)\n", @@ -697,7 +715,8 @@ "\n", "print(\"Status of postequilibration:\", rdata_reduced[\"posteq_status\"])\n", "print(\n", - " \"Number of steps employed in postequilibration:\", rdata_reduced[\"posteq_numsteps\"]\n", + " \"Number of steps employed in postequilibration:\",\n", + " rdata_reduced[\"posteq_numsteps\"],\n", ")\n", "print(\n", " \"Number of backward steps employed in postequilibration:\",\n", @@ -733,7 +752,9 @@ ], "source": [ "# Call adjoint postequilibration with model with singular Jacobian\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "model.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "solver = model.getSolver()\n", "solver.setNewtonMaxSteps(10)\n", "solver.setSensitivityMethod(amici.SensitivityMethod.adjoint)\n", @@ -741,9 +762,12 @@ "rdata = amici.runAmiciSimulation(model, solver, edata)\n", "\n", "print(\"Status of postequilibration:\", rdata[\"posteq_status\"])\n", - "print(\"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"])\n", "print(\n", - " \"Number of backward steps employed in postequilibration:\", rdata[\"posteq_numstepsB\"]\n", + " \"Number of steps employed in postequilibration:\", rdata[\"posteq_numsteps\"]\n", + ")\n", + "print(\n", + " \"Number of backward steps employed in postequilibration:\",\n", + " rdata[\"posteq_numstepsB\"],\n", ")\n", "print(\"Computed gradient:\", rdata[\"sllh\"])" ] @@ -890,7 +914,9 @@ "edata.setTimepoints(np.array([0.0, 0.1, 1.0]))\n", "\n", "# create the solver object and run the simulation, singular Jacobian, enforce Newton solver for sensitivities\n", - "model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "model.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "solver = model.getSolver()\n", "solver.setNewtonMaxSteps(10)\n", "solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", @@ -1086,7 +1112,9 @@ "solver_reduced = model_reduced.getSolver()\n", "solver_reduced.setNewtonMaxSteps(10)\n", "solver_reduced.setSensitivityMethod(amici.SensitivityMethod.adjoint)\n", - "solver_reduced.setSensitivityMethodPreequilibration(amici.SensitivityMethod.adjoint)\n", + "solver_reduced.setSensitivityMethodPreequilibration(\n", + " amici.SensitivityMethod.adjoint\n", + ")\n", "solver_reduced.setSensitivityOrder(amici.SensitivityOrder.first)\n", "rdata_reduced = amici.runAmiciSimulation(model_reduced, solver_reduced, edata)\n", "\n", @@ -1193,14 +1221,18 @@ "solver_reduced.setAbsoluteToleranceSteadyState(1e-3)\n", "solver_reduced.setRelativeToleranceSteadyStateSensi(1e-2)\n", "solver_reduced.setAbsoluteToleranceSteadyStateSensi(1e-3)\n", - "rdata_reduced_lax = amici.runAmiciSimulation(model_reduced, solver_reduced, edata)\n", + "rdata_reduced_lax = amici.runAmiciSimulation(\n", + " model_reduced, solver_reduced, edata\n", + ")\n", "\n", "# run with strict tolerances\n", "solver_reduced.setRelativeToleranceSteadyState(1e-12)\n", "solver_reduced.setAbsoluteToleranceSteadyState(1e-16)\n", "solver_reduced.setRelativeToleranceSteadyStateSensi(1e-12)\n", "solver_reduced.setAbsoluteToleranceSteadyStateSensi(1e-16)\n", - "rdata_reduced_strict = amici.runAmiciSimulation(model_reduced, solver_reduced, edata)\n", + "rdata_reduced_strict = amici.runAmiciSimulation(\n", + " model_reduced, solver_reduced, edata\n", + ")\n", "\n", "# compare ODE outputs\n", "print(\"\\nODE solver steps, which were necessary to reach steady state:\")\n", diff --git a/python/examples/example_errors.ipynb b/python/examples/example_errors.ipynb index 3ff62df851..5e07803d96 100644 --- a/python/examples/example_errors.ipynb +++ b/python/examples/example_errors.ipynb @@ -75,7 +75,9 @@ "outputs": [], "source": [ "petab_problem = benchmark_models_petab.get_problem(\"Fujita_SciSignal2010\")\n", - "amici_model = import_petab_problem(petab_problem, verbose=False, force_compile=False)\n", + "amici_model = import_petab_problem(\n", + " petab_problem, verbose=False, force_compile=False\n", + ")\n", "\n", "np.random.seed(2991)\n", "problem_parameters = dict(\n", @@ -91,9 +93,12 @@ " scaled_parameters=True,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"Status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\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", @@ -144,7 +149,8 @@ ")\n", "\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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 succesfully.\")\n", @@ -164,7 +170,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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 succesfully.\")" @@ -204,11 +211,12 @@ " scaled_parameters=True,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"Status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", - " \"AMICI_TOO_MUCH_WORK\"\n", - "]" + "assert [\n", + " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + "] == [\"AMICI_TOO_MUCH_WORK\"]" ] }, { @@ -307,10 +315,13 @@ " scaled_parameters=True,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"Status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", - "assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\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", @@ -385,7 +396,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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 succesfully.\")" @@ -409,7 +421,9 @@ "outputs": [], "source": [ "petab_problem = benchmark_models_petab.get_problem(\"Weber_BMC2015\")\n", - "amici_model = import_petab_problem(petab_problem, verbose=False, force_compile=False)\n", + "amici_model = import_petab_problem(\n", + " petab_problem, verbose=False, force_compile=False\n", + ")\n", "\n", "np.random.seed(4)\n", "problem_parameters = dict(\n", @@ -425,9 +439,12 @@ " scaled_parameters=True,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"Status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", + "assert [\n", + " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + "] == [\n", " \"AMICI_ERROR\",\n", " \"AMICI_NOT_RUN\",\n", "]" @@ -478,7 +495,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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])" ] @@ -517,11 +535,12 @@ " scaled_parameters=True,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"Status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", - "assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", - " \"AMICI_FIRST_RHSFUNC_ERR\"\n", - "]" + "assert [\n", + " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + "] == [\"AMICI_FIRST_RHSFUNC_ERR\"]" ] }, { @@ -614,7 +633,9 @@ "unscaled_parameter = dict(\n", " zip(\n", " amici_model.getParameterIds(),\n", - " starmap(amici.getUnscaledParameter, zip(edata.parameters, edata.pscale)),\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\")))" @@ -684,7 +705,9 @@ "amici_solver = amici_model.getSolver()\n", "amici_solver.setSensitivityMethod(amici.SensitivityMethod.forward)\n", "amici_solver.setSensitivityOrder(amici.SensitivityOrder.first)\n", - "amici_model.setSteadyStateSensitivityMode(amici.SteadyStateSensitivityMode.newtonOnly)\n", + "amici_model.setSteadyStateSensitivityMode(\n", + " amici.SteadyStateSensitivityMode.newtonOnly\n", + ")\n", "\n", "np.random.seed(2020)\n", "problem_parameters = dict(\n", @@ -701,14 +724,15 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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 [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\n", - " \"AMICI_ERROR\"\n", - " ]" + " assert [\n", + " amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]\n", + " ] == [\"AMICI_ERROR\"]" ] }, { @@ -771,7 +795,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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])" ] @@ -809,7 +834,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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])" ] @@ -856,12 +882,15 @@ ")\n", "\n", "print(\n", - " \"Status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\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 [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\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", @@ -923,7 +952,9 @@ "source": [ "# Reduce relative tolerance for integration\n", "amici_solver = amici_model.getSolver()\n", - "amici_solver.setRelativeTolerance(1 / 100 * amici_solver.getRelativeTolerance())\n", + "amici_solver.setRelativeTolerance(\n", + " 1 / 100 * amici_solver.getRelativeTolerance()\n", + ")\n", "\n", "res = simulate_petab(\n", " petab_problem=petab_problem,\n", @@ -933,7 +964,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", @@ -959,7 +991,8 @@ " print(f\"Relaxing tolerances by factor {10 ** log10_relaxation_factor}\")\n", " amici_solver = amici_model.getSolver()\n", " amici_solver.setRelativeToleranceSteadyState(\n", - " amici_solver.getRelativeToleranceSteadyState() * 10**log10_relaxation_factor\n", + " amici_solver.getRelativeToleranceSteadyState()\n", + " * 10**log10_relaxation_factor\n", " )\n", "\n", " res = simulate_petab(\n", @@ -978,7 +1011,8 @@ " print(\"-> Failed.\\n\")\n", "\n", "print(\n", - " \"status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", @@ -1018,7 +1052,8 @@ " solver=amici_solver,\n", ")\n", "print(\n", - " \"status:\", [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]]\n", + " \"status:\",\n", + " [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]],\n", ")\n", "\n", "rdata = res[RDATAS][0]\n", @@ -1028,7 +1063,9 @@ "print(f\"{rdata.preeq_numsteps=}\")\n", "# hard to reproduce on GHA\n", "if os.getenv(\"GITHUB_ACTIONS\") is None:\n", - " assert [amici.simulation_status_to_str(rdata.status) for rdata in res[RDATAS]] == [\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", diff --git a/python/examples/example_jax/ExampleJax.ipynb b/python/examples/example_jax/ExampleJax.ipynb index 1a2f5da3d7..efda5b458e 100644 --- a/python/examples/example_jax/ExampleJax.ipynb +++ b/python/examples/example_jax/ExampleJax.ipynb @@ -264,7 +264,9 @@ "source": [ "from amici.petab_import import import_petab_problem\n", "\n", - "amici_model = import_petab_problem(petab_problem, force_compile=True, verbose=False)" + "amici_model = import_petab_problem(\n", + " petab_problem, force_compile=True, verbose=False\n", + ")" ] }, { @@ -329,7 +331,9 @@ "\n", "def amici_hcb_sllh(parameters: jnp.array):\n", " sllh = amici_hcb_base(parameters)[\"sllh\"]\n", - " return jnp.asarray(tuple(sllh[par_id] for par_id in petab_problem.x_free_ids))" + " return jnp.asarray(\n", + " tuple(sllh[par_id] for par_id in petab_problem.x_free_ids)\n", + " )" ] }, { @@ -459,7 +463,9 @@ "metadata": {}, "outputs": [], "source": [ - "llh_jax, sllh_jax = jax_objective_with_parameter_transform(scaled_parameters_np)" + "llh_jax, sllh_jax = jax_objective_with_parameter_transform(\n", + " scaled_parameters_np\n", + ")" ] }, { @@ -667,7 +673,8 @@ "grad_jax = np.asarray(sllh_jax)\n", "rel_diff = (grad_amici - grad_jax) / grad_jax\n", "pd.DataFrame(\n", - " index=r[\"sllh\"].keys(), data=dict(amici=grad_amici, jax=grad_jax, rel_diff=rel_diff)\n", + " index=r[\"sllh\"].keys(),\n", + " data=dict(amici=grad_amici, jax=grad_jax, rel_diff=rel_diff),\n", ")" ] }, @@ -688,7 +695,9 @@ "outputs": [], "source": [ "jax.config.update(\"jax_enable_x64\", True)\n", - "llh_jax, sllh_jax = jax_objective_with_parameter_transform(scaled_parameters_np)" + "llh_jax, sllh_jax = jax_objective_with_parameter_transform(\n", + " scaled_parameters_np\n", + ")" ] }, { @@ -877,7 +886,8 @@ "grad_jax = np.asarray(sllh_jax)\n", "rel_diff = (grad_amici - grad_jax) / grad_jax\n", "pd.DataFrame(\n", - " index=r[\"sllh\"].keys(), data=dict(amici=grad_amici, jax=grad_jax, rel_diff=rel_diff)\n", + " index=r[\"sllh\"].keys(),\n", + " data=dict(amici=grad_amici, jax=grad_jax, rel_diff=rel_diff),\n", ")" ] } diff --git a/python/examples/example_large_models/example_performance_optimization.ipynb b/python/examples/example_large_models/example_performance_optimization.ipynb index 46d4064b66..31a9fc1729 100644 --- a/python/examples/example_large_models/example_performance_optimization.ipynb +++ b/python/examples/example_large_models/example_performance_optimization.ipynb @@ -221,7 +221,9 @@ " plt.show()\n", "\n", " import_times = df.sort_values(\"nprocs\")[\"time\"].values\n", - " percent_change = (import_times[0] - min(import_times[1:])) / import_times[0] * 100\n", + " percent_change = (\n", + " (import_times[0] - min(import_times[1:])) / import_times[0] * 100\n", + " )\n", " if percent_change > 0:\n", " print(f\"Import time decreased by up to ~{percent_change:.0f}%.\")\n", " else:\n", @@ -283,7 +285,10 @@ "source": [ "figsize(8, 4)\n", "compilation_time_s = [3022.453, 289.518]\n", - "labels = [\"g++ (Ubuntu 12.2.0-3ubuntu1) 12.2.0\", \"Ubuntu clang version 15.0.2-1\"]\n", + "labels = [\n", + " \"g++ (Ubuntu 12.2.0-3ubuntu1) 12.2.0\",\n", + " \"Ubuntu clang version 15.0.2-1\",\n", + "]\n", "plt.bar(labels, compilation_time_s)\n", "plt.ylim(ymin=0)\n", "plt.title(\"Choice of compiler - FröhlichGer2022\")\n", diff --git a/python/examples/example_petab/petab.ipynb b/python/examples/example_petab/petab.ipynb index 83483232f6..689d793f56 100644 --- a/python/examples/example_petab/petab.ipynb +++ b/python/examples/example_petab/petab.ipynb @@ -339,7 +339,10 @@ " for x_id, x_val in zip(petab_problem.x_ids, petab_problem.x_nominal_scaled)\n", "}\n", "simulate_petab(\n", - " petab_problem, amici_model, problem_parameters=parameters, scaled_parameters=True\n", + " petab_problem,\n", + " amici_model,\n", + " problem_parameters=parameters,\n", + " scaled_parameters=True,\n", ")" ] }, diff --git a/python/examples/example_presimulation/createModelPresimulation.py b/python/examples/example_presimulation/createModelPresimulation.py index fed24d151f..1db454a6b5 100644 --- a/python/examples/example_presimulation/createModelPresimulation.py +++ b/python/examples/example_presimulation/createModelPresimulation.py @@ -49,7 +49,8 @@ Rule( "PROT_dephospho", - prot(phospho="p", drug=None, kin=None) >> prot(phospho="u", drug=None, kin=None), + prot(phospho="p", drug=None, kin=None) + >> prot(phospho="u", drug=None, kin=None), Parameter("kdephospho_prot", 0.1), ) diff --git a/python/examples/example_splines/ExampleSplines.ipynb b/python/examples/example_splines/ExampleSplines.ipynb index ef882155c4..d376ba91e5 100644 --- a/python/examples/example_splines/ExampleSplines.ipynb +++ b/python/examples/example_splines/ExampleSplines.ipynb @@ -64,7 +64,11 @@ " build_dir = os.path.join(BUILD_PATH, model_name)\n", " rmtree(build_dir, ignore_errors=True)\n", " return _simulate(\n", - " sbml_model, parameters, build_dir=build_dir, model_name=model_name, **kwargs\n", + " sbml_model,\n", + " parameters,\n", + " build_dir=build_dir,\n", + " model_name=model_name,\n", + " **kwargs\n", " )\n", "\n", "\n", @@ -724,7 +728,10 @@ " evaluate_at=amici.sbml_utils.amici_time_symbol,\n", " nodes=amici.splines.UniformGrid(0, 1, number_of_nodes=3),\n", " values_at_nodes=[-2, 1, -1],\n", - " extrapolate=(None, \"constant\"), # no extrapolation required on the left side\n", + " extrapolate=(\n", + " None,\n", + " \"constant\",\n", + " ), # no extrapolation required on the left side\n", ")" ] }, @@ -1144,7 +1151,15 @@ "outputs": [], "source": [ "nruns = 6 # number of replicates\n", - "num_nodes = [5, 10, 15, 20, 25, 30, 40] # benchmark model import for these node numbers\n", + "num_nodes = [\n", + " 5,\n", + " 10,\n", + " 15,\n", + " 20,\n", + " 25,\n", + " 30,\n", + " 40,\n", + "] # benchmark model import for these node numbers\n", "amici_only_nodes = [\n", " 50,\n", " 75,\n", @@ -1203,13 +1218,15 @@ " if n in num_nodes:\n", " with tempfile.TemporaryDirectory() as tmpdir:\n", " t0 = time.perf_counter_ns()\n", - " amici.SbmlImporter(sbml_model, discard_annotations=True).sbml2amici(\n", - " \"benchmark\", tmpdir\n", - " )\n", + " amici.SbmlImporter(\n", + " sbml_model, discard_annotations=True\n", + " ).sbml2amici(\"benchmark\", tmpdir)\n", " dt = time.perf_counter_ns() - t0\n", " timings_piecewise.append(dt / 1e9)\n", " # Append benchmark data to dataframe\n", - " df_amici = pd.DataFrame(dict(num_nodes=n, time=timings_amici, use_annotations=True))\n", + " df_amici = pd.DataFrame(\n", + " dict(num_nodes=n, time=timings_amici, use_annotations=True)\n", + " )\n", " df_piecewise = pd.DataFrame(\n", " dict(num_nodes=n, time=timings_piecewise, use_annotations=False)\n", " )\n", @@ -1219,7 +1236,9 @@ " )\n", " else:\n", " df = pd.concat(\n", - " [df, df_amici, df_piecewise], ignore_index=True, verify_integrity=True\n", + " [df, df_amici, df_piecewise],\n", + " ignore_index=True,\n", + " verify_integrity=True,\n", " )" ] }, @@ -1262,7 +1281,9 @@ "ax.set_ylabel(\"model import time (s)\")\n", "ax.set_xlabel(\"number of spline nodes\")\n", "ax.set_yscale(\"log\")\n", - "ax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0f}\"))\n", + "ax.yaxis.set_major_formatter(\n", + " mpl.ticker.FuncFormatter(lambda x, pos: f\"{x:.0f}\")\n", + ")\n", "ax.xaxis.set_ticks(\n", " [\n", " 10,\n", diff --git a/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb b/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb index 05e9bd8e4b..8e3ee6db10 100644 --- a/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb +++ b/python/examples/example_splines_swameye/ExampleSplinesSwameye2003.ipynb @@ -217,7 +217,9 @@ " os.path.join(\"Swameye_PNAS2003\", \"swameye2003_model.xml\")\n", ")\n", "sbml_model = sbml_doc.getModel()\n", - "spline.add_to_sbml_model(sbml_model, auto_add=True, y_nominal=0.1, y_constant=True)" + "spline.add_to_sbml_model(\n", + " sbml_model, auto_add=True, y_nominal=0.1, y_constant=True\n", + ")" ] }, { @@ -402,7 +404,9 @@ "source": [ "# Load existing results if available\n", "if os.path.exists(f\"{name}.h5\"):\n", - " pypesto_result = pypesto.store.read_result(f\"{name}.h5\", problem=pypesto_problem)\n", + " pypesto_result = pypesto.store.read_result(\n", + " f\"{name}.h5\", problem=pypesto_problem\n", + " )\n", "else:\n", " pypesto_result = None\n", "# Overwrite\n", @@ -531,19 +535,25 @@ "\n", "def simulate_pEpoR(x=None, **kwargs):\n", " problem, rdata = _simulate(x, **kwargs)\n", - " assert problem.objective.amici_model.getObservableIds()[0].startswith(\"pEpoR\")\n", + " assert problem.objective.amici_model.getObservableIds()[0].startswith(\n", + " \"pEpoR\"\n", + " )\n", " return rdata[\"t\"], rdata[\"y\"][:, 0]\n", "\n", "\n", "def simulate_pSTAT5(x=None, **kwargs):\n", " problem, rdata = _simulate(x, **kwargs)\n", - " assert problem.objective.amici_model.getObservableIds()[1].startswith(\"pSTAT5\")\n", + " assert problem.objective.amici_model.getObservableIds()[1].startswith(\n", + " \"pSTAT5\"\n", + " )\n", " return rdata[\"t\"], rdata[\"y\"][:, 1]\n", "\n", "\n", "def simulate_tSTAT5(x=None, **kwargs):\n", " problem, rdata = _simulate(x, **kwargs)\n", - " assert problem.objective.amici_model.getObservableIds()[-1].startswith(\"tSTAT5\")\n", + " assert problem.objective.amici_model.getObservableIds()[-1].startswith(\n", + " \"tSTAT5\"\n", + " )\n", " return rdata[\"t\"], rdata[\"y\"][:, -1]\n", "\n", "\n", @@ -551,9 +561,15 @@ "df_measurements = petab.measurements.get_measurement_df(\n", " os.path.join(\"Swameye_PNAS2003\", \"swameye2003_measurements.tsv\")\n", ")\n", - "df_pEpoR = df_measurements[df_measurements[\"observableId\"].str.startswith(\"pEpoR\")]\n", - "df_pSTAT5 = df_measurements[df_measurements[\"observableId\"].str.startswith(\"pSTAT5\")]\n", - "df_tSTAT5 = df_measurements[df_measurements[\"observableId\"].str.startswith(\"tSTAT5\")]" + "df_pEpoR = df_measurements[\n", + " df_measurements[\"observableId\"].str.startswith(\"pEpoR\")\n", + "]\n", + "df_pSTAT5 = df_measurements[\n", + " df_measurements[\"observableId\"].str.startswith(\"pSTAT5\")\n", + "]\n", + "df_tSTAT5 = df_measurements[\n", + " df_measurements[\"observableId\"].str.startswith(\"tSTAT5\")\n", + "]" ] }, { @@ -766,7 +782,9 @@ "source": [ "# Create spline for pEpoR\n", "nodes = [0, 2.5, 5.0, 7.5, 10.0, 12.5, 15.0, 17.5, 20, 25, 30, 35, 40, 50, 60]\n", - "values_at_nodes = [sp.Symbol(f\"pEpoR_t{str(t).replace('.', '_dot_')}\") for t in nodes]\n", + "values_at_nodes = [\n", + " sp.Symbol(f\"pEpoR_t{str(t).replace('.', '_dot_')}\") for t in nodes\n", + "]\n", "spline = amici.splines.CubicHermiteSpline(\n", " sbml_id=\"pEpoR\",\n", " evaluate_at=amici.sbml_utils.amici_time_symbol,\n", @@ -858,7 +876,9 @@ " os.path.join(\"Swameye_PNAS2003\", \"swameye2003_model.xml\")\n", ")\n", "sbml_model = sbml_doc.getModel()\n", - "spline.add_to_sbml_model(sbml_model, auto_add=True, y_nominal=0.1, y_constant=True)" + "spline.add_to_sbml_model(\n", + " sbml_model, auto_add=True, y_nominal=0.1, y_constant=True\n", + ")" ] }, { @@ -1033,7 +1053,9 @@ " name = f\"Swameye_PNAS2003_15nodes_FD_reg{regstrength}\"\n", " pypesto_problem.fix_parameters(\n", " pypesto_problem.x_names.index(\"regularization_strength\"),\n", - " np.log10(regstrength), # parameter is specified as log10 scale in PEtab\n", + " np.log10(\n", + " regstrength\n", + " ), # parameter is specified as log10 scale in PEtab\n", " )\n", " regproblem = copy.deepcopy(pypesto_problem)\n", "\n", @@ -1051,7 +1073,9 @@ " new_ids = [str(i) for i in range(n_starts)]\n", " else:\n", " last_id = max(int(i) for i in regresult.optimize_result.id)\n", - " new_ids = [str(i) for i in range(last_id + 1, last_id + n_starts + 1)]\n", + " new_ids = [\n", + " str(i) for i in range(last_id + 1, last_id + n_starts + 1)\n", + " ]\n", " regresult = pypesto.optimize.minimize(\n", " regproblem,\n", " n_starts=n_starts,\n", @@ -1103,7 +1127,9 @@ "stats = []\n", "for regstrength in regstrengths:\n", " t, pEpoR = simulate_pEpoR(\n", - " N=None, problem=regproblems[regstrength], result=regresults[regstrength]\n", + " N=None,\n", + " problem=regproblems[regstrength],\n", + " result=regresults[regstrength],\n", " )\n", " assert np.array_equal(df_pEpoR[\"time\"], t[:-1])\n", " pEpoR = pEpoR[:-1]\n", @@ -1138,8 +1164,12 @@ ], "source": [ "# Visualize the results of the multistarts for a chosen regularization strength\n", - "ax = pypesto.visualize.waterfall(regresults[chosen_regstrength], size=[6.5, 3.5])\n", - "ax.set_title(f\"Waterfall plot (regularization strength = {chosen_regstrength})\")\n", + "ax = pypesto.visualize.waterfall(\n", + " regresults[chosen_regstrength], size=[6.5, 3.5]\n", + ")\n", + "ax.set_title(\n", + " f\"Waterfall plot (regularization strength = {chosen_regstrength})\"\n", + ")\n", "ax.set_ylim(ax.get_ylim()[0], 100);" ] }, @@ -1171,7 +1201,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -1234,7 +1266,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -1291,7 +1325,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -1343,7 +1379,8 @@ "# Plot ML fit for pEpoR (single regularization strength with noise model)\n", "fig, ax = plt.subplots(figsize=(6.5, 3.5))\n", "t, pEpoR = simulate_pEpoR(\n", - " problem=regproblems[chosen_regstrength], result=regresults[chosen_regstrength]\n", + " problem=regproblems[chosen_regstrength],\n", + " result=regresults[chosen_regstrength],\n", ")\n", "sigma_pEpoR = 0.0274 + 0.1 * pEpoR\n", "ax.fill_between(\n", @@ -1448,9 +1485,12 @@ "source": [ "# Create spline for pEpoR\n", "nodes = [0, 5, 10, 20, 60]\n", - "values_at_nodes = [sp.Symbol(f\"pEpoR_t{str(t).replace('.', '_dot_')}\") for t in nodes]\n", + "values_at_nodes = [\n", + " sp.Symbol(f\"pEpoR_t{str(t).replace('.', '_dot_')}\") for t in nodes\n", + "]\n", "derivatives_at_nodes = [\n", - " sp.Symbol(f\"derivative_pEpoR_t{str(t).replace('.', '_dot_')}\") for t in nodes[:-1]\n", + " sp.Symbol(f\"derivative_pEpoR_t{str(t).replace('.', '_dot_')}\")\n", + " for t in nodes[:-1]\n", "]\n", "spline = amici.splines.CubicHermiteSpline(\n", " sbml_id=\"pEpoR\",\n", @@ -1535,7 +1575,9 @@ " os.path.join(\"Swameye_PNAS2003\", \"swameye2003_model.xml\")\n", ")\n", "sbml_model = sbml_doc.getModel()\n", - "spline.add_to_sbml_model(sbml_model, auto_add=True, y_nominal=0.1, y_constant=True)" + "spline.add_to_sbml_model(\n", + " sbml_model, auto_add=True, y_nominal=0.1, y_constant=True\n", + ")" ] }, { @@ -1734,7 +1776,9 @@ " name = f\"Swameye_PNAS2003_5nodes_reg{regstrength}\"\n", " pypesto_problem.fix_parameters(\n", " pypesto_problem.x_names.index(\"regularization_strength\"),\n", - " np.log10(regstrength), # parameter is specified as log10 scale in PEtab\n", + " np.log10(\n", + " regstrength\n", + " ), # parameter is specified as log10 scale in PEtab\n", " )\n", " regproblem = copy.deepcopy(pypesto_problem)\n", "\n", @@ -1752,7 +1796,9 @@ " new_ids = [str(i) for i in range(n_starts)]\n", " else:\n", " last_id = max(int(i) for i in regresult.optimize_result.id)\n", - " new_ids = [str(i) for i in range(last_id + 1, last_id + n_starts + 1)]\n", + " new_ids = [\n", + " str(i) for i in range(last_id + 1, last_id + n_starts + 1)\n", + " ]\n", " regresult = pypesto.optimize.minimize(\n", " regproblem,\n", " n_starts=n_starts,\n", @@ -1802,7 +1848,9 @@ "stats = []\n", "for regstrength in regstrengths:\n", " t, pEpoR = simulate_pEpoR(\n", - " N=None, problem=regproblems[regstrength], result=regresults[regstrength]\n", + " N=None,\n", + " problem=regproblems[regstrength],\n", + " result=regresults[regstrength],\n", " )\n", " assert np.array_equal(df_pEpoR[\"time\"], t[:-1])\n", " pEpoR = pEpoR[:-1]\n", @@ -1837,8 +1885,12 @@ ], "source": [ "# Visualize the results of the multistarts for a chosen regularization strength\n", - "ax = pypesto.visualize.waterfall(regresults[chosen_regstrength], size=[6.5, 3.5])\n", - "ax.set_title(f\"Waterfall plot (regularization strength = {chosen_regstrength})\")\n", + "ax = pypesto.visualize.waterfall(\n", + " regresults[chosen_regstrength], size=[6.5, 3.5]\n", + ")\n", + "ax.set_title(\n", + " f\"Waterfall plot (regularization strength = {chosen_regstrength})\"\n", + ")\n", "ax.set_ylim(ax.get_ylim()[0], 100);" ] }, @@ -1870,7 +1922,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -1932,7 +1986,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -1989,7 +2045,9 @@ " )\n", " if regstrength == chosen_regstrength:\n", " kwargs = dict(\n", - " color=\"black\", label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\", zorder=2\n", + " color=\"black\",\n", + " label=f\"$\\\\mathbf{{\\\\lambda = {regstrength}}}$\",\n", + " zorder=2,\n", " )\n", " else:\n", " kwargs = dict(label=f\"$\\\\lambda = {regstrength}$\", alpha=0.5)\n", @@ -2041,7 +2099,8 @@ "# Plot ML fit for pEpoR (single regularization strength with noise model)\n", "fig, ax = plt.subplots(figsize=(6.5, 3.5))\n", "t, pEpoR = simulate_pEpoR(\n", - " problem=regproblems[chosen_regstrength], result=regresults[chosen_regstrength]\n", + " problem=regproblems[chosen_regstrength],\n", + " result=regresults[chosen_regstrength],\n", ")\n", "sigma_pEpoR = 0.0274 + 0.1 * pEpoR\n", "ax.fill_between(\n", diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index bd50894ec3..502174fe15 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -1570,7 +1570,9 @@ "solver.setSensitivityMethod(\n", " amici.SensitivityMethod.forward\n", ") # forward sensitivity analysis\n", - "solver.setSensitivityOrder(amici.SensitivityOrder.first) # first-order sensitivities\n", + "solver.setSensitivityOrder(\n", + " amici.SensitivityOrder.first\n", + ") # first-order sensitivities\n", "\n", "rdata = amici.runAmiciSimulation(model, solver)\n", "\n", @@ -1758,7 +1760,9 @@ "for ip in range(model.np()):\n", " plist = [ip]\n", " p = p_orig.copy()\n", - " err_norm = check_grad(func, grad, p[plist], \"llh\", p, [ip], epsilon=epsilon)\n", + " err_norm = check_grad(\n", + " func, grad, p[plist], \"llh\", p, [ip], epsilon=epsilon\n", + " )\n", " print(\"sllh: p[%d]: |error|_2: %f\" % (ip, err_norm))\n", "\n", "print()\n", @@ -1779,7 +1783,9 @@ "for ip in range(model.np()):\n", " plist = [ip]\n", " p = p_orig.copy()\n", - " err_norm = check_grad(func, grad, p[plist], \"sigmay\", p, [ip], epsilon=epsilon)\n", + " err_norm = check_grad(\n", + " func, grad, p[plist], \"sigmay\", p, [ip], epsilon=epsilon\n", + " )\n", " print(\"ssigmay: p[%d]: |error|_2: %f\" % (ip, err_norm))" ] }, @@ -1798,7 +1804,9 @@ "solver.setSensitivityMethod(\n", " amici.SensitivityMethod.forward\n", ") # forward sensitivity analysis\n", - "solver.setSensitivityOrder(amici.SensitivityOrder.first) # first-order sensitivities\n", + "solver.setSensitivityOrder(\n", + " amici.SensitivityOrder.first\n", + ") # first-order sensitivities\n", "model.requireSensitivitiesForAllParameters()\n", "solver.setRelativeTolerance(1e-12)\n", "rdata = amici.runAmiciSimulation(model, solver, edata)\n", @@ -1824,7 +1832,9 @@ " for ip in range(4):\n", " fd_approx = fd(model.getParameters(), ip, eps, symbol=symbol)\n", "\n", - " axes[ip, 0].plot(edata.getTimepoints(), rdata[f\"s{symbol}\"][:, ip, :], \"r-\")\n", + " axes[ip, 0].plot(\n", + " edata.getTimepoints(), rdata[f\"s{symbol}\"][:, ip, :], \"r-\"\n", + " )\n", " axes[ip, 0].plot(edata.getTimepoints(), fd_approx, \"k--\")\n", " axes[ip, 0].set_ylabel(f\"sensitivity {symbol}\")\n", " axes[ip, 0].set_xlabel(\"time\")\n", diff --git a/python/sdist/amici/__init__.py b/python/sdist/amici/__init__.py index 7160acb475..bcb7387fbf 100644 --- a/python/sdist/amici/__init__.py +++ b/python/sdist/amici/__init__.py @@ -66,9 +66,9 @@ def _imported_from_setup() -> bool: # requires the AMICI extension during its installation, but seems # unlikely... frame_path = os.path.realpath(os.path.expanduser(frame.filename)) - if frame_path == os.path.join(package_root, "setup.py") or frame_path.endswith( - f"{sep}setuptools{sep}build_meta.py" - ): + if frame_path == os.path.join( + package_root, "setup.py" + ) or frame_path.endswith(f"{sep}setuptools{sep}build_meta.py"): return True return False @@ -203,6 +203,8 @@ def _get_default_argument(func: Callable, arg: str) -> Any: import inspect signature = inspect.signature(func) - if (default := signature.parameters[arg].default) is not inspect.Parameter.empty: + if ( + default := signature.parameters[arg].default + ) is not inspect.Parameter.empty: return default raise ValueError(f"No default value for argument {arg} of {func}.") diff --git a/python/sdist/amici/__main__.py b/python/sdist/amici/__main__.py index b8fbc77c0f..165f5d9516 100644 --- a/python/sdist/amici/__main__.py +++ b/python/sdist/amici/__main__.py @@ -21,7 +21,9 @@ def print_info(): if hdf5_enabled: features.append("HDF5") - print(f"AMICI ({sys.platform}) version {__version__} ({','.join(features)})") + print( + f"AMICI ({sys.platform}) version {__version__} ({','.join(features)})" + ) if __name__ == "__main__": diff --git a/python/sdist/amici/conserved_quantities_demartino.py b/python/sdist/amici/conserved_quantities_demartino.py index 03e3730b06..5b04fa1479 100644 --- a/python/sdist/amici/conserved_quantities_demartino.py +++ b/python/sdist/amici/conserved_quantities_demartino.py @@ -60,7 +60,9 @@ def compute_moiety_conservation_laws( if not done: # construct interaction matrix - J, J2, fields = _fill(stoichiometric_list, engaged_species, num_species) + J, J2, fields = _fill( + stoichiometric_list, engaged_species, num_species + ) # seed random number generator if rng_seed is not False: @@ -142,14 +144,23 @@ def log(*args, **kwargs): if not engaged_species_idxs: continue log( - f"Moiety number {i + 1} engages {len(engaged_species_idxs)} " "species:" + f"Moiety number {i + 1} engages {len(engaged_species_idxs)} " + "species:" ) - for species_idx, coefficient in zip(engaged_species_idxs, coefficients): - name = species_names[species_idx] if species_names else species_idx + for species_idx, coefficient in zip( + engaged_species_idxs, coefficients + ): + name = ( + species_names[species_idx] + if species_names + else species_idx + ) log(f"\t{name}\t{coefficient}") -def _qsort(k: int, km: int, order: MutableSequence[int], pivots: Sequence[int]) -> None: +def _qsort( + k: int, km: int, order: MutableSequence[int], pivots: Sequence[int] +) -> None: """Quicksort Recursive implementation of the quicksort algorithm @@ -233,7 +244,9 @@ def _kernel( matrix2[i].append(1) order: List[int] = list(range(num_species)) - pivots = [matrix[i][0] if len(matrix[i]) else _MAX for i in range(num_species)] + pivots = [ + matrix[i][0] if len(matrix[i]) else _MAX for i in range(num_species) + ] done = False while not done: @@ -253,7 +266,10 @@ def _kernel( for i in range(len(matrix[order[j + 1]])): min2 = min( min2, - abs(matrix2[order[j + 1]][0] / matrix2[order[j + 1]][i]), + abs( + matrix2[order[j + 1]][0] + / matrix2[order[j + 1]][i] + ), ) if min2 > min1: @@ -293,7 +309,9 @@ def _kernel( kernel_dim = 0 for i in range(num_species): - done = all(matrix[i][j] >= num_reactions for j in range(len(matrix[i]))) + done = all( + matrix[i][j] >= num_reactions for j in range(len(matrix[i])) + ) if done and len(matrix[i]): for j in range(len(matrix[i])): RSolutions[kernel_dim].append(matrix[i][j] - num_reactions) @@ -334,7 +352,8 @@ def _kernel( assert int_kernel_dim <= kernel_dim assert len(cls_species_idxs) == len(cls_coefficients), ( - "Inconsistent number of conserved quantities in coefficients and " "species" + "Inconsistent number of conserved quantities in coefficients and " + "species" ) return ( kernel_dim, @@ -474,7 +493,10 @@ def _is_linearly_dependent( for i in range(len(matrix[order[j + 1]])): min2 = min( min2, - abs(matrix2[order[j + 1]][0] / matrix2[order[j + 1]][i]), + abs( + matrix2[order[j + 1]][0] + / matrix2[order[j + 1]][i] + ), ) if min2 > min1: # swap @@ -556,7 +578,9 @@ def _monte_carlo( considered otherwise the algorithm retries Monte Carlo up to max_iter """ dim = len(matched) - num = [int(2 * random.uniform(0, 1)) if len(J[i]) else 0 for i in range(dim)] + num = [ + int(2 * random.uniform(0, 1)) if len(J[i]) else 0 for i in range(dim) + ] numtot = sum(num) def compute_h(): @@ -728,7 +752,10 @@ def _relax( for i in range(len(matrix[order[j + 1]])): min2 = min( min2, - abs(matrix2[order[j + 1]][0] / matrix2[order[j + 1]][i]), + abs( + matrix2[order[j + 1]][0] + / matrix2[order[j + 1]][i] + ), ) if min2 > min1: # swap @@ -787,7 +814,9 @@ def _relax( row_k[matrix[j][a]] -= matrix2[j][a] * matrix2[k][i] # filter matrix[k] = [ - row_idx for row_idx, row_val in enumerate(row_k) if row_val != 0 + row_idx + for row_idx, row_val in enumerate(row_k) + if row_val != 0 ] matrix2[k] = [row_val for row_val in row_k if row_val != 0] diff --git a/python/sdist/amici/custom_commands.py b/python/sdist/amici/custom_commands.py index 2e69800fc7..d54060a009 100644 --- a/python/sdist/amici/custom_commands.py +++ b/python/sdist/amici/custom_commands.py @@ -171,7 +171,8 @@ def build_extension(self, ext: CMakeExtension) -> None: build_dir = self.build_lib if self.inplace == 0 else os.getcwd() build_dir = Path(build_dir).absolute().as_posix() ext.cmake_configure_options = [ - x.replace("${build_dir}", build_dir) for x in ext.cmake_configure_options + x.replace("${build_dir}", build_dir) + for x in ext.cmake_configure_options ] super().build_extension(ext) diff --git a/python/sdist/amici/cxxcodeprinter.py b/python/sdist/amici/cxxcodeprinter.py index e043590db6..e6e377b331 100644 --- a/python/sdist/amici/cxxcodeprinter.py +++ b/python/sdist/amici/cxxcodeprinter.py @@ -68,7 +68,9 @@ def doprint(self, expr: sp.Expr, assign_to: Optional[str] = None) -> str: def _print_min_max(self, expr, cpp_fun: str, sympy_fun): # C++ doesn't like mixing int and double for arguments for min/max, # therefore, we just always convert to float - arg0 = sp.Float(expr.args[0]) if expr.args[0].is_number else expr.args[0] + arg0 = ( + sp.Float(expr.args[0]) if expr.args[0].is_number else expr.args[0] + ) if len(expr.args) == 1: return self._print(arg0) return "%s%s(%s, %s)" % ( @@ -108,7 +110,8 @@ def _get_sym_lines_array( C++ code as list of lines """ return [ - " " * indent_level + f"{variable}[{index}] = " f"{self.doprint(math)};" + " " * indent_level + f"{variable}[{index}] = " + f"{self.doprint(math)};" for index, math in enumerate(equations) if math not in [0, 0.0] ] @@ -150,7 +153,9 @@ def format_regular_line(symbol, math, index): if self.extract_cse: # Extract common subexpressions cse_sym_prefix = "__amici_cse_" - symbol_generator = numbered_symbols(cls=sp.Symbol, prefix=cse_sym_prefix) + symbol_generator = numbered_symbols( + cls=sp.Symbol, prefix=cse_sym_prefix + ) replacements, reduced_exprs = sp.cse( equations, symbols=symbol_generator, @@ -166,7 +171,9 @@ def format_regular_line(symbol, math, index): sorted_symbols = toposort( { identifier: { - s for s in definition.free_symbols if s in expr_dict + s + for s in definition.free_symbols + if s in expr_dict } for (identifier, definition) in expr_dict.items() } @@ -182,7 +189,9 @@ def format_line(symbol: sp.Symbol): f"= {self.doprint(math)};" ) elif math not in [0, 0.0]: - return format_regular_line(symbol, math, symbol_to_idx[symbol]) + return format_regular_line( + symbol, math, symbol_to_idx[symbol] + ) return [ line @@ -251,7 +260,9 @@ def csc_matrix( symbol_row_vals.append(row) idx += 1 - symbol_name = f"d{rownames[row].name}" f"_d{colnames[col].name}" + symbol_name = ( + f"d{rownames[row].name}" f"_d{colnames[col].name}" + ) if identifier: symbol_name += f"_{identifier}" symbol_list.append(symbol_name) diff --git a/python/sdist/amici/de_export.py b/python/sdist/amici/de_export.py index d89ddd9879..e91e5f4c59 100644 --- a/python/sdist/amici/de_export.py +++ b/python/sdist/amici/de_export.py @@ -70,9 +70,13 @@ # Template for model simulation main.cpp file CXX_MAIN_TEMPLATE_FILE = os.path.join(amiciSrcPath, "main.template.cpp") # Template for model/swig/CMakeLists.txt -SWIG_CMAKE_TEMPLATE_FILE = os.path.join(amiciSwigPath, "CMakeLists_model.cmake") +SWIG_CMAKE_TEMPLATE_FILE = os.path.join( + amiciSwigPath, "CMakeLists_model.cmake" +) # Template for model/CMakeLists.txt -MODEL_CMAKE_TEMPLATE_FILE = os.path.join(amiciSrcPath, "CMakeLists.template.cmake") +MODEL_CMAKE_TEMPLATE_FILE = os.path.join( + amiciSrcPath, "CMakeLists.template.cmake" +) IDENTIFIER_PATTERN = re.compile(r"^[a-zA-Z_]\w*$") DERIVATIVE_PATTERN = re.compile(r"^d(x_rdata|xdot|\w+?)d(\w+?)(?:_explicit)?$") @@ -285,7 +289,8 @@ def arguments(self, ode: bool = True) -> str: " const realtype *k, const int ip", ), "sigmaz": _FunctionInfo( - "realtype *sigmaz, const realtype t, const realtype *p, " "const realtype *k", + "realtype *sigmaz, const realtype t, const realtype *p, " + "const realtype *k", ), "sroot": _FunctionInfo( "realtype *stau, const realtype t, const realtype *x, " @@ -326,7 +331,8 @@ def arguments(self, ode: bool = True) -> str: assume_pow_positivity=True, ), "x0": _FunctionInfo( - "realtype *x0, const realtype t, const realtype *p, " "const realtype *k" + "realtype *x0, const realtype t, const realtype *p, " + "const realtype *k" ), "x0_fixedParameters": _FunctionInfo( "realtype *x0_fixedParameters, const realtype t, " @@ -1015,15 +1021,21 @@ def transform_dxdt_to_concentration(species_id, dxdt): # we need to flatten out assignments in the compartment in # order to ensure that we catch all species dependencies - v = smart_subs_dict(v, si.symbols[SymbolId.EXPRESSION], "value") + v = smart_subs_dict( + v, si.symbols[SymbolId.EXPRESSION], "value" + ) dv_dt = v.diff(amici_time_symbol) # we may end up with a time derivative of the compartment # volume due to parameter rate rules comp_rate_vars = [ - p for p in v.free_symbols if p in si.symbols[SymbolId.SPECIES] + p + for p in v.free_symbols + if p in si.symbols[SymbolId.SPECIES] ] for var in comp_rate_vars: - dv_dt += v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] + dv_dt += ( + v.diff(var) * si.symbols[SymbolId.SPECIES][var]["dt"] + ) dv_dx = v.diff(species_id) xdot = (dxdt - dv_dt * species_id) / (dv_dx * species_id + v) return xdot @@ -1042,7 +1054,9 @@ def transform_dxdt_to_concentration(species_id, dxdt): return dxdt / v # create dynamics without respecting conservation laws first - dxdt = smart_multiply(si.stoichiometric_matrix, MutableDenseMatrix(fluxes)) + dxdt = smart_multiply( + si.stoichiometric_matrix, MutableDenseMatrix(fluxes) + ) for ix, ((species_id, species), formula) in enumerate( zip(symbols[SymbolId.SPECIES].items(), dxdt) ): @@ -1052,7 +1066,9 @@ def transform_dxdt_to_concentration(species_id, dxdt): if species["amount"]: species["dt"] = formula else: - species["dt"] = transform_dxdt_to_concentration(species_id, formula) + species["dt"] = transform_dxdt_to_concentration( + species_id, formula + ) # create all basic components of the DE model and add them. for symbol_name in symbols: @@ -1100,11 +1116,14 @@ def transform_dxdt_to_concentration(species_id, dxdt): # fill in 'self._sym' based on prototypes and components in ode_model self.generate_basic_variables() self._has_quadratic_nllh = all( - llh["dist"] in ["normal", "lin-normal", "log-normal", "log10-normal"] + llh["dist"] + in ["normal", "lin-normal", "log-normal", "log10-normal"] for llh in si.symbols[SymbolId.LLHY].values() ) - self._process_sbml_rate_of(symbols) # substitute SBML-rateOf constructs + self._process_sbml_rate_of( + symbols + ) # substitute SBML-rateOf constructs def _process_sbml_rate_of(self, symbols) -> None: """Substitute any SBML-rateOf constructs in the model equations""" @@ -1144,7 +1163,10 @@ def get_rate(symbol: sp.Symbol): for i_state in range(len(self.eq("x0"))): if rate_ofs := self._eqs["x0"][i_state].find(rate_of_func): self._eqs["x0"][i_state] = self._eqs["x0"][i_state].subs( - {rate_of: get_rate(rate_of.args[0]) for rate_of in rate_ofs} + { + rate_of: get_rate(rate_of.args[0]) + for rate_of in rate_ofs + } ) for component in chain( @@ -1171,7 +1193,10 @@ def get_rate(symbol: sp.Symbol): component.set_val( component.get_val().subs( - {rate_of: get_rate(rate_of.args[0]) for rate_of in rate_ofs} + { + rate_of: get_rate(rate_of.args[0]) + for rate_of in rate_ofs + } ) ) @@ -1267,14 +1292,21 @@ def add_conservation_law( )[0] except StopIteration: raise ValueError( - f"Specified state {state} was not found in the " f"model states." + f"Specified state {state} was not found in the " + f"model states." ) state_id = self._differential_states[ix].get_id() # \sum_{i≠j}(a_i * x_i)/a_j target_expression = ( - sp.Add(*(c_i * x_i for x_i, c_i in coefficients.items() if x_i != state)) + sp.Add( + *( + c_i * x_i + for x_i, c_i in coefficients.items() + if x_i != state + ) + ) / coefficients[state] ) @@ -1471,7 +1503,9 @@ def sparseeq(self, name) -> sp.Matrix: self._generate_sparse_symbol(name) return self._sparseeqs[name] - def colptrs(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: + def colptrs( + self, name: str + ) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the column pointers for a sparsified symbolic variable. @@ -1488,7 +1522,9 @@ def colptrs(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: self._generate_sparse_symbol(name) return self._colptrs[name] - def rowvals(self, name: str) -> Union[List[sp.Number], List[List[sp.Number]]]: + def rowvals( + self, name: str + ) -> Union[List[sp.Number], List[List[sp.Number]]]: """ Returns (and constructs if necessary) the row values for a sparsified symbolic variable. @@ -1556,7 +1592,9 @@ def _generate_symbol(self, name: str) -> None: """ if name in self._variable_prototype: components = self._variable_prototype[name]() - self._syms[name] = sp.Matrix([comp.get_id() for comp in components]) + self._syms[name] = sp.Matrix( + [comp.get_id() for comp in components] + ) if name == "y": self._syms["my"] = sp.Matrix( [comp.get_measurement_symbol() for comp in components] @@ -1740,7 +1778,9 @@ def get_appearance_counts(self, idxs: List[int]) -> List[int]: return [ free_symbols_dt.count(str(self._differential_states[idx].get_id())) - + free_symbols_expr.count(str(self._differential_states[idx].get_id())) + + free_symbols_expr.count( + str(self._differential_states[idx].get_id()) + ) for idx in idxs ] @@ -1825,7 +1865,9 @@ def _compute_equation(self, name: str) -> None: time_symbol = sp.Matrix([amici_time_symbol]) if name in self._equation_prototype: - self._equation_from_components(name, self._equation_prototype[name]()) + self._equation_from_components( + name, self._equation_prototype[name]() + ) elif name in self._total_derivative_prototypes: args = self._total_derivative_prototypes[name] @@ -1915,7 +1957,9 @@ def _compute_equation(self, name: str) -> None: if any(sym in eq.free_symbols for sym in k) ] eq = self.eq("x0") - self._eqs[name] = sp.Matrix([eq[ix] for ix in self._x0_fixedParameters_idx]) + self._eqs[name] = sp.Matrix( + [eq[ix] for ix in self._x0_fixedParameters_idx] + ) elif name == "dtotal_cldx_rdata": x_rdata = self.sym("x_rdata") @@ -1928,7 +1972,9 @@ def _compute_equation(self, name: str) -> None: elif name == "dtcldx": # this is always zero - self._eqs[name] = sp.zeros(self.num_cons_law(), self.num_states_solver()) + self._eqs[name] = sp.zeros( + self.num_cons_law(), self.num_states_solver() + ) elif name == "dtcldp": # force symbols @@ -1953,15 +1999,21 @@ def _compute_equation(self, name: str) -> None: elif name == "dx_rdatadp": if self.num_cons_law(): - self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("p")) + self._eqs[name] = smart_jacobian( + self.eq("x_rdata"), self.sym("p") + ) else: # so far, dx_rdatadp is only required for sx_rdata # in case of no conservation laws, C++ code will directly use # sx, we don't need this - self._eqs[name] = sp.zeros(self.num_states_rdata(), self.num_par()) + self._eqs[name] = sp.zeros( + self.num_states_rdata(), self.num_par() + ) elif name == "dx_rdatadtcl": - self._eqs[name] = smart_jacobian(self.eq("x_rdata"), self.sym("tcl")) + self._eqs[name] = smart_jacobian( + self.eq("x_rdata"), self.sym("tcl") + ) elif name == "dxdotdx_explicit": # force symbols @@ -2024,7 +2076,9 @@ def _compute_equation(self, name: str) -> None: self._eqs[name] = event_eqs elif name == "z": - event_observables = [sp.zeros(self.num_eventobs(), 1) for _ in self._events] + event_observables = [ + sp.zeros(self.num_eventobs(), 1) for _ in self._events + ] event_ids = [e.get_id() for e in self._events] # TODO: get rid of this stupid 1-based indexing as soon as we can # the matlab interface @@ -2032,7 +2086,9 @@ def _compute_equation(self, name: str) -> None: event_ids.index(event_obs.get_event()) + 1 for event_obs in self._event_observables ] - for (iz, ie), event_obs in zip(enumerate(z2event), self._event_observables): + for (iz, ie), event_obs in zip( + enumerate(z2event), self._event_observables + ): event_observables[ie - 1][iz] = event_obs.get_val() self._eqs[name] = event_observables @@ -2050,7 +2106,10 @@ def _compute_equation(self, name: str) -> None: ] if name == "dzdx": for ie in range(self.num_events()): - dtaudx = -self.eq("drootdx")[ie, :] / self.eq("drootdt_total")[ie] + dtaudx = ( + -self.eq("drootdx")[ie, :] + / self.eq("drootdt_total")[ie] + ) for iz in range(self.num_eventobs()): if ie != self._z2event[iz] - 1: continue @@ -2121,7 +2180,9 @@ def _compute_equation(self, name: str) -> None: ) # finish chain rule for the state variables - tmp_eq += smart_multiply(self.eq("ddeltaxdx")[ie], tmp_dxdp) + tmp_eq += smart_multiply( + self.eq("ddeltaxdx")[ie], tmp_dxdp + ) event_eqs.append(tmp_eq) @@ -2140,7 +2201,9 @@ def _compute_equation(self, name: str) -> None: # that we need to reverse the order here for cl in reversed(self._conservation_laws) ] - ).col_join(smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x)) + ).col_join( + smart_jacobian(self.eq("w")[self.num_cons_law() :, :], x) + ) elif match_deriv: self._derivative(match_deriv[1], match_deriv[2], name) @@ -2214,7 +2277,10 @@ def _derivative(self, eq: str, var: str, name: str = None) -> None: and cv not in self._lock_total_derivative and var != cv and min(self.sym(cv).shape) - and ((eq, var) not in ignore_chainrule or ignore_chainrule[(eq, var)] != cv) + and ( + (eq, var) not in ignore_chainrule + or ignore_chainrule[(eq, var)] != cv + ) ] if len(chainvars): self._lock_total_derivative += chainvars @@ -2317,9 +2383,14 @@ def _total_derivative( dxdz = self.sym_or_eq(name, dxdz_name) # Save time for large models if one multiplicand is zero, # which is not checked for by sympy - if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix(dxdz): + if not smart_is_zero_matrix(dydx) and not smart_is_zero_matrix( + dxdz + ): dydx_times_dxdz = smart_multiply(dydx, dxdz) - if dxdz.shape[1] == 1 and self._eqs[name].shape[1] != dxdz.shape[1]: + if ( + dxdz.shape[1] == 1 + and self._eqs[name].shape[1] != dxdz.shape[1] + ): for iz in range(self._eqs[name].shape[1]): self._eqs[name][:, iz] += dydx_times_dxdz else: @@ -2345,7 +2416,9 @@ def sym_or_eq(self, name: str, varname: str) -> sp.Matrix: # within a column may differ from the initialization of symbols here, # so those are not safe to use. Not removing them from signature as # this would break backwards compatibility. - if var_in_function_signature(name, varname, self.is_ode()) and varname not in [ + if var_in_function_signature( + name, varname, self.is_ode() + ) and varname not in [ "dwdx", "dwdp", ]: @@ -2469,7 +2542,8 @@ def state_has_fixed_parameter_initial_condition(self, ix: int) -> bool: if not isinstance(ic, sp.Basic): return False return any( - fp in (c.get_id() for c in self._constants) for fp in ic.free_symbols + fp in (c.get_id() for c in self._constants) + for fp in ic.free_symbols ) def state_has_conservation_law(self, ix: int) -> bool: @@ -2786,7 +2860,9 @@ def __init__( # include/amici/model.h for details) self.model: DEModel = de_model self.model._code_printer.known_functions.update( - splines.spline_user_functions(self.model.splines, self._get_index("p")) + splines.spline_user_functions( + self.model.splines, self._get_index("p") + ) ) # To only generate a subset of functions, apply subselection here @@ -2802,7 +2878,9 @@ def generate_model_code(self) -> None: Generates the native C++ code for the loaded model and a Matlab script that can be run to compile a mex file from the C++ code """ - with _monkeypatched(sp.Pow, "_eval_derivative", _custom_pow_eval_derivative): + with _monkeypatched( + sp.Pow, "_eval_derivative", _custom_pow_eval_derivative + ): self._prepare_model_folder() self._generate_c_code() self._generate_m_code() @@ -2864,7 +2942,9 @@ def _generate_c_code(self) -> None: self._write_swig_files() self._write_module_setup() - shutil.copy(CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp")) + shutil.copy( + CXX_MAIN_TEMPLATE_FILE, os.path.join(self.model_path, "main.cpp") + ) def _compile_c_code( self, @@ -2942,8 +3022,10 @@ def _generate_m_code(self) -> None: lines = [ "% This compile script was automatically created from" " Python SBML import.", - "% If mex compiler is set up within MATLAB, it can be run" " from MATLAB ", - "% in order to compile a mex-file from the Python" " generated C++ files.", + "% If mex compiler is set up within MATLAB, it can be run" + " from MATLAB ", + "% in order to compile a mex-file from the Python" + " generated C++ files.", "", f"modelName = '{self.model_name}';", "amimodel.compileAndLinkModel(modelName, '', [], [], [], []);", @@ -2975,7 +3057,10 @@ def _get_index(self, name: str) -> Dict[sp.Symbol, int]: else: raise ValueError(f"Unknown symbolic array: {name}") - return {strip_pysb(symbol).name: index for index, symbol in enumerate(symbols)} + return { + strip_pysb(symbol).name: index + for index, symbol in enumerate(symbols) + } def _write_index_files(self, name: str) -> None: """ @@ -3030,7 +3115,8 @@ def _write_function_file(self, function: str) -> None: if function in sparse_functions: equations = self.model.sparseeq(function) elif ( - not self.allow_reinit_fixpar_initcond and function == "sx0_fixedParameters" + not self.allow_reinit_fixpar_initcond + and function == "sx0_fixedParameters" ): # Not required. Will create empty function body. equations = sp.Matrix() @@ -3096,14 +3182,21 @@ def _write_function_file(self, function: str) -> None: if iszero and not ( (sym == "y" and "Jy" in function) - or (sym == "w" and "xdot" in function and len(self.model.sym(sym))) + or ( + sym == "w" + and "xdot" in function + and len(self.model.sym(sym)) + ) ): continue lines.append(f'#include "{sym}.h"') # include return symbols - if function in self.model.sym_names() and function not in non_unique_id_symbols: + if ( + function in self.model.sym_names() + and function not in non_unique_id_symbols + ): lines.append(f'#include "{function}.h"') lines.extend( @@ -3255,7 +3348,9 @@ def _generate_function_index( return lines - def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: + def _get_function_body( + self, function: str, equations: sp.Matrix + ) -> List[str]: """ Generate C++ code for body of function ``function``. @@ -3294,7 +3389,9 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: + str(len(self.model._x0_fixedParameters_idx)) + "> _x0_fixedParameters_idxs = {", " " - + ", ".join(str(x) for x in self.model._x0_fixedParameters_idx), + + ", ".join( + str(x) for x in self.model._x0_fixedParameters_idx + ), " };", "", # Set all parameters that are to be reset to 0, so that the @@ -3331,7 +3428,9 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: lines.extend(get_switch_statement("ip", cases, 1)) elif function == "x0_fixedParameters": - for index, formula in zip(self.model._x0_fixedParameters_idx, equations): + for index, formula in zip( + self.model._x0_fixedParameters_idx, equations + ): lines.append( f" if(std::find(reinitialization_state_idxs.cbegin(), " f"reinitialization_state_idxs.cend(), {index}) != " @@ -3365,7 +3464,10 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: outer_cases[ie] = copy.copy(inner_lines) lines.extend(get_switch_statement("ie", outer_cases, 1)) - elif function in sensi_functions and equations.shape[1] == self.model.num_par(): + elif ( + function in sensi_functions + and equations.shape[1] == self.model.num_par() + ): cases = { ipar: self.model._code_printer._get_sym_lines_array( equations[:, ipar], function, 0 @@ -3398,7 +3500,8 @@ def _get_function_body(self, function: str, equations: sp.Matrix) -> List[str]: lines.extend(get_switch_statement(iterator, cases, 1)) elif ( - function in self.model.sym_names() and function not in non_unique_id_symbols + function in self.model.sym_names() + and function not in non_unique_id_symbols ): if function in sparse_functions: symbols = list(map(sp.Symbol, self.model.sparsesym(function))) @@ -3425,19 +3528,21 @@ def _get_create_splines_body(self): body = ["return {"] for ispl, spline in enumerate(self.model.splines): if isinstance(spline.nodes, splines.UniformGrid): - nodes = f"{ind8}{{{spline.nodes.start}, {spline.nodes.stop}}}, " + nodes = ( + f"{ind8}{{{spline.nodes.start}, {spline.nodes.stop}}}, " + ) else: nodes = f"{ind8}{{{', '.join(map(str, spline.nodes))}}}, " # vector with the node values - values = f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " + values = ( + f"{ind8}{{{', '.join(map(str, spline.values_at_nodes))}}}, " + ) # vector with the slopes if spline.derivatives_by_fd: slopes = f"{ind8}{{}}," else: - slopes = ( - f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," - ) + slopes = f"{ind8}{{{', '.join(map(str, spline.derivatives_at_nodes))}}}," body.extend( [ @@ -3460,7 +3565,8 @@ def _get_create_splines_body(self): body.append(ind8 + bc_to_cpp[bc]) except KeyError: raise ValueError( - f"Unknown boundary condition '{bc}' " "found in spline object" + f"Unknown boundary condition '{bc}' " + "found in spline object" ) extrapolate_to_cpp = { None: "SplineExtrapolation::noExtrapolation, ", @@ -3474,12 +3580,15 @@ def _get_create_splines_body(self): body.append(ind8 + extrapolate_to_cpp[extr]) except KeyError: raise ValueError( - f"Unknown extrapolation '{extr}' " "found in spline object" + f"Unknown extrapolation '{extr}' " + "found in spline object" ) line = ind8 line += "true, " if spline.derivatives_by_fd else "false, " line += ( - "true, " if isinstance(spline.nodes, splines.UniformGrid) else "false, " + "true, " + if isinstance(spline.nodes, splines.UniformGrid) + else "false, " ) line += "true" if spline.logarithmic_parametrization else "false" body.append(line) @@ -3558,10 +3667,12 @@ def _write_model_header_cpp(self) -> None: "NK": self.model.num_const(), "O2MODE": "amici::SecondOrderMode::none", # using code printer ensures proper handling of nan/inf - "PARAMETERS": self.model._code_printer.doprint(self.model.val("p"))[1:-1], - "FIXED_PARAMETERS": self.model._code_printer.doprint(self.model.val("k"))[ - 1:-1 - ], + "PARAMETERS": self.model._code_printer.doprint( + self.model.val("p") + )[1:-1], + "FIXED_PARAMETERS": self.model._code_printer.doprint( + self.model.val("k") + )[1:-1], "PARAMETER_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "p" ), @@ -3576,12 +3687,16 @@ def _write_model_header_cpp(self) -> None: ), "OBSERVABLE_TRAFO_INITIALIZER_LIST": "\n".join( f"ObservableScaling::{trafo.value}, // y[{idx}]" - for idx, trafo in enumerate(self.model.get_observable_transformations()) + for idx, trafo in enumerate( + self.model.get_observable_transformations() + ) ), "EXPRESSION_NAMES_INITIALIZER_LIST": self._get_symbol_name_initializer_list( "w" ), - "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list("p"), + "PARAMETER_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( + "p" + ), "STATE_IDS_INITIALIZER_LIST": self._get_symbol_id_initializer_list( "x_rdata" ), @@ -3662,16 +3777,22 @@ def _write_model_header_cpp(self) -> None: indexfield, nobody=True, ) - tpl_data[f"{func_name.upper()}_{indexfield.upper()}_DEF"] = "" + tpl_data[ + f"{func_name.upper()}_{indexfield.upper()}_DEF" + ] = "" tpl_data[ f"{func_name.upper()}_{indexfield.upper()}_IMPL" ] = impl continue - tpl_data[f"{func_name.upper()}_DEF"] = get_function_extern_declaration( + tpl_data[ + f"{func_name.upper()}_DEF" + ] = get_function_extern_declaration( func_name, self.model_name, self.model.is_ode() ) - tpl_data[f"{func_name.upper()}_IMPL"] = get_model_override_implementation( + tpl_data[ + f"{func_name.upper()}_IMPL" + ] = get_model_override_implementation( func_name, self.model_name, self.model.is_ode() ) if func_name in sparse_functions: @@ -3902,7 +4023,9 @@ def get_function_extern_declaration(fun: str, name: str, ode: bool) -> str: return f"extern {f.return_type} {fun}_{name}({f.arguments(ode)});" -def get_sunindex_extern_declaration(fun: str, name: str, indextype: str) -> str: +def get_sunindex_extern_declaration( + fun: str, name: str, indextype: str +) -> str: """ Constructs the function declaration for an index function of a given function diff --git a/python/sdist/amici/de_model.py b/python/sdist/amici/de_model.py index c5363511e7..77d9013ad2 100644 --- a/python/sdist/amici/de_model.py +++ b/python/sdist/amici/de_model.py @@ -67,7 +67,8 @@ def __init__( hasattr(identifier, "name") and identifier.name in RESERVED_SYMBOLS ): raise ValueError( - f'Cannot add model quantity with name "{name}", ' f"please rename." + f'Cannot add model quantity with name "{name}", ' + f"please rename." ) self._identifier: sp.Symbol = identifier @@ -301,7 +302,9 @@ class DifferentialState(State): """ - def __init__(self, identifier: sp.Symbol, name: str, init: sp.Expr, dt: sp.Expr): + def __init__( + self, identifier: sp.Symbol, name: str, init: sp.Expr, dt: sp.Expr + ): """ Create a new State instance. Extends :meth:`ModelQuantity.__init__` by ``dt`` @@ -335,7 +338,8 @@ def set_conservation_law(self, law: ConservationLaw) -> None: """ if not isinstance(law, ConservationLaw): raise TypeError( - f"conservation law must have type ConservationLaw" f", was {type(law)}" + f"conservation law must have type ConservationLaw" + f", was {type(law)}" ) self._conservation_law = law @@ -425,13 +429,17 @@ def __init__( def get_measurement_symbol(self) -> sp.Symbol: if self._measurement_symbol is None: - self._measurement_symbol = generate_measurement_symbol(self.get_id()) + self._measurement_symbol = generate_measurement_symbol( + self.get_id() + ) return self._measurement_symbol def get_regularization_symbol(self) -> sp.Symbol: if self._regularization_symbol is None: - self._regularization_symbol = generate_regularization_symbol(self.get_id()) + self._regularization_symbol = generate_regularization_symbol( + self.get_id() + ) return self._regularization_symbol @@ -556,7 +564,9 @@ class Parameter(ModelQuantity): sensitivities may be computed, abbreviated by ``p``. """ - def __init__(self, identifier: sp.Symbol, name: str, value: numbers.Number): + def __init__( + self, identifier: sp.Symbol, name: str, value: numbers.Number + ): """ Create a new Expression instance. @@ -579,7 +589,9 @@ class Constant(ModelQuantity): sensitivities cannot be computed, abbreviated by ``k``. """ - def __init__(self, identifier: sp.Symbol, name: str, value: numbers.Number): + def __init__( + self, identifier: sp.Symbol, name: str, value: numbers.Number + ): """ Create a new Expression instance. diff --git a/python/sdist/amici/gradient_check.py b/python/sdist/amici/gradient_check.py index c0992ba213..27e2d671d3 100644 --- a/python/sdist/amici/gradient_check.py +++ b/python/sdist/amici/gradient_check.py @@ -193,7 +193,8 @@ def check_derivatives( fields.append("x") leastsquares_applicable = ( - solver.getSensitivityMethod() == SensitivityMethod.forward and edata is not None + solver.getSensitivityMethod() == SensitivityMethod.forward + and edata is not None ) if ( @@ -333,4 +334,6 @@ def _check_results( if type(result) is float: result = np.array(result) - _check_close(result=result, expected=expected, atol=atol, rtol=rtol, field=field) + _check_close( + result=result, expected=expected, atol=atol, rtol=rtol, field=field + ) diff --git a/python/sdist/amici/import_utils.py b/python/sdist/amici/import_utils.py index 7299b89beb..77a2add60b 100644 --- a/python/sdist/amici/import_utils.py +++ b/python/sdist/amici/import_utils.py @@ -45,14 +45,18 @@ def __init__(self, data): s = "Circular dependencies exist among these items: {{{}}}".format( ", ".join( "{!r}:{!r}".format(key, value) - for key, value in sorted({str(k): v for k, v in data.items()}.items()) + for key, value in sorted( + {str(k): v for k, v in data.items()}.items() + ) ) ) super(CircularDependencyError, self).__init__(s) self.data = data -setattr(sys.modules["toposort"], "CircularDependencyError", CircularDependencyError) +setattr( + sys.modules["toposort"], "CircularDependencyError", CircularDependencyError +) annotation_namespace = "https://github.com/AMICI-dev/AMICI" @@ -215,7 +219,8 @@ def noise_distribution_to_cost_function( y_string = "log(2*{sigma}*{m}) + Abs(log({y}) - log({m})) / {sigma}" elif noise_distribution == "log10-laplace": y_string = ( - "log(2*{sigma}*{m}*log(10)) " "+ Abs(log({y}, 10) - log({m}, 10)) / {sigma}" + "log(2*{sigma}*{m}*log(10)) " + "+ Abs(log({y}, 10) - log({m}, 10)) / {sigma}" ) elif noise_distribution in ["binomial", "lin-binomial"]: # Binomial noise model parameterized via success probability p @@ -236,7 +241,9 @@ def noise_distribution_to_cost_function( f"- {{m}} * log({{sigma}})" ) else: - raise ValueError(f"Cost identifier {noise_distribution} not recognized.") + raise ValueError( + f"Cost identifier {noise_distribution} not recognized." + ) def nllh_y_string(str_symbol): y, m, sigma = _get_str_symbol_identifiers(str_symbol) @@ -278,7 +285,8 @@ def smart_subs_dict( Substituted symbolic expression """ s = [ - (eid, expr[field] if field is not None else expr) for eid, expr in subs.items() + (eid, expr[field] if field is not None else expr) + for eid, expr in subs.items() ] if reverse: s.reverse() @@ -309,7 +317,9 @@ def smart_subs(element: sp.Expr, old: sp.Symbol, new: sp.Expr) -> sp.Expr: return element.subs(old, new) if element.has(old) else element -def toposort_symbols(symbols: SymbolDef, field: Optional[str] = None) -> SymbolDef: +def toposort_symbols( + symbols: SymbolDef, field: Optional[str] = None +) -> SymbolDef: """ Topologically sort symbol definitions according to their interdependency @@ -386,7 +396,9 @@ def _parse_special_functions(sym: sp.Expr, toplevel: bool = True) -> sp.Expr: if sym.__class__.__name__ in fun_mappings: return fun_mappings[sym.__class__.__name__](*args) - elif sym.__class__.__name__ == "piecewise" or isinstance(sym, sp.Piecewise): + elif sym.__class__.__name__ == "piecewise" or isinstance( + sym, sp.Piecewise + ): if isinstance(sym, sp.Piecewise): # this is sympy piecewise, can't be nested denested_args = args @@ -438,7 +450,9 @@ def _denest_piecewise( # piece was picked previous_was_picked = sp.false # recursively denest those first - for sub_coeff, sub_cond in grouper(_denest_piecewise(cond.args), 2, True): + for sub_coeff, sub_cond in grouper( + _denest_piecewise(cond.args), 2, True + ): # flatten the individual pieces pick_this = sp.And(sp.Not(previous_was_picked), sub_cond) if sub_coeff == sp.true: @@ -519,7 +533,9 @@ def _parse_heaviside_trigger(trigger: sp.Expr) -> sp.Expr: # or(x,y) = not(and(not(x),not(y)) if isinstance(trigger, sp.Or): - return 1 - sp.Mul(*[1 - _parse_heaviside_trigger(arg) for arg in trigger.args]) + return 1 - sp.Mul( + *[1 - _parse_heaviside_trigger(arg) for arg in trigger.args] + ) if isinstance(trigger, sp.And): return sp.Mul(*[_parse_heaviside_trigger(arg) for arg in trigger.args]) @@ -530,7 +546,9 @@ def _parse_heaviside_trigger(trigger: sp.Expr) -> sp.Expr: ) -def grouper(iterable: Iterable, n: int, fillvalue: Any = None) -> Iterable[Tuple[Any]]: +def grouper( + iterable: Iterable, n: int, fillvalue: Any = None +) -> Iterable[Tuple[Any]]: """ Collect data into fixed-length chunks or blocks @@ -662,7 +680,9 @@ def generate_regularization_symbol(observable_id: Union[str, sp.Symbol]): return symbol_with_assumptions(f"r{observable_id}") -def generate_flux_symbol(reaction_index: int, name: Optional[str] = None) -> sp.Symbol: +def generate_flux_symbol( + reaction_index: int, name: Optional[str] = None +) -> sp.Symbol: """ Generate identifier symbol for a reaction flux. This function will always return the same unique python object for a diff --git a/python/sdist/amici/logging.py b/python/sdist/amici/logging.py index 5f548de7a1..2648fc5b28 100644 --- a/python/sdist/amici/logging.py +++ b/python/sdist/amici/logging.py @@ -166,7 +166,8 @@ def get_logger( _setup_logger(**kwargs) elif kwargs: warnings.warn( - "AMICI logger already exists, ignoring keyword " "arguments to setup_logger" + "AMICI logger already exists, ignoring keyword " + "arguments to setup_logger" ) logger = logging.getLogger(logger_name) @@ -193,7 +194,8 @@ def decorator_timer(func): def wrapper_timer(*args, **kwargs): # append pluses to indicate recursion level recursion_level = sum( - frame.function == "wrapper_timer" and frame.filename == __file__ + frame.function == "wrapper_timer" + and frame.filename == __file__ for frame in getouterframes(currentframe(), context=0) ) diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py index 91ca6449f6..23ebfdbbc4 100644 --- a/python/sdist/amici/numpy.py +++ b/python/sdist/amici/numpy.py @@ -220,7 +220,8 @@ def __init__(self, rdata: Union[ReturnDataPtr, ReturnData]): """ if not isinstance(rdata, (ReturnDataPtr, ReturnData)): raise TypeError( - f"Unsupported pointer {type(rdata)}, must be" f"amici.ExpDataPtr!" + f"Unsupported pointer {type(rdata)}, must be" + f"amici.ExpDataPtr!" ) self._field_dimensions = { "ts": [rdata.nt], @@ -293,7 +294,9 @@ def __getitem__( return super().__getitem__(item) - def by_id(self, entity_id: str, field: str = None, model: Model = None) -> np.array: + def by_id( + self, entity_id: str, field: str = None, model: Model = None + ) -> np.array: """ Get the value of a given field for a named entity. @@ -311,11 +314,17 @@ def by_id(self, entity_id: str, field: str = None, model: Model = None) -> np.ar if field in {"x", "x0", "x_ss", "sx", "sx0", "sx_ss"}: ids = (model and model.getStateIds()) or self._swigptr.state_ids elif field in {"w"}: - ids = (model and model.getExpressionIds()) or self._swigptr.expression_ids + ids = ( + model and model.getExpressionIds() + ) or self._swigptr.expression_ids elif field in {"y", "sy", "sigmay"}: - ids = (model and model.getObservableIds()) or self._swigptr.observable_ids + ids = ( + model and model.getObservableIds() + ) or self._swigptr.observable_ids elif field in {"sllh"}: - ids = (model and model.getParameterIds()) or self._swigptr.parameter_ids + ids = ( + model and model.getParameterIds() + ) or self._swigptr.parameter_ids else: raise NotImplementedError( f"Subsetting {field} by ID is not implemented or not possible." @@ -348,7 +357,8 @@ def __init__(self, edata: Union[ExpDataPtr, ExpData]): """ if not isinstance(edata, (ExpDataPtr, ExpData)): raise TypeError( - f"Unsupported pointer {type(edata)}, must be" f"amici.ExpDataPtr!" + f"Unsupported pointer {type(edata)}, must be" + f"amici.ExpDataPtr!" ) self._field_dimensions = { # observables "observedData": [edata.nt(), edata.nytrue()], @@ -411,7 +421,9 @@ def _entity_type_from_id( return symbol else: if entity_id in getattr( - rdata if isinstance(rdata, amici.ReturnData) else rdata._swigptr, + rdata + if isinstance(rdata, amici.ReturnData) + else rdata._swigptr, f"{entity_type.lower()}_ids", ): return symbol diff --git a/python/sdist/amici/pandas.py b/python/sdist/amici/pandas.py index dd240242af..8a2eb5049d 100644 --- a/python/sdist/amici/pandas.py +++ b/python/sdist/amici/pandas.py @@ -107,7 +107,9 @@ def getDataObservablesAsDataFrame( _get_names_or_ids(model, "Observable", by_id=by_id) ): datadict[obs] = npdata["observedData"][i_time, i_obs] - datadict[obs + "_std"] = npdata["observedDataStdDev"][i_time, i_obs] + datadict[obs + "_std"] = npdata["observedDataStdDev"][ + i_time, i_obs + ] # add conditions _fill_conditions_dict(datadict, model, edata, by_id=by_id) @@ -396,12 +398,16 @@ def _fill_conditions_dict( datadict[par] = model.getFixedParameters()[i_par] if len(edata.fixedParametersPreequilibration): - datadict[par + "_preeq"] = edata.fixedParametersPreequilibration[i_par] + datadict[par + "_preeq"] = edata.fixedParametersPreequilibration[ + i_par + ] else: datadict[par + "_preeq"] = np.nan if len(edata.fixedParametersPresimulation): - datadict[par + "_presim"] = edata.fixedParametersPresimulation[i_par] + datadict[par + "_presim"] = edata.fixedParametersPresimulation[ + i_par + ] else: datadict[par + "_presim"] = np.nan return datadict @@ -526,7 +532,9 @@ def _get_expression_cols(model: AmiciModel, by_id: bool) -> List[str]: ) -def _get_names_or_ids(model: AmiciModel, variable: str, by_id: bool) -> List[str]: +def _get_names_or_ids( + model: AmiciModel, variable: str, by_id: bool +) -> List[str]: """ Obtains a unique list of identifiers for the specified variable. First tries model.getVariableNames and then uses model.getVariableIds. @@ -674,22 +682,33 @@ def constructEdataFromDataFrame( ) # fill in preequilibration parameters - if any([overwrite_preeq[key] != condition[key] for key in overwrite_preeq]): - edata.fixedParametersPreequilibration = _get_specialized_fixed_parameters( - model, condition, overwrite_preeq, by_id=by_id + if any( + [overwrite_preeq[key] != condition[key] for key in overwrite_preeq] + ): + edata.fixedParametersPreequilibration = ( + _get_specialized_fixed_parameters( + model, condition, overwrite_preeq, by_id=by_id + ) ) elif len(overwrite_preeq): - edata.fixedParametersPreequilibration = copy.deepcopy(edata.fixedParameters) + edata.fixedParametersPreequilibration = copy.deepcopy( + edata.fixedParameters + ) # fill in presimulation parameters if any( - [overwrite_presim[key] != condition[key] for key in overwrite_presim.keys()] + [ + overwrite_presim[key] != condition[key] + for key in overwrite_presim.keys() + ] ): edata.fixedParametersPresimulation = _get_specialized_fixed_parameters( model, condition, overwrite_presim, by_id=by_id ) elif len(overwrite_presim.keys()): - edata.fixedParametersPresimulation = copy.deepcopy(edata.fixedParameters) + edata.fixedParametersPresimulation = copy.deepcopy( + edata.fixedParameters + ) # fill in presimulation time if "t_presim" in condition.keys(): @@ -739,7 +758,9 @@ def getEdataFromDataFrame( # aggregate features that define a condition # fixed parameters - condition_parameters = _get_names_or_ids(model, "FixedParameter", by_id=by_id) + condition_parameters = _get_names_or_ids( + model, "FixedParameter", by_id=by_id + ) # preeq and presim parameters for par in _get_names_or_ids(model, "FixedParameter", by_id=by_id): if par + "_preeq" in df.columns: @@ -758,7 +779,9 @@ def getEdataFromDataFrame( selected = np.ones((len(df),), dtype=bool) for par_label, par in row.items(): if math.isnan(par): - selected = selected & np.isnan(df[par_label].astype(float).values) + selected = selected & np.isnan( + df[par_label].astype(float).values + ) else: selected = selected & (df[par_label] == par) edata_df = df[selected] diff --git a/python/sdist/amici/parameter_mapping.py b/python/sdist/amici/parameter_mapping.py index 9f4d3b24dd..f1cf75a150 100644 --- a/python/sdist/amici/parameter_mapping.py +++ b/python/sdist/amici/parameter_mapping.py @@ -126,7 +126,9 @@ class ParameterMapping(Sequence): List of parameter mappings for specific conditions. """ - def __init__(self, parameter_mappings: List[ParameterMappingForCondition] = None): + def __init__( + self, parameter_mappings: List[ParameterMappingForCondition] = None + ): super().__init__() if parameter_mappings is None: parameter_mappings = [] @@ -146,7 +148,9 @@ def __getitem__( def __len__(self): return len(self.parameter_mappings) - def append(self, parameter_mapping_for_condition: ParameterMappingForCondition): + def append( + self, parameter_mapping_for_condition: ParameterMappingForCondition + ): """Append a condition specific parameter mapping.""" self.parameter_mappings.append(parameter_mapping_for_condition) @@ -188,7 +192,8 @@ def fill_in_parameters( set(problem_parameters.keys()) - parameter_mapping.free_symbols ): warnings.warn( - "The following problem parameters were not used: " + str(unused_parameters), + "The following problem parameters were not used: " + + str(unused_parameters), RuntimeWarning, ) @@ -262,10 +267,12 @@ def _get_par(model_par, value, mapping): return value map_preeq_fix = { - key: _get_par(key, val, map_preeq_fix) for key, val in map_preeq_fix.items() + key: _get_par(key, val, map_preeq_fix) + for key, val in map_preeq_fix.items() } map_sim_fix = { - key: _get_par(key, val, map_sim_fix) for key, val in map_sim_fix.items() + key: _get_par(key, val, map_sim_fix) + for key, val in map_sim_fix.items() } map_sim_var = { key: _get_par(key, val, dict(map_sim_fix, **map_sim_var)) @@ -289,7 +296,9 @@ def _get_par(model_par, value, mapping): # variable parameters and parameter scale # parameter list from mapping dict - parameters = [map_sim_var[par_id] for par_id in amici_model.getParameterIds()] + parameters = [ + map_sim_var[par_id] for par_id in amici_model.getParameterIds() + ] # scales list from mapping dict scales = [ @@ -317,7 +326,8 @@ def _get_par(model_par, value, mapping): # fixed parameters preequilibration if map_preeq_fix: fixed_pars_preeq = [ - map_preeq_fix[par_id] for par_id in amici_model.getFixedParameterIds() + map_preeq_fix[par_id] + for par_id in amici_model.getFixedParameterIds() ] edata.fixedParametersPreequilibration = fixed_pars_preeq @@ -325,7 +335,8 @@ def _get_par(model_par, value, mapping): # fixed parameters simulation if map_sim_fix: fixed_pars_sim = [ - map_sim_fix[par_id] for par_id in amici_model.getFixedParameterIds() + map_sim_fix[par_id] + for par_id in amici_model.getFixedParameterIds() ] edata.fixedParameters = fixed_pars_sim @@ -370,11 +381,14 @@ def scale_parameter(value: numbers.Number, petab_scale: str) -> numbers.Number: if petab_scale == LOG: return np.log(value) raise ValueError( - f"Unknown parameter scale {petab_scale}. " f"Must be from {(LIN, LOG, LOG10)}" + f"Unknown parameter scale {petab_scale}. " + f"Must be from {(LIN, LOG, LOG10)}" ) -def unscale_parameter(value: numbers.Number, petab_scale: str) -> numbers.Number: +def unscale_parameter( + value: numbers.Number, petab_scale: str +) -> numbers.Number: """Bring parameter from scale to linear scale. :param value: @@ -392,7 +406,8 @@ def unscale_parameter(value: numbers.Number, petab_scale: str) -> numbers.Number if petab_scale == LOG: return np.exp(value) raise ValueError( - f"Unknown parameter scale {petab_scale}. " f"Must be from {(LIN, LOG, LOG10)}" + f"Unknown parameter scale {petab_scale}. " + f"Must be from {(LIN, LOG, LOG10)}" ) diff --git a/python/sdist/amici/petab_import.py b/python/sdist/amici/petab_import.py index eb994a4304..23fe4394f0 100644 --- a/python/sdist/amici/petab_import.py +++ b/python/sdist/amici/petab_import.py @@ -306,7 +306,9 @@ def import_petab_problem( if petab_problem.mapping_df is not None: # It's partially supported. Remove at your own risk... - raise NotImplementedError("PEtab v2.0.0 mapping tables are not yet supported.") + raise NotImplementedError( + "PEtab v2.0.0 mapping tables are not yet supported." + ) model_name = model_name or petab_problem.model.model_id @@ -366,7 +368,9 @@ def import_petab_problem( model = model_module.getModel() check_model(amici_model=model, petab_problem=petab_problem) - logger.info(f"Successfully loaded model {model_name} " f"from {model_output_dir}.") + logger.info( + f"Successfully loaded model {model_name} " f"from {model_output_dir}." + ) return model @@ -383,7 +387,9 @@ def check_model( amici_ids = amici_ids_free | set(amici_model.getFixedParameterIds()) petab_ids_free = set( - petab_problem.parameter_df.loc[petab_problem.parameter_df[ESTIMATE] == 1].index + petab_problem.parameter_df.loc[ + petab_problem.parameter_df[ESTIMATE] == 1 + ].index ) amici_ids_free_required = petab_ids_free.intersection(amici_ids) @@ -429,7 +435,9 @@ def _create_model_name(folder: Union[str, Path]) -> str: return os.path.split(os.path.normpath(folder))[-1] -def _can_import_model(model_name: str, model_output_dir: Union[str, Path]) -> bool: +def _can_import_model( + model_name: str, model_output_dir: Union[str, Path] +) -> bool: """ Check whether a module of that name can already be imported. """ @@ -555,7 +563,8 @@ def import_model_sbml( if petab_problem.observable_df is None: raise NotImplementedError( - "PEtab import without observables table " "is currently not supported." + "PEtab import without observables table " + "is currently not supported." ) assert isinstance(petab_problem.model, SbmlModel) @@ -596,8 +605,10 @@ def import_model_sbml( ) sbml_model = sbml_importer.sbml - allow_n_noise_pars = not petab.lint.observable_table_has_nontrivial_noise_formula( - petab_problem.observable_df + allow_n_noise_pars = ( + not petab.lint.observable_table_has_nontrivial_noise_formula( + petab_problem.observable_df + ) ) if ( petab_problem.measurement_df is not None @@ -632,7 +643,9 @@ def import_model_sbml( # so we add any output parameters to the SBML model. # this should be changed to something more elegant # - formulas = chain((val["formula"] for val in observables.values()), sigmas.values()) + formulas = chain( + (val["formula"] for val in observables.values()), sigmas.values() + ) output_parameters = OrderedDict() for formula in formulas: # we want reproducible parameter ordering upon repeated import @@ -649,10 +662,13 @@ def import_model_sbml( ): output_parameters[sym] = None logger.debug( - "Adding output parameters to model: " f"{list(output_parameters.keys())}" + "Adding output parameters to model: " + f"{list(output_parameters.keys())}" ) output_parameter_defaults = output_parameter_defaults or {} - if extra_pars := (set(output_parameter_defaults) - set(output_parameters.keys())): + if extra_pars := ( + set(output_parameter_defaults) - set(output_parameters.keys()) + ): raise ValueError( f"Default output parameter values were given for {extra_pars}, " "but they those are not output parameters." @@ -691,7 +707,8 @@ def import_model_sbml( # Can only reset parameters after preequilibration if they are fixed. fixed_parameters.append(PREEQ_INDICATOR_ID) logger.debug( - "Adding preequilibration indicator " f"constant {PREEQ_INDICATOR_ID}" + "Adding preequilibration indicator " + f"constant {PREEQ_INDICATOR_ID}" ) logger.debug(f"Adding initial assignments for {initial_states.keys()}") for assignee_id in initial_states: @@ -773,7 +790,9 @@ def import_model_sbml( def get_observation_model( observable_df: pd.DataFrame, -) -> Tuple[Dict[str, Dict[str, str]], Dict[str, str], Dict[str, Union[str, float]]]: +) -> Tuple[ + Dict[str, Dict[str, str]], Dict[str, str], Dict[str, Union[str, float]] +]: """ Get observables, sigmas, and noise distributions from PEtab observation table in a format suitable for @@ -805,7 +824,9 @@ def get_observation_model( # 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) + sp.sympify(observable["formula"], locals=_clash): sp.Symbol( + observable_id + ) for observable_id, observable in observables.items() } for observable_id, formula in sigmas.items(): @@ -871,7 +892,9 @@ def show_model_info(sbml_model: "libsbml.Model"): """Log some model quantities""" logger.info(f"Species: {len(sbml_model.getListOfSpecies())}") - logger.info("Global parameters: " + str(len(sbml_model.getListOfParameters()))) + logger.info( + "Global parameters: " + str(len(sbml_model.getListOfParameters())) + ) logger.info(f"Reactions: {len(sbml_model.getListOfReactions())}") diff --git a/python/sdist/amici/petab_import_pysb.py b/python/sdist/amici/petab_import_pysb.py index d10f6c92f1..8036d1358d 100644 --- a/python/sdist/amici/petab_import_pysb.py +++ b/python/sdist/amici/petab_import_pysb.py @@ -23,7 +23,9 @@ logger = get_logger(__name__, logging.WARNING) -def _add_observation_model(pysb_model: pysb.Model, petab_problem: petab.Problem): +def _add_observation_model( + pysb_model: pysb.Model, petab_problem: petab.Problem +): """Extend PySB model by observation model as defined in the PEtab observables table""" @@ -65,7 +67,9 @@ def _add_observation_model(pysb_model: pysb.Model, petab_problem: petab.Problem) local_syms[sigma_id] = sigma_expr -def _add_initialization_variables(pysb_model: pysb.Model, petab_problem: petab.Problem): +def _add_initialization_variables( + pysb_model: pysb.Model, petab_problem: petab.Problem +): """Add initialization variables to the PySB model to support initial conditions specified in the PEtab condition table. @@ -92,7 +96,8 @@ def _add_initialization_variables(pysb_model: pysb.Model, petab_problem: petab.P # Can only reset parameters after preequilibration if they are fixed. fixed_parameters.append(PREEQ_INDICATOR_ID) logger.debug( - "Adding preequilibration indicator constant " f"{PREEQ_INDICATOR_ID}" + "Adding preequilibration indicator constant " + f"{PREEQ_INDICATOR_ID}" ) logger.debug(f"Adding initial assignments for {initial_states.keys()}") @@ -131,7 +136,9 @@ def _add_initialization_variables(pysb_model: pysb.Model, petab_problem: petab.P pysb_model.add_component(formula) for initial in pysb_model.initials: - if match_complex_pattern(initial.pattern, species_pattern, exact=True): + if match_complex_pattern( + initial.pattern, species_pattern, exact=True + ): logger.debug( "The PySB model has an initial defined for species " f"{assignee_id}, but this species also has an initial " @@ -231,7 +238,9 @@ def import_model_pysb( petab_noise_distributions_to_amici, ) - constant_parameters = get_fixed_parameters(petab_problem) + fixed_parameters + constant_parameters = ( + get_fixed_parameters(petab_problem) + fixed_parameters + ) if petab_problem.observable_df is None: observables = None @@ -246,7 +255,9 @@ def import_model_pysb( sigmas = {obs_id: f"{obs_id}_sigma" for obs_id in observables} - noise_distrs = petab_noise_distributions_to_amici(petab_problem.observable_df) + noise_distrs = petab_noise_distributions_to_amici( + petab_problem.observable_df + ) from amici.pysb_import import pysb2amici diff --git a/python/sdist/amici/petab_objective.py b/python/sdist/amici/petab_objective.py index 0656619227..87409ee446 100644 --- a/python/sdist/amici/petab_objective.py +++ b/python/sdist/amici/petab_objective.py @@ -144,7 +144,11 @@ def simulate_petab( # number of amici simulations will be number of unique # (preequilibrationConditionId, simulationConditionId) pairs. # Can be optimized by checking for identical condition vectors. - if simulation_conditions is None and parameter_mapping is None and edatas is None: + if ( + simulation_conditions is None + and parameter_mapping is None + and edatas is None + ): simulation_conditions = ( petab_problem.get_simulation_conditions_from_measurement_df() ) @@ -262,7 +266,8 @@ def aggregate_sllh( if petab_scale and petab_problem is None: raise ValueError( - "Please provide the PEtab problem, when using " "`petab_scale=True`." + "Please provide the PEtab problem, when using " + "`petab_scale=True`." ) # Check for issues in all condition simulation results. @@ -280,7 +285,9 @@ def aggregate_sllh( for condition_parameter_mapping, edata, rdata in zip( parameter_mapping, edatas, rdatas ): - for sllh_parameter_index, condition_parameter_sllh in enumerate(rdata.sllh): + for sllh_parameter_index, condition_parameter_sllh in enumerate( + rdata.sllh + ): # Get PEtab parameter ID # Use ExpData if it provides a parameter list, else default to # Model. @@ -301,9 +308,11 @@ def aggregate_sllh( if petab_scale: # `ParameterMappingForCondition` objects provide the scale in # terms of `petab.C` constants already, not AMICI equivalents. - model_parameter_scale = condition_parameter_mapping.scale_map_sim_var[ - model_parameter_id - ] + model_parameter_scale = ( + condition_parameter_mapping.scale_map_sim_var[ + model_parameter_id + ] + ) petab_parameter_scale = petab_problem.parameter_df.loc[ petab_parameter_id, PARAMETER_SCALE ] @@ -362,7 +371,9 @@ def rescale_sensitivity( scale[(LOG10, LOG)] = lambda s: scale[(LIN, LOG)](scale[(LOG10, LIN)](s)) if (old_scale, new_scale) not in scale: - raise NotImplementedError(f"Old scale: {old_scale}. New scale: {new_scale}.") + raise NotImplementedError( + f"Old scale: {old_scale}. New scale: {new_scale}." + ) return scale[(old_scale, new_scale)](sensitivity) @@ -497,14 +508,18 @@ def create_parameter_mapping( if parameter_mapping_kwargs is None: parameter_mapping_kwargs = {} - prelim_parameter_mapping = petab.get_optimization_to_simulation_parameter_mapping( - condition_df=petab_problem.condition_df, - measurement_df=petab_problem.measurement_df, - parameter_df=petab_problem.parameter_df, - observable_df=petab_problem.observable_df, - mapping_df=petab_problem.mapping_df, - model=petab_problem.model, - **dict(default_parameter_mapping_kwargs, **parameter_mapping_kwargs), + prelim_parameter_mapping = ( + petab.get_optimization_to_simulation_parameter_mapping( + condition_df=petab_problem.condition_df, + measurement_df=petab_problem.measurement_df, + parameter_df=petab_problem.parameter_df, + observable_df=petab_problem.observable_df, + mapping_df=petab_problem.mapping_df, + model=petab_problem.model, + **dict( + default_parameter_mapping_kwargs, **parameter_mapping_kwargs + ), + ) ) parameter_mapping = ParameterMapping() @@ -539,12 +554,21 @@ def _get_initial_state_sbml( else initial_assignment ) elif type_code == libsbml.SBML_PARAMETER: - value = element.getValue() if initial_assignment is None else initial_assignment + value = ( + element.getValue() + if initial_assignment is None + else initial_assignment + ) elif type_code == libsbml.SBML_COMPARTMENT: - value = element.getSize() if initial_assignment is None else initial_assignment + value = ( + element.getSize() + if initial_assignment is None + else initial_assignment + ) else: raise NotImplementedError( - f"Don't know what how to handle {element_id} in " "condition table." + f"Don't know what how to handle {element_id} in " + "condition table." ) return value @@ -560,7 +584,9 @@ def _get_initial_state_pysb( ( initial.value for initial in petab_problem.model.model.initials - if match_complex_pattern(initial.pattern, species_pattern, exact=True) + if match_complex_pattern( + initial.pattern, species_pattern, exact=True + ) ), 0.0, ) @@ -617,9 +643,9 @@ def _set_initial_state( scale_map[init_par_id] = petab.LIN else: # parametric initial state - scale_map[init_par_id] = petab_problem.parameter_df[PARAMETER_SCALE].get( - value, petab.LIN - ) + scale_map[init_par_id] = petab_problem.parameter_df[ + PARAMETER_SCALE + ].get(value, petab.LIN) def create_parameter_mapping_for_condition( @@ -657,13 +683,17 @@ def create_parameter_mapping_for_condition( condition_map_sim ) != len(condition_scale_map_sim): raise AssertionError( - "Number of parameters and number of parameter " "scales do not match." + "Number of parameters and number of parameter " + "scales do not match." ) - if len(condition_map_preeq) and len(condition_map_preeq) != len(condition_map_sim): + if len(condition_map_preeq) and len(condition_map_preeq) != len( + condition_map_sim + ): logger.debug(f"Preequilibration parameter map: {condition_map_preeq}") logger.debug(f"Simulation parameter map: {condition_map_sim}") raise AssertionError( - "Number of parameters for preequilbration " "and simulation do not match." + "Number of parameters for preequilbration " + "and simulation do not match." ) ########################################################################## @@ -745,7 +775,9 @@ def create_parameter_mapping_for_condition( ( condition_scale_map_preeq_var, condition_scale_map_preeq_fix, - ) = _subset_dict(condition_scale_map_preeq, variable_par_ids, fixed_par_ids) + ) = _subset_dict( + condition_scale_map_preeq, variable_par_ids, fixed_par_ids + ) condition_map_sim_var, condition_map_sim_fix = _subset_dict( condition_map_sim, variable_par_ids, fixed_par_ids @@ -755,9 +787,13 @@ def create_parameter_mapping_for_condition( condition_scale_map_sim, variable_par_ids, fixed_par_ids ) - logger.debug("Fixed parameters preequilibration: " f"{condition_map_preeq_fix}") + logger.debug( + "Fixed parameters preequilibration: " f"{condition_map_preeq_fix}" + ) logger.debug("Fixed parameters simulation: " f"{condition_map_sim_fix}") - logger.debug("Variable parameters preequilibration: " f"{condition_map_preeq_var}") + logger.debug( + "Variable parameters preequilibration: " f"{condition_map_preeq_var}" + ) logger.debug("Variable parameters simulation: " f"{condition_map_sim_var}") petab.merge_preeq_and_sim_pars_condition( @@ -879,7 +915,10 @@ def create_edata_for_condition( states_in_condition_table = get_states_in_condition_table( petab_problem, condition=condition ) - if condition.get(PREEQUILIBRATION_CONDITION_ID) and states_in_condition_table: + if ( + condition.get(PREEQUILIBRATION_CONDITION_ID) + and states_in_condition_table + ): state_ids = amici_model.getStateIds() state_idx_reinitalization = [ state_ids.index(s) @@ -898,7 +937,9 @@ def create_edata_for_condition( # timepoints # find replicate numbers of time points - timepoints_w_reps = _get_timepoints_with_replicates(df_for_condition=measurement_df) + timepoints_w_reps = _get_timepoints_with_replicates( + df_for_condition=measurement_df + ) edata.setTimepoints(timepoints_w_reps) ########################################################################## @@ -951,7 +992,9 @@ def _get_timepoints_with_replicates( timepoints_w_reps = [] for time in timepoints: # subselect for time - df_for_time = df_for_condition[df_for_condition.time.astype(float) == time] + df_for_time = df_for_condition[ + df_for_condition.time.astype(float) == time + ] # rep number is maximum over rep numbers for observables n_reps = max(df_for_time.groupby([OBSERVABLE_ID, TIME]).size()) # append time point n_rep times @@ -984,7 +1027,9 @@ def _get_measurements_and_sigmas( arrays for measurement and sigmas """ # prepare measurement matrix - y = np.full(shape=(len(timepoints_w_reps), len(observable_ids)), fill_value=np.nan) + y = np.full( + shape=(len(timepoints_w_reps), len(observable_ids)), fill_value=np.nan + ) # prepare sigma matrix sigma_y = y.copy() @@ -1013,10 +1058,12 @@ def _get_measurements_and_sigmas( y[time_ix_for_obs_ix[observable_ix], observable_ix] = measurement[ MEASUREMENT ] - if isinstance(measurement.get(NOISE_PARAMETERS, None), numbers.Number): - sigma_y[time_ix_for_obs_ix[observable_ix], observable_ix] = measurement[ - NOISE_PARAMETERS - ] + if isinstance( + measurement.get(NOISE_PARAMETERS, None), numbers.Number + ): + sigma_y[ + time_ix_for_obs_ix[observable_ix], observable_ix + ] = measurement[NOISE_PARAMETERS] return y, sigma_y @@ -1054,7 +1101,9 @@ def rdatas_to_measurement_df( t = list(rdata.ts) # extract rows for condition - cur_measurement_df = petab.get_rows_for_condition(measurement_df, condition) + cur_measurement_df = petab.get_rows_for_condition( + measurement_df, condition + ) # iterate over entries for the given condition # note: this way we only generate a dataframe entry for every diff --git a/python/sdist/amici/petab_simulate.py b/python/sdist/amici/petab_simulate.py index e7d66a71f4..32c1ef8955 100644 --- a/python/sdist/amici/petab_simulate.py +++ b/python/sdist/amici/petab_simulate.py @@ -53,7 +53,10 @@ def simulate_without_noise(self, **kwargs) -> pd.DataFrame: in the Simulator constructor (including the PEtab problem). """ if AMICI_MODEL in {*kwargs, *dir(self)} and ( - any(k in kwargs for k in inspect.signature(import_petab_problem).parameters) + any( + k in kwargs + for k in inspect.signature(import_petab_problem).parameters + ) ): print( "Arguments related to the PEtab import are unused if " diff --git a/python/sdist/amici/petab_util.py b/python/sdist/amici/petab_util.py index a9666d84ac..9108b108bc 100644 --- a/python/sdist/amici/petab_util.py +++ b/python/sdist/amici/petab_util.py @@ -27,7 +27,9 @@ def get_states_in_condition_table( raise NotImplementedError() species_check_funs = { - MODEL_TYPE_SBML: lambda x: _element_is_sbml_state(petab_problem.sbml_model, x), + MODEL_TYPE_SBML: lambda x: _element_is_sbml_state( + petab_problem.sbml_model, x + ), MODEL_TYPE_PYSB: lambda x: _element_is_pysb_pattern( petab_problem.model.model, x ), @@ -36,7 +38,9 @@ def get_states_in_condition_table( resolve_mapping(petab_problem.mapping_df, col): (None, None) if condition is None else ( - petab_problem.condition_df.loc[condition[SIMULATION_CONDITION_ID], col], + petab_problem.condition_df.loc[ + condition[SIMULATION_CONDITION_ID], col + ], petab_problem.condition_df.loc[ condition[PREEQUILIBRATION_CONDITION_ID], col ] @@ -60,7 +64,9 @@ def get_states_in_condition_table( pysb.bng.generate_equations(petab_problem.model.model) try: - spm = pysb.pattern.SpeciesPatternMatcher(model=petab_problem.model.model) + spm = pysb.pattern.SpeciesPatternMatcher( + model=petab_problem.model.model + ) except NotImplementedError as e: raise NotImplementedError( "Requires https://github.com/pysb/pysb/pull/570. " diff --git a/python/sdist/amici/pysb_import.py b/python/sdist/amici/pysb_import.py index 14bcc1d790..469df23a2c 100644 --- a/python/sdist/amici/pysb_import.py +++ b/python/sdist/amici/pysb_import.py @@ -265,8 +265,12 @@ def ode_model_from_pysb_importer( _process_pysb_parameters(model, ode, constant_parameters) if compute_conservation_laws: _process_pysb_conservation_laws(model, ode) - _process_pysb_observables(model, ode, observables, sigmas, noise_distributions) - _process_pysb_expressions(model, ode, observables, sigmas, noise_distributions) + _process_pysb_observables( + model, ode, observables, sigmas, noise_distributions + ) + _process_pysb_expressions( + model, ode, observables, sigmas, noise_distributions + ) ode._has_quadratic_nllh = not noise_distributions or all( noise_distr in ["normal", "lin-normal", "log-normal", "log10-normal"] for noise_distr in noise_distributions.values() @@ -392,7 +396,9 @@ def _process_pysb_species(pysb_model: pysb.Model, ode_model: DEModel) -> None: for ix, specie in enumerate(pysb_model.species): init = sp.sympify("0.0") for ic in pysb_model.odes.model.initials: - if pysb.pattern.match_complex_pattern(ic.pattern, specie, exact=True): + if pysb.pattern.match_complex_pattern( + ic.pattern, specie, exact=True + ): # we don't want to allow expressions in initial conditions if ic.value in pysb_model.expressions: init = pysb_model.expressions[ic.value.name].expand_expr() @@ -400,7 +406,9 @@ def _process_pysb_species(pysb_model: pysb.Model, ode_model: DEModel) -> None: init = ic.value ode_model.add_component( - DifferentialState(sp.Symbol(f"__s{ix}"), f"{specie}", init, xdot[ix]) + DifferentialState( + sp.Symbol(f"__s{ix}"), f"{specie}", init, xdot[ix] + ) ) logger.debug(f"Finished Processing PySB species ") @@ -474,7 +482,8 @@ def _process_pysb_expressions( include_derived=True ) | pysb_model.expressions_dynamic(include_derived=True): if any( - isinstance(symbol, pysb.Tag) for symbol in expr.expand_expr().free_symbols + isinstance(symbol, pysb.Tag) + for symbol in expr.expand_expr().free_symbols ): # we only need explicit instantiations of expressions with tags, # which are defined in the derived expressions. The abstract @@ -531,11 +540,15 @@ def _add_expression( :param ode_model: see :py:func:`_process_pysb_expressions` """ - ode_model.add_component(Expression(sym, name, _parse_special_functions(expr))) + ode_model.add_component( + Expression(sym, name, _parse_special_functions(expr)) + ) if name in observables: noise_dist = ( - noise_distributions.get(name, "normal") if noise_distributions else "normal" + noise_distributions.get(name, "normal") + if noise_distributions + else "normal" ) y = sp.Symbol(f"{name}") @@ -543,7 +556,9 @@ def _add_expression( obs = Observable(y, name, sym, transformation=trafo) ode_model.add_component(obs) - sigma_name, sigma_value = _get_sigma_name_and_value(pysb_model, name, sigmas) + sigma_name, sigma_value = _get_sigma_name_and_value( + pysb_model, name, sigmas + ) sigma = sp.Symbol(sigma_name) ode_model.add_component(SigmaY(sigma, f"{sigma_name}", sigma_value)) @@ -552,10 +567,14 @@ def _add_expression( my = generate_measurement_symbol(obs.get_id()) cost_fun_expr = sp.sympify( cost_fun_str, - locals=dict(zip(_get_str_symbol_identifiers(name), (y, my, sigma))), + locals=dict( + zip(_get_str_symbol_identifiers(name), (y, my, sigma)) + ), ) ode_model.add_component( - LogLikelihoodY(sp.Symbol(f"llh_{name}"), f"llh_{name}", cost_fun_expr) + LogLikelihoodY( + sp.Symbol(f"llh_{name}"), f"llh_{name}", cost_fun_expr + ) ) @@ -585,7 +604,9 @@ def _get_sigma_name_and_value( sigma_name = sigmas[obs_name] try: # find corresponding Expression instance - sigma_expr = next(x for x in pysb_model.expressions if x.name == sigma_name) + sigma_expr = next( + x for x in pysb_model.expressions if x.name == sigma_name + ) except StopIteration: raise ValueError( f"value of sigma {obs_name} is not a " f"valid expression." @@ -643,7 +664,9 @@ def _process_pysb_observables( @log_execution_time("computing PySB conservation laws", logger) -def _process_pysb_conservation_laws(pysb_model: pysb.Model, ode_model: DEModel) -> None: +def _process_pysb_conservation_laws( + pysb_model: pysb.Model, ode_model: DEModel +) -> None: """ Removes species according to conservation laws to ensure that the jacobian has full rank @@ -657,7 +680,9 @@ def _process_pysb_conservation_laws(pysb_model: pysb.Model, ode_model: DEModel) monomers_without_conservation_law = set() for rule in pysb_model.rules: - monomers_without_conservation_law |= _get_unconserved_monomers(rule, pysb_model) + monomers_without_conservation_law |= _get_unconserved_monomers( + rule, pysb_model + ) monomers_without_conservation_law |= ( _compute_monomers_with_fixed_initial_conditions(pysb_model) @@ -731,7 +756,9 @@ def _generate_cl_prototypes( """ cl_prototypes = dict() - _compute_possible_indices(cl_prototypes, pysb_model, ode_model, excluded_monomers) + _compute_possible_indices( + cl_prototypes, pysb_model, ode_model, excluded_monomers + ) _compute_dependency_idx(cl_prototypes) _compute_target_index(cl_prototypes, ode_model) @@ -839,7 +866,9 @@ def _compute_dependency_idx(cl_prototypes: CL_Prototype) -> None: prototype_j["dependency_idx"][idx] |= {monomer_i} -def _compute_target_index(cl_prototypes: CL_Prototype, ode_model: DEModel) -> None: +def _compute_target_index( + cl_prototypes: CL_Prototype, ode_model: DEModel +) -> None: """ Computes the target index for every monomer @@ -899,7 +928,9 @@ def _compute_target_index(cl_prototypes: CL_Prototype, ode_model: DEModel) -> No # multimers has a low upper bound and the species count does not # vary too much across conservation laws, this approximation # should be fine - prototype["fillin"] = prototype["appearance_count"] * prototype["species_count"] + prototype["fillin"] = ( + prototype["appearance_count"] * prototype["species_count"] + ) # we might end up with the same index for multiple monomers, so loop until # we have a set of unique target indices @@ -950,7 +981,9 @@ def _cl_has_cycle(monomer: str, cl_prototypes: CL_Prototype) -> bool: root = monomer return any( _is_in_cycle(connecting_monomer, cl_prototypes, visited, root) - for connecting_monomer in prototype["dependency_idx"][prototype["target_index"]] + for connecting_monomer in prototype["dependency_idx"][ + prototype["target_index"] + ] ) @@ -994,7 +1027,9 @@ def _is_in_cycle( return any( _is_in_cycle(connecting_monomer, cl_prototypes, visited, root) - for connecting_monomer in prototype["dependency_idx"][prototype["target_index"]] + for connecting_monomer in prototype["dependency_idx"][ + prototype["target_index"] + ] ) @@ -1011,9 +1046,9 @@ def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None: target_indices = _get_target_indices(cl_prototypes) for monomer, prototype in cl_prototypes.items(): - if target_indices.count(prototype["target_index"]) > 1 or _cl_has_cycle( - monomer, cl_prototypes - ): + if target_indices.count( + prototype["target_index"] + ) > 1 or _cl_has_cycle(monomer, cl_prototypes): # compute how much fillin the next best target_index would yield # we exclude already existing target indices to avoid that @@ -1022,7 +1057,9 @@ def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None: # solution but prevents infinite loops for target_index in list(set(target_indices)): try: - local_idx = prototype["possible_indices"].index(target_index) + local_idx = prototype["possible_indices"].index( + target_index + ) except ValueError: local_idx = None @@ -1037,13 +1074,16 @@ def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None: idx = np.argmin(prototype["appearance_counts"]) prototype["local_index"] = idx - prototype["alternate_target_index"] = prototype["possible_indices"][idx] - prototype["alternate_appearance_count"] = prototype["appearance_counts"][ - idx - ] + prototype["alternate_target_index"] = prototype[ + "possible_indices" + ][idx] + prototype["alternate_appearance_count"] = prototype[ + "appearance_counts" + ][idx] prototype["alternate_fillin"] = ( - prototype["alternate_appearance_count"] * prototype["species_count"] + prototype["alternate_appearance_count"] + * prototype["species_count"] ) prototype["diff_fillin"] = ( @@ -1052,13 +1092,18 @@ def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None: else: prototype["diff_fillin"] = -1 - if all(prototype["diff_fillin"] == -1 for prototype in cl_prototypes.values()): + if all( + prototype["diff_fillin"] == -1 for prototype in cl_prototypes.values() + ): raise RuntimeError( - "Could not compute a valid set of conservation " "laws for this model!" + "Could not compute a valid set of conservation " + "laws for this model!" ) # this puts prototypes with high diff_fillin last - cl_prototypes = sorted(cl_prototypes.items(), key=lambda kv: kv[1]["diff_fillin"]) + cl_prototypes = sorted( + cl_prototypes.items(), key=lambda kv: kv[1]["diff_fillin"] + ) cl_prototypes = {proto[0]: proto[1] for proto in cl_prototypes} for monomer in cl_prototypes: @@ -1072,12 +1117,15 @@ def _greedy_target_index_update(cl_prototypes: CL_Prototype) -> None: # are recomputed on the fly) if prototype["diff_fillin"] > -1 and ( - _get_target_indices(cl_prototypes).count(prototype["target_index"]) > 1 + _get_target_indices(cl_prototypes).count(prototype["target_index"]) + > 1 or _cl_has_cycle(monomer, cl_prototypes) ): prototype["fillin"] = prototype["alternate_fillin"] prototype["target_index"] = prototype["alternate_target_index"] - prototype["appearance_count"] = prototype["alternate_appearance_count"] + prototype["appearance_count"] = prototype[ + "alternate_appearance_count" + ] del prototype["possible_indices"][prototype["local_index"]] del prototype["appearance_counts"][prototype["local_index"]] @@ -1176,9 +1224,12 @@ def _flatten_conservation_laws( for cl in conservation_laws: # only update if we changed something if any( - _apply_conseration_law_sub(cl, sub) for sub in conservation_law_subs + _apply_conseration_law_sub(cl, sub) + for sub in conservation_law_subs ): - conservation_law_subs = _get_conservation_law_subs(conservation_laws) + conservation_law_subs = _get_conservation_law_subs( + conservation_laws + ) def _apply_conseration_law_sub( @@ -1289,7 +1340,9 @@ def has_fixed_parameter_ic( ( ic for ic, condition in enumerate(pysb_model.initials) - if pysb.pattern.match_complex_pattern(condition[0], specie, exact=True) + if pysb.pattern.match_complex_pattern( + condition[0], specie, exact=True + ) ), None, ) @@ -1322,7 +1375,9 @@ def extract_monomers( ] -def _get_unconserved_monomers(rule: pysb.Rule, pysb_model: pysb.Model) -> Set[str]: +def _get_unconserved_monomers( + rule: pysb.Rule, pysb_model: pysb.Model +) -> Set[str]: """ Constructs the set of monomer names for which the specified rule changes the stoichiometry of the monomer in the specified model. @@ -1338,11 +1393,16 @@ def _get_unconserved_monomers(rule: pysb.Rule, pysb_model: pysb.Model) -> Set[st """ unconserved_monomers = set() - if not rule.delete_molecules and len(rule.product_pattern.complex_patterns) == 0: + if ( + not rule.delete_molecules + and len(rule.product_pattern.complex_patterns) == 0 + ): # if delete_molecules is not True but we have a degradation rule, # we have to actually go through the reactions that are created by # the rule - for reaction in [r for r in pysb_model.reactions if rule.name in r["rule"]]: + for reaction in [ + r for r in pysb_model.reactions if rule.name in r["rule"] + ]: unconserved_monomers |= _get_changed_stoichiometries( [pysb_model.species[ix] for ix in reaction["reactants"]], [pysb_model.species[ix] for ix in reaction["products"]], @@ -1395,7 +1455,9 @@ def pysb_model_from_path(pysb_model_file: Union[str, Path]) -> pysb.Model: :return: The pysb Model instance """ - pysb_model_module_name = os.path.splitext(os.path.split(pysb_model_file)[-1])[0] + pysb_model_module_name = os.path.splitext( + os.path.split(pysb_model_file)[-1] + )[0] import importlib.util diff --git a/python/sdist/amici/sbml_import.py b/python/sdist/amici/sbml_import.py index 30cedca361..f6c46aa4a3 100644 --- a/python/sdist/amici/sbml_import.py +++ b/python/sdist/amici/sbml_import.py @@ -211,7 +211,9 @@ def _process_document(self) -> None: """ # Ensure we got a valid SBML model, otherwise further processing # might lead to undefined results - log_execution_time("validating SBML", logger)(self.sbml_doc.validateSBML)() + log_execution_time("validating SBML", logger)( + self.sbml_doc.validateSBML + )() _check_lib_sbml_errors(self.sbml_doc, self.show_sbml_warnings) # Flatten "comp" model? Do that before any other converters are run @@ -251,7 +253,9 @@ def _process_document(self) -> None: self.sbml_doc.convert )(convert_config) - convert_config = sbml.SBMLLocalParameterConverter().getDefaultProperties() + convert_config = ( + sbml.SBMLLocalParameterConverter().getDefaultProperties() + ) log_execution_time("converting SBML local parameters", logger)( self.sbml_doc.convert )(convert_config) @@ -471,7 +475,9 @@ def _build_ode_model( See :py:func:`sbml2amici` for parameters. """ - constant_parameters = list(constant_parameters) if constant_parameters else [] + constant_parameters = ( + list(constant_parameters) if constant_parameters else [] + ) hardcode_symbols = set(hardcode_symbols) if hardcode_symbols else {} if invalid := (set(constant_parameters) & set(hardcode_symbols)): @@ -494,7 +500,9 @@ def _build_ode_model( self._reset_symbols() self.sbml_parser_settings.setParseLog( - sbml.L3P_PARSE_LOG_AS_LOG10 if log_as_log10 else sbml.L3P_PARSE_LOG_AS_LN + sbml.L3P_PARSE_LOG_AS_LOG10 + if log_as_log10 + else sbml.L3P_PARSE_LOG_AS_LN ) self._process_sbml( constant_parameters=constant_parameters, @@ -531,7 +539,9 @@ def _build_ode_model( simplify=simplify, cache_simplify=cache_simplify, ) - ode_model.import_from_sbml_importer(self, compute_cls=compute_conservation_laws) + ode_model.import_from_sbml_importer( + self, compute_cls=compute_conservation_laws + ) return ode_model @log_execution_time("importing SBML", logger) @@ -583,7 +593,10 @@ def check_support(self) -> None: # the "required" attribute is only available in SBML Level 3 for i_plugin in range(self.sbml.getNumPlugins()): plugin = self.sbml.getPlugin(i_plugin) - if self.sbml_doc.getPkgRequired(plugin.getPackageName()) is False: + if ( + self.sbml_doc.getPkgRequired(plugin.getPackageName()) + is False + ): # if not "required", this has no impact on model # simulation, and we can safely ignore it @@ -599,7 +612,9 @@ def check_support(self) -> None: raise SBMLException( "The following fbc extension elements are " "currently not supported: " - + ", ".join(list(map(str, plugin.getListOfAllElements()))) + + ", ".join( + list(map(str, plugin.getListOfAllElements())) + ) ) continue @@ -671,7 +686,8 @@ def check_event_support(self) -> None: trigger_sbml = event.getTrigger() if trigger_sbml is None: logger.warning( - f"Event {event_id} trigger has no trigger, " "so will be skipped." + f"Event {event_id} trigger has no trigger, " + "so will be skipped." ) continue if trigger_sbml.getMath() is None: @@ -698,7 +714,9 @@ def _gather_locals(self, hardcode_symbols: Sequence[str] = None) -> None: self._gather_base_locals(hardcode_symbols=hardcode_symbols) self._gather_dependent_locals() - def _gather_base_locals(self, hardcode_symbols: Sequence[str] = None) -> None: + def _gather_base_locals( + self, hardcode_symbols: Sequence[str] = None + ) -> None: """ Populate self.local_symbols with pure symbol definitions that do not depend on any other symbol. @@ -743,7 +761,10 @@ def _gather_base_locals(self, hardcode_symbols: Sequence[str] = None) -> None: for x_ref in _get_list_of_species_references(self.sbml): if not x_ref.isSetId(): continue - if x_ref.isSetStoichiometry() and not self.is_assignment_rule_target(x_ref): + if ( + x_ref.isSetStoichiometry() + and not self.is_assignment_rule_target(x_ref) + ): value = sp.Float(x_ref.getStoichiometry()) else: value = _get_identifier_symbol(x_ref) @@ -822,7 +843,9 @@ def _process_species(self) -> None: Get species information from SBML model. """ if self.sbml.isSetConversionFactor(): - conversion_factor = symbol_with_assumptions(self.sbml.getConversionFactor()) + conversion_factor = symbol_with_assumptions( + self.sbml.getConversionFactor() + ) else: conversion_factor = 1 @@ -834,7 +857,9 @@ def _process_species(self) -> None: "compartment": _get_species_compartment_symbol(s), "constant": s.getConstant() or s.getBoundaryCondition(), "amount": s.getHasOnlySubstanceUnits(), - "conversion_factor": symbol_with_assumptions(s.getConversionFactor()) + "conversion_factor": symbol_with_assumptions( + s.getConversionFactor() + ) if s.isSetConversionFactor() else conversion_factor, "index": len(self.symbols[SymbolId.SPECIES]), @@ -859,7 +884,9 @@ def _process_species_initial(self): # targets to have InitialAssignments. species = self.symbols[SymbolId.SPECIES].get(species_id, None) - ia_initial = self._get_element_initial_assignment(species_variable.getId()) + ia_initial = self._get_element_initial_assignment( + species_variable.getId() + ) if ia_initial is not None: initial = ia_initial if species: @@ -872,7 +899,9 @@ def _process_species_initial(self): all_rateof_dummies.append(rateof_dummies) # don't assign this since they need to stay in order - sorted_species = toposort_symbols(self.symbols[SymbolId.SPECIES], "init") + sorted_species = toposort_symbols( + self.symbols[SymbolId.SPECIES], "init" + ) for species, rateof_dummies in zip( self.symbols[SymbolId.SPECIES].values(), all_rateof_dummies ): @@ -965,7 +994,9 @@ def add_d_dt( variable0 = smart_subs(variable0, species_id, species["init"]) for species in self.symbols[SymbolId.SPECIES].values(): - species["init"] = smart_subs(species["init"], variable, variable0) + species["init"] = smart_subs( + species["init"], variable, variable0 + ) # add compartment/parameter species self.symbols[SymbolId.SPECIES][variable] = { @@ -1032,7 +1063,8 @@ def _process_parameters( ] for parameter in fixed_parameters: if ( - self._get_element_initial_assignment(parameter.getId()) is not None + self._get_element_initial_assignment(parameter.getId()) + is not None or self.is_assignment_rule_target(parameter) or self.is_rate_rule_target(parameter) ): @@ -1073,8 +1105,12 @@ def _process_parameters( for par in self.sbml.getListOfParameters(): if ( ia := self._get_element_initial_assignment(par.getId()) - ) is not None and ia.find(sp.core.function.UndefinedFunction("rateOf")): - self.symbols[SymbolId.EXPRESSION][_get_identifier_symbol(par)] = { + ) is not None and ia.find( + sp.core.function.UndefinedFunction("rateOf") + ): + self.symbols[SymbolId.EXPRESSION][ + _get_identifier_symbol(par) + ] = { "name": par.getName() if par.isSetName() else par.getId(), "value": ia, } @@ -1127,9 +1163,9 @@ def _process_reactions(self): # rate of change in species concentration) now occurs # in the `dx_dt` method in "de_export.py", which also # accounts for possibly variable compartments. - self.stoichiometric_matrix[species["index"], reaction_index] += ( - sign * stoichiometry * species["conversion_factor"] - ) + self.stoichiometric_matrix[ + species["index"], reaction_index + ] += (sign * stoichiometry * species["conversion_factor"]) if reaction.isSetId(): sym_math = self._local_symbols[reaction.getId()] else: @@ -1205,9 +1241,9 @@ def _process_rule_algebraic(self, rule: sbml.AlgebraicRule): continue # and there must also not be a rate rule or assignment # rule for it - if self.is_assignment_rule_target(sbml_var) or self.is_rate_rule_target( + if self.is_assignment_rule_target( sbml_var - ): + ) or self.is_rate_rule_target(sbml_var): continue # Furthermore, if the entity is a Species object, its value # must not be determined by reactions, which means that it @@ -1225,7 +1261,11 @@ def _process_rule_algebraic(self, rule: sbml.AlgebraicRule): :, ] ) - if is_species and not is_boundary_condition and is_involved_in_reaction: + if ( + is_species + and not is_boundary_condition + and is_involved_in_reaction + ): continue free_variables.add(symbol) @@ -1275,14 +1315,22 @@ def _process_rule_algebraic(self, rule: sbml.AlgebraicRule): symbol["init"] = sp.Float(symbol.pop("value")) # if not a species, add a zeros row to the stoichiometric # matrix - if (isinstance(symbol["init"], float) and np.isnan(symbol["init"])) or ( - isinstance(symbol["init"], sp.Number) and symbol["init"] == sp.nan + if ( + isinstance(symbol["init"], float) + and np.isnan(symbol["init"]) + ) or ( + isinstance(symbol["init"], sp.Number) + and symbol["init"] == sp.nan ): # placeholder, needs to be determined in IC calculation symbol["init"] = sp.Float(0.0) - self.stoichiometric_matrix = self.stoichiometric_matrix.row_insert( - self.stoichiometric_matrix.shape[0], - sp.SparseMatrix([[0] * self.stoichiometric_matrix.shape[1]]), + self.stoichiometric_matrix = ( + self.stoichiometric_matrix.row_insert( + self.stoichiometric_matrix.shape[0], + sp.SparseMatrix( + [[0] * self.stoichiometric_matrix.shape[1]] + ), + ) ) elif var_ix != self.stoichiometric_matrix.shape[0] - 1: # if not the last col, move it to the end @@ -1358,7 +1406,9 @@ def _convert_event_assignment_parameter_targets_to_species(self): This is for the convenience of only implementing event assignments for "species". """ - parameter_targets = _collect_event_assignment_parameter_targets(self.sbml) + parameter_targets = _collect_event_assignment_parameter_targets( + self.sbml + ) for parameter_target in parameter_targets: # Parameter rate rules already exist as species. if parameter_target in self.symbols[SymbolId.SPECIES]: @@ -1379,7 +1429,9 @@ def _convert_event_assignment_parameter_targets_to_species(self): "Unexpected error. The parameter target of an " "event assignment was processed twice." ) - parameter_def = self.symbols[symbol_id].pop(parameter_target) + parameter_def = self.symbols[symbol_id].pop( + parameter_target + ) if parameter_def is None: # this happens for parameters that have initial assignments # or are assignment rule targets @@ -1387,7 +1439,9 @@ def _convert_event_assignment_parameter_targets_to_species(self): ia_init = self._get_element_initial_assignment(par.getId()) parameter_def = { "name": par.getName() if par.isSetName() else par.getId(), - "value": sp.Float(par.getValue()) if ia_init is None else ia_init, + "value": sp.Float(par.getValue()) + if ia_init is None + else ia_init, } # Fixed parameters are added as species such that they can be # targets of events. @@ -1428,9 +1482,9 @@ def get_empty_bolus_value() -> sp.Float: # Species has a compartment "compartment" in species_def ): - concentration_species_by_compartment[species_def["compartment"]].append( - species - ) + concentration_species_by_compartment[ + species_def["compartment"] + ].append(species) for ievent, event in enumerate(events): # get the event id (which is optional unfortunately) @@ -1453,7 +1507,9 @@ def get_empty_bolus_value() -> sp.Float: event_assignments = event.getListOfEventAssignments() compartment_event_assignments = set() for event_assignment in event_assignments: - variable_sym = symbol_with_assumptions(event_assignment.getVariable()) + variable_sym = symbol_with_assumptions( + event_assignment.getVariable() + ) if event_assignment.getMath() is None: # Ignore event assignments with no change in value. continue @@ -1518,10 +1574,14 @@ def get_empty_bolus_value() -> sp.Float: for index in range(len(bolus)): if bolus[index] != get_empty_bolus_value(): bolus[index] -= state_vector[index] - bolus[index] = bolus[index].subs(get_empty_bolus_value(), sp.Float(0.0)) + bolus[index] = bolus[index].subs( + get_empty_bolus_value(), sp.Float(0.0) + ) initial_value = ( - trigger_sbml.getInitialValue() if trigger_sbml is not None else True + trigger_sbml.getInitialValue() + if trigger_sbml is not None + else True ) if self.symbols[SymbolId.ALGEBRAIC_EQUATION] and not initial_value: # in principle this could be implemented, requires running @@ -1567,7 +1627,9 @@ def _process_observables( See :py:func:`sbml2amici`. """ - _validate_observables(observables, sigmas, noise_distributions, events=False) + _validate_observables( + observables, sigmas, noise_distributions, events=False + ) # add user-provided observables or make all species, and compartments # with assignment rules, observable @@ -1589,7 +1651,9 @@ def _process_observables( # check for nesting of observables (unsupported) observable_syms = set(self.symbols[SymbolId.OBSERVABLE].keys()) for obs in self.symbols[SymbolId.OBSERVABLE].values(): - if any(sym in observable_syms for sym in obs["value"].free_symbols): + if any( + sym in observable_syms for sym in obs["value"].free_symbols + ): raise ValueError( "Nested observables are not supported, " f"but observable `{obs['name']} = {obs['value']}` " @@ -1598,7 +1662,9 @@ def _process_observables( elif observables is None: self._generate_default_observables() - _check_symbol_nesting(self.symbols[SymbolId.OBSERVABLE], "eventObservable") + _check_symbol_nesting( + self.symbols[SymbolId.OBSERVABLE], "eventObservable" + ) self._process_log_likelihood(sigmas, noise_distributions) @@ -1636,7 +1702,10 @@ def _process_event_observables( for obs, definition in event_observables.items(): self.add_local_symbol(obs, symbol_with_assumptions(obs)) # check corresponding event exists - if sp.Symbol(definition["event"]) not in self.symbols[SymbolId.EVENT]: + if ( + sp.Symbol(definition["event"]) + not in self.symbols[SymbolId.EVENT] + ): raise ValueError( "Could not find an event with the event identifier " f'{definition["event"]} for the event observable with name' @@ -1769,7 +1838,9 @@ def _process_log_likelihood( self.symbols[sigma_symbol] = { symbol_with_assumptions(f"sigma_{obs_id}"): { "name": f'sigma_{obs["name"]}', - "value": self._sympy_from_sbml_math(sigmas.get(str(obs_id), "1.0")), + "value": self._sympy_from_sbml_math( + sigmas.get(str(obs_id), "1.0") + ), } for obs_id, obs in self.symbols[obs_symbol].items() } @@ -1820,7 +1891,9 @@ def _process_initial_assignments(self): continue sym_math = self._make_initial( - smart_subs_dict(sym_math, self.symbols[SymbolId.EXPRESSION], "value") + smart_subs_dict( + sym_math, self.symbols[SymbolId.EXPRESSION], "value" + ) ) self.initial_assignments[_get_identifier_symbol(ia)] = sym_math @@ -1883,7 +1956,9 @@ def _make_initial( if "init" in species: sym_math = smart_subs(sym_math, species_id, species["init"]) - sym_math = smart_subs(sym_math, self._local_symbols["time"], sp.Float(0)) + sym_math = smart_subs( + sym_math, self._local_symbols["time"], sp.Float(0) + ) sym_math = _dummy_to_rateof(sym_math, rateof_to_dummy) @@ -1919,7 +1994,9 @@ def process_conservation_laws(self, ode_model) -> None: # add algebraic variables to species_solver as they were ignored above ndifferential = len(ode_model._differential_states) nalgebraic = len(ode_model._algebraic_states) - species_solver.extend(list(range(ndifferential, ndifferential + nalgebraic))) + species_solver.extend( + list(range(ndifferential, ndifferential + nalgebraic)) + ) # Check, whether species_solver is empty now. As currently, AMICI # cannot handle ODEs without species, CLs must be switched off in this @@ -1929,7 +2006,9 @@ def process_conservation_laws(self, ode_model) -> None: species_solver = list(range(ode_model.num_states_rdata())) # prune out species from stoichiometry and - self.stoichiometric_matrix = self.stoichiometric_matrix[species_solver, :] + self.stoichiometric_matrix = self.stoichiometric_matrix[ + species_solver, : + ] # add the found CLs to the ode_model for cl in conservation_laws: @@ -1953,7 +2032,9 @@ def _get_conservation_laws_demartino( compute_moiety_conservation_laws, ) - sm = self.stoichiometric_matrix[: len(self.symbols[SymbolId.SPECIES]), :] + sm = self.stoichiometric_matrix[ + : len(self.symbols[SymbolId.SPECIES]), : + ] try: stoichiometric_list = [float(entry) for entry in sm.T.flat()] @@ -1976,7 +2057,9 @@ def _get_conservation_laws_demartino( stoichiometric_list, *sm.shape, rng_seed=32, - species_names=[str(x.get_id()) for x in ode_model._differential_states], + species_names=[ + str(x.get_id()) for x in ode_model._differential_states + ], ) # Sparsify conserved quantities @@ -1988,7 +2071,9 @@ def _get_conservation_laws_demartino( # `A * x0 = total_cl` and bring it to reduced row echelon form. The # pivot species are the ones to be eliminated. The resulting state # expressions are sparse and void of any circular dependencies. - A = sp.zeros(len(cls_coefficients), len(ode_model._differential_states)) + A = sp.zeros( + len(cls_coefficients), len(ode_model._differential_states) + ) for i_cl, (cl, coefficients) in enumerate( zip(cls_state_idxs, cls_coefficients) ): @@ -2025,7 +2110,9 @@ def _get_conservation_laws_rref( try: S = np.asarray( - self.stoichiometric_matrix[: len(self.symbols[SymbolId.SPECIES]), :], + self.stoichiometric_matrix[ + : len(self.symbols[SymbolId.SPECIES]), : + ], dtype=float, ) except TypeError: @@ -2221,7 +2308,8 @@ def _replace_in_all_expressions( if old not in self.symbols[symbol]: continue self.symbols[symbol] = { - smart_subs(k, old, new): v for k, v in self.symbols[symbol].items() + smart_subs(k, old, new): v + for k, v in self.symbols[symbol].items() } for symbol in [ @@ -2261,7 +2349,9 @@ def _replace_in_all_expressions( **self.symbols[SymbolId.SPECIES], **self.symbols[SymbolId.ALGEBRAIC_STATE], }.values(): - state["init"] = smart_subs(state["init"], old, self._make_initial(new)) + state["init"] = smart_subs( + state["init"], old, self._make_initial(new) + ) if "dt" in state: state["dt"] = smart_subs(state["dt"], old, new) @@ -2346,10 +2436,14 @@ def _sympy_from_sbml_math( if isinstance(formula, sp.Expr): formula = _parse_special_functions_sbml(formula) - _check_unsupported_functions_sbml(formula, expression_type=ele_name) + _check_unsupported_functions_sbml( + formula, expression_type=ele_name + ) return formula - def _get_element_initial_assignment(self, element_id: str) -> Union[sp.Expr, None]: + def _get_element_initial_assignment( + self, element_id: str + ) -> Union[sp.Expr, None]: """ Extract value of sbml variable according to its initial assignment @@ -2444,7 +2538,8 @@ def _check_lib_sbml_errors( error = sbml_doc.getError(i_error) # we ignore any info messages for now if error.getSeverity() >= sbml.LIBSBML_SEV_ERROR or ( - show_warnings and error.getSeverity() >= sbml.LIBSBML_SEV_WARNING + show_warnings + and error.getSeverity() >= sbml.LIBSBML_SEV_WARNING ): logger.error( f"libSBML {error.getCategoryAsString()} " @@ -2453,7 +2548,9 @@ def _check_lib_sbml_errors( ) if num_error + num_fatal: - raise SBMLException("SBML Document failed to load (see error messages above)") + raise SBMLException( + "SBML Document failed to load (see error messages above)" + ) def _parse_event_trigger(trigger: sp.Expr) -> sp.Expr: @@ -2704,7 +2801,9 @@ def _check_unsupported_functions_sbml( raise SBMLException(str(err)) -def _parse_special_functions_sbml(sym: sp.Expr, toplevel: bool = True) -> sp.Expr: +def _parse_special_functions_sbml( + sym: sp.Expr, toplevel: bool = True +) -> sp.Expr: try: return _parse_special_functions(sym, toplevel) except RuntimeError as err: diff --git a/python/sdist/amici/sbml_utils.py b/python/sdist/amici/sbml_utils.py index cce2a6c4fa..66c9d01bbc 100644 --- a/python/sdist/amici/sbml_utils.py +++ b/python/sdist/amici/sbml_utils.py @@ -161,7 +161,9 @@ def add_species( ) compartment_id = compartments[0].getId() elif not model.getCompartment(compartment_id): - raise SbmlMissingComponentIdError(f"No compartment with ID {compartment_id}.") + raise SbmlMissingComponentIdError( + f"No compartment with ID {compartment_id}." + ) sp = model.createSpecies() if sp.setIdAttribute(species_id) != libsbml.LIBSBML_OPERATION_SUCCESS: @@ -532,6 +534,8 @@ def _parse_logical_operators( return math_str if " xor(" in math_str or " Xor(" in math_str: - raise SBMLException("Xor is currently not supported as logical " "operation.") + raise SBMLException( + "Xor is currently not supported as logical " "operation." + ) return (math_str.replace("&&", "&")).replace("||", "|") diff --git a/python/sdist/amici/splines.py b/python/sdist/amici/splines.py index 4d7f0a318a..fdb0912045 100644 --- a/python/sdist/amici/splines.py +++ b/python/sdist/amici/splines.py @@ -122,19 +122,25 @@ def __init__( stop = sp.nsimplify(sp.sympify(stop)) if step is None: if number_of_nodes is None: - raise ValueError("One of step/number_of_nodes must be specified!") + raise ValueError( + "One of step/number_of_nodes must be specified!" + ) if not isinstance(number_of_nodes, Integral): raise TypeError("Length must be an integer!") if number_of_nodes < 2: raise ValueError("Length must be at least 2!") step = (stop - start) / (number_of_nodes - 1) elif number_of_nodes is not None: - raise ValueError("Only one of step/number_of_nodes can be specified!") + raise ValueError( + "Only one of step/number_of_nodes can be specified!" + ) else: step = sp.nsimplify(sp.sympify(step)) if start > stop: - raise ValueError(f"Start point {start} greater than stop point {stop}!") + raise ValueError( + f"Start point {start} greater than stop point {stop}!" + ) if step <= 0: raise ValueError(f"Step size {step} must be strictly positive!") @@ -191,7 +197,8 @@ def __array__(self, dtype=None) -> np.ndarray: def __repr__(self) -> str: return ( - f"UniformGrid(start={self.start}, stop={self.stop}, " f"step={self.step})" + f"UniformGrid(start={self.start}, stop={self.stop}, " + f"step={self.step})" ) @@ -311,7 +318,10 @@ def __init__( if not isinstance(evaluate_at, sp.Basic): # It may still be e.g. a list! raise ValueError(f"Invalid evaluate_at = {evaluate_at}!") - if evaluate_at != amici_time_symbol and evaluate_at != sbml_time_symbol: + if ( + evaluate_at != amici_time_symbol + and evaluate_at != sbml_time_symbol + ): logger.warning( "At the moment AMICI only supports evaluate_at = (model time). " "Annotations with correct piecewise MathML formulas " @@ -321,7 +331,9 @@ def __init__( if not isinstance(nodes, UniformGrid): nodes = np.asarray([sympify_noeval(x) for x in nodes]) - values_at_nodes = np.asarray([sympify_noeval(y) for y in values_at_nodes]) + values_at_nodes = np.asarray( + [sympify_noeval(y) for y in values_at_nodes] + ) if len(nodes) != len(values_at_nodes): raise ValueError( @@ -343,7 +355,10 @@ def __init__( ) bc, extrapolate = self._normalize_bc_and_extrapolate(bc, extrapolate) - if bc == ("periodic", "periodic") and values_at_nodes[0] != values_at_nodes[-1]: + if ( + bc == ("periodic", "periodic") + and values_at_nodes[0] != values_at_nodes[-1] + ): raise ValueError( "If the spline is to be periodic, " "the first and last elements of values_at_nodes must be equal!" @@ -564,15 +579,21 @@ def check_if_valid(self, importer: sbml_import.SbmlImporter) -> None: fixed_parameters: List[sp.Symbol] = list( importer.symbols[SymbolId.FIXED_PARAMETER].keys() ) - species: List[sp.Symbol] = list(importer.symbols[SymbolId.SPECIES].keys()) + species: List[sp.Symbol] = list( + importer.symbols[SymbolId.SPECIES].keys() + ) for x in self.nodes: if not x.free_symbols.issubset(fixed_parameters): - raise ValueError("nodes should only depend on constant parameters!") + raise ValueError( + "nodes should only depend on constant parameters!" + ) for y in self.values_at_nodes: if y.free_symbols.intersection(species): - raise ValueError("values_at_nodes should not depend on model species!") + raise ValueError( + "values_at_nodes should not depend on model species!" + ) fixed_parameters_values = [ importer.symbols[SymbolId.FIXED_PARAMETER][fp]["value"] @@ -585,7 +606,9 @@ def check_if_valid(self, importer: sbml_import.SbmlImporter) -> None: if not np.all(np.diff(nodes_values) >= 0): raise ValueError("nodes should be strictly increasing!") - def poly(self, i: Integral, *, x: Union[Real, sp.Basic] = None) -> sp.Basic: + def poly( + self, i: Integral, *, x: Union[Real, sp.Basic] = None + ) -> sp.Basic: """ Get the polynomial interpolant on the ``(nodes[i], nodes[i+1])`` interval. The polynomial is written in Horner form with respect to the scaled @@ -633,7 +656,9 @@ def poly_variable(self, x: Union[Real, sp.Basic], i: Integral) -> sp.Basic: return self._poly_variable(x, i) @abstractmethod - def _poly_variable(self, x: Union[Real, sp.Basic], i: Integral) -> sp.Basic: + def _poly_variable( + self, x: Union[Real, sp.Basic], i: Integral + ) -> sp.Basic: """This function (and not poly_variable) should be implemented by the subclasses""" raise NotImplementedError() @@ -786,13 +811,17 @@ def _formula( x = self._to_base_interval(x) extr_left, extr_right = None, None else: - extr_left, extr_right = self._extrapolation_formulas(x, extrapolate) + extr_left, extr_right = self._extrapolation_formulas( + x, extrapolate + ) if extr_left is not None: pieces.append((extr_left, x < self.nodes[0])) for i in range(len(self.nodes) - 2): - pieces.append((self.segment_formula(i, x=x), x < self.nodes[i + 1])) + pieces.append( + (self.segment_formula(i, x=x), x < self.nodes[i + 1]) + ) if extr_right is not None: pieces.append((self.segment_formula(-1, x=x), x < self.nodes[-1])) @@ -828,7 +857,9 @@ def _to_base_interval( """For periodic splines, maps the real point `x` to the reference period.""" if self.bc != ("periodic", "periodic"): - raise ValueError("_to_base_interval makes no sense with non-periodic bc") + raise ValueError( + "_to_base_interval makes no sense with non-periodic bc" + ) xA = self.nodes[0] xB = self.nodes[-1] @@ -871,7 +902,9 @@ def squared_L2_norm_of_curvature(self) -> sp.Basic: integral = sp.sympify(0) for i in range(len(self.nodes) - 1): formula = self.poly(i, x=x).diff(x, 2) ** 2 - integral += sp.integrate(formula, (x, self.nodes[i], self.nodes[i + 1])) + integral += sp.integrate( + formula, (x, self.nodes[i], self.nodes[i + 1]) + ) return sp.simplify(integral) def integrate( @@ -902,7 +935,9 @@ def integrate( return formula.integrate((x, z0, z1)) if k0 + 1 == k1: - return formula.integrate((x, z0, xB)) + formula.integrate((x, xA, z1)) + return formula.integrate((x, z0, xB)) + formula.integrate( + (x, xA, z1) + ) return ( formula.integrate((x, z0, xB)) @@ -979,7 +1014,9 @@ def _annotation_children(self) -> Dict[str, Union[str, List[str]]]: assert amici_time_symbol not in x.free_symbols children["spline_grid"] = [sbml_mathml(x) for x in self.nodes] - children["spline_values"] = [sbml_mathml(y) for y in self.values_at_nodes] + children["spline_values"] = [ + sbml_mathml(y) for y in self.values_at_nodes + ] return children @@ -1042,10 +1079,12 @@ def add_to_sbml_model( # Autoadd parameters if auto_add is True or auto_add == "spline": - if not model.getParameter(str(self.sbml_id)) and not model.getSpecies( + if not model.getParameter( str(self.sbml_id) - ): - add_parameter(model, self.sbml_id, constant=False, units=y_units) + ) and not model.getSpecies(str(self.sbml_id)): + add_parameter( + model, self.sbml_id, constant=False, units=y_units + ) if auto_add is True: if isinstance(x_nominal, collections.abc.Sequence): @@ -1056,7 +1095,9 @@ def add_to_sbml_model( ) for i in range(len(x_nominal) - 1): if x[i] >= x[i + 1]: - raise ValueError("x_nominal must be strictly increasing!") + raise ValueError( + "x_nominal must be strictly increasing!" + ) elif x_nominal is None: x_nominal = len(self.nodes) * [None] else: @@ -1065,7 +1106,9 @@ def add_to_sbml_model( raise TypeError("x_nominal must be a Sequence!") for _x, _val in zip(self.nodes, x_nominal): if _x.is_Symbol and not model.getParameter(_x.name): - add_parameter(model, _x.name, value=_val, units=x_units) + add_parameter( + model, _x.name, value=_val, units=x_units + ) if isinstance(y_nominal, collections.abc.Sequence): if len(y_nominal) != len(self.values_at_nodes): @@ -1119,7 +1162,9 @@ def add_to_sbml_model( k = sp.Piecewise((3, sp.cos(s) < 0), (1, True)) formula = x0 + T * (sp.atan(sp.tan(s)) / (2 * sp.pi) + k / 4) assert amici_time_symbol not in formula.free_symbols - par = add_parameter(model, parameter_id, constant=False, units=x_units) + par = add_parameter( + model, parameter_id, constant=False, units=x_units + ) retcode = par.setAnnotation( f'' ) @@ -1127,13 +1172,17 @@ def add_to_sbml_model( raise SbmlAnnotationError("Could not set SBML annotation!") add_assignment_rule(model, parameter_id, formula) - def _replace_in_all_expressions(self, old: sp.Symbol, new: sp.Symbol) -> None: + def _replace_in_all_expressions( + self, old: sp.Symbol, new: sp.Symbol + ) -> None: if self.sbml_id == old: self._sbml_id = new self._x = self.evaluate_at.subs(old, new) if not isinstance(self.nodes, UniformGrid): self._nodes = [x.subs(old, new) for x in self.nodes] - self._values_at_nodes = [y.subs(old, new) for y in self.values_at_nodes] + self._values_at_nodes = [ + y.subs(old, new) for y in self.values_at_nodes + ] @staticmethod def is_spline(rule: libsbml.AssignmentRule) -> bool: @@ -1179,7 +1228,9 @@ def from_annotation( must be hard-coded into this function here (at the moment). """ if annotation.tag != f"{{{annotation_namespace}}}spline": - raise ValueError("The given annotation is not an AMICI spline annotation!") + raise ValueError( + "The given annotation is not an AMICI spline annotation!" + ) attributes = {} for key, value in annotation.items(): @@ -1215,14 +1266,17 @@ def from_annotation( if attributes["spline_method"] == "cubic_hermite": cls = CubicHermiteSpline else: - raise ValueError(f"Unknown spline method {attributes['spline_method']}!") + raise ValueError( + f"Unknown spline method {attributes['spline_method']}!" + ) del attributes["spline_method"] kwargs = cls._from_annotation(attributes, children) if attributes: raise ValueError( - "Unprocessed attributes in spline annotation!\n" + str(attributes) + "Unprocessed attributes in spline annotation!\n" + + str(attributes) ) if children: @@ -1293,7 +1347,9 @@ def _from_annotation( ) if "spline_values" not in children: - raise ValueError("Required spline annotation 'spline_values' missing!") + raise ValueError( + "Required spline annotation 'spline_values' missing!" + ) kwargs["values_at_nodes"] = children.pop("spline_values") return kwargs @@ -1312,7 +1368,9 @@ def _parameters(self) -> Set[sp.Symbol]: parameters.update(y.free_symbols) return parameters - def ode_model_symbol(self, importer: sbml_import.SbmlImporter) -> sp.Function: + def ode_model_symbol( + self, importer: sbml_import.SbmlImporter + ) -> sp.Function: """ Returns the `sympy` object to be used by :py:class:`amici.de_export.ODEModel`. @@ -1412,7 +1470,9 @@ def plot( nodes = np.asarray(self.nodes) xlim = (float(nodes[0]), float(nodes[-1])) nodes = np.linspace(*xlim, npoints) - ax.plot(nodes, [float(self.evaluate(x).subs(parameters)) for x in nodes]) + ax.plot( + nodes, [float(self.evaluate(x).subs(parameters)) for x in nodes] + ) ax.plot( self.nodes, [float(y.subs(parameters)) for y in self.values_at_nodes], @@ -1517,7 +1577,9 @@ def __init__( if not isinstance(nodes, UniformGrid): nodes = np.asarray([sympify_noeval(x) for x in nodes]) - values_at_nodes = np.asarray([sympify_noeval(y) for y in values_at_nodes]) + values_at_nodes = np.asarray( + [sympify_noeval(y) for y in values_at_nodes] + ) if len(nodes) != len(values_at_nodes): # NB this would be checked in AbstractSpline.__init__() @@ -1529,13 +1591,19 @@ def __init__( ) bc, extrapolate = self._normalize_bc_and_extrapolate(bc, extrapolate) - if bc[0] == "zeroderivative+natural" or bc[1] == "zeroderivative+natural": + if ( + bc[0] == "zeroderivative+natural" + or bc[1] == "zeroderivative+natural" + ): raise ValueError( - "zeroderivative+natural bc not supported by " "CubicHermiteSplines!" + "zeroderivative+natural bc not supported by " + "CubicHermiteSplines!" ) if derivatives_at_nodes is None: - derivatives_at_nodes = _finite_differences(nodes, values_at_nodes, bc) + derivatives_at_nodes = _finite_differences( + nodes, values_at_nodes, bc + ) self._derivatives_by_fd = True else: derivatives_at_nodes = np.asarray( @@ -1620,7 +1688,9 @@ def d_scaled(self, i: Integral) -> sp.Expr: return self.derivatives_at_nodes[i] / self.values_at_nodes[i] return self.derivatives_at_nodes[i] - def _poly_variable(self, x: Union[Real, sp.Basic], i: Integral) -> sp.Basic: + def _poly_variable( + self, x: Union[Real, sp.Basic], i: Integral + ) -> sp.Basic: assert 0 <= i < len(self.nodes) - 1 dx = self.nodes[i + 1] - self.nodes[i] with evaluate(False): @@ -1662,7 +1732,9 @@ def _parameters(self) -> Set[sp.Symbol]: parameters.update(d.free_symbols) return parameters - def _replace_in_all_expressions(self, old: sp.Symbol, new: sp.Symbol) -> None: + def _replace_in_all_expressions( + self, old: sp.Symbol, new: sp.Symbol + ) -> None: super()._replace_in_all_expressions(old, new) self._derivatives_at_nodes = [ d.subs(old, new) for d in self.derivatives_at_nodes @@ -1697,7 +1769,9 @@ def __str__(self) -> str: return s + " [" + ", ".join(cmps) + "]" -def _finite_differences(xx: np.ndarray, yy: np.ndarray, bc: NormalizedBC) -> np.ndarray: +def _finite_differences( + xx: np.ndarray, yy: np.ndarray, bc: NormalizedBC +) -> np.ndarray: dd = [] if bc[0] == "periodic": @@ -1736,7 +1810,9 @@ def _finite_differences(xx: np.ndarray, yy: np.ndarray, bc: NormalizedBC) -> np. "At least 3 nodes are needed " "for computing finite differences with natural bc!" ) - fd = _natural_fd(yy[-1], xx[-2] - xx[-1], yy[-2], xx[-3] - xx[-2], yy[-3]) + fd = _natural_fd( + yy[-1], xx[-2] - xx[-1], yy[-2], xx[-3] - xx[-2], yy[-3] + ) else: fd = _onesided_fd(yy[-2], yy[-1], xx[-1] - xx[-2]) dd.append(fd) diff --git a/python/sdist/amici/swig.py b/python/sdist/amici/swig.py index 49cbe711de..ef75646389 100644 --- a/python/sdist/amici/swig.py +++ b/python/sdist/amici/swig.py @@ -16,21 +16,29 @@ class TypeHintFixer(ast.NodeTransformer): "size_t": ast.Name("int"), "bool": ast.Name("bool"), "std::unique_ptr< amici::Solver >": ast.Constant("Solver"), - "amici::InternalSensitivityMethod": ast.Constant("InternalSensitivityMethod"), + "amici::InternalSensitivityMethod": ast.Constant( + "InternalSensitivityMethod" + ), "amici::InterpolationType": ast.Constant("InterpolationType"), "amici::LinearMultistepMethod": ast.Constant("LinearMultistepMethod"), "amici::LinearSolver": ast.Constant("LinearSolver"), "amici::Model *": ast.Constant("Model"), "amici::Model const *": ast.Constant("Model"), - "amici::NewtonDampingFactorMode": ast.Constant("NewtonDampingFactorMode"), - "amici::NonlinearSolverIteration": ast.Constant("NonlinearSolverIteration"), + "amici::NewtonDampingFactorMode": ast.Constant( + "NewtonDampingFactorMode" + ), + "amici::NonlinearSolverIteration": ast.Constant( + "NonlinearSolverIteration" + ), "amici::ObservableScaling": ast.Constant("ObservableScaling"), "amici::ParameterScaling": ast.Constant("ParameterScaling"), "amici::RDataReporting": ast.Constant("RDataReporting"), "amici::SensitivityMethod": ast.Constant("SensitivityMethod"), "amici::SensitivityOrder": ast.Constant("SensitivityOrder"), "amici::Solver *": ast.Constant("Solver"), - "amici::SteadyStateSensitivityMode": ast.Constant("SteadyStateSensitivityMode"), + "amici::SteadyStateSensitivityMode": ast.Constant( + "SteadyStateSensitivityMode" + ), "amici::realtype": ast.Name("float"), "DoubleVector": ast.Constant("Sequence[float]"), "IntVector": ast.Name("Sequence[int]"), diff --git a/python/sdist/amici/swig_wrappers.py b/python/sdist/amici/swig_wrappers.py index 6f859b0956..f56f3bd5d2 100644 --- a/python/sdist/amici/swig_wrappers.py +++ b/python/sdist/amici/swig_wrappers.py @@ -107,7 +107,8 @@ def runAmiciSimulation( """ if ( model.ne > 0 - and solver.getSensitivityMethod() == amici_swig.SensitivityMethod.adjoint + and solver.getSensitivityMethod() + == amici_swig.SensitivityMethod.adjoint and solver.getSensitivityOrder() == amici_swig.SensitivityOrder.first ): warnings.warn( @@ -168,7 +169,8 @@ def runAmiciSimulations( """ if ( model.ne > 0 - and solver.getSensitivityMethod() == amici_swig.SensitivityMethod.adjoint + and solver.getSensitivityMethod() + == amici_swig.SensitivityMethod.adjoint and solver.getSensitivityOrder() == amici_swig.SensitivityOrder.first ): warnings.warn( @@ -311,7 +313,9 @@ def _log_simulation(rdata: amici_swig.ReturnData): ) -def _ids_and_names_to_rdata(rdata: amici_swig.ReturnData, model: amici_swig.Model): +def _ids_and_names_to_rdata( + rdata: amici_swig.ReturnData, model: amici_swig.Model +): """Copy entity IDs and names from a Model to ReturnData.""" for entity_type in ( "State", diff --git a/python/sdist/setup.py b/python/sdist/setup.py index 2c30eaa7ef..4d65634cc2 100755 --- a/python/sdist/setup.py +++ b/python/sdist/setup.py @@ -130,7 +130,9 @@ def get_extensions(): source_dir="amici", cmake_configure_options=[ *global_cmake_configure_options, - "-Werror=dev" if "GITHUB_ACTIONS" in os.environ else "-Wno-error=dev", + "-Werror=dev" + if "GITHUB_ACTIONS" in os.environ + else "-Wno-error=dev", "-DAMICI_PYTHON_BUILD_EXT_ONLY=ON", f"-DPython3_EXECUTABLE={Path(sys.executable).as_posix()}", ], diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 1c07ddacac..9ab64b91d7 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -41,7 +41,9 @@ def sbml_example_presimulation_module(): constant_parameters=constant_parameters, ) - yield amici.import_model_module(module_name=module_name, module_path=outdir) + yield amici.import_model_module( + module_name=module_name, module_path=outdir + ) @pytest.fixture(scope="session") diff --git a/python/tests/pysb_test_models/bngwiki_egfr_simple_deletemolecules.py b/python/tests/pysb_test_models/bngwiki_egfr_simple_deletemolecules.py index 4723d3ac36..767c239c5d 100644 --- a/python/tests/pysb_test_models/bngwiki_egfr_simple_deletemolecules.py +++ b/python/tests/pysb_test_models/bngwiki_egfr_simple_deletemolecules.py @@ -76,7 +76,9 @@ ) # Transphosphorylation of EGFR by RTK -Rule("egfr_transphos", EGFR(CR1=ANY, Y1068="U") >> EGFR(CR1=ANY, Y1068="P"), kp3) +Rule( + "egfr_transphos", EGFR(CR1=ANY, Y1068="U") >> EGFR(CR1=ANY, Y1068="P"), kp3 +) # Dephosphorylation Rule("egfr_dephos", EGFR(Y1068="P") >> EGFR(Y1068="U"), km3) diff --git a/python/tests/splines_utils.py b/python/tests/splines_utils.py index 3c61e089c9..0746207ddb 100644 --- a/python/tests/splines_utils.py +++ b/python/tests/splines_utils.py @@ -267,7 +267,9 @@ def create_petab_problem( problem.to_files( model_file=os.path.join(folder, f"{model_name}_model.xml"), condition_file=os.path.join(folder, f"{model_name}_conditions.tsv"), - measurement_file=os.path.join(folder, f"{model_name}_measurements.tsv"), + measurement_file=os.path.join( + folder, f"{model_name}_measurements.tsv" + ), parameter_file=os.path.join(folder, f"{model_name}_parameters.tsv"), observable_file=os.path.join(folder, f"{model_name}_observables.tsv"), yaml_file=os.path.join(folder, f"{model_name}.yaml"), @@ -370,10 +372,14 @@ def simulate_splines( ) if petab_problem is None and amici_model is not None: - raise ValueError("if amici_model is given, petab_problem must be given too") + raise ValueError( + "if amici_model is given, petab_problem must be given too" + ) if petab_problem is not None and initial_values is None: - raise ValueError("if petab_problem is given, initial_values must be given too") + raise ValueError( + "if petab_problem is given, initial_values must be given too" + ) if petab_problem is None: # Create PEtab problem @@ -467,14 +473,18 @@ def simulate_splines( ) -def compute_ground_truth(splines, initial_values, times, params_true, params_sorted): +def compute_ground_truth( + splines, initial_values, times, params_true, params_sorted +): x_true_sym = sp.Matrix( [ integrate_spline(spline, None, times, iv) for (spline, iv) in zip(splines, initial_values) ] ).transpose() - groundtruth = {"x_true": np.asarray(x_true_sym.subs(params_true), dtype=float)} + groundtruth = { + "x_true": np.asarray(x_true_sym.subs(params_true), dtype=float) + } sx_by_state = [ x_true_sym[:, i].jacobian(params_sorted).subs(params_true) for i in range(x_true_sym.shape[1]) @@ -572,7 +582,9 @@ def check_splines( # Sort splines/ics/parameters as in the AMICI model splines = [splines[species_to_index(name)] for name in state_ids] - initial_values = [initial_values[species_to_index(name)] for name in state_ids] + initial_values = [ + initial_values[species_to_index(name)] for name in state_ids + ] def param_by_name(id): for p in params_true.keys(): @@ -672,7 +684,9 @@ def param_by_name(id): ) elif debug == "print": sx_err_abs = abs(rdata["sx"] - sx_true) - sx_err_rel = np.where(sx_err_abs == 0, 0, sx_err_abs / abs(sx_true)) + sx_err_rel = np.where( + sx_err_abs == 0, 0, sx_err_abs / abs(sx_true) + ) print(f"sx_atol={sx_atol} sx_rtol={sx_rtol}") print("sx_err_abs:") print(np.squeeze(sx_err_abs)) @@ -701,7 +715,9 @@ def param_by_name(id): if sllh_atol is None: sllh_atol = np.finfo(float).eps sllh_err_abs = abs(sllh).max() - if (sllh_err_abs > sllh_atol and debug is not True) or debug == "print": + if ( + sllh_err_abs > sllh_atol and debug is not True + ) or debug == "print": print(f"sllh_atol={sllh_atol}") print(f"sllh_err_abs = {sllh_err_abs}") if not debug: @@ -710,7 +726,11 @@ def param_by_name(id): assert sllh is None # Try different parameter lists - if not skip_sensitivity and (not use_adjoint) and parameter_lists is not None: + if ( + not skip_sensitivity + and (not use_adjoint) + and parameter_lists is not None + ): for plist in parameter_lists: amici_model.setParameterList(plist) amici_model.setTimepoints(rdata.t) diff --git a/python/tests/test_compare_conservation_laws_sbml.py b/python/tests/test_compare_conservation_laws_sbml.py index fa01533583..4d6a453b52 100644 --- a/python/tests/test_compare_conservation_laws_sbml.py +++ b/python/tests/test_compare_conservation_laws_sbml.py @@ -28,7 +28,9 @@ def edata_fixture(): 2, 0, 0, - np.array([0.0, 0.0, 0.0, 1.0, 2.0, 2.0, 4.0, float("inf"), float("inf")]), + np.array( + [0.0, 0.0, 0.0, 1.0, 2.0, 2.0, 4.0, float("inf"), float("inf")] + ), ) edata_full.setObservedData([3.14] * 18) edata_full.fixedParameters = np.array([1.0, 2.0]) @@ -129,7 +131,9 @@ def test_compare_conservation_laws_sbml(models, edata_fixture): assert model_without_cl.nx_rdata == model_with_cl.nx_rdata assert model_with_cl.nx_solver < model_without_cl.nx_solver assert len(model_with_cl.getStateIdsSolver()) == model_with_cl.nx_solver - assert len(model_without_cl.getStateIdsSolver()) == model_without_cl.nx_solver + assert ( + len(model_without_cl.getStateIdsSolver()) == model_without_cl.nx_solver + ) # ----- compare simulations wo edata, sensi = 0, states ------------------ # run simulations @@ -258,9 +262,15 @@ def test_adjoint_pre_and_post_equilibration(models, edata_fixture): ) # assert all are close - assert_allclose(rff_cl["sllh"], rfa_cl["sllh"], rtol=1.0e-5, atol=1.0e-8) - assert_allclose(rfa_cl["sllh"], raa_cl["sllh"], rtol=1.0e-5, atol=1.0e-8) - assert_allclose(raa_cl["sllh"], rff_cl["sllh"], rtol=1.0e-5, atol=1.0e-8) + assert_allclose( + rff_cl["sllh"], rfa_cl["sllh"], rtol=1.0e-5, atol=1.0e-8 + ) + assert_allclose( + rfa_cl["sllh"], raa_cl["sllh"], rtol=1.0e-5, atol=1.0e-8 + ) + assert_allclose( + raa_cl["sllh"], rff_cl["sllh"], rtol=1.0e-5, atol=1.0e-8 + ) # compare fully adjoint approach to simulation with singular # Jacobian diff --git a/python/tests/test_conserved_quantities_demartino.py b/python/tests/test_conserved_quantities_demartino.py index 58b6f2e9a6..339743cc4e 100644 --- a/python/tests/test_conserved_quantities_demartino.py +++ b/python/tests/test_conserved_quantities_demartino.py @@ -167,7 +167,8 @@ def data_demartino2014(): S = [ int(item) for sl in [ - entry.decode("ascii").strip().split("\t") for entry in data.readlines() + entry.decode("ascii").strip().split("\t") + for entry in data.readlines() ] for item in sl ] @@ -177,7 +178,9 @@ def data_demartino2014(): r"https://github.com/AMICI-dev/AMICI/files/11430970/test-ecoli-met.txt", timeout=10, ) - row_names = [entry.decode("ascii").strip() for entry in io.BytesIO(response.read())] + row_names = [ + entry.decode("ascii").strip() for entry in io.BytesIO(response.read()) + ] return S, row_names @@ -194,7 +197,9 @@ def test_kernel_demartino2014(data_demartino2014, quiet=True): ), "Unexpected dimension of stoichiometric matrix" # Expected number of metabolites per conservation law found after kernel() - expected_num_species = [53] + [2] * 11 + [6] + [3] * 2 + [2] * 15 + [3] + [2] * 5 + expected_num_species = ( + [53] + [2] * 11 + [6] + [3] * 2 + [2] * 15 + [3] + [2] * 5 + ) ( kernel_dim, @@ -222,7 +227,9 @@ def test_kernel_demartino2014(data_demartino2014, quiet=True): assert ( engaged_species == demartino2014_kernel_engaged_species ), "Wrong engaged metabolites reported" - assert len(conserved_moieties) == 128, "Wrong number of conserved moieties reported" + assert ( + len(conserved_moieties) == 128 + ), "Wrong number of conserved moieties reported" # Assert that each conserved moiety has the correct number of metabolites for i in range(int_kernel_dim - 2): @@ -770,7 +777,9 @@ def test_fill_demartino2014(data_demartino2014): assert not any(fields[len(ref_for_fields) :]) -def compute_moiety_conservation_laws_demartino2014(data_demartino2014, quiet=False): +def compute_moiety_conservation_laws_demartino2014( + data_demartino2014, quiet=False +): """Compute conserved quantities for De Martino's published results for E. coli network""" stoichiometric_list, row_names = data_demartino2014 @@ -799,7 +808,9 @@ def compute_moiety_conservation_laws_demartino2014(data_demartino2014, quiet=Fal def test_compute_moiety_conservation_laws_demartino2014(data_demartino2014): """Invoke test case and benchmarking for De Martino's published results for E. coli network""" - compute_moiety_conservation_laws_demartino2014(data_demartino2014, quiet=False) + compute_moiety_conservation_laws_demartino2014( + data_demartino2014, quiet=False + ) @skip_on_valgrind @@ -830,7 +841,9 @@ def test_compute_moiety_conservation_laws_simple(): stoichiometric_matrix = sp.Matrix( [[-1.0, 1.0], [-1.0, 1.0], [1.0, -1.0], [1.0, -1.0]] ) - stoichiometric_list = [float(entry) for entry in stoichiometric_matrix.T.flat()] + stoichiometric_list = [ + float(entry) for entry in stoichiometric_matrix.T.flat() + ] num_tries = 1000 found_all_n_times = 0 diff --git a/python/tests/test_edata.py b/python/tests/test_edata.py index c2d2ea470b..9c4d9b9edc 100644 --- a/python/tests/test_edata.py +++ b/python/tests/test_edata.py @@ -16,7 +16,9 @@ def test_edata_sensi_unscaling(model_units_module): sx0 = (3, 3, 3, 3) - parameter_scales_log10 = [amici.ParameterScaling.log10.value] * len(parameters0) + parameter_scales_log10 = [amici.ParameterScaling.log10.value] * len( + parameters0 + ) amici_parameter_scales_log10 = amici.parameterScalingFromIntVector( parameter_scales_log10 ) diff --git a/python/tests/test_events.py b/python/tests/test_events.py index 30a65c4043..d2a177bded 100644 --- a/python/tests/test_events.py +++ b/python/tests/test_events.py @@ -15,8 +15,12 @@ @pytest.fixture( params=[ pytest.param("events_plus_heavisides", marks=skip_on_valgrind), - pytest.param("piecewise_plus_event_simple_case", marks=skip_on_valgrind), - pytest.param("piecewise_plus_event_semi_complicated", marks=skip_on_valgrind), + pytest.param( + "piecewise_plus_event_simple_case", marks=skip_on_valgrind + ), + pytest.param( + "piecewise_plus_event_semi_complicated", marks=skip_on_valgrind + ), pytest.param( "piecewise_plus_event_trigger_depends_on_state", marks=skip_on_valgrind, @@ -73,7 +77,9 @@ def get_model_definition(model_name): if model_name == "event_state_dep_ddeltax_dtpx": return model_definition_event_state_dep_ddeltax_dtpx() - raise NotImplementedError(f"Model with name {model_name} is not implemented.") + raise NotImplementedError( + f"Model with name {model_name} is not implemented." + ) def model_definition_events_plus_heavisides(): @@ -290,7 +296,9 @@ def x_expected(t, k1, k2, inflow_1, decay_1, decay_2, bolus): def get_early_x(t): # compute dynamics before event - x_1 = equil * (1 - np.exp(-decay_1 * t)) + k1 * np.exp(-decay_1 * t) + x_1 = equil * (1 - np.exp(-decay_1 * t)) + k1 * np.exp( + -decay_1 * t + ) x_2 = k2 * np.exp(-decay_2 * t) return np.array([[x_1], [x_2]]) @@ -304,9 +312,9 @@ def get_early_x(t): # compute dynamics after event inhom = np.exp(decay_1 * event_time) * tau_x1 - x_1 = equil * (1 - np.exp(decay_1 * (event_time - t))) + inhom * np.exp( - -decay_1 * t - ) + x_1 = equil * ( + 1 - np.exp(decay_1 * (event_time - t)) + ) + inhom * np.exp(-decay_1 * t) x_2 = tau_x2 * np.exp(decay_2 * event_time) * np.exp(-decay_2 * t) x = np.array([[x_1], [x_2]]) @@ -463,7 +471,8 @@ def x_expected(t, x_1_0, alpha, beta, gamma, delta): else: # after third event triggered x = ( - ((x_1_0 + alpha) * alpha + (beta - alpha)) * delta + (gamma - beta) + ((x_1_0 + alpha) * alpha + (beta - alpha)) * delta + + (gamma - beta) ) ** 2 * 2 + (t - gamma) return np.array((x,)) @@ -546,7 +555,9 @@ def x_expected(t, x_1_0, x_2_0, alpha, beta, gamma, delta, eta): x_1 = x_1_heaviside_1 * np.exp(delta * (t - heaviside_1)) else: x_1_heaviside_1 = gamma * np.exp(-(heaviside_1 - t_event_1)) - x_1_at_event_2 = x_1_heaviside_1 * np.exp(delta * (t_event_2 - heaviside_1)) + x_1_at_event_2 = x_1_heaviside_1 * np.exp( + delta * (t_event_2 - heaviside_1) + ) x_2_at_event_2 = x_2_0 * np.exp(-eta * t_event_2) x1_after_event_2 = x_1_at_event_2 + x_2_at_event_2 x_1 = x1_after_event_2 * np.exp(-(t - t_event_2)) @@ -671,8 +682,12 @@ def sx_expected(t, parameters): def test_models(model): amici_model, parameters, timepoints, x_expected, sx_expected = model - result_expected_x = np.array([x_expected(t, **parameters) for t in timepoints]) - result_expected_sx = np.array([sx_expected(t, parameters) for t in timepoints]) + result_expected_x = np.array( + [x_expected(t, **parameters) for t in timepoints] + ) + result_expected_sx = np.array( + [sx_expected(t, parameters) for t in timepoints] + ) # assert correctness of trajectories check_trajectories_without_sensitivities(amici_model, result_expected_x) diff --git a/python/tests/test_hdf5.py b/python/tests/test_hdf5.py index c47d8653eb..232f22be8c 100644 --- a/python/tests/test_hdf5.py +++ b/python/tests/test_hdf5.py @@ -32,7 +32,9 @@ def _modify_solver_attrs(solver): getattr(solver, attr)(cval) -@pytest.mark.skipif(not amici.hdf5_enabled, reason="AMICI was compiled without HDF5") +@pytest.mark.skipif( + not amici.hdf5_enabled, reason="AMICI was compiled without HDF5" +) def test_solver_hdf5_roundtrip(sbml_example_presimulation_module): """TestCase class for AMICI HDF5 I/O""" diff --git a/python/tests/test_heavisides.py b/python/tests/test_heavisides.py index 26e9e7e96a..c3bea26a0c 100644 --- a/python/tests/test_heavisides.py +++ b/python/tests/test_heavisides.py @@ -53,8 +53,12 @@ def model(request): def test_models(model): amici_model, parameters, timepoints, x_expected, sx_expected = model - result_expected_x = np.array([x_expected(t, **parameters) for t in timepoints]) - result_expected_sx = np.array([sx_expected(t, **parameters) for t in timepoints]) + result_expected_x = np.array( + [x_expected(t, **parameters) for t in timepoints] + ) + result_expected_sx = np.array( + [sx_expected(t, **parameters) for t in timepoints] + ) # Does the AMICI simulation match the analytical solution? check_trajectories_without_sensitivities(amici_model, result_expected_x) @@ -71,7 +75,9 @@ def get_model_definition(model_name): elif model_name == "piecewise_many_conditions": return model_definition_piecewise_many_conditions() else: - raise NotImplementedError(f"Model with name {model_name} is not implemented.") + raise NotImplementedError( + f"Model with name {model_name} is not implemented." + ) def model_definition_state_and_parameter_dependent_heavisides(): @@ -136,8 +142,12 @@ def sx_expected(t, alpha, beta, gamma, delta, eta, zeta): sx_1_zeta = np.exp(alpha * t) else: # Never trust Wolfram Alpha... - sx_1_alpha = zeta * tau_1 * np.exp(alpha * tau_1 - beta * (t - tau_1)) - sx_1_beta = zeta * (tau_1 - t) * np.exp(alpha * tau_1 - beta * (t - tau_1)) + sx_1_alpha = ( + zeta * tau_1 * np.exp(alpha * tau_1 - beta * (t - tau_1)) + ) + sx_1_beta = ( + zeta * (tau_1 - t) * np.exp(alpha * tau_1 - beta * (t - tau_1)) + ) sx_1_gamma = ( zeta * (alpha + beta) diff --git a/python/tests/test_misc.py b/python/tests/test_misc.py index 299d5d15b1..5a88fda6f8 100644 --- a/python/tests/test_misc.py +++ b/python/tests/test_misc.py @@ -119,7 +119,9 @@ def test_monkeypatch(): assert (t**n).diff(t).subs(vals) is sp.nan # check that we can monkeypatch it out - with _monkeypatched(sp.Pow, "_eval_derivative", _custom_pow_eval_derivative): + with _monkeypatched( + sp.Pow, "_eval_derivative", _custom_pow_eval_derivative + ): assert (t**n).diff(t).subs(vals) is not sp.nan # check that the monkeypatch is transient diff --git a/python/tests/test_observable_events.py b/python/tests/test_observable_events.py index 28e298a312..2887308ff6 100644 --- a/python/tests/test_observable_events.py +++ b/python/tests/test_observable_events.py @@ -62,7 +62,9 @@ def model_neuron_def(): } } - event_observables = {"z1": {"name": "z1", "event": "event_1", "formula": "time"}} + event_observables = { + "z1": {"name": "z1", "event": "event_1", "formula": "time"} + } return ( initial_assignments, parameters, diff --git a/python/tests/test_pandas.py b/python/tests/test_pandas.py index e904fce7cc..21c58bcaff 100644 --- a/python/tests/test_pandas.py +++ b/python/tests/test_pandas.py @@ -49,8 +49,12 @@ def test_pandas_import_export(sbml_example_presimulation_module, case): assert case[fp] == getattr(edata_reconstructed[0], fp) else: - assert model.getFixedParameters() == getattr(edata_reconstructed[0], fp) + assert model.getFixedParameters() == getattr( + edata_reconstructed[0], fp + ) - assert model.getFixedParameters() == getattr(edata_reconstructed[0], fp) + assert model.getFixedParameters() == getattr( + edata_reconstructed[0], fp + ) assert getattr(edata[0], fp) == case[fp] diff --git a/python/tests/test_parameter_mapping.py b/python/tests/test_parameter_mapping.py index f74cb7e666..ae66e23f53 100644 --- a/python/tests/test_parameter_mapping.py +++ b/python/tests/test_parameter_mapping.py @@ -37,9 +37,16 @@ def test_parameter_mapping_for_condition_default_args(): expected_scale_map_preeq_fix = {"sim_par2": "lin"} expected_scale_map_sim_fix = {"sim_par2": "lin"} - assert par_map_for_condition.scale_map_sim_var == expected_scale_map_sim_var - assert par_map_for_condition.scale_map_preeq_fix == expected_scale_map_preeq_fix - assert par_map_for_condition.scale_map_sim_fix == expected_scale_map_sim_fix + assert ( + par_map_for_condition.scale_map_sim_var == expected_scale_map_sim_var + ) + assert ( + par_map_for_condition.scale_map_preeq_fix + == expected_scale_map_preeq_fix + ) + assert ( + par_map_for_condition.scale_map_sim_fix == expected_scale_map_sim_fix + ) @skip_on_valgrind diff --git a/python/tests/test_petab_import.py b/python/tests/test_petab_import.py index f6db30f18a..fc978b76ea 100644 --- a/python/tests/test_petab_import.py +++ b/python/tests/test_petab_import.py @@ -60,7 +60,9 @@ def test_get_fixed_parameters(simple_sbml_model): ) ) parameter_df = petab.get_parameter_df( - pd.DataFrame({petab.PARAMETER_ID: ["p3", "p4"], petab.ESTIMATE: [0, 1]}) + pd.DataFrame( + {petab.PARAMETER_ID: ["p3", "p4"], petab.ESTIMATE: [0, 1]} + ) ) print(condition_df) print(parameter_df) @@ -122,7 +124,9 @@ def test_default_output_parameters(simple_sbml_model): ) assert ( 1.0 - == sbml_importer.sbml.getParameter("observableParameter1_obs1").getValue() + == sbml_importer.sbml.getParameter( + "observableParameter1_obs1" + ).getValue() ) with pytest.raises(ValueError): diff --git a/python/tests/test_petab_objective.py b/python/tests/test_petab_objective.py index 1b2436ceab..e31e693d11 100755 --- a/python/tests/test_petab_objective.py +++ b/python/tests/test_petab_objective.py @@ -50,7 +50,9 @@ def test_simulate_petab_sensitivities(lotka_volterra): for scaled_gradients in [True, False]: _problem_parameters = problem_parameters.copy() if scaled_parameters: - _problem_parameters = petab_problem.scale_parameters(problem_parameters) + _problem_parameters = petab_problem.scale_parameters( + problem_parameters + ) results[(scaled_parameters, scaled_gradients)] = pd.Series( amici.petab_objective.simulate_petab( petab_problem=petab_problem, diff --git a/python/tests/test_petab_simulate.py b/python/tests/test_petab_simulate.py index 385f98e05e..febea5fd50 100644 --- a/python/tests/test_petab_simulate.py +++ b/python/tests/test_petab_simulate.py @@ -53,7 +53,9 @@ def test_subset_call(petab_problem): simulator0 = PetabSimulator(petab_problem) assert not (Path(model_output_dir) / model_name).is_dir() - simulator0.simulate(model_name=model_name, model_output_dir=model_output_dir) + simulator0.simulate( + model_name=model_name, model_output_dir=model_output_dir + ) # Model name is handled correctly assert simulator0.amici_model.getName() == model_name # Check model output directory is created, by diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index 3407c7aea0..a42bc6354d 100644 --- a/python/tests/test_preequilibration.py +++ b/python/tests/test_preequilibration.py @@ -142,7 +142,9 @@ def test_manual_preequilibration(preeq_fixture): rdata_sim[variable], atol=1e-6, rtol=1e-6, - err_msg=str(dict(pscale=pscale, plist=plist, variable=variable)), + err_msg=str( + dict(pscale=pscale, plist=plist, variable=variable) + ), ) @@ -386,7 +388,9 @@ def test_equilibration_methods_with_adjoints(preeq_fixture): atol=1e-6, rtol=1e-6, err_msg=str( - dict(variable=variable, setting1=setting1, setting2=setting2) + dict( + variable=variable, setting1=setting1, setting2=setting2 + ) ), ) @@ -535,11 +539,15 @@ def test_steadystate_computation_mode(preeq_fixture): assert rdatas[mode]["status"] == amici.AMICI_SUCCESS assert np.all( - rdatas[amici.SteadyStateComputationMode.integrationOnly]["preeq_status"][0] + rdatas[amici.SteadyStateComputationMode.integrationOnly][ + "preeq_status" + ][0] == [0, 1, 0] ) assert ( - rdatas[amici.SteadyStateComputationMode.integrationOnly]["preeq_numsteps"][0][0] + rdatas[amici.SteadyStateComputationMode.integrationOnly][ + "preeq_numsteps" + ][0][0] == 0 ) @@ -548,7 +556,10 @@ def test_steadystate_computation_mode(preeq_fixture): == [1, 0, 0] ) assert ( - rdatas[amici.SteadyStateComputationMode.newtonOnly]["preeq_numsteps"][0][0] > 0 + rdatas[amici.SteadyStateComputationMode.newtonOnly]["preeq_numsteps"][ + 0 + ][0] + > 0 ) # assert correct results @@ -575,7 +586,9 @@ def test_simulation_errors(preeq_fixture): ) = preeq_fixture solver.setSensitivityOrder(amici.SensitivityOrder.first) - solver.setSensitivityMethodPreequilibration(amici.SensitivityMethod.forward) + solver.setSensitivityMethodPreequilibration( + amici.SensitivityMethod.forward + ) model.setSteadyStateSensitivityMode( amici.SteadyStateSensitivityMode.integrationOnly ) @@ -606,10 +619,16 @@ def test_simulation_errors(preeq_fixture): rdata = amici.runAmiciSimulation(model, solver, e) assert rdata["status"] != amici.AMICI_SUCCESS assert rdata._swigptr.messages[0].severity == amici.LogSeverity_debug - assert rdata._swigptr.messages[0].identifier == "CVODES:CVode:RHSFUNC_FAIL" + assert ( + rdata._swigptr.messages[0].identifier + == "CVODES:CVode:RHSFUNC_FAIL" + ) assert rdata._swigptr.messages[1].severity == amici.LogSeverity_debug assert rdata._swigptr.messages[1].identifier == "EQUILIBRATION_FAILURE" - assert "exceedingly long simulation time" in rdata._swigptr.messages[1].message + assert ( + "exceedingly long simulation time" + in rdata._swigptr.messages[1].message + ) assert rdata._swigptr.messages[2].severity == amici.LogSeverity_error assert rdata._swigptr.messages[2].identifier == "OTHER" assert rdata._swigptr.messages[3].severity == amici.LogSeverity_debug diff --git a/python/tests/test_pregenerated_models.py b/python/tests/test_pregenerated_models.py index 23f6ff8969..5a110cdfc2 100755 --- a/python/tests/test_pregenerated_models.py +++ b/python/tests/test_pregenerated_models.py @@ -97,7 +97,10 @@ def test_pregenerated_model(sub_test, case): verify_simulation_opts["atol"] = 1e-5 verify_simulation_opts["rtol"] = 1e-2 - if model_name.startswith("model_robertson") and case == "sensiforwardSPBCG": + if ( + model_name.startswith("model_robertson") + and case == "sensiforwardSPBCG" + ): verify_simulation_opts["atol"] = 1e-3 verify_simulation_opts["rtol"] = 1e-3 @@ -115,7 +118,9 @@ def test_pregenerated_model(sub_test, case): if ( edata and model_name != "model_neuron_o2" - and not (model_name == "model_robertson" and case == "sensiforwardSPBCG") + and not ( + model_name == "model_robertson" and case == "sensiforwardSPBCG" + ) ): if isinstance(edata, amici.amici.ExpData): edatas = [edata, edata] @@ -224,14 +229,18 @@ def verify_simulation_results( subfields = expected_results["diagnosis"].keys() else: - attrs = [field for field in fields if field in expected_results.attrs.keys()] + attrs = [ + field for field in fields if field in expected_results.attrs.keys() + ] if "diagnosis" in expected_results.keys(): subfields = [ field for field in fields if field in expected_results["diagnosis"].keys() ] - fields = [field for field in fields if field in expected_results.keys()] + fields = [ + field for field in fields if field in expected_results.keys() + ] if expected_results.attrs["status"][0] != 0: assert rdata["status"] == expected_results.attrs["status"][0] @@ -272,4 +281,6 @@ def verify_simulation_results( ) for attr in attrs: - _check_results(rdata, attr, expected_results.attrs[attr], atol=atol, rtol=rtol) + _check_results( + rdata, attr, expected_results.attrs[attr], atol=atol, rtol=rtol + ) diff --git a/python/tests/test_pysb.py b/python/tests/test_pysb.py index 81273d783a..52ca3a320f 100644 --- a/python/tests/test_pysb.py +++ b/python/tests/test_pysb.py @@ -141,7 +141,9 @@ def test_compare_to_pysb_simulation(example): with amici.add_path(os.path.dirname(pysb.examples.__file__)): with amici.add_path( - os.path.join(os.path.dirname(__file__), "..", "tests", "pysb_test_models") + os.path.join( + os.path.dirname(__file__), "..", "tests", "pysb_test_models" + ) ): # load example pysb.SelfExporter.cleanup() # reset pysb @@ -185,7 +187,9 @@ def test_compare_to_pysb_simulation(example): observables=list(pysb_model.observables.keys()), ) - amici_model_module = amici.import_model_module(pysb_model.name, outdir) + amici_model_module = amici.import_model_module( + pysb_model.name, outdir + ) model_pysb = amici_model_module.getModel() model_pysb.setTimepoints(tspan) @@ -196,7 +200,9 @@ def test_compare_to_pysb_simulation(example): rdata = amici.runAmiciSimulation(model_pysb, solver) # check agreement of species simulations - assert np.isclose(rdata["x"], pysb_simres.species, 1e-4, 1e-4).all() + assert np.isclose( + rdata["x"], pysb_simres.species, 1e-4, 1e-4 + ).all() if example not in [ "fricker_2010_apoptosis", @@ -375,4 +381,6 @@ def test_energy(): solver.setRelativeTolerance(1e-14) solver.setAbsoluteTolerance(1e-14) - check_derivatives(amici_model, solver, epsilon=1e-4, rtol=1e-2, atol=1e-2) + check_derivatives( + amici_model, solver, epsilon=1e-4, rtol=1e-2, atol=1e-2 + ) diff --git a/python/tests/test_rdata.py b/python/tests/test_rdata.py index cbfc6dc7a9..29ea401932 100644 --- a/python/tests/test_rdata.py +++ b/python/tests/test_rdata.py @@ -23,7 +23,9 @@ def test_rdata_by_id(rdata_by_id_fixture): assert_array_equal(rdata.by_id(model.getStateIds()[1]), rdata.x[:, 1]) assert_array_equal(rdata.by_id(model.getStateIds()[1], "x"), rdata.x[:, 1]) - assert_array_equal(rdata.by_id(model.getStateIds()[1], "x", model), rdata.x[:, 1]) + assert_array_equal( + rdata.by_id(model.getStateIds()[1], "x", model), rdata.x[:, 1] + ) assert_array_equal( rdata.by_id(model.getObservableIds()[0], "y", model), rdata.y[:, 0] diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index c0e377aa29..7c4a67c0a2 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -138,7 +138,9 @@ def observable_dependent_error_model(simple_sbml_model): "observable_s1_scaled": "0.02 * observable_s1_scaled", }, ) - yield amici.import_model_module(module_name=model_name, module_path=tmpdir) + yield amici.import_model_module( + module_name=model_name, module_path=tmpdir + ) @skip_on_valgrind @@ -159,7 +161,9 @@ def test_sbml2amici_observable_dependent_error( rtol=1.0e-5, atol=1.0e-8, ) - assert_allclose(rdata.sigmay[:, 1], 0.02 * rdata.y[:, 1], rtol=1.0e-5, atol=1.0e-8) + assert_allclose( + rdata.sigmay[:, 1], 0.02 * rdata.y[:, 1], rtol=1.0e-5, atol=1.0e-8 + ) edata = amici.ExpData(rdata, 1.0, 0.0) edata.setObservedDataStdDev(np.nan) @@ -204,7 +208,9 @@ def model_steadystate_module(): observables = amici.assignmentRules2observables( sbml_importer.sbml, - filter_function=lambda variable: variable.getId().startswith("observable_") + filter_function=lambda variable: variable.getId().startswith( + "observable_" + ) and not variable.getId().endswith("_sigma"), ) @@ -218,7 +224,9 @@ def model_steadystate_module(): sigmas={"observable_x1withsigma": "observable_x1withsigma_sigma"}, ) - yield amici.import_model_module(module_name=module_name, module_path=outdir) + yield amici.import_model_module( + module_name=module_name, module_path=outdir + ) @pytest.fixture(scope="session") @@ -231,7 +239,9 @@ def model_units_module(): with TemporaryDirectory() as outdir: sbml_importer.sbml2amici(model_name=module_name, output_dir=outdir) - yield amici.import_model_module(module_name=module_name, module_path=outdir) + yield amici.import_model_module( + module_name=module_name, module_path=outdir + ) def test_presimulation(sbml_example_presimulation_module): @@ -331,7 +341,9 @@ def test_steadystate_simulation(model_steadystate_module): solver.setRelativeTolerance(1e-12) solver.setAbsoluteTolerance(1e-12) - check_derivatives(model, solver, edata[0], atol=1e-3, rtol=1e-3, epsilon=1e-4) + check_derivatives( + model, solver, edata[0], atol=1e-3, rtol=1e-3, epsilon=1e-4 + ) # Run some additional tests which need a working Model, # but don't need precomputed expectations. @@ -414,7 +426,9 @@ def model_test_likelihoods(): noise_distributions=noise_distributions, ) - yield amici.import_model_module(module_name=module_name, module_path=outdir) + yield amici.import_model_module( + module_name=module_name, module_path=outdir + ) @skip_on_valgrind @@ -520,7 +534,9 @@ def test_units(model_units_module): @skip_on_valgrind -@pytest.mark.skipif(os.name == "nt", reason="Avoid `CERTIFICATE_VERIFY_FAILED` error") +@pytest.mark.skipif( + os.name == "nt", reason="Avoid `CERTIFICATE_VERIFY_FAILED` error" +) def test_sympy_exp_monkeypatch(): """ This model contains a removeable discontinuity at t=0 that requires @@ -565,7 +581,9 @@ def test_sympy_exp_monkeypatch(): # print sensitivity-related results assert rdata["status"] == amici.AMICI_SUCCESS - check_derivatives(model, solver, None, atol=1e-2, rtol=1e-2, epsilon=1e-3) + check_derivatives( + model, solver, None, atol=1e-2, rtol=1e-2, epsilon=1e-3 + ) def normal_nllh(m, y, sigma): @@ -602,7 +620,8 @@ def log_laplace_nllh(m, y, sigma): def log10_laplace_nllh(m, y, sigma): return sum( - np.log(2 * sigma * m * np.log(10)) + np.abs(np.log10(y) - np.log10(m)) / sigma + np.log(2 * sigma * m * np.log(10)) + + np.abs(np.log10(y) - np.log10(m)) / sigma ) diff --git a/python/tests/test_sbml_import_special_functions.py b/python/tests/test_sbml_import_special_functions.py index 2042bf2fba..9d8f447511 100644 --- a/python/tests/test_sbml_import_special_functions.py +++ b/python/tests/test_sbml_import_special_functions.py @@ -51,7 +51,9 @@ def model_special_likelihoods(): noise_distributions=noise_distributions, ) - yield amici.import_model_module(module_name=module_name, module_path=outdir) + yield amici.import_model_module( + module_name=module_name, module_path=outdir + ) @skip_on_valgrind @@ -210,9 +212,13 @@ def test_rateof(): assert_approx_equal(rdata["xdot"][i_S1], rdata["xdot"][i_p2]) assert_array_almost_equal_nulp(rdata.by_id("S3"), t, 10) - assert_array_almost_equal_nulp(rdata.by_id("S2"), 2 * rdata.by_id("S3")) + assert_array_almost_equal_nulp( + rdata.by_id("S2"), 2 * rdata.by_id("S3") + ) assert_array_almost_equal_nulp( rdata.by_id("S4")[1:], 0.5 * np.diff(rdata.by_id("S3")), 10 ) assert_array_almost_equal_nulp(rdata.by_id("p3"), 0) - assert_array_almost_equal_nulp(rdata.by_id("p2"), 1 + rdata.by_id("S1")) + assert_array_almost_equal_nulp( + rdata.by_id("p2"), 1 + rdata.by_id("S1") + ) diff --git a/python/tests/test_splines.py b/python/tests/test_splines.py index 93e8063a5d..a7fe01e84a 100644 --- a/python/tests/test_splines.py +++ b/python/tests/test_splines.py @@ -61,7 +61,9 @@ def test_multiple_splines(**kwargs): tols5 = (tols5, tols5, tols5) tols = [] - for t0, t1, t2, t3, t4, t5 in zip(tols0, tols1, tols2, tols3, tols4, tols5): + for t0, t1, t2, t3, t4, t5 in zip( + tols0, tols1, tols2, tols3, tols4, tols5 + ): keys = set().union( t0.keys(), t1.keys(), t2.keys(), t3.keys(), t4.keys(), t5.keys() ) diff --git a/python/tests/test_splines_python.py b/python/tests/test_splines_python.py index 5d7ac15e0c..4c4de5ccfc 100644 --- a/python/tests/test_splines_python.py +++ b/python/tests/test_splines_python.py @@ -217,7 +217,9 @@ def check_gradient(spline, t, params, params_values, expected, rel_tol=1e-9): value = spline.evaluate(t) subs = {pname: pvalue for (pname, pvalue) in zip(params, params_values)} for p, exp in zip(params, expected): - assert math.isclose(float(value.diff(p).subs(subs)), exp, rel_tol=rel_tol) + assert math.isclose( + float(value.diff(p).subs(subs)), exp, rel_tol=rel_tol + ) @skip_on_valgrind @@ -240,7 +242,9 @@ def test_SplineUniformSensitivity(): rel_tol=1e-5, ) check_gradient(spline, 1.0 / 3, params, params_values, [0.0, 0.0, 5.0]) - check_gradient(spline, 0.50, params, params_values, [0.1875, -0.125, 2.625]) + check_gradient( + spline, 0.50, params, params_values, [0.1875, -0.125, 2.625] + ) check_gradient(spline, 2.0 / 3, params, params_values, [0.0, 0.0, 0.0]) check_gradient( spline, @@ -275,7 +279,9 @@ def test_SplineNonUniformSensitivity(): check_gradient(spline, 0.10, params, params_values, [0.0, 0.0, 5.0]) check_gradient(spline, 0.30, params, params_values, [-0.45, -0.3, 3.6]) check_gradient(spline, 0.50, params, params_values, [0.0, 0.0, 0.0]) - check_gradient(spline, 0.75, params, params_values, [-2.625, 0.4375, 0.921875]) + check_gradient( + spline, 0.75, params, params_values, [-2.625, 0.4375, 0.921875] + ) check_gradient(spline, 1.00, params, params_values, [-6.0, 1.0, 3.0]) diff --git a/python/tests/test_splines_short.py b/python/tests/test_splines_short.py index 37df5f5db9..59e54a3279 100644 --- a/python/tests/test_splines_short.py +++ b/python/tests/test_splines_short.py @@ -99,7 +99,9 @@ def test_splines_plist(): # Real spline #3 xx = UniformGrid(0, 5, number_of_nodes=6) p1, p2, p3, p4, p5 = sp.symbols("p1 p2 p3 p4 p5") - yy = np.asarray([p1 + p2, p2 * p3, p4, sp.cos(p1 + p3), p4 * sp.log(p1), p3]) + yy = np.asarray( + [p1 + p2, p2 * p3, p4, sp.cos(p1 + p3), p4 * sp.log(p1), p3] + ) dd = np.asarray([-0.75, -0.875, p5, 0.125, 1.15057181, 0.0]) params = {p1: 1.0, p2: 0.5, p3: 1.5, p4: -0.25, p5: -0.5} # print([y.subs(params).evalf() for y in yy]) diff --git a/python/tests/test_swig_interface.py b/python/tests/test_swig_interface.py index bfb7c74129..b1afa0b76b 100644 --- a/python/tests/test_swig_interface.py +++ b/python/tests/test_swig_interface.py @@ -49,7 +49,9 @@ def test_copy_constructors(pysb_example_presimulation_module): obj_clone = obj.clone() - assert get_val(obj, attr) == get_val(obj_clone, attr), f"{obj} - {attr}" + assert get_val(obj, attr) == get_val( + obj_clone, attr + ), f"{obj} - {attr}" # `None` values are skipped in `test_model_instance_settings`. @@ -165,7 +167,9 @@ def test_model_instance_settings(pysb_example_presimulation_module): # The new model has the default settings. model_default_settings = amici.get_model_settings(model) for name in model_instance_settings: - if (name == "InitialStates" and not model.hasCustomInitialStates()) or ( + if ( + name == "InitialStates" and not model.hasCustomInitialStates() + ) or ( name == ( "getInitialStateSensitivities", @@ -177,7 +181,8 @@ def test_model_instance_settings(pysb_example_presimulation_module): assert model_default_settings[name] == [] else: assert ( - model_default_settings[name] == model_instance_settings[name][i_default] + model_default_settings[name] + == model_instance_settings[name][i_default] ), name # The grouped setter method works. @@ -224,7 +229,9 @@ def test_interdependent_settings(pysb_example_presimulation_module): # Some values need to be transformed to be tested in Python # (e.g. SWIG objects). Default transformer is no transformation # (the identity function). - getter_transformers = {setting: (lambda x: x) for setting in original_settings} + getter_transformers = { + setting: (lambda x: x) for setting in original_settings + } getter_transformers.update( { # Convert from SWIG object. @@ -314,7 +321,9 @@ def test_unhandled_settings(pysb_example_presimulation_module): name for names in model_instance_settings for name in ( - names if isinstance(names, tuple) else (f"get{names}", f"set{names}") + names + if isinstance(names, tuple) + else (f"get{names}", f"set{names}") ) ] @@ -378,7 +387,9 @@ def test_model_instance_settings_custom_x0(pysb_example_presimulation_module): assert not model.hasCustomInitialStateSensitivities() settings = amici.get_model_settings(model) model.setInitialStates(model.getInitialStates()) - model.setUnscaledInitialStateSensitivities(model.getInitialStateSensitivities()) + model.setUnscaledInitialStateSensitivities( + model.getInitialStateSensitivities() + ) amici.set_model_settings(model, settings) assert not model.hasCustomInitialStates() assert not model.hasCustomInitialStateSensitivities() diff --git a/python/tests/util.py b/python/tests/util.py index 14f514c997..dde10eb454 100644 --- a/python/tests/util.py +++ b/python/tests/util.py @@ -31,7 +31,9 @@ def create_amici_model(sbml_model, model_name, **kwargs) -> AmiciModel: else tempfile.mkdtemp() ) - sbml_importer.sbml2amici(model_name=model_name, output_dir=output_dir, **kwargs) + sbml_importer.sbml2amici( + model_name=model_name, output_dir=output_dir, **kwargs + ) model_module = import_model_module(model_name, output_dir) return model_module.getModel() @@ -111,7 +113,9 @@ def create_event_assignment(target, assignment): create_event_assignment(event_target, event_assignment) else: - create_event_assignment(event_def["target"], event_def["assignment"]) + create_event_assignment( + event_def["target"], event_def["assignment"] + ) if to_file: libsbml.writeSBMLToFile(document, to_file) @@ -133,7 +137,9 @@ def check_trajectories_without_sensitivities( solver.setAbsoluteTolerance(1e-15) solver.setRelativeTolerance(1e-12) rdata = runAmiciSimulation(amici_model, solver=solver) - _check_close(rdata["x"], result_expected_x, field="x", rtol=5e-9, atol=1e-13) + _check_close( + rdata["x"], result_expected_x, field="x", rtol=5e-9, atol=1e-13 + ) def check_trajectories_with_forward_sensitivities( @@ -153,5 +159,9 @@ def check_trajectories_with_forward_sensitivities( solver.setAbsoluteToleranceFSA(1e-15) solver.setRelativeToleranceFSA(1e-13) rdata = runAmiciSimulation(amici_model, solver=solver) - _check_close(rdata["x"], result_expected_x, field="x", rtol=1e-10, atol=1e-12) - _check_close(rdata["sx"], result_expected_sx, field="sx", rtol=1e-7, atol=1e-9) + _check_close( + rdata["x"], result_expected_x, field="x", rtol=1e-10, atol=1e-12 + ) + _check_close( + rdata["sx"], result_expected_sx, field="sx", rtol=1e-7, atol=1e-9 + ) diff --git a/tests/benchmark-models/evaluate_benchmark.py b/tests/benchmark-models/evaluate_benchmark.py index 749b70e3e5..0c6e2e4122 100644 --- a/tests/benchmark-models/evaluate_benchmark.py +++ b/tests/benchmark-models/evaluate_benchmark.py @@ -29,12 +29,15 @@ ratios = ( pd.concat( - [df[sensi] / df["t_sim"].values for sensi in ["t_fwd", "t_adj"]] + [df.np], + [df[sensi] / df["t_sim"].values for sensi in ["t_fwd", "t_adj"]] + + [df.np], axis=1, ) .reset_index() .melt(id_vars=["index", "np"]) - .rename(columns={"index": "model", "variable": "sensitivity", "value": "ratio"}) + .rename( + columns={"index": "model", "variable": "sensitivity", "value": "ratio"} + ) ) ratios["sensitivity"] = ratios["sensitivity"].replace( {"t_fwd": "forward", "t_adj": "adjoint"} diff --git a/tests/benchmark-models/test_petab_benchmark.py b/tests/benchmark-models/test_petab_benchmark.py index 936b3a5ab5..0b3c6d80e0 100755 --- a/tests/benchmark-models/test_petab_benchmark.py +++ b/tests/benchmark-models/test_petab_benchmark.py @@ -19,7 +19,9 @@ RTOL: float = 1e-2 benchmark_path = ( - Path(__file__).parent.parent.parent / "Benchmark-Models-PEtab" / "Benchmark-Models" + Path(__file__).parent.parent.parent + / "Benchmark-Models-PEtab" + / "Benchmark-Models" ) # reuse compiled models from test_benchmark_collection.sh benchmark_outdir = Path(__file__).parent.parent.parent / "test_bmc" @@ -67,11 +69,15 @@ def test_benchmark_gradient(model, scale): # only fail on linear scale pytest.skip() - petab_problem = petab.Problem.from_yaml(benchmark_path / model / (model + ".yaml")) + petab_problem = petab.Problem.from_yaml( + benchmark_path / model / (model + ".yaml") + ) petab.flatten_timepoint_specific_output_overrides(petab_problem) # Only compute gradient for estimated parameters. - parameter_df_free = petab_problem.parameter_df.loc[petab_problem.x_free_ids] + parameter_df_free = petab_problem.parameter_df.loc[ + petab_problem.x_free_ids + ] parameter_ids = list(parameter_df_free.index) # Setup AMICI objects. diff --git a/tests/benchmark-models/test_petab_model.py b/tests/benchmark-models/test_petab_model.py index 0135faa5bf..cf85147535 100755 --- a/tests/benchmark-models/test_petab_model.py +++ b/tests/benchmark-models/test_petab_model.py @@ -91,7 +91,8 @@ def parse_cli_args(): "-o", "--simulation-file", dest="simulation_file", - help="File to write simulation result to, in PEtab" "measurement table format.", + help="File to write simulation result to, in PEtab" + "measurement table format.", ) return parser.parse_args() @@ -167,10 +168,14 @@ def main(): times["np"] = sum(problem.parameter_df[petab.ESTIMATE]) - pd.Series(times).to_csv(f"./tests/benchmark-models/{args.model_name}_benchmark.csv") + pd.Series(times).to_csv( + f"./tests/benchmark-models/{args.model_name}_benchmark.csv" + ) for rdata in rdatas: - assert rdata.status == amici.AMICI_SUCCESS, f"Simulation failed for {rdata.id}" + assert ( + rdata.status == amici.AMICI_SUCCESS + ), f"Simulation failed for {rdata.id}" # create simulation PEtab table sim_df = rdatas_to_measurement_df( @@ -217,7 +222,8 @@ def main(): if np.isclose(llh, ref_llh, rtol=rtol, atol=atol): logger.info( - f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}." + tolstr + f"Computed llh {llh:.4e} matches reference {ref_llh:.4e}." + + tolstr ) else: logger.error( diff --git a/tests/conftest.py b/tests/conftest.py index 4d2f5521ff..9e90400518 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -48,7 +48,9 @@ def get_all_semantic_case_ids(): suite""" pattern = re.compile(r"\d{5}") return sorted( - str(x.name) for x in SBML_SEMANTIC_CASES_DIR.iterdir() if pattern.match(x.name) + str(x.name) + for x in SBML_SEMANTIC_CASES_DIR.iterdir() + if pattern.match(x.name) ) @@ -78,7 +80,9 @@ def pytest_generate_tests(metafunc): def pytest_sessionfinish(session, exitstatus): """Process test results""" global passed_ids - terminalreporter = session.config.pluginmanager.get_plugin("terminalreporter") + terminalreporter = session.config.pluginmanager.get_plugin( + "terminalreporter" + ) terminalreporter.ensure_newline() # parse test names to get passed case IDs (don't know any better way to # access fixture values) @@ -100,9 +104,14 @@ def write_passed_tags(passed_ids, out=sys.stdout): passed_component_tags |= cur_component_tags passed_test_tags |= cur_test_tags - out.write("\nAt least one test with the following component tags has " "passed:\n") + out.write( + "\nAt least one test with the following component tags has " + "passed:\n" + ) out.write(" " + "\n ".join(sorted(passed_component_tags))) - out.write("\n\nAt least one test with the following test tags has " "passed:\n") + out.write( + "\n\nAt least one test with the following test tags has " "passed:\n" + ) out.write(" " + "\n ".join(sorted(passed_test_tags))) @@ -132,7 +141,9 @@ def get_tags_for_test(test_id: str) -> Tuple[Set[str], Set[str]]: test_tags = set() for line in f: if line.startswith("testTags:"): - test_tags = set(re.split(r"[ ,:]", line[len("testTags:") :].strip())) + test_tags = set( + re.split(r"[ ,:]", line[len("testTags:") :].strip()) + ) test_tags.discard("") if line.startswith("componentTags:"): component_tags = set( diff --git a/tests/generateTestConfig/example.py b/tests/generateTestConfig/example.py index 12e5d2e53c..963fdf64ce 100644 --- a/tests/generateTestConfig/example.py +++ b/tests/generateTestConfig/example.py @@ -18,7 +18,9 @@ def dict2hdf5(object, dictionary): dtype = "f8" else: dtype = " List[int]: def pytest_addoption(parser): """Add pytest CLI options""" parser.addoption("--petab-cases", help="Test cases to run") - parser.addoption("--only-pysb", help="Run only PySB tests", action="store_true") + parser.addoption( + "--only-pysb", help="Run only PySB tests", action="store_true" + ) parser.addoption( "--only-sbml", help="Run only SBML tests", @@ -44,7 +46,10 @@ def pytest_generate_tests(metafunc): """Parameterize tests""" # Run for all PEtab test suite cases - if "case" in metafunc.fixturenames and "model_type" in metafunc.fixturenames: + if ( + "case" in metafunc.fixturenames + and "model_type" in metafunc.fixturenames + ): # Get CLI option cases = metafunc.config.getoption("--petab-cases") if cases: @@ -59,7 +64,9 @@ def pytest_generate_tests(metafunc): (case, "sbml", version) for version in ("v1.0.0", "v2.0.0") for case in ( - test_numbers if test_numbers else get_cases("sbml", version=version) + test_numbers + if test_numbers + else get_cases("sbml", version=version) ) ] elif metafunc.config.getoption("--only-pysb"): diff --git a/tests/petab_test_suite/test_petab_suite.py b/tests/petab_test_suite/test_petab_suite.py index 5bd33c1b6c..35ee3adcfc 100755 --- a/tests/petab_test_suite/test_petab_suite.py +++ b/tests/petab_test_suite/test_petab_suite.py @@ -57,7 +57,9 @@ def _test_case(case, model_type, version): # compile amici model if case.startswith("0006"): petab.flatten_timepoint_specific_output_overrides(problem) - model_name = f"petab_{model_type}_test_case_{case}" f"_{version.replace('.', '_')}" + model_name = ( + f"petab_{model_type}_test_case_{case}" f"_{version.replace('.', '_')}" + ) model_output_dir = f"amici_models/{model_name}" model = import_petab_problem( petab_problem=problem, @@ -79,9 +81,13 @@ def _test_case(case, model_type, version): rdatas = ret["rdatas"] chi2 = sum(rdata["chi2"] for rdata in rdatas) llh = ret["llh"] - simulation_df = rdatas_to_measurement_df(rdatas, model, problem.measurement_df) + simulation_df = rdatas_to_measurement_df( + rdatas, model, problem.measurement_df + ) petab.check_measurement_df(simulation_df, problem.observable_df) - simulation_df = simulation_df.rename(columns={petab.MEASUREMENT: petab.SIMULATION}) + simulation_df = simulation_df.rename( + columns={petab.MEASUREMENT: petab.SIMULATION} + ) simulation_df[petab.TIME] = simulation_df[petab.TIME].astype(int) solution = petabtests.load_solution(case, model_type, version=version) gt_chi2 = solution[petabtests.CHI2] @@ -118,13 +124,17 @@ def _test_case(case, model_type, version): ): logger.log( logging.DEBUG, - f"x_ss: {model.getStateIds()} " f"{[rdata.x_ss for rdata in rdatas]}", + f"x_ss: {model.getStateIds()} " + f"{[rdata.x_ss for rdata in rdatas]}", + ) + logger.log( + logging.ERROR, f"Expected simulations:\n{gt_simulation_dfs}" ) - logger.log(logging.ERROR, f"Expected simulations:\n{gt_simulation_dfs}") logger.log(logging.ERROR, f"Actual simulations:\n{simulation_df}") logger.log( logging.DEBUG if chi2s_match else logging.ERROR, - f"CHI2: simulated: {chi2}, expected: {gt_chi2}," f" match = {chi2s_match}", + f"CHI2: simulated: {chi2}, expected: {gt_chi2}," + f" match = {chi2s_match}", ) logger.log( logging.DEBUG if simulations_match else logging.ERROR, @@ -135,7 +145,9 @@ def _test_case(case, model_type, version): if not all([llhs_match, simulations_match]) or not chi2s_match: logger.error(f"Case {case} failed.") - raise AssertionError(f"Case {case}: Test results do not match " "expectations") + raise AssertionError( + f"Case {case}: Test results do not match " "expectations" + ) logger.info(f"Case {case} passed.") diff --git a/tests/testSBMLSuite.py b/tests/testSBMLSuite.py index 25f000ca55..cfad477ac4 100755 --- a/tests/testSBMLSuite.py +++ b/tests/testSBMLSuite.py @@ -44,7 +44,9 @@ def sbml_test_dir(): sys.path = old_path -def test_sbml_testsuite_case(test_number, result_path, sbml_semantic_cases_dir): +def test_sbml_testsuite_case( + test_number, result_path, sbml_semantic_cases_dir +): test_id = format_test_id(test_number) model_dir = None @@ -92,7 +94,9 @@ def test_sbml_testsuite_case(test_number, result_path, sbml_semantic_cases_dir): raise RuntimeError("Simulation failed unexpectedly") # verify - simulated = verify_results(settings, rdata, results, wrapper, model, atol, rtol) + simulated = verify_results( + settings, rdata, results, wrapper, model, atol, rtol + ) # record results write_result_file(simulated, test_id, result_path) @@ -117,7 +121,10 @@ def verify_results(settings, rdata, expected, wrapper, model, atol, rtol): # collect states simulated = pd.DataFrame( rdata["y"], - columns=[obs["name"] for obs in wrapper.symbols[SymbolId.OBSERVABLE].values()], + columns=[ + obs["name"] + for obs in wrapper.symbols[SymbolId.OBSERVABLE].values() + ], ) simulated["time"] = rdata["ts"] # collect parameters @@ -137,7 +144,10 @@ def verify_results(settings, rdata, expected, wrapper, model, atol, rtol): # hasOnlySubstanceUnits="true", but then request results as concentrations. requested_concentrations = [ s - for s in settings["concentration"].replace(" ", "").replace("\n", "").split(",") + for s in settings["concentration"] + .replace(" ", "") + .replace("\n", "") + .split(",") if s ] # We only need to convert species that have only substance units @@ -147,7 +157,8 @@ def verify_results(settings, rdata, expected, wrapper, model, atol, rtol): **wrapper.symbols[SymbolId.SPECIES], **wrapper.symbols[SymbolId.ALGEBRAIC_STATE], }.items() - if str(state_id) in requested_concentrations and state.get("amount", False) + if str(state_id) in requested_concentrations + and state.get("amount", False) ] amounts_to_concentrations( concentration_species, wrapper, simulated, requested_concentrations @@ -220,7 +231,9 @@ def concentrations_to_amounts( # Species with OnlySubstanceUnits don't have to be converted as long # as we don't request concentrations for them. Only applies when # called from amounts_to_concentrations. - if (is_amt and species not in requested_concentrations) or comp is None: + if ( + is_amt and species not in requested_concentrations + ) or comp is None: continue simulated.loc[:, species] *= simulated.loc[ @@ -228,7 +241,9 @@ def concentrations_to_amounts( ] -def write_result_file(simulated: pd.DataFrame, test_id: str, result_path: Path): +def write_result_file( + simulated: pd.DataFrame, test_id: str, result_path: Path +): """ Create test result file for upload to http://raterule.caltech.edu/Facilities/Database @@ -245,10 +260,14 @@ def get_amount_and_variables(settings): """Read amount and species from settings file""" # species for which results are expected as amounts - amount_species = settings["amount"].replace(" ", "").replace("\n", "").split(",") + amount_species = ( + settings["amount"].replace(" ", "").replace("\n", "").split(",") + ) # IDs of all variables for which results are expected/provided - variables = settings["variables"].replace(" ", "").replace("\n", "").split(",") + variables = ( + settings["variables"].replace(" ", "").replace("\n", "").split(",") + ) return amount_species, variables