Add Model.simulate

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
 

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,17 +58,17 @@
   },
   {
    "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",
     "# instantiate model\n",
     "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",

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
 

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),
+        )

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,15 +595,15 @@ 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
     while min(edata.get_observed_data()) < 0:
         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")

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"}
 ]
 
 %}