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

# ## A simple `A <-> B` reaction between 2 species
# with 1st-order kinetics in both directions, taken to equilibrium
# 
# See also the experiment _"1D/reactions/reaction_1"_ ; this is the "single-compartment" version of it.
# 
# 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.modules.reactions.reaction_dynamics import ReactionDynamics

import numpy as np
import plotly.express as px
from src.modules.visualization.graphic_log import GraphicLog


# In[3]:


# Initialize the HTML logging (for the graphics)
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")


# # Initialize the System

# In[4]:


# Initialize the reaction
chem_data = chem(names=["A", "B"])

# Reaction A <-> B , with 1st-order kinetics in both directions
chem_data.add_reaction(reactants=["A"], products=["B"], 
                       forward_rate=3., reverse_rate=2.)

print("Number of reactions: ", chem_data.number_of_reactions())


# In[5]:


chem_data.describe_reactions()


# In[6]:


# Send a plot of the network of reactions to the HTML log file
graph_data = chem_data.prepare_graph_network()
GraphicLog.export_plot(graph_data, "vue_cytoscape_1")


# # Start the simulation

# In[7]:


dynamics = ReactionDynamics(chem_data=chem_data)


# In[8]:


# Initial concentrations of all the chemicals, in index order
dynamics.set_conc([10., 50.])


# In[9]:


dynamics.describe_state()


# In[10]:


dynamics.get_history()


# ## Start the reaction

# In[11]:


# First step of reaction
dynamics.single_compartment_react(initial_step=0.1, n_steps=1,
                                  snapshots={"initial_caption": "first reaction step"})


# In[12]:


dynamics.get_history()


# In[13]:


# Numerous more fixed steps
dynamics.single_compartment_react(initial_step=0.1, n_steps=10,
                                  snapshots={"initial_caption": "2nd reaction step",
                                             "final_caption": "last reaction step"})


# #### NOTE: for demonstration purposes, we're using FIXED time steps...  Typically, one would use the option for automated variable time steps (see experiment `react 2`)

# In[14]:


dynamics.get_history()


# ### Check the final equilibrium

# In[15]:


dynamics.get_system_conc()


# NOTE: Consistent with the 3/2 ratio of forward/reverse rates (and the 1:1 stoichiometry, and the 1st order reactions), the systems settles in the following equilibrium:
# 
# [A] = 23.99316406
#  
# [B] = 36.00683594
# 

# In[16]:


# Verify that the reaction has reached equilibrium
dynamics.is_in_equilibrium()


# ### Note that, because of the high initial concentration of B relative to A, the overall reaction has proceeded **in reverse**

# ## Plots of changes of concentration with time

# In[17]:


dynamics.plot_curves(colors=['blue', 'orange'])


# ### Note the raggedness of the left-side (early times) of the curves.  
# ### In experiment `react_2` this simulation gets repeated with an adaptive variable time resolution that takes smaller steps at the beginning, when the reaction is proceeding faster   
# ### By contrast, here we used FIXED steps (see below), which is generally a bad approach

# In[18]:


dynamics.plot_curves(colors=['blue', 'orange'], show_intervals=True)


# In[19]:


df = dynamics.get_history()


# In[20]:


df


# ## Now investigate A_dot, i.e. d[A]/dt

# In[21]:


A = list(df.A)


# In[22]:


A


# In[23]:


len(A)


# In[24]:


A_dot = np.gradient(A, 0.1)      # 0.1 is the constant step size


# In[25]:


A_dot


# In[26]:


df['A_dot'] = A_dot     # Add a column to the dataframe


# In[27]:


df


# In[28]:


fig = px.line(data_frame=dynamics.get_history(), x="SYSTEM TIME", y=["A", "A_dot"], 
              title="Changes in concentration of A with time (blue) , and changes in its rate of change (brown)",
              color_discrete_sequence = ['navy', 'brown'],
              labels={"value":"concentration (blue) /<br> concentration change per unit time (brown)", "variable":"Chemical"})
fig.show()


# ### At t=0, [A]=10 and [A] has a high rate of change (70);  
# ### as the system approaches equilibrium, [A] approaches a value of 24, and its rate of change decays to zero.
# 
# The curves are jagged because of limitations of numerically estimating derivatives, as well as the large time steps taken (especially in the early times, when there's a lot of change.)

# In[ ]: