From 6c2d9cca9f2340bd8591c76fbceb71c2e4d2e671 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Mon, 6 Oct 2025 10:04:38 +0200 Subject: [PATCH] Add `Model.simulate` Add `Model.simulate()` as a convenience function to run simulations without having to create a `Solver` object explicitly. This is a wrapper for both `amici.run_simulation` and `amici.run_simulations`, depending on the type of the `edata` argument. It also supports passing some `Solver` options as keyword arguments. --- CHANGELOG.md | 8 +- .../getting_started/GettingStarted.ipynb | 87 +++++----- python/sdist/amici/__init__.py | 2 +- python/sdist/amici/swig_wrappers.py | 77 ++++++++- python/tests/test_sbml_import.py | 15 +- swig/amici.i | 5 +- swig/model.i | 149 ++++++++++++++++++ 7 files changed, 277 insertions(+), 66 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 83a419fd06..836df019ea 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,6 +34,8 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni * For a more consistent API, all function names are now snake_case instead of camelCase. * `Model.getSolver` has been renamed to `Model.create_solver`. +* `amici.runAmiciSimulation` and `amici.runAmiciSimulations` have been renamed + to `amici.run_simulation` and `amici.run_simulations`. * The following deprecated functionality has been removed: * The complete MATLAB interface has been removed. * `NonlinearSolverIteration::functional` has been removed, @@ -54,7 +56,11 @@ See also our [versioning policy](https://amici.readthedocs.io/en/latest/versioni `DataArray`s include the identifiers and are often more convenient than the plain numpy arrays. This allows for easy subselection and plotting of the results, and conversion to DataFrames. - +* `Model.simulate()` has been added as a convenience function to run + simulations without having to create a `Solver` object explicitly. + This is a wrapper for both `amici.run_simulation` and + `amici.run_simulations`, depending on the type of the `edata` argument. + It also supports passing some `Solver` options as keyword arguments. ## v0.X Series diff --git a/doc/examples/getting_started/GettingStarted.ipynb b/doc/examples/getting_started/GettingStarted.ipynb index c58e7f2ddc..1e01d26e81 100644 --- a/doc/examples/getting_started/GettingStarted.ipynb +++ b/doc/examples/getting_started/GettingStarted.ipynb @@ -21,14 +21,14 @@ }, { "cell_type": "code", - "execution_count": 1, "metadata": {}, - "outputs": [], "source": [ "import amici\n", "\n", "sbml_importer = amici.SbmlImporter(\"model_steadystate_scaled.xml\")" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -39,14 +39,14 @@ }, { "cell_type": "code", - "execution_count": 2, "metadata": {}, - "outputs": [], "source": [ "model_name = \"model_steadystate\"\n", "model_dir = \"model_dir\"\n", "sbml_importer.sbml2amici(model_name, model_dir)" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -58,9 +58,7 @@ }, { "cell_type": "code", - "execution_count": 3, "metadata": {}, - "outputs": [], "source": [ "# load the model module\n", "model_module = amici.import_model_module(model_name, model_dir)\n", @@ -68,7 +66,9 @@ "model = model_module.get_model()\n", "# instantiate solver\n", "solver = model.create_solver()" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -77,10 +77,10 @@ }, { "cell_type": "code", - "execution_count": 4, "metadata": {}, + "source": "model.set_parameter_by_name(\"p1\", 1e-3)", "outputs": [], - "source": "model.set_parameter_by_name(\"p1\", 1e-3)" + "execution_count": null }, { "cell_type": "markdown", @@ -89,10 +89,10 @@ }, { "cell_type": "code", - "execution_count": 5, "metadata": {}, + "source": "solver.set_absolute_tolerance(1e-10)", "outputs": [], - "source": "solver.set_absolute_tolerance(1e-10)" + "execution_count": null }, { "cell_type": "markdown", @@ -104,18 +104,18 @@ { "cell_type": "markdown", "metadata": {}, - "source": "Model simulations can be executed using the [amici.run_simulation](https://amici.readthedocs.io/en/latest/generated/amici.html#amici.run_simulation) routine. By default, the model does not contain any timepoints for which the model is to be simulated. Here we define a simulation timecourse with two timepoints at `0` and `1` and then run the simulation." + "source": "Model simulations can be executed using the [Model.simulate](https://amici.readthedocs.io/en/latest/generated/amici.amici.html#amici.amici.Model.simulate) method (or, alternatively, [amici.run_simulation](https://amici.readthedocs.io/en/latest/generated/amici.html#amici.run_simulation)). By default, the model does not contain any output timepoints for which the model is to be simulated. Here we define a simulation timecourse with two output timepoints at `0` and `1` and then run the simulation." }, { "cell_type": "code", - "execution_count": 6, "metadata": {}, - "outputs": [], "source": [ "# set timepoints\n", "model.set_timepoints([0, 1])\n", - "rdata = amici.run_simulation(model, solver)" - ] + "rdata = model.simulate(solver=solver)" + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -126,24 +126,12 @@ }, { "cell_type": "code", - "execution_count": 7, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "array([[0.1 , 0.4 , 0.7 ],\n", - " [0.98208413, 0.51167992, 0.10633388]])" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], + "metadata": {}, "source": [ "rdata.x" - ] + ], + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", @@ -152,21 +140,22 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "('x1', 'x2', 'x3')" - ] - }, - "execution_count": 8, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": "model.get_state_names()" + "metadata": {}, + "source": "model.get_state_names()", + "outputs": [], + "execution_count": null + }, + { + "metadata": {}, + "cell_type": "markdown", + "source": "For convenience, most results stored in `ReturnData` can also be retrieved as [xarray.DataArray](https://docs.xarray.dev/en/stable/index.html) objects that already include the respective row and column names. This can be accessed via the `xr` attribute of `ReturnData`. Here, we access the model state `x` as `DataArray` object to convert it to a `pandas.DataFrame`:" + }, + { + "metadata": {}, + "cell_type": "code", + "source": "rdata.xr.x.to_pandas()", + "outputs": [], + "execution_count": null }, { "cell_type": "markdown", diff --git a/python/sdist/amici/__init__.py b/python/sdist/amici/__init__.py index 4a6a1f9ecb..c4bfa62b2d 100644 --- a/python/sdist/amici/__init__.py +++ b/python/sdist/amici/__init__.py @@ -150,7 +150,7 @@ def get_model(self) -> amici.Model: """Create a model instance.""" ... - AmiciModel = Union[amici.Model, amici.ModelPtr] + AmiciModel = amici.Model | amici.ModelPtr else: ModelModule = ModuleType diff --git a/python/sdist/amici/swig_wrappers.py b/python/sdist/amici/swig_wrappers.py index 126c7bfd0c..e388285f0d 100644 --- a/python/sdist/amici/swig_wrappers.py +++ b/python/sdist/amici/swig_wrappers.py @@ -1,8 +1,11 @@ """Convenience wrappers for the swig interface""" +from __future__ import annotations import logging import warnings from typing import Any +from collections.abc import Sequence +import contextlib import amici import amici.amici as amici_swig @@ -12,8 +15,11 @@ AmiciExpDataVector, AmiciModel, AmiciSolver, + SensitivityMethod, + SensitivityOrder, + Solver, ) -from . import numpy +from . import numpy, ReturnDataView from .logging import get_logger logger = get_logger(__name__, log_level=logging.DEBUG) @@ -33,7 +39,7 @@ def run_simulation( model: AmiciModel, solver: AmiciSolver, edata: AmiciExpData | None = None, -) -> "numpy.ReturnDataView": +) -> ReturnDataView: """ Convenience wrapper around :py:func:`amici.amici.run_simulation` (generated by swig) @@ -79,7 +85,7 @@ def run_simulations( edata_list: AmiciExpDataVector, failfast: bool = True, num_threads: int = 1, -) -> list["numpy.ReturnDataView"]: +) -> list[ReturnDataView]: """ Convenience wrapper for loops of amici.runAmiciSimulation @@ -267,5 +273,70 @@ def _ids_and_names_to_rdata( f"{entity_type.lower()}_{name_or_id.lower()}", names_or_ids, ) + rdata.state_ids_solver = model.get_state_ids_solver() rdata.state_names_solver = model.get_state_names_solver() + + +@contextlib.contextmanager +def _solver_settings(solver, sensi_method=None, sensi_order=None): + """Context manager to temporarily apply solver settings.""" + old_method = old_order = None + + if sensi_method is not None: + old_method = solver.get_sensitivity_method() + if isinstance(sensi_method, str): + sensi_method = SensitivityMethod[sensi_method] + solver.set_sensitivity_method(sensi_method) + + if sensi_order is not None: + old_order = solver.get_sensitivity_order() + if isinstance(sensi_order, str): + sensi_order = SensitivityOrder[sensi_order] + solver.set_sensitivity_order(sensi_order) + + try: + yield solver + finally: + if old_method is not None: + solver.set_sensitivity_method(old_method) + if old_order is not None: + solver.set_sensitivity_order(old_order) + + +def _Model__simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpData | AmiciExpDataVector | None = None, + failfast: bool = True, + num_threads: int = 1, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> ReturnDataView | list[ReturnDataView]: + """ + For use in `swig/model.i` to avoid code duplication in subclasses. + + Keep in sync with `Model.simulate` and `ModelPtr.simulate`. + + """ + if solver is None: + solver = self.create_solver() + + with _solver_settings( + solver=solver, sensi_method=sensi_method, sensi_order=sensi_order + ): + if isinstance(edata, Sequence): + return run_simulations( + model=_get_ptr(self), + solver=_get_ptr(solver), + edata_list=edata, + failfast=failfast, + num_threads=num_threads, + ) + + return run_simulation( + model=_get_ptr(self), + solver=_get_ptr(solver), + edata=_get_ptr(edata), + ) diff --git a/python/tests/test_sbml_import.py b/python/tests/test_sbml_import.py index 178ae762ea..fa52232d4f 100644 --- a/python/tests/test_sbml_import.py +++ b/python/tests/test_sbml_import.py @@ -118,10 +118,7 @@ def test_nosensi(tempdir): model = model_module.get_model() model.set_timepoints(np.linspace(0, 60, 61)) - solver = model.create_solver() - solver.set_sensitivity_order(amici.SensitivityOrder.first) - solver.set_sensitivity_method(amici.SensitivityMethod.forward) - rdata = amici.run_simulation(model, solver) + rdata = model.simulate(sensi_order="first", sensi_method="forward") assert rdata.status == amici.AMICI_ERROR @@ -598,7 +595,7 @@ def test_likelihoods(model_test_likelihoods): # run model once to create an edata - rdata = amici.run_simulation(model, solver) + rdata = model.simulate(solver=solver) sigmas = rdata["y"].max(axis=0) * 0.05 edata = amici.ExpData(rdata, sigmas, []) # just make all observables positive since some are logarithmic @@ -606,7 +603,7 @@ def test_likelihoods(model_test_likelihoods): edata = amici.ExpData(rdata, sigmas, []) # and now run for real and also compute likelihood values - rdata = amici.run_simulations(model, solver, [edata])[0] + rdata = model.simulate(solver=solver, edata=[edata])[0] # check if the values make overall sense assert np.isfinite(rdata["llh"]) @@ -1054,8 +1051,7 @@ def test_regression_2700(tempdir): model_module = import_model_module(model_name, tempdir) model = model_module.get_model() model.set_timepoints([0, 1, 2]) - solver = model.create_solver() - rdata = amici.run_simulation(model, solver) + rdata = model.simulate() assert np.all(rdata.by_id("pp") == [1, 1, 1]) @@ -1093,8 +1089,7 @@ def test_heaviside_init_values_and_bool_to_float_conversion(tempdir): model = model_module.get_model() model.set_timepoints([0, 1, 2]) - solver = model.create_solver() - rdata = amici.run_simulation(model, solver) + rdata = model.simulate() assert np.all(rdata.by_id("a") == np.array([2, 2, 2])), rdata.by_id("a") assert np.all(rdata.by_id("b") == np.array([0, 1, 1])), rdata.by_id("b") diff --git a/swig/amici.i b/swig/amici.i index a956033045..7b7929c078 100644 --- a/swig/amici.i +++ b/swig/amici.i @@ -373,11 +373,12 @@ if sys.platform == 'win32': // import additional types for typehints // also import np for use in __repr__ functions %pythonbegin %{ -from typing import TYPE_CHECKING, Iterable, Union +from typing import TYPE_CHECKING, Iterable, Union, overload from collections.abc import Sequence import numpy as np if TYPE_CHECKING: import numpy + from .numpy import ReturnDataView %} %pythoncode %{ @@ -418,7 +419,7 @@ __all__ = [ x for x in dir(sys.modules[__name__]) if not x.startswith('_') - and x not in {"np", "sys", "os", "numpy", "IntEnum", "enum", "pi", "TYPE_CHECKING", "Iterable", "Sequence", "Path"} + and x not in {"np", "sys", "os", "numpy", "IntEnum", "enum", "pi", "TYPE_CHECKING", "Iterable", "Sequence", "Path", "Union", "overload"} ] %} diff --git a/swig/model.i b/swig/model.i index 4a7436e0ec..b3699aaa26 100644 --- a/swig/model.i +++ b/swig/model.i @@ -113,6 +113,80 @@ using namespace amici; %pythoncode %{ def __deepcopy__(self, memo): return self.clone() + +@overload +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpData | None = None, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> ReturnDataView: ... + + +@overload +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpDataVector | None = None, + failfast: bool = True, + num_threads: int = 1, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> list[ReturnDataView]: ... + + +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpData | AmiciExpDataVector | None = None, + failfast: bool = True, + num_threads: int = 1, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> ReturnDataView | list[ReturnDataView]: + """Simulate model with given solver and experimental data. + + :param solver: + Solver to use for simulation. Defaults to :meth:`Model.get_solver`. + :param edata: + Experimental data to use for simulation. + A single :class:`ExpData` instance or a sequence of such instances. + If `None`, no experimental data is used and the model is simulated + as is. + :param sensi_method: + Sensitivity method to use for simulation. + If `None`, the solver's current sensitivity method is used. + :param sensi_order: + Sensitivity order to use for simulation. + If `None`, the solvers's current sensitivity order is used. + :param failfast: + Whether to stop simulations on first failure. + Only relevant if `edata` is a sequence of :class:`ExpData` instances. + :param num_threads: + Number of threads to use for simulation. + Only relevant if AMICI was compiled with OpenMP support and if `edata` + is a sequence of :class:`ExpData` instances. + :return: + A single :class:`ReturnDataView` instance containing the simulation + results if `edata` is a single :class:`ExpData` instance or `None`. + If `edata` is a sequence of :class:`ExpData` instances, a list of + :class:`ReturnDataView` instances is returned. + """ + from .swig_wrappers import _Model__simulate + + return _Model__simulate( + self, + solver=solver, + edata=edata, + failfast=failfast, + num_threads=num_threads, + sensi_method=sensi_method, + sensi_order=sensi_order, + ) %} }; @@ -120,6 +194,81 @@ def __deepcopy__(self, memo): %pythoncode %{ def __deepcopy__(self, memo): return self.clone() + + +@overload +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpData | None = None, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> ReturnDataView: ... + + +@overload +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpDataVector | None = None, + failfast: bool = True, + num_threads: int = 1, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> list[ReturnDataView]: ... + + +def simulate( + self: AmiciModel, + *, + solver: Solver | None = None, + edata: AmiciExpData | AmiciExpDataVector | None = None, + failfast: bool = True, + num_threads: int = 1, + sensi_method: SensitivityMethod | str = None, + sensi_order: SensitivityOrder | str = None, +) -> ReturnDataView | list[ReturnDataView]: + """Simulate model with given solver and experimental data. + + :param solver: + Solver to use for simulation. Defaults to :meth:`Model.get_solver`. + :param edata: + Experimental data to use for simulation. + A single :class:`ExpData` instance or a sequence of such instances. + If `None`, no experimental data is used and the model is simulated + as is. + :param sensi_method: + Sensitivity method to use for simulation. + If `None`, the solver's current sensitivity method is used. + :param sensi_order: + Sensitivity order to use for simulation. + If `None`, the solvers's current sensitivity order is used. + :param failfast: + Whether to stop simulations on first failure. + Only relevant if `edata` is a sequence of :class:`ExpData` instances. + :param num_threads: + Number of threads to use for simulation. + Only relevant if AMICI was compiled with OpenMP support and if `edata` + is a sequence of :class:`ExpData` instances. + :return: + A single :class:`ReturnDataView` instance containing the simulation + results if `edata` is a single :class:`ExpData` instance or `None`. + If `edata` is a sequence of :class:`ExpData` instances, a list of + :class:`ReturnDataView` instances is returned. + """ + from .swig_wrappers import _Model__simulate + + return _Model__simulate( + self, + solver=solver, + edata=edata, + failfast=failfast, + num_threads=num_threads, + sensi_method=sensi_method, + sensi_order=sensi_order, + ) %} };