-
Notifications
You must be signed in to change notification settings - Fork 0
/
utils.py
105 lines (73 loc) · 2.78 KB
/
utils.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
from abc import abstractmethod
from typing import Tuple
from numpy import array, ma
class StochasticAlgorithm:
@abstractmethod
def step(self) -> float:
pass
@abstractmethod
def get_state(self):
pass
@abstractmethod
def update_molecule(self, molecule: str, qnt: int):
pass
class OdeAlgorithm:
@abstractmethod
def solve(self, precision: float, end_time: float):
pass
class Reaction:
"""
Class that simulates a chemical reaction, each instance contains the reactants,
the products and a kinetic constant.
"""
def __init__(self, name: str, reactants: dict, products: dict, kinetic: float):
self.name = name
self.reactants: dict = reactants
self.products: dict = products
self.kinetic: float = kinetic
def show_reaction(self) -> str:
reaction = ""
for r, l in self.reactants.items():
reaction += str(l) + str(r) + " "
reaction += " --> "
for r, l in self.products.items():
reaction += str(l) + str(r) + " "
return reaction + " (" + str(self.kinetic) + ")"
def solveOde(model_: OdeAlgorithm, end_time: float = 100, precision: float = 0.1):
return model_.solve(precision, end_time)
def simulation(model_: StochasticAlgorithm, events: list, end_time: float):
current_time: float = 0
hist, times = [], []
while current_time < end_time:
# perform one step of computation
dt = model_.step()
if dt == -1:
current_time = end_time
else:
current_time += dt
# updating the current time and history
times.append(current_time)
hist.append(list(model_.get_state()))
# check if the event must occur
# the list of events is sorted by the fist element (time) and when a specific event
# do not occur because is too early then we ignore all the rest because they have higher time values.
for time, molecule, qnt in events:
if time <= current_time:
model_.update_molecule(molecule, qnt)
events.pop(0)
else:
break
return array(hist), times
def avg_simulations(list_simulations: list, list_times: list) -> Tuple:
"""
Takes a list of simulation and return a sort of means
"""
n_molecules = list_simulations[0].shape[1]
lens = [len(i) for i in list_simulations]
simulations = ma.empty((max(lens), len(list_simulations), n_molecules))
times = ma.empty((max(lens), len(list_times)))
simulations.mask, times.mask = True, True
for idx, (s, t) in enumerate(zip(list_simulations, list_times)):
simulations[:len(s), idx] = s
times[:len(t), idx] = t
return simulations.mean(axis=1), times.mean(axis=-1)