From 5da8e21015299e20ad033d42d965585fd4891546 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 27 Jul 2023 11:30:08 +0200 Subject: [PATCH 01/10] Evaluate and plot symbolic expressions based on simulation results WIP --- python/sdist/amici/plotting.py | 51 +++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index da718c1ec7..aa93bd7011 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -3,15 +3,18 @@ -------- Plotting related functions """ -from typing import Iterable, Optional +from typing import Iterable, Optional, Sequence, Union import matplotlib.pyplot as plt import pandas as pd import seaborn as sns +import sympy as sp from matplotlib.axes import Axes from . import Model, ReturnDataView +StrOrExpr = Union[str, sp.Expr] + def plot_state_trajectories( rdata: ReturnDataView, @@ -115,3 +118,49 @@ def plot_jacobian(rdata: ReturnDataView): # backwards compatibility plotStateTrajectories = plot_state_trajectories plotObservableTrajectories = plot_observable_trajectories + + +def evaluate(expr: StrOrExpr, rdata: ReturnDataView) -> "numpy.array": + """Evaluate a symbolic expression based on the given simulation outputs. + + :param expr: + A symbolic expression, e.g. a sympy expression or a string that can be sympified. + Can include state variable, expression, and observable IDs, depending on whether + the respective data is available in the simulation results. + Parameters are not yet supported. + :param rdata: + The simulation results. + + :return: + The evaluated expression for the simulation output timepoints. + """ + from sympy.utilities.lambdify import lambdify + + if isinstance(expr, str): + expr = sp.sympify(expr) + + arg_names = list(sorted(expr.free_symbols, key=lambda x: x.name)) + func = lambdify(arg_names, expr, "numpy") + args = [rdata.by_id(arg.name) for arg in arg_names] + return func(*args) + + +def plot_expressions( + exprs: Union[Sequence[StrOrExpr], StrOrExpr], rdata: ReturnDataView +) -> None: + """Plot the given expressions evaluated on the given simulation outputs. + + :param exprs: + A symbolic expression, e.g. a sympy expression or a string that can be sympified. + Can include state variable, expression, and observable IDs, depending on whether + the respective data is available in the simulation results. + Parameters are not yet supported. + :param rdata: + The simulation results. + """ + if not isinstance(exprs, Sequence): + exprs = [exprs] + + for expr in exprs: + plt.plot(rdata.t, evaluate(expr, rdata), label=str(expr)) + plt.legend() From 509ba9e27455de345310c98436ba50fc30899032 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Thu, 27 Jul 2023 12:25:28 +0200 Subject: [PATCH 02/10] .. --- python/sdist/amici/plotting.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index aa93bd7011..e99be7d969 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -6,6 +6,7 @@ from typing import Iterable, Optional, Sequence, Union import matplotlib.pyplot as plt +import numpy as np import pandas as pd import seaborn as sns import sympy as sp @@ -120,7 +121,7 @@ def plot_jacobian(rdata: ReturnDataView): plotObservableTrajectories = plot_observable_trajectories -def evaluate(expr: StrOrExpr, rdata: ReturnDataView) -> "numpy.array": +def evaluate(expr: StrOrExpr, rdata: ReturnDataView) -> np.array: """Evaluate a symbolic expression based on the given simulation outputs. :param expr: From f96640d31f74661e0222c7844b958dba82678c3c Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Sep 2023 21:57:34 +0200 Subject: [PATCH 03/10] refactor --- python/sdist/amici/numpy.py | 28 ++++++++++++++++++++++++++++ python/sdist/amici/plotting.py | 30 +----------------------------- 2 files changed, 29 insertions(+), 29 deletions(-) diff --git a/python/sdist/amici/numpy.py b/python/sdist/amici/numpy.py index 23ebfdbbc4..d9b34b6447 100644 --- a/python/sdist/amici/numpy.py +++ b/python/sdist/amici/numpy.py @@ -10,9 +10,12 @@ import amici import numpy as np +import sympy as sp from . import ExpData, ExpDataPtr, Model, ReturnData, ReturnDataPtr +StrOrExpr = Union[str, sp.Expr] + class SwigPtrView(collections.abc.Mapping): """ @@ -429,3 +432,28 @@ def _entity_type_from_id( return symbol raise KeyError(f"Unknown symbol {entity_id}.") + + +def evaluate(expr: StrOrExpr, rdata: ReturnDataView) -> np.array: + """Evaluate a symbolic expression based on the given simulation outputs. + + :param expr: + A symbolic expression, e.g. a sympy expression or a string that can be sympified. + Can include state variable, expression, and observable IDs, depending on whether + the respective data is available in the simulation results. + Parameters are not yet supported. + :param rdata: + The simulation results. + + :return: + The evaluated expression for the simulation output timepoints. + """ + from sympy.utilities.lambdify import lambdify + + if isinstance(expr, str): + expr = sp.sympify(expr) + + arg_names = list(sorted(expr.free_symbols, key=lambda x: x.name)) + func = lambdify(arg_names, expr, "numpy") + args = [rdata.by_id(arg.name) for arg in arg_names] + return func(*args) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index e99be7d969..330e74edea 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -6,15 +6,12 @@ from typing import Iterable, Optional, Sequence, Union import matplotlib.pyplot as plt -import numpy as np import pandas as pd import seaborn as sns -import sympy as sp from matplotlib.axes import Axes from . import Model, ReturnDataView - -StrOrExpr = Union[str, sp.Expr] +from .numpy import StrOrExpr, evaluate def plot_state_trajectories( @@ -121,31 +118,6 @@ def plot_jacobian(rdata: ReturnDataView): plotObservableTrajectories = plot_observable_trajectories -def evaluate(expr: StrOrExpr, rdata: ReturnDataView) -> np.array: - """Evaluate a symbolic expression based on the given simulation outputs. - - :param expr: - A symbolic expression, e.g. a sympy expression or a string that can be sympified. - Can include state variable, expression, and observable IDs, depending on whether - the respective data is available in the simulation results. - Parameters are not yet supported. - :param rdata: - The simulation results. - - :return: - The evaluated expression for the simulation output timepoints. - """ - from sympy.utilities.lambdify import lambdify - - if isinstance(expr, str): - expr = sp.sympify(expr) - - arg_names = list(sorted(expr.free_symbols, key=lambda x: x.name)) - func = lambdify(arg_names, expr, "numpy") - args = [rdata.by_id(arg.name) for arg in arg_names] - return func(*args) - - def plot_expressions( exprs: Union[Sequence[StrOrExpr], StrOrExpr], rdata: ReturnDataView ) -> None: From 5d69d170deb1f9f06f3efe961f9953f53361a290 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Sep 2023 22:16:25 +0200 Subject: [PATCH 04/10] Add test --- python/tests/test_rdata.py | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/python/tests/test_rdata.py b/python/tests/test_rdata.py index 29ea401932..ac7659f363 100644 --- a/python/tests/test_rdata.py +++ b/python/tests/test_rdata.py @@ -2,7 +2,8 @@ import amici import numpy as np import pytest -from numpy.testing import assert_array_equal +from amici.numpy import evaluate +from numpy.testing import assert_almost_equal, assert_array_equal @pytest.fixture(scope="session") @@ -39,3 +40,27 @@ def test_rdata_by_id(rdata_by_id_fixture): assert_array_equal( rdata.by_id(model.getStateIds()[1], "sx", model), rdata.sx[:, :, 1] ) + + +def test_evaluate(rdata_by_id_fixture): + # get IDs of model components + model, rdata = rdata_by_id_fixture + expr0_id = model.getExpressionIds()[0] + state1_id = model.getStateIds()[1] + observable0_id = model.getObservableIds()[0] + + # ensure `evaluate` works for atoms + expr0 = rdata.by_id(expr0_id) + assert_array_equal(expr0, evaluate(expr0_id, rdata=rdata)) + + state1 = rdata.by_id(state1_id) + assert_array_equal(state1, evaluate(state1_id, rdata=rdata)) + + observable0 = rdata.by_id(observable0_id) + assert_array_equal(observable0, evaluate(observable0_id, rdata=rdata)) + + # ensure `evaluate` works for expressions + assert_almost_equal( + expr0 + state1 * observable0, + evaluate(f"{expr0_id} + {state1_id} * {observable0_id}", rdata=rdata), + ) From 8d6a15dd303b44b4d2df77c85df8400f1bcc568f Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Tue, 19 Sep 2023 23:03:45 +0200 Subject: [PATCH 05/10] Example --- .../ExampleSteadystate.ipynb | 26 +++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index 502174fe15..809e4c73dc 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -920,10 +920,32 @@ "source": [ "import amici.plotting\n", "\n", - "amici.plotting.plotStateTrajectories(rdata, model=None)\n", - "amici.plotting.plotObservableTrajectories(rdata, model=None)" + "amici.plotting.plot_state_trajectories(rdata, model=None)\n", + "amici.plotting.plot_observable_trajectories(rdata, model=None)" ] }, + { + "cell_type": "markdown", + "source": [ + "We can also evaluate symbolic expressions of model quantities using `amici.numpy.evaluate`, or directly plot the results using `amici.plotting.plot_expressions`:" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [ + "amici.plotting.plot_expressions(\n", + " \"observable_x1 + observable_x2 + observable_x3\", rdata=rdata\n", + ")" + ], + "metadata": { + "collapsed": false + } + }, { "cell_type": "markdown", "metadata": {}, From 929bd12a1d650f7aeff5a3d436cef4485899ce4a Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Sep 2023 10:57:27 +0200 Subject: [PATCH 06/10] fixup --- python/sdist/amici/plotting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index 330e74edea..17fe4d1c5a 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -131,7 +131,7 @@ def plot_expressions( :param rdata: The simulation results. """ - if not isinstance(exprs, Sequence): + if not isinstance(exprs, Sequence) or isinstance(exprs, str): exprs = [exprs] for expr in exprs: From c2f3f7dc0e0cfdd19a31f4e04c3d7668533889d8 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Sep 2023 14:12:59 +0200 Subject: [PATCH 07/10] Always execute --- python/examples/example_steadystate/ExampleSteadystate.ipynb | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index 809e4c73dc..dd60305bef 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -2045,7 +2045,10 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.7" + "version": "3.7.7", + "nbsphinx": { + "execute": "always" + } }, "toc": { "base_numbering": 1, From b21c4e0ac149d5d4001129e80cf3234865ea4aef Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Sep 2023 14:36:23 +0200 Subject: [PATCH 08/10] Always execute --- .../examples/example_steadystate/ExampleSteadystate.ipynb | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index dd60305bef..f3e88e85dd 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -2045,10 +2045,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.7", - "nbsphinx": { - "execute": "always" - } + "version": "3.7.7" }, "toc": { "base_numbering": 1, @@ -2062,6 +2059,9 @@ "toc_position": {}, "toc_section_display": true, "toc_window_display": false + }, + "nbsphinx": { + "execute": "always" } }, "nbformat": 4, From 5abd04564e502ccfc8e7d0ea23ae6db81c53afcd Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Sep 2023 16:19:36 +0200 Subject: [PATCH 09/10] xlabel --- python/sdist/amici/plotting.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/sdist/amici/plotting.py b/python/sdist/amici/plotting.py index 17fe4d1c5a..bd1f3a8ba1 100644 --- a/python/sdist/amici/plotting.py +++ b/python/sdist/amici/plotting.py @@ -136,4 +136,6 @@ def plot_expressions( for expr in exprs: plt.plot(rdata.t, evaluate(expr, rdata), label=str(expr)) + plt.legend() + plt.gca().set_xlabel("$t$") From 68503b0800f5c1446e9b1f4592dbe763df936221 Mon Sep 17 00:00:00 2001 From: Daniel Weindl Date: Wed, 20 Sep 2023 16:23:12 +0200 Subject: [PATCH 10/10] fixup --- python/examples/example_steadystate/ExampleSteadystate.ipynb | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/examples/example_steadystate/ExampleSteadystate.ipynb b/python/examples/example_steadystate/ExampleSteadystate.ipynb index 4c28e5dfed..b57ed522aa 100644 --- a/python/examples/example_steadystate/ExampleSteadystate.ipynb +++ b/python/examples/example_steadystate/ExampleSteadystate.ipynb @@ -354,7 +354,7 @@ "source": [ "### Importing the module and loading the model\n", "\n", - "If everything went well, we need to add the previously selected model output directory to our PYTHON_PATH and are then ready to load newly generated model:" + "If everything went well, we can now import the newly generated Python module containing our model:" ] }, { @@ -392,7 +392,7 @@ "source": [ "model = model_module.getModel()\n", "\n", - "print(\"Model name:\", model.getName())\n", + "print(\"Model name: \", model.getName())\n", "print(\"Model parameters:\", model.getParameterIds())\n", "print(\"Model outputs: \", model.getObservableIds())\n", "print(\"Model states: \", model.getStateIds())"