Make get_y return a function and use a clever default for t_span

ROADMAP.rst

@@ -1,7 +1,6 @@
 Initial code idea:
 
 y0 = [1.0, 0.0]
-t_span = [0.0, 100.0]
 freeenergies = [0.0, 20.0, -10.0]
 
 model = parse("A -> A‡ -> B")
@@ -13,8 +12,8 @@ model = parse("A -> A‡ -> B")
 k = eyring(model.B @ freeenergies)
 
 dydt = get_dydt(model, k)
-t, y, r = get_y(dydt, t_span, y0)
-y[:, -1]
+y, r = get_y(dydt, y0)
+y(y.t_max)
 
 In some advanced tutorial, I might show how to create a polymerization model
 (Pn + P -> Pn+1‡ -> Pn+1) using data extrapolated from data on small to

docs/examples/hydrogen-abstraction-microkinetics.py

@@ -31,9 +31,9 @@
 y0_HCl = 1 / (np.sum(model.compounds["HCl"].atommasses) * 1e3)
 y0 = [y0_CH4, y0_Cl, 0.0, 0.0, y0_HCl]
 print(y0)
-t_span = [0, 3e-3]
+
 dydt = api.get_dydt(model.scheme, k_eck)
-t, y, r = api.get_y(dydt, y0=y0, t_span=t_span, method="Radau")
+t, y, r = api.get_y(dydt, y0=y0, method="Radau")
 
 print(model.scheme.compounds)
 print(y[:, -1])

docsrc/examples/hydrogen-abstraction-drc.py

@@ -28,13 +28,11 @@
 y0_HCl = 1 / (np.sum(model.compounds["HCl"].atommasses) * 1e3)
 y0 = [y0_CH4, y0_Cl, 0.0, 0.0, y0_HCl]
 print(y0)
-t_span = [0, 3e-3]
 
 t, drc = api.get_drc(
     model.scheme,
     model.compounds,
     y0,
-    t_span,
     scale="M-1 s-1",
     method="BDF",
     num=500,

docsrc/examples/hydrogen-abstraction-microkinetics.py

@@ -31,17 +31,19 @@
 y0_HCl = 1 / (np.sum(model.compounds["HCl"].atommasses) * 1e3)
 y0 = [y0_CH4, y0_Cl, 0.0, 0.0, y0_HCl]
 print(y0)
-t_span = [0, 3e-3]
+
 dydt = api.get_dydt(model.scheme, k_eck)
-t, y, r = api.get_y(dydt, y0=y0, t_span=t_span, method="Radau")
+y, r = api.get_y(dydt, y0=y0, method="Radau")
 
 print(model.scheme.compounds)
-print(y[:, -1])
+print(y(y.t_max))
+
+t = np.linspace(y.t_min, 5e-3)
 
 fig, ax = plt.subplots()
 for i, name in enumerate(model.scheme.compounds):
     if not api.is_transition_state(name):
-        ax.plot(1e3 * t, 1e3 * y[i], label=f"{name}")
+        ax.plot(1e3 * t, 1e3 * y(t)[i], label=f"{name}")
 
 ax.set_ylabel("Concentration [mM]")
 ax.set_xlabel("Time [ms]")

docsrc/notebooks/#3 Correcting for systematic error in computational microkinetics.ipynb

@@ -51,7 +51,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.9"
+   "version": "3.8.5"
   }
  },
  "nbformat": 4,

docsrc/tutorials/1-rate_constant.rst

@@ -29,11 +29,12 @@ Simple concurrent first order reactions using the rates above would be:
 ...     A -> AC‡ -> C
 ... """)
 >>> dydt = simulate.get_dydt(scheme, k)
->>> t, y, r = simulate.get_y(dydt, y0=[1.0, 0.0, 0.0, 0.0, 0.0], method="Radau")
+>>> y, r = simulate.get_y(dydt, y0=[1.0, 0.0, 0.0, 0.0, 0.0], method="Radau")
+>>> t = np.linspace(y.t_min, 5.0)
 >>> plt.clf()
 >>> for i, compound in enumerate(scheme.compounds):
 ...    if not compound.endswith("‡"):
-...        plt.plot(t, y[i], label=compound)
+...        plt.plot(t, y(t)[i], label=compound)
 [...]
 >>> plt.legend()
 <...>

docsrc/tutorials/2-equilibrium-constants.rst

@@ -21,9 +21,9 @@ So let's check it:
 (['Cd2p', 'MeNH2', '[Cd(MeNH2)4]2p'],
  ['Cd2p + 4 MeNH2 -> [Cd(MeNH2)4]2p', '[Cd(MeNH2)4]2p -> Cd2p + 4 MeNH2'])
 >>> dydt = simulate.get_dydt(scheme, [K, 1.])
->>> t, y, r = simulate.get_y(dydt, y0=[0., 0., 1.])
->>> y[:, -1]
+>>> y, r = simulate.get_y(dydt, y0=[0., 0., 1.])
+>>> y(y.t_max)
 array([0.01608807, 0.06435228, 0.98391193])
->>> Kobs = y[:, -1][2] / (y[:, -1][0] *  y[:, -1][1]**4)
+>>> Kobs = y(y.t_max)[2] / (y(y.t_max)[0] *  y(y.t_max)[1]**4)
 >>> np.log10(Kobs)
 6.55

docsrc/tutorials/4-curtin-hammett.rst

@@ -24,7 +24,7 @@ The selectivity is defined as the ratio between both products:
 >>> k
 array([1.445e+13, 6.212e+12, 2.905e+05, 2.310e+04])
 
->>> t, y, r = simulate.get_y(simulate.get_dydt(scheme, k), y0=[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], method="BDF")
->>> S = y[:, -1][scheme.compounds.index("P1")] / y[:, -1][scheme.compounds.index("P2")]
+>>> y, r = simulate.get_y(simulate.get_dydt(scheme, k), y0=[1.0, 0.0, 0.0, 0.0, 0.0, 0.0], method="BDF")
+>>> S = y(y.t_max)[scheme.compounds.index("P1")] / y(y.t_max)[scheme.compounds.index("P2")]
 >>> S
 5.408

docsrc/tutorials/degree-of-rate-control/degree-of-rate-control.rst

@@ -31,13 +31,14 @@ Now the graph:
 
 >>> import matplotlib.pyplot as plt
 >>> y0 = [1., 1000., 0, 1., 0., 0., 0., 0.]
->>> t, y, r = simulate.get_y(dydt, y0, [0.0, 30 * 60.0], method="BDF")
+>>> y, r = simulate.get_y(dydt, y0, [0.0, 30 * 60.0], method="BDF")
+>>> t = np.linspace(y.t_min, 500.0)
 >>> plt.clf()
 >>> for i, compound in enumerate(scheme.compounds):
 ...    if not compound.endswith("‡"):
-...        plt.plot(t, y[i], label=compound)
+...        plt.plot(t, y(t)[i], label=compound)
 [...]
->>> drc = [derivative(lambda k: simulate.get_dydt(scheme, k)(0.0, y)[-1], np.array(k), 1e-4) for y in y.T]
+>>> drc = [derivative(lambda k: simulate.get_dydt(scheme, k)(0.0, y)[-1], np.array(k), 1e-4) for y in y(t).T]
 >>> plt.plot(t, drc, label="DRC")
 [...]
 >>> plt.legend()

docsrc/tutorials/simulation/0-basic-reaction-simulation.rst

@@ -56,17 +56,20 @@ We are going to simulate 10 seconds, starting with an initial concentration of
 1 molar of A (the concentration units evidently depend on the units of the
 reaction rate constant).
 
->>> t, y, r = simulate.get_y(dydt, y0=[1., 0.], method="Radau")
+>>> y, r = simulate.get_y(dydt, y0=[1., 0.], method="Radau")
+
+>>> import numpy as np
+>>> t = np.linspace(y.t_min, 5.0)
 
 The simulation data is stored in `t` (points in time) and `y` (concentrations).
 In order to generate the graph shown at the beginning of the present tutorial,
 we do the following:
 
 >>> import matplotlib.pyplot as plt
 >>> fig, ax = plt.subplots()
->>> ax.plot(t, y[0], label="A")
+>>> ax.plot(t, y(t)[0], label="A")
 [...]
->>> ax.plot(t, y[1], label="B")
+>>> ax.plot(t, y(t)[1], label="B")
 [...]
 >>> ax.legend()
 <...>
@@ -80,5 +83,5 @@ Text(...)
 We can see that the reaction went to full completion by checking the final
 concentrations:
 
->>> y[:, -1]
+>>> y(y.t_max)
 array([0.0000, 1.0000])

docsrc/tutorials/simulation/1-basic-reaction-simulation.rst

@@ -145,16 +145,18 @@ enzyme-substrate complex yet.
 Let's now do a one minute simulation with ``get_y`` (methods Radau or BDF are
 recommended for likely stiff equations such as those):
 
->>> t = 60.0  # seconds
->>> t, y, r = simulate.get_y(dydt, y0, [0.0, t], method="Radau")
+>>> y, r = simulate.get_y(dydt, y0, method="Radau")
+
+>>> import numpy as np
+>>> t = np.linspace(y.t_min, 60.0)  # seconds
 
 We can graph concentrations over time with ``t`` and ``y``:
 
 >>> import matplotlib.pyplot as plt
 >>> plt.clf()
 >>> for i, compound in enumerate(scheme.compounds):
 ...    if not compound.endswith("‡"):
-...        plt.plot(t, y[i], label=compound)
+...        plt.plot(t, y(t)[i], label=compound)
 [...]
 >>> plt.legend()
 <...>
@@ -174,7 +176,7 @@ Text(...)
 The simulation time was enough to convert all substrate into products and
 regenerate the initial enzyme molecules:
 
->>> y[:, -1]
+>>> y(y.t_max)
 array([0.05, 0.00, 0.00, 0.00, 1.00])
 
 Getting rates back
@@ -183,7 +185,7 @@ Getting rates back
 From the solution above, we can get the rate of production formation over time:
 
 >>> import numpy as np
->>> dy = np.array([dydt(t, y) for t, y in zip(t, y.T)]).T
+>>> dy = np.array([dydt(t, y) for t, y in zip(t, y(t).T)]).T
 >>> plt.clf()
 >>> plt.plot(t, dy[scheme.compounds.index("P")], label="P")
 [...]
@@ -201,7 +203,7 @@ Text(...)
 
 Furthermore, we can get the turnover frequency (TOF) as:
 
->>> total_enzyme = y[scheme.compounds.index("E"), :] + y[scheme.compounds.index("ES"), :]
+>>> total_enzyme = y(t)[scheme.compounds.index("E"), :] + y(t)[scheme.compounds.index("ES"), :]
 >>> tof = dy[scheme.compounds.index("P")] / total_enzyme
 >>> plt.clf()
 >>> plt.plot(t, tof, label="P")

overreact/api.py

@@ -262,13 +262,14 @@ def get_freeenergies(
     return enthalpies - temperature * entropies + _np.asanyarray(bias)
 
 
+# TODO(schneiderfelipe): this should probably be deprecated but it is very
+# good as a starting point for sensitivity analyses in general.
 def get_drc(
     scheme,
     compounds,
     y0,
-    t_span=(0.0, 10.0),
+    t_span=None,
     method="Radau",
-    num=50,
     qrrho=True,
     scale="l mol-1 s-1",
     temperature=298.15,
@@ -281,10 +282,9 @@ def get_drc(
     --------
     >>> model = parse_model("data/tanaka1996/UMP2/6-311G(2df,2pd)/model.jk")
     """
-    i = 2  # Gi
     x0 = _np.zeros(len(scheme.compounds))
 
-    def func(x, full_output=False):
+    def func(t, x=0.0, i=-1):
         bias = _np.copy(x0)
         bias[i] += x
 
@@ -301,18 +301,18 @@ def func(x, full_output=False):
             # molecularity=molecularity,
             # volume=volume,
         )
-        t, y, r = get_y(
-            get_dydt(scheme, k), y0=y0, t_span=t_span, method=method, num=num
-        )
+        _, r = get_y(get_dydt(scheme, k), y0=y0, t_span=t_span, method=method)
 
-        if full_output:
-            return t, y, _np.log(r)
-        return _np.log(r)
+        return _np.log(r(t))
 
-    drc = (
-        -constants.R * temperature * _derivative(func, x0=0.0, dx=dx, n=1, order=order)
-    )
-    return func(0.0, full_output=True)[0], drc
+    def drc(t, i=-1):
+        return (
+            -constants.R
+            * temperature
+            * _derivative(lambda x: func(t, x, i), x0=0.0, dx=dx, n=1, order=order)
+        )
+
+    return drc
 
 
 def get_k(