From 3da923f4d955791978a6be07711f0a7f9b609aa2 Mon Sep 17 00:00:00 2001 From: "Felipe S. S. Schneider" Date: Sat, 28 Nov 2020 12:33:41 -0300 Subject: [PATCH] Allow the user to choose a simulation time --- .../0-basic-reaction-simulation.rst | 2 +- overreact/cli.py | 14 +++++++---- overreact/simulate.py | 24 ++++++++++++------- tests/test_simulate.py | 8 +++---- 4 files changed, 30 insertions(+), 18 deletions(-) diff --git a/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst b/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst index 39693a45..17a68eae 100644 --- a/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst +++ b/docsrc/tutorials/simulation/0-basic-reaction-simulation.rst @@ -84,4 +84,4 @@ We can see that the reaction went to full completion by checking the final concentrations: >>> y(y.t_max) -array([0.0000, 1.0000]) +array([0.00, 1.00]) diff --git a/overreact/cli.py b/overreact/cli.py index c7b23d3e..2fc80e1b 100644 --- a/overreact/cli.py +++ b/overreact/cli.py @@ -474,7 +474,9 @@ def _yield_kinetics(self): # TODO(schneiderfelipe): the following is inefficient but probably OK y0[self.model.scheme.compounds.index(name)] = quantity - y, r = api.get_y(dydt, y0=y0) + y, r = api.get_y( + dydt, y0=y0, method=self.method, rtol=self.rtol, atol=self.atol + ) conc_table = Table( Column("no", justify="right"), Column("compound", justify="left"), @@ -492,9 +494,13 @@ def _yield_kinetics(self): ) yield conc_table - # TODO(schneiderfelipe): we can get a max time now based on the - # changes through time: stop when the graph gets boring. - t = np.linspace(y.t_min, y.t_max) + eps = 1e-4 + t_max = y.t_max + if np.any(r(y.t_min) > eps): + while np.all(r(t_max) < eps): + t_max *= 0.9 + + t = np.linspace(y.t_min, t_max) if self.plot: for i, name in enumerate(self.model.scheme.compounds): if not core.is_transition_state(name): diff --git a/overreact/simulate.py b/overreact/simulate.py index 79eddc4c..ea70fd62 100644 --- a/overreact/simulate.py +++ b/overreact/simulate.py @@ -27,7 +27,9 @@ jnp = np -def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11): +def get_y( + dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11, max_time=5 * 60 +): """Simulate a reaction scheme from its rate function. This uses scipy's ``solve_ivp`` under the hood. @@ -49,6 +51,9 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11): normally unsuited. "Radau", "BDF" or "LSODA" are good choices. rtol, atol : array-like See `scipy.integrade.solve_ivp` for details. + max_time : float, optional + If `t_span` is not given, an interval will be estimated, but it can't + be larger than this parameter. Returns ------- @@ -71,10 +76,10 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11): used to produce a suitable vector of timepoints for, e.g., plotting: >>> y.t_min, y.t_max - (0.0, 10.0) + (0.0, 5.0) >>> t = np.linspace(y.t_min, y.t_max) >>> t - array([ 0. , 0.20408163, ..., 9.79591837, 10. ]) + array([..., 0.20408163, ..., 4.89795918, ...]) Both `y` and `r` can be used to check concentrations and rates in any point in time. In particular, both are vectorized: @@ -90,26 +95,27 @@ def get_y(dydt, y0, t_span=None, method="Radau", rtol=1e-5, atol=1e-11): y0 = np.asanyarray(y0) if t_span is None: - n_halflives = 10.0 + n_halflives = 5.0 # ensure < 5% remaining material in the worst case halflife_estimate = 1.0 if hasattr(dydt, "k"): halflife_estimate = ( np.max( [ - np.max(y0) / 2.0, # zeroth-order half-life - np.log(2.0), # first-order half-life - 1.0 / np.min(y0[np.nonzero(y0)]), # second-order half-life + np.max(y0) / 2.0, # zeroth-order halflife + np.log(2.0), # first-order halflife + 1.0 / np.min(y0[np.nonzero(y0)]), # second-order halflife ] ) / np.min(dydt.k) ) + logger.info(f"largest halflife guess = {halflife_estimate} s") t_span = [ 0.0, - n_halflives * halflife_estimate, + min(n_halflives * halflife_estimate, max_time), ] - logger.info(f"simulation time span = {t_span} s") + logger.info(f"simulation time span = {t_span} s") jac = None if hasattr(dydt, "jac"): diff --git a/tests/test_simulate.py b/tests/test_simulate.py index 45234358..f82e3db8 100644 --- a/tests/test_simulate.py +++ b/tests/test_simulate.py @@ -44,7 +44,7 @@ def test_get_y_propagates_reaction_automatically(): y, r = simulate.get_y(simulate.get_dydt(scheme, [1.0, 1.0]), y0=y0) assert y.t_min == 0.0 - assert y.t_max == 1000.0 + assert y.t_max == 300.0 assert y(y.t_min) == pytest.approx(y0) assert y(y.t_max) == pytest.approx( [1.668212890625, 0.6728515625, 0.341787109375], 9e-5 @@ -86,11 +86,11 @@ def test_get_y_conservation_in_equilibria(): t = np.linspace(y.t_min, y.t_max, num=100) assert y.t_min == 0.0 - assert y.t_max == 10.0 + assert y.t_max == 5.0 assert y(y.t_min) == pytest.approx(y0) - assert y(y.t_max) == pytest.approx([0.5, 0.5], 3e-5) + assert y(y.t_max) == pytest.approx([0.5, 0.5], 5e-5) assert r(y.t_min) == pytest.approx([-1, 1]) - assert r(y.t_max) == pytest.approx([0.0, 0.0], abs=3e-5) + assert r(y.t_max) == pytest.approx([0.0, 0.0], abs=5e-5) assert y.t_min == t[0] assert y.t_max == t[-1]