#!/usr/bin/env python
# coding: utf-8

# ## 2 COUPLED reactions:  `A + B <-> C`  and  `C + D <-> E` , 
# ### with 1st-order kinetics for each species, taken to equilibrium
# 
# Both reactions are stronger in their respective forward rates.  For the most part, "C" is produced by the 1st reaction, and consumed by the 2nd one
# 
# Diffusion not applicable (just 1 bin)
# 
# LAST REVISED: July 14, 2023

# In[1]:


import set_path      # Importing this module will add the project's home directory to sys.path


# In[2]:


from experiments.get_notebook_info import get_notebook_basename

from src.modules.chemicals.chem_data import ChemData as chem
from src.life_1D.bio_sim_1d import BioSim1D

import plotly.express as px
from src.modules.html_log.html_log import HtmlLog as log
from src.modules.visualization.graphic_log import GraphicLog


# In[3]:


# Initialize the HTML logging
log_file = get_notebook_basename() + ".log.htm"    # Use the notebook base filename for the log file

# Set up the use of some specified graphic (Vue) components
GraphicLog.config(filename=log_file,
                  components=["vue_cytoscape_1"],
                  extra_js="https://cdnjs.cloudflare.com/ajax/libs/cytoscape/3.21.2/cytoscape.umd.js")


# In[4]:


# Initialize the system
chem_data = chem(names=["A", "B", "C", "D", "E"])     # NOTE: Diffusion not applicable (just 1 bin)

# Specify the reactions


# Reactions A + B <-> C  and  C + D <-> E , with 1st-order kinetics for each species
chem_data.add_reaction(reactants=["A", "B"], products=["C"], forward_rate=5., reverse_rate=2.)
chem_data.add_reaction(reactants=["C", "D"], products=["E"], forward_rate=8., reverse_rate=4.)

bio = BioSim1D(n_bins=1, chem_data=chem_data)

bio.set_all_uniform_concentrations( [3., 5., 1., 0.4, 0.1] )

bio.describe_state()


# In[5]:


# Save the state of the concentrations of all species at bin 0
bio.add_snapshot(bio.bin_snapshot(bin_address = 0), caption="Initial state")
bio.get_history()


# In[6]:


chem_data.describe_reactions()


# In[7]:


# Send a header and a plot to the HTML log file
log.write("2 COUPLED reactions:  A + B <-> C  and  C + D <-> E",
          style=log.h2)
# Send the plot to the HTML log file
graph_data = chem_data.prepare_graph_network()
GraphicLog.export_plot(graph_data, "vue_cytoscape_1")


# ### <a name="sec_2_first_step"></a>First step

# In[8]:


# First step
bio.react(time_step=0.01, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()


# 1 bins and 5 species:
#  [[2.27 ]
#  [4.27 ]
#  [1.702]
#  [0.372]
#  [0.128]]

# In[9]:


bio.get_history()


# In[10]:


# Identical 2nd step
bio.react(time_step=0.01, n_steps=1, snapshots={"sample_bin": 0})
bio.describe_state()


# 1 bins and 5 species:
#  [[1.819395  ]
#  [3.819395  ]
#  [2.10707348]
#  [0.32646848]
#  [0.17353152]]

# In[11]:


bio.get_history()


# ### <a name="sec_2"></a>Numerous more steps

# In[12]:


# Numerous more identical steps, to equilibrium
bio.react(time_step=0.01, n_steps=200, snapshots={"sample_bin": 0, "frequency": 10})
bio.describe_state()


# 1 bins and 5 species:
#  [[0.50508029]
#  [2.50508029]
#  [3.16316668]
#  [0.06824696]
#  [0.43175304]]

# ### <a name="sec_2_equilibrium"></a>Equilibrium

# In[13]:


# Verify that each reaction has reached equilibrium
bio.reaction_dynamics.is_in_equilibrium(rxn_index=0, conc=bio.bin_snapshot(bin_address = 0))


# In[14]:


bio.reaction_dynamics.is_in_equilibrium(rxn_index=1, conc=bio.bin_snapshot(bin_address = 0))


# In[15]:


# Do a consistent check with the equilibrium concentrations:

A_eq = bio.bin_concentration(0, 0)
B_eq = bio.bin_concentration(0, 1)
C_eq = bio.bin_concentration(0, 2)
D_eq = bio.bin_concentration(0, 3)
E_eq = bio.bin_concentration(0, 4)

Rf0 = chem_data.get_forward_rate(0)
Rb0 = chem_data.get_reverse_rate(0)

Rf1 = chem_data.get_forward_rate(1)
Rb1 = chem_data.get_reverse_rate(1)

equil = -(Rf0 * A_eq * B_eq - Rf1 * C_eq * D_eq) + (Rb0 * C_eq - Rb1 * E_eq)

print("\nAt equilibrium: ", equil, " (this should be close to 0 at equilibrium)")


# In[16]:


bio.get_history()


# # Plots of changes of concentration with time

# In[17]:


fig = px.line(data_frame=bio.get_history(), x="SYSTEM TIME", y=["A", "B", "C", "D", "E"], 
              title="2 COUPLED reactions:  A + B <-> C  and  C + D <-> E . Changes in concentrations",
              color_discrete_sequence = ['navy', 'cyan', 'red', 'orange', 'green'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()


# A and B get consumed.  
# C gets produced by the 1st reaction more quickly than consumed by the 2nd one.
# D gets consumed, while E gets produced.

# In[ ]: