abstochkin.utils

  1"""
  2Some utility functions on generating random number streams, measuring a
  3function's runtime, calculating a statistic measuring the goodness of fit
  4when comparing time series data, and performing unit conversion of
  5kinetic parameters.
  6"""
  7
  8#  Copyright (c) 2024-2025, Alex Plakantonakis.
  9#
 10#  This program is free software: you can redistribute it and/or modify
 11#  it under the terms of the GNU General Public License as published by
 12#  the Free Software Foundation, either version 3 of the License, or
 13#  (at your option) any later version.
 14#
 15#  This program is distributed in the hope that it will be useful,
 16#  but WITHOUT ANY WARRANTY; without even the implied warranty of
 17#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 18#  GNU General Public License for more details.
 19#
 20#  You should have received a copy of the GNU General Public License
 21#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 22
 23import functools
 24from time import perf_counter
 25
 26import numpy as np
 27from numpy import random
 28from scipy.constants import N_A
 29
 30
 31def rng_streams(n: int, random_state: int):
 32    """
 33    Generate independent streams of random numbers spawned from the
 34    same initial seed.
 35
 36    Parameters
 37    ----------
 38    n : int
 39        number of generators/streams to generate.
 40    random_state : int, optional
 41        initial seed from which n new seeds are spawned.
 42
 43    Returns
 44    -------
 45    list[PCG64DXSM Generator objects]
 46        List of `n` generators of independent streams of random numbers
 47
 48    Notes
 49    -----
 50    See https://numpy.org/doc/stable/reference/random/parallel.html for more info.
 51
 52    On PCG64DXSM:
 53        - https://numpy.org/doc/stable/reference/random/upgrading-pcg64.html#upgrading-pcg64
 54        - https://numpy.org/doc/stable/reference/random/bit_generators/pcg64dxsm.html
 55
 56    Examples
 57    --------
 58    >>> rng = rng_streams(5)  # make 5 random number generators
 59    >>> a = rng[1].integers(1, 10, 100)
 60    """
 61    ss = random.SeedSequence(random_state)
 62    seeds = ss.spawn(n)
 63
 64    # return [random.default_rng(seed) for seed in seeds]  # PCG64 generator objects
 65
 66    # Return PCG64DXSM generator objects
 67    return [random.Generator(random.PCG64DXSM(seed)) for seed in seeds]
 68
 69
 70def measure_runtime(fcn):
 71    """ Decorator for measuring the duration of a function's execution. """
 72    @functools.wraps(fcn)
 73    def inner(*args, **kwargs):
 74        start = perf_counter()
 75        fcn(*args, **kwargs)
 76        stop = perf_counter()
 77        duration = stop - start
 78        if duration < 60:
 79            msg = f"Simulation Runtime: {duration:.3f} sec"
 80        elif 60 <= duration < 3600:
 81            msg = f"Simulation Runtime: {duration / 60:.3f} min"
 82        else:
 83            msg = f"Simulation Runtime: {duration / 3600:.3f} hr"
 84        print(msg)
 85
 86    return inner
 87
 88
 89def r_squared(actual: np.array, theoretical: np.array) -> float:
 90    """
 91    Compute the coefficient of determination, $R^2$.
 92
 93    In the case of comparing the average AbStochKin-simulated species
 94    trajectory to its deterministic trajectory. Since the latter is only
 95    meaningful for a homogeneous population, R² should be
 96    close to `1` for a simulated homogeneous process.
 97    For a heterogeneous process, it can be interpreted as how close
 98    the simulated trajectory is to the deterministic trajectory of a
 99    *homogeneous* process. In this case, $R^2$ would not be expected
100    to be close to $1$ and the importance of looking at this metric
101    is questionable.
102
103    Parameters
104    ----------
105    actual : numpy.array
106        Actual data obtained through a simulation.
107    theoretical : numpy.array
108        Theoretical data to compare the actual data to.
109
110    Returns
111    -------
112    float
113        The coefficient of determination, $R^2$.
114    """
115    # sst: total sum of squares for simulation avg trajectory
116    sst = np.nansum((actual - np.nanmean(actual)) ** 2)
117    # ssr: sum of square residuals (data vs deterministic prediction)
118    ssr = np.nansum((actual - theoretical) ** 2)
119
120    return 1 - ssr / sst if sst != 0 else np.nan
121
122
123def macro_to_micro(macro_val: float | int | list[float | int, ...] | tuple[float | int, float | int],
124                   volume: float | int,
125                   order: int = 0,
126                   *,
127                   inverse: bool = False) -> float | list[float, ...] | tuple[float, float]:
128    """
129    Convert a kinetic parameter value from macroscopic to microscopic form.
130
131    The ABK algorithm uses microscopic kinetic constants, thus necessitating
132    the conversion of any molar quantities to their microscopic counterpart.
133    For a kinetic parameter, the microscopic form is interpreted as the number
134    of transition events per second (or whatever the time unit may be).
135    For a molar quantity, its microscopic form is the number of particles in
136    the given volume.
137
138    Parameters
139    ----------
140    macro_val : float or int
141        The value of the parameter to be converted, expressed in terms of
142        molar quantities.
143    volume : float or int
144        The volume, in liters, in which the process that the given parameter
145        value is a descriptor of.
146    order : int, default: 0
147        The order of the process whose kinetic parameter is to be converted.
148        The default value of 0 is for parameters (such as Km or K50) whose
149        units are molarity.
150    inverse : bool, default: False
151        Perform the inverse of this operation. That is, convert
152        from microscopic to macroscopic form.
153
154    Returns
155    --------
156    float
157        A kinetic parameter is returned in units of reciprocal seconds.
158        A molar quantity is returned as the number of particles in the
159        given volume.
160
161    Notes
162    -----
163    - A kinetic parameter for a 1st order process will remain unchanged
164    because its units are already reciprocal seconds.
165
166    Reference
167    ---------
168    Plakantonakis, Alex. “Agent-based Kinetics: A Nonspatial Stochastic Method
169    for Simulating the Dynamics of Heterogeneous Populations.”
170    OSF Preprints, 26 July 2019. Web. Section 2.1.
171    """
172    assert volume > 0, "The volume has to be a positive quantity."
173    assert order >= 0, "The process order cannot be negative."
174
175    denom = (N_A * volume) ** (order - 1) if not inverse else 1 / (N_A * volume) ** (order - 1)
176
177    if isinstance(macro_val, (list, tuple)):
178        assert all([True if val >= 0 else False for val in macro_val]), \
179            "The parameter values cannot be negative."
180        micro_vals = [val / denom for val in macro_val]
181        return micro_vals if isinstance(macro_val, list) else tuple(micro_vals)
182    else:
183        assert macro_val >= 0, "The parameter value cannot be negative."
184        return macro_val / denom