abstochkin.base

  1""" Base class, AbStochKin, for initializing and storing all data for
  2performing stochastic simulations using the Agent-based Kinetics
  3method. A simulation project can be initialized and run as shown in
  4the examples below.
  5
  6Example
  7-------
  8>>> from abstochkin import AbStochKin
  9>>> sim = AbStochKin()
 10>>> sim.add_process_from_str('A -> ', 0.2)  # degradation process
 11>>> sim.simulate(p0={'A': 100},t_max=20,plot_backend=plotly)
 12>>> # All data for the above simulation is stored in `sim.sims[0]`.
 13>>>
 14>>> # Now set up a new simulation without actually running it.
 15>>> sim.simulate(p0={'A': 10},t_max=10,n=50,run=False,plot_backend=plotly)
 16>>> # All data for the new simulation is stored in `sim.sims[1]`.
 17>>> # The simulation can then be manually run using methods
 18>>> # documented in the class `Simulation`.
 19
 20"""
 21#  Copyright (c) 2024-2025, Alex Plakantonakis.
 22#
 23#  This program is free software: you can redistribute it and/or modify
 24#  it under the terms of the GNU General Public License as published by
 25#  the Free Software Foundation, either version 3 of the License, or
 26#  (at your option) any later version.
 27#
 28#  This program is distributed in the hope that it will be useful,
 29#  but WITHOUT ANY WARRANTY; without even the implied warranty of
 30#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 31#  GNU General Public License for more details.
 32#
 33#  You should have received a copy of the GNU General Public License
 34#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 35
 36import re
 37from ast import literal_eval
 38from concurrent.futures import ProcessPoolExecutor
 39from typing import Any, Literal
 40
 41from .het_calcs import get_het_processes
 42from .process import Process, MichaelisMentenProcess, ReversibleProcess, \
 43    RegulatedProcess, RegulatedMichaelisMentenProcess
 44from .simulation import Simulation
 45
 46
 47class AbStochKin:
 48    """ Base class for Agent-based Kinetics (AbStochKin) simulator.
 49
 50    Attributes
 51    ----------
 52    volume : float, default : None, optional
 53        The volume *in liters* of the compartment in which the processes
 54        are taking place.
 55    volume_unit : str, default : 'L', optional
 56        A string of the volume unit. The default value is 'L' for liters.
 57    time_unit : str, default : 'sec', optional
 58        A string of the time unit to be used for describing the kinetics
 59        of the given processes.
 60    processes : list
 61        A list of the processes that the AbStochKin object has.
 62    het_processes : list
 63        A list of the processes where population heterogeneity in one
 64        of the parameters is to be modeled. This list is a subset of
 65        the `processes` attribute.
 66    sims : list
 67        A list of all simulations performed for the given set of processes.
 68        Each member of the list is an object of the `Simulation` class and
 69        contains all data for that simulation.
 70    """
 71
 72    def __init__(self,
 73                 volume: float = None,
 74                 volume_unit: str = 'L',
 75                 time_unit: str = 'sec'):
 76        self.volume = volume
 77        self.volume_unit = volume_unit
 78        self.time_unit = time_unit
 79
 80        self.processes = list()
 81        self.het_processes = list()
 82        self.sims = list()
 83
 84    def add_processes_from_file(self, filename: str):
 85        """ Add a batch of processes from a text file. """
 86        with open(filename) as f:
 87            lines = f.readlines()
 88
 89        for line in lines:
 90            self.extract_process_from_str(line)
 91
 92    def extract_process_from_str(self, process_str: str):
 93        """
 94        Extract a process and all of its specified parameters from a string.
 95
 96        This functions parses a string specifying all values and parameters
 97        needed to define a process. It then creates a Process object
 98        based on the extracted data.
 99        """
100        process_str = process_str.replace(' ', '').replace('"', '\'')
101
102        proc_str = process_str.split(',')[0]
103
104        # Extract string parameters first
105        patt_str_params = r"(\w+)=('[\w\s,]+')"
106        str_params = re.findall(patt_str_params, process_str)
107
108        process_str_remain = re.sub(r"\w+='[\w\s,]+'", '', process_str)
109
110        patt_num_params = r"(\w+)=([\[\(,.\s\d\)\]]+)"
111        num_params = re.findall(patt_num_params, process_str_remain)
112
113        all_params = list()
114        for name, val_str in num_params + str_params:
115            while val_str.endswith(','):
116                val_str = val_str[:-1]
117            all_params.append((name, literal_eval(val_str)))
118
119        self.add_process_from_str(proc_str, **dict(all_params))
120
121    def add_process_from_str(self,
122                             process_str: str,
123                             /,
124                             k: float | int | list[float | int, ...] | tuple[float | int, float | int],
125                             **kwargs):
126        """
127        Add a process by specifying a string: 'reactants -> products'.
128        Additional arguments determine if a specialized process
129        (such as a reversible, regulated, or Michaelis-Menten process)
130        is to be defined.
131        """
132        kwargs.setdefault('volume', self.volume)
133
134        if '<->' in process_str or 'k_rev' in kwargs:  # reversible process
135            self.processes.append(ReversibleProcess.from_string(process_str, k, **kwargs))
136        elif 'catalyst' in kwargs and 'Km' in kwargs and 'regulating_species' not in kwargs:
137            self.processes.append(MichaelisMentenProcess.from_string(process_str, k, **kwargs))
138        elif 'regulating_species' in kwargs and 'alpha' in kwargs and 'K50' in kwargs and 'nH' in kwargs:
139            if 'catalyst' in kwargs and 'Km' in kwargs:  # Regulated MM process
140                self.processes.append(RegulatedMichaelisMentenProcess.from_string(process_str,
141                                                                                  k, **kwargs))
142            else:  # Regulated process
143                self.processes.append(RegulatedProcess.from_string(process_str, k, **kwargs))
144        else:  # simple unidirectional process
145            self.processes.append(Process.from_string(process_str, k, **kwargs))
146
147    def add_process(self,
148                    /,
149                    reactants: dict,
150                    products: dict,
151                    k: float | int | list[float | int, ...] | tuple[float | int, float | int],
152                    **kwargs):
153        """
154        Add a process by using a dictionary for the reactants and products.
155        Additional arguments determine if a specialized process
156        (such as a reversible, regulated, or Michaelis-Menten process)
157        is to be defined.
158        """
159        kwargs.setdefault('volume', self.volume)
160
161        if 'k_rev' in kwargs:  # reversible process
162            self.processes.append(ReversibleProcess(reactants, products, k, **kwargs))
163        elif 'catalyst' in kwargs and 'Km' in kwargs and 'regulating_species' not in kwargs:
164            self.processes.append(MichaelisMentenProcess(reactants, products, k, **kwargs))
165        elif 'regulating_species' in kwargs and 'alpha' in kwargs and 'K50' in kwargs and 'nH' in kwargs:
166            if 'catalyst' in kwargs and 'Km' in kwargs:  # Regulated MM process
167                self.processes.append(
168                    RegulatedMichaelisMentenProcess(reactants, products, k, **kwargs))
169            else:  # Regulated process
170                self.processes.append(RegulatedProcess(reactants, products, k, **kwargs))
171        else:  # simple unidirectional process
172            self.processes.append(Process(reactants, products, k, **kwargs))
173
174    def del_process_from_str(self,
175                             process_str: str,
176                             /,
177                             k: float | int | list[float | int, ...] | tuple[float | int],
178                             **kwargs):
179        """ Delete a process by specifying a string: 'reactants -> products'. """
180        kwargs.setdefault('volume', self.volume)
181
182        try:
183            if '<->' in process_str or 'k_rev' in kwargs:  # reversible process
184                self.processes.remove(ReversibleProcess.from_string(process_str, k, **kwargs))
185            elif 'catalyst' in kwargs and 'Km' in kwargs and 'regulating_species' not in kwargs:
186                self.processes.remove(MichaelisMentenProcess.from_string(process_str, k, **kwargs))
187            elif 'regulating_species' in kwargs and 'alpha' in kwargs and 'K50' in kwargs and 'nH' in kwargs:
188                if 'catalyst' in kwargs and 'Km' in kwargs:  # Regulated MM process
189                    self.processes.remove(
190                        RegulatedMichaelisMentenProcess.from_string(process_str, k, **kwargs))
191                else:  # Regulated process
192                    self.processes.remove(RegulatedProcess.from_string(process_str, k, **kwargs))
193            else:  # simple unidirectional process
194                self.processes.remove(Process.from_string(process_str, k, **kwargs))
195        except ValueError:
196            print("Process to be removed was not found.")
197        else:
198            print(f"Removed: {process_str}, k = {k}, kwargs = {kwargs}")
199
200    def del_process(self,
201                    /,
202                    reactants: dict,
203                    products: dict,
204                    k: float | int | list[float | int, ...] | tuple[float | int, float | int],
205                    **kwargs):
206        """ Delete a process by using a dictionary for the reactants and products. """
207        kwargs.setdefault('volume', self.volume)
208
209        try:
210            if 'k_rev' in kwargs:  # reversible process
211                self.processes.remove(ReversibleProcess(reactants, products, k, **kwargs))
212            elif 'catalyst' in kwargs and 'Km' in kwargs and 'regulating_species' not in kwargs:
213                self.processes.remove(MichaelisMentenProcess(reactants, products, k, **kwargs))
214            elif 'regulating_species' in kwargs and 'alpha' in kwargs and 'K50' in kwargs and 'nH' in kwargs:
215                if 'catalyst' in kwargs and 'Km' in kwargs:  # Regulated MM process
216                    self.processes.remove(
217                        RegulatedMichaelisMentenProcess(reactants, products, k, **kwargs))
218                else:  # Regulated process
219                    self.processes.remove(RegulatedProcess(reactants, products, k, **kwargs))
220            else:  # simple unidirectional process
221                self.processes.remove(Process(reactants, products, k, **kwargs))
222        except ValueError:
223            print("Process to be removed was not found.")
224        else:
225            lhs, rhs = Process(reactants, products, k)._reconstruct_string()
226            print(f"Removed: " + " -> ".join([lhs, rhs]) + f", k = {k}, kwargs = {kwargs}")
227
228    def simulate(self,
229                 /,
230                 p0: dict,
231                 t_max: float | int,
232                 dt: float = 0.01,
233                 n: int = 100,
234                 *,
235                 random_seed: int = 19,
236                 solve_odes: bool = True,
237                 ode_method: str = 'RK45',
238                 run: bool = True,
239                 show_plots: bool = True,
240                 plot_backend: Literal['matplotlib', 'plotly'] = 'matplotlib',
241                 multithreading: bool = True,
242                 max_agents_by_species: dict = None,
243                 max_agents_multiplier: int = 2,
244                 _return_simulation: bool = False):
245        """
246        Start an AbStochKin simulation by creating an instance of the class
247        `Simulation`. The resulting object is appended to the list
248        in the class attribute `AbStochKin.sims`.
249
250        Parameters
251        ----------
252        p0 : dict[str: int]
253            Dictionary specifying the initial population sizes of all
254            species in the given processes.
255        t_max : float or int
256            Numerical value of the end of simulated time in the units
257            specified in the class attribute `AbStochKin.time_unit`.
258        dt : float, default: 0.1, optional
259            The duration of the time interval that the simulation's
260            algorithm considers. The current implementation only
261            supports a fixed time step interval whose value is `dt`.
262        n : int, default: 100, optional
263            The number of repetitions of the simulation to be performed.
264        random_seed : int, default: 19, optional
265            A number used to seed the random number generator.
266        solve_odes : bool, default: True, optional
267            Specify whether to numerically solve the system of
268            ODEs defined from the given set of processes.
269        ode_method : str, default: RK45, optional
270            Available ODE methods: RK45, RK23, DOP853, Radau, BDF, LSODA.
271        run : bool, default: True, optional
272            Specify whether to run an AbStochKin simulation.
273        show_plots : bool, default: True, optional
274            Specify whether to graph the results of the AbStochKin simulation.
275        plot_backend : str, default: 'matplotlib', optional
276            `Matplotlib` and `Plotly` are currently supported.
277        multithreading : bool, default: True, optional
278            Specify whether to parallelize the simulation
279            using multithreading. If `False`, the ensemble
280            of simulations is run sequentially.
281        max_agents_by_species : None or dict, default: dict
282            Specification of the maximum number of agents that each
283            species should have when running the simulation.
284            If `None`, that a default approach will
285            be taken by the class `Simulation` and the number
286            for each species will be automatically determined
287            (see method `Simulation._setup_runtime_data()` for details).
288            The entries in the dictionary should be
289            `species name (string): number (int)`.
290        max_agents_multiplier : float or int, default: 2
291            This parameter is used to calculate the maximum number of
292            agents of each species that the simulation engine allocates
293            memory for. This be determined by multiplying the maximum value
294            of the ODE time trajectory for this species by the
295            multiplier value specified here.
296        _return_simulation : bool
297            Determines if the `self.simulate` method returns a `Simulation`
298            object or appends it to the list `self.sims`.
299            Returning a `Simulation` object is needed when calling the method
300            `simulate_series_in_parallel`.
301        """
302
303        if max_agents_by_species is None:
304            max_agents_by_species = dict()
305
306        self.het_processes = get_het_processes(self.processes)
307
308        sim = Simulation(p0,
309                         t_max,
310                         dt,
311                         n,
312                         self.processes,
313                         random_state=random_seed,
314                         do_solve_ODEs=solve_odes,
315                         ODE_method=ode_method,
316                         do_run=run,
317                         show_graphs=show_plots,
318                         graph_backend=plot_backend,
319                         use_multithreading=multithreading,
320                         max_agents=max_agents_by_species,
321                         max_agents_multiplier=max_agents_multiplier,
322                         time_unit=self.time_unit)
323
324        if _return_simulation:
325            assert run, "Must run individual simulations if a series of " \
326                        "simulations is to be run with multiprocessing."
327
328            # Set un-pickleable objects to None for data serialization to work
329            sim.algo_sequence = None
330            sim.progress_bar = None
331
332            return sim
333        else:
334            self.sims.append(sim)
335
336    def simulate_series_in_parallel(self,
337                                    series_kwargs: list[dict[str, Any], ...],
338                                    *,
339                                    max_workers: int = None):
340        """
341        Perform a series of simulations in parallel by initializing
342        separate processes. Each process runs a simulation and appends
343        a `Simulation` object in the list `self.sims`.
344
345        Parameters
346        ----------
347        series_kwargs : list of dict
348            A list containing dictionaries of the keyword arguments for
349            performing each simulation in the series. The number of elements
350            in the list is the number of simulations that will be run.
351        max_workers : int, default: None
352            The maximum number of processes to be used for performing
353            the given series of simulations. If None, then as many worker
354            processes will be created as the machine has processors.
355
356        Examples
357        --------
358        - Run a series of simulations by varying the initial population size of A.
359        >>>  from abstochkin import AbStochKin
360        >>>
361        >>>  sim = AbStochKin()
362        >>>  sim.add_process_from_str("A -> B", 0.3, catalyst='E', Km=10)
363        >>>  series_kwargs = [{"p0": {'A': i, 'B': 0, 'E': 10}, "t_max": 10} for i in range(40, 51)]
364        >>>  sim.simulate_series_in_parallel(series_kwargs)
365        """
366        extra_opts = {"show_plots": False, "_return_simulation": True}
367        with ProcessPoolExecutor(max_workers=max_workers) as executor:
368            futures = [executor.submit(self.simulate, **(kwargs | extra_opts)) for kwargs in
369                       series_kwargs]
370            for future in futures:
371                self.sims.append(future.result())