from scipy import integrate
import numpy as np
from matplotlib import cm
from bokeh.io import output_notebook, show, output_file, push_notebook
from bokeh.layouts import column, row
from bokeh.models import ColumnDataSource, Slider, TextInput, ColorBar, LinearColorMapper, ContinuousColorMapper
from bokeh.plotting import curdoc, figure, show
from bokeh.application import Application
from bokeh.application.handlers import FunctionHandler
from bokeh.palettes import brewer, plasma
from bokeh.transform import transform
from IPython.display import display, Markdown
def diff(x, concs, ks,reacts,matrix):
dt_list = []
for val in (list(matrix.values())):
tot = 0
for key,value in val.items():
react = reacts[key]
tot += value*eval(react)
dt_list.append(tot)
for i in range(len(concs)-len(dt_list)):
dt_list.append(0)
return dt_list
class ODE_systems:
def __init__(self, reactions, matrix,ks,concs,to_run = None):
self.reactions = reactions
self.matrix = matrix
self.ks = ks
self.concs = concs
if to_run != None:
self.to_run = to_run
else:
self.to_run = None
def run_ODEs(self):
to_ins = (list(self.matrix.values())[0])
sol = integrate.solve_ivp(diff, [0,time], self.concs, t_eval = timespan, args=[self.ks,self.reactions,self.matrix],method='BDF', rtol=1e-8, atol=1e-10)
#sol = integrate.solve_ivp(diff, [0,time], self.concs, t_eval = timespan, args=[self.ks,self.reactions,self.matrix],method='BDF', rtol=1e-8, atol=1e-10)
#output_notebook()
return sol
def get_matrices(self):
systems_list = []
if self.to_run != None:
for run in self.to_run:
new_matrix, new_ks, new_reactions, new_concs = preprocess_matrix(self.matrix,self.ks,self.reactions,self.concs,run)
new_reactions = preprocess_reactions(new_reactions,run)
print(f' reaction matrix: {dict(new_matrix)}')
print(f' rate constants: {list(new_ks)}')
print(f' reactions: {dict(new_reactions)}')
print(f' initial conditions: {list(new_concs)}')
print(f' number of constants: {len(new_ks)}')
print('')
systems_list.append([dict(new_matrix), list(new_ks), dict(new_reactions), list(new_concs), len(new_ks)])
else:
print(f' reaction matrix: {self.matrix}')
print(f' rate constants: {self.ks}')
print(f' reactions: {self.reactions}')
print(f' initial conditions: {self.concs}')
print('')
systems_list.append([self.matrix, self.ks, self.reactions, self.concs])
return systems_list
# Constructing kinetic sheme here: please define the reaction constants ks, initial concentrations cons, reaction matrix and the coefficients matrix
ks = [0.03, 30000, 19000, 3300] #define all rate constants in a list
concs = [0,250,0,250,0,0] #define initial conditions
#define the reactions that can happen in a dictionary, write them down as strings, this is more modular
comps = {
0: "e",
1: "TCE",
2: "TCE*-",
3: "O2",
4: "O2*",
5: "TCEO2*"
}
reactions = {
1:"ks[0]",
2:"ks[1]*concs[0]*concs[1]",
3:"ks[2]*concs[0]*concs[3]",
4:"ks[3]*concs[3]*concs[2]"
}
matrix = {
"e": {1:1, 2:-1, 3:-1},#0
"TCE": {2:-1},#1
"TCE*": {2:1,4:-1},#2
"O2": {3:-1,4:-1},#3
"O2*": {3:1},#4
"TCEO2*": {4:1},#5
}
#define for each species which reactions apply, and with which stoichiometric constant, so{"dXdt": {reaction:stoichiometric_constant}, etc.}
#if more reactions apply it would look like {"dEdt": {1:-1, 2:1}}; for structuring of more complex system look at test_6 at the bottom
# Enable the output of Bokeh plots in the notebook
output_notebook()
#Define the timespan, resolution of the plots and a number of steps for ki variation sliders
t_res = 30
#timespan = range(0,200,t_res)
timespan = np.linspace(0, 550, t_res)
time = timespan[-1]
num_steps = 10
# Initial data
x = ODE_systems(reactions,matrix,ks,concs).run_ODEs().t
y = ODE_systems(reactions,matrix,ks,concs).run_ODEs().y
# Define the function to create the Bokeh document
def modify_doc(doc):
# Initial data
#x = ODE_systems(reactions,matrix,ks,concs).run_ODEs().t
#y = ODE_systems(reactions,matrix,ks,concs).run_ODEs().y
sources = []
plot = figure(title=f"Kinetic profiles", x_axis_label='Time, a.u.', y_axis_label='Concentration, a.u.')
for i in range(len(concs)):
source = ColumnDataSource(data=dict(x=x, y=ODE_systems(reactions,matrix,ks,concs).run_ODEs().y[i]))
sources.append(source)
plot.line('x', 'y', source=source, line_width=4, color = brewer['Set1'][len(concs)][i], legend_label = list(matrix)[i])
# Create sliders and text inputs dynamically based on the number of parameters
sliders = []
start_inputs = []
end_inputs = []
for i in range(len(ks)):
slider = Slider(start=0, end=ks[i]*2, value=ks[i], step=(ks[i]*2)/num_steps, title=f"k {i+1}")
sliders.append(slider)
start_input = TextInput(value=str(slider.start), title="Start")
end_input = TextInput(value=str(slider.end), title="End")
start_inputs.append(start_input)
end_inputs.append(end_input)
# Define callback function to update the plot
def update(attr, old, new):
ks_new = [slider.value for slider in sliders]
new_y = ODE_systems(reactions,matrix,ks_new,concs).run_ODEs().y
for i in range(len(concs)):
sources[i].data = dict(x=x, y=new_y[i])
if handle:
push_notebook(handle=handle)
# Define callbacks to update sliders based on text input values
def update_slider_range(attr, old, new, slider, start_input, end_input):
start = float(start_input.value)
end = float(end_input.value)
step = (end - start) / num_steps
slider.start = start
slider.end = end
slider.step = step
slider.value = max(slider.start, min(slider.value, slider.end)) # Ensure current value is within new range
# Attach the callback to each slider and text input
for slider, start_input, end_input in zip(sliders, start_inputs, end_inputs):
slider.on_change('value', update)
start_input.on_change('value', lambda attr, old, new, s=slider, si=start_input, ei=end_input: update_slider_range(attr, old, new, s, si, ei))
end_input.on_change('value', lambda attr, old, new, s=slider, si=start_input, ei=end_input: update_slider_range(attr, old, new, s, si, ei))
# Arrange layout
slider_rows = [row(slider, start_input, end_input) for slider, start_input, end_input in zip(sliders, start_inputs, end_inputs)]
layout = column(plot, *slider_rows)
# Add the layout to the document
doc.add_root(layout)
# Create an application handler with the function
handler = FunctionHandler(modify_doc)
app = Application(handler)
# Show the plot with the Bokeh server application and get the notebook handle
handle = show(app, notebook_url="localhost:8888", notebook_handle=True)
def generate_latex(reactions, matrix, comps):
# Replace concs[i] with component names in reactions
def replace_concs(reaction, comps):
for idx, name in comps.items():
reaction = reaction.replace(f"concs[{idx}]", name)
return reaction
# Generate LaTeX for reactions
latex_string = "## Reactions\n\n"
latex_string += "\\[\n\\begin{align*}\n"
for i, reaction in reactions.items():
reaction = replace_concs(reaction, comps)
latex_string += f"r_{i}: & \\quad {reaction} \\\\\n"
latex_string += "\\end{align*}\n\\]\n\n"
# Generate LaTeX for matrix
latex_string += "## Kinetic equations\n\n"
latex_string += "\\[\n\\begin{align*}\n"
for eq, terms in matrix.items():
eq_lhs = eq.replace('_', '\\_')
eq_rhs = " + ".join(f"{coef}({replace_concs(reactions[term], comps)})" for term, coef in terms.items())
eq_rhs = eq_rhs.replace('+ -', '-')
eq_rhs = eq_rhs.replace('1(', '')
eq_rhs = eq_rhs.replace('(-1', '-')
eq_rhs = eq_rhs.replace(')', '')
latex_string += f"\\frac{{d{eq_lhs}}}{{dt}} & : {eq_rhs} \\\\\n"
latex_string += "\\end{align*}\n\\]\n"
return latex_string
# Generate the LaTeX code
latex_code = generate_latex(reactions, matrix, comps)
# Display the LaTeX code in a new markdown cell
display(Markdown(latex_code))
[ \begin{align*} r_1: & \quad ks[0] \\ r_2: & \quad ks[1]*e*TCE \\ r_3: & \quad ks[2]*e*O2 \\ r_4: & \quad ks[3]*O2*TCE*- \\ \end{align*} ]
[ \begin{align*} \frac{de}{dt} & : ks[0] -ks[1]*e*TCE -ks[2]*e*O2 \\ \frac{dTCE}{dt} & : -ks[1]*e*TCE \\ \frac{dTCE*}{dt} & : ks[1]*e*TCE -ks[3]*O2*TCE*- \\ \frac{dO2}{dt} & : -ks[2]*e*O2 -ks[3]*O2*TCE*- \\ \frac{dO2*}{dt} & : ks[2]*e*O2 \\ \frac{dTCEO2*}{dt} & : ks[3]*O2*TCE*- \\ \end{align*} ]
#import numpy as np
#from matplotlib import cm
#from bokeh.plotting import figure, show
#from bokeh.models import ColorBar, LinearColorMapper, ContinuousColorMapper
#from bokeh.io import output_notebook
#from bokeh.transform import transform
# Display Bokeh plots inline in Jupyter Notebook
#output_notebook()
#Define the timespan, resolution of the plots and a number of steps for ki variation sliders
t_res = 10
#timespan = range(0,200,t_res)
timespan = np.linspace(0, 550, t_res)
time = timespan[-1]
# Generate a matrix (2D array)
concO2 = np.linspace(0, 250, 300)
concTCE = np.linspace(0, 500, 200)
#X, Y = np.meshgrid(concO2, concTCE)
#for i in range (len(concO2))
Z = np.zeros((len(concO2), len(concTCE)))
for i in range (len(concO2)):
for j in range (len(concTCE)):
concs[3]=concO2[i]
concs[1]=concTCE[j]
Z[i,j]=np.max(ODE_systems(reactions,matrix,ks,concs).run_ODEs().y[5])
# Define the colormap
colormap = cm.get_cmap("plasma")
# Normalize the matrix values for the colormap
min_val = Z.min()
max_val = Z.max()
norm = cm.colors.Normalize(vmin=min_val, vmax=max_val)
# Convert normalized values to colors
colors = colormap(norm(Z))
# Convert the colormap to Bokeh-compatible format
colors = (colors[:, :, :3] * 255).astype(np.uint8)
colors_hex = ["#%02x%02x%02x" % tuple(color) for color in colors.reshape(-1, 3)]
# Create the Bokeh plot
p = figure(title="Contour plot with colormap", match_aspect=True, toolbar_location='right')
mapper = LinearColorMapper(palette="Plasma256", low=min_val, high=max_val)
# Add the image data to the plot
p.image(image=[Z], x=0, y=0, dw=250, dh=500, palette="Plasma256")
# Add the color bar
color_bar = ColorBar(color_mapper=mapper, location=(0, 0), label_standoff=12)
p.add_layout(color_bar, 'right')
# Show the plot
show(p)
# Data export section
#file_path = 'C:\EAU\Data\Calculations and modeling\Kinetics\TCE S oxidation'
file_path = 'O2_TCE_contour_ext'
#os.makedirs(os.path.dirname(file_path), exist_ok=True)
np.savetxt(file_path, Z, delimiter=',')