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

# ## 2 COUPLED reactions of different speeds, forming a "cascade":  
# ### `A <-> B` (fast) and `B <-> C` (slow)
# All 1st order. Taken to equilibrium. Both reactions are mostly forward.
# The concentration of the intermediate product B manifests 1 oscillation (transient "overshoot")
# 
# (Adaptive variable time resolution is used, with extensive diagnostics.)
# 
# LAST REVISED: Feb. 5, 2023

# ## Bathtub analogy:
# A is initially full, while B and C are empty.  
# Tubs are progressively lower (reactions are mostly forward.)  
# A BIG pipe connects A and B: fast kinetics.  A small pipe connects B and C: slow kinetics. 
# 
# INTUITION: B, unable to quickly drain into C while at the same time being blasted by a hefty inflow from A,  
# will experience a transient surge, in excess of its final equilibrium level.
# 
# * [Compare with the final reaction plot (the red line is B)](#cascade_1_plot)

# ![2 Coupled Reactions](../../docs/2_coupled_reactions.png)

# In[1]:


# Extend the sys.path variable, to contain the project's root directory
import set_path
set_path.add_ancestor_dir_to_syspath(2)  # The number of levels to go up 
                                         # to reach the project's home, from the folder containing this notebook


# In[2]:


from experiments.get_notebook_info import get_notebook_basename

from src.modules.reactions.reaction_data import ReactionData 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
# Specify the chemicals and the reactions

# In[4]:


# Specify the chemicals
chem_data = chem(names=["A", "B", "C"])

# Reaction A <-> B (fast)
chem_data.add_reaction(reactants=["A"], products=["B"],
                       forward_rate=64., reverse_rate=8.) 

# Reaction B <-> C (slow)
chem_data.add_reaction(reactants=["B"], products=["C"],
                       forward_rate=12., 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(reaction_data=chem_data)


# ### Set the initial concentrations of all the chemicals, in their index order

# In[8]:


dynamics.set_conc([50., 0, 0.], snapshot=True)


# In[9]:


dynamics.describe_state()


# In[10]:


dynamics.get_history()


# ## Run the reaction

# In[11]:


dynamics.set_diagnostics()       # To save diagnostic information about the call to single_compartment_react()
#dynamics.verbose_list = [1, 2, 3]      # Uncomment for detailed run information (meant for debugging the adaptive variable time step)

# The changes of concentrations vary very rapidly early on; 
# so, we'll be using the dynamic_substeps option to increase time resolution,
# as long as the reaction remains "fast" (based on a threshold of % change, as specified by fast_threshold)
dynamics.single_compartment_react(time_step=0.02, reaction_duration=0.4,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  dynamic_substeps=10, rel_fast_threshold=190)


# ### Note: the argument  _dynamic_step=10_  splits the time steps in 10 for any reactions that are "fast-changing" (as determined using the given value of  _fast_threshold_ )

# In[12]:


df = dynamics.get_history()
df


# In[13]:


dynamics.explain_time_advance()


# In[14]:


# Expand the early part
df.loc[:59]


# In[15]:


# Expand the last part
df.loc[60:]


# ### Notice:
# * the reaction proceeds in smaller steps in the earlier times (until t=0.160, in line 80), when the concentrations are changing much more rapidly 
# 
# * between lines 30 and 80 (time 0.060 to 0.160), only rection #1 (B <-> C) is regarded as fast-changing (based on the fast_threshold we specified in the _simulation run_); at earlier times, both reactions were regarded as fast-changing
# 
# * "fast-changing" and "slow-changing" is NOT the same thing as "fast" and "slow" reaction kinetics.  For example, reaction #1, though it has much slower kinetics than reaction #0, involves large concentration changes
# 
# * after step 80, both reactions are regarded as slow-changing, and no more intermediate steps are used

# ### Check the final equilibrium

# In[16]:


# Verify that all the reactions have reached equilibrium
dynamics.is_in_equilibrium(tolerance=4)


# ## <a name="cascade_1_plot"></a>Plots of changes of concentration with time

# In[17]:


dynamics.plot_curves(title="Coupled reactions A <-> B and B <-> C",
                     colors=['blue', 'red', 'green'])


# In[18]:


dynamics.curve_intersection(t_start=0, t_end=0.05, var1="A", var2="B")


# In[19]:


dynamics.curve_intersection(t_start=0, t_end=0.05, var1="A", var2="C")


# In[20]:


dynamics.curve_intersection(t_start=0.05, t_end=0.1, var1="B", var2="C")


# In[ ]:





# # EVERYTHING BELOW IS FOR DIAGNOSTIC INSIGHT
# 
# ### Perform some verification

# In[21]:


# Verify that the stoichiometry is respected at every reaction step/substep (NOTE: it requires earlier activation of saving diagnostic data)
dynamics.stoichiometry_checker_entire_run()


# ### Take a peek at the diagnostic data saved during the earlier reaction simulation

# In[22]:


# This table contains a subset of a typical system history
dynamics.diagnostic_data_baselines.get()


# In[23]:


# Concentration increments due to reaction 0 (A <-> B)
# Note that [C] is not affected
dynamics.get_diagnostic_data(rxn_index=0)


# In[24]:


# Concentration increments due to reaction 1 (B <-> C)
# Note that [A] is not affected
dynamics.get_diagnostic_data(rxn_index=1)


# In[25]:


# Expand the last part of the above table
dynamics.get_diagnostic_data(rxn_index=1).loc[60:]


# In[ ]:





# ### Provide a detailed explanation of all the steps/substeps of the reactions, from the saved diagnostic data

# In[26]:


# dynamics.explain_reactions()     # Uncomment if desired


# In[ ]:





# # Re-run with very small constanst steps

# We'll use constant steps of size 0.0005, which is 1/4 of the smallest steps (the "substep" size) previously used in the variable-step run

# In[27]:


dynamics2 = ReactionDynamics(reaction_data=chem_data)


# In[28]:


dynamics2.set_conc([50., 0, 0.], snapshot=True)


# In[29]:


dynamics2.single_compartment_react(time_step=0.0005, reaction_duration=0.4,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  )      


# In[30]:


fig = px.line(data_frame=dynamics2.get_history(), x="SYSTEM TIME", y=["A", "B", "C"], 
              title="Coupled reactions A <-> B and B <-> C (run with constant steps)",
              color_discrete_sequence = ['blue', 'red', 'green'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()


# In[31]:


dynamics2.curve_intersection(t_start=0, t_end=0.05, var1="A", var2="B")


# In[32]:


dynamics2.curve_intersection(t_start=0, t_end=0.05, var1="A", var2="C")


# In[33]:


dynamics2.curve_intersection(t_start=0.05, t_end=0.1, var1="B", var2="C")


# In[34]:


df2 = dynamics2.get_history()
df2


# ## Notice that we now did 801 steps - vs. the 93 of the earlier variable-resolution run!

# ### Let's compare some entries with the coarser previous variable-time run

# **At time t=0.002:**

# In[35]:


df.loc[1]


# In[36]:


df2.loc[4]   # High-precision result from fine fixed-resolution run


# In[37]:


old = df.loc[1][['SYSTEM TIME', 'A', 'B', 'C']].to_numpy().astype(float)
old


# In[38]:


new = df2.loc[4][['SYSTEM TIME', 'A', 'B', 'C']].to_numpy().astype(float)
new


# In[39]:


new - old


# In[ ]:





# **At time t=0.032**, when [A] and [C] are almost equal:

# In[40]:


df.loc[16]


# In[41]:


df2.loc[64]   # High-precision result from fine fixed-resolution run


# In[ ]:





# **At time t=0.14**:

# In[42]:


df.loc[70]


# In[43]:


df2.loc[280]   # High-precision result from fine fixed-resolution run


# In[ ]:





# **At time t=0.26**:

# In[44]:


df.loc[85]


# In[45]:


df2.loc[520]   # High-precision result from fine fixed-resolution run


# In[ ]:





# ## Let's compare the plots of [B] from the earlier (variable-step) run, and the latest (high-precision fixed-step) one

# In[46]:


# Earlier run
fig1 = px.line(data_frame=dynamics.get_history(), x="SYSTEM TIME", y=["B"], 
              color_discrete_sequence = ['red'],
              labels={"value":"concentration", "variable":"Variable-step run"})
fig1.show()


# In[47]:


# Latest run (High-precision result from fine fixed-resolution run)
fig2 = px.line(data_frame=dynamics2.get_history(), x="SYSTEM TIME", y=["B"], 
              color_discrete_sequence = ['orange'],
              labels={"value":"concentration", "variable":"Fine fixed-step run"})
fig2.show()


# In[48]:


import plotly.graph_objects as go

all_fig = go.Figure(data=fig1.data + fig2.data)    # Note that the + is concatenating lists
all_fig.update_layout(title="The 2 runs contrasted")
all_fig['data'][0]['name']="B (variable steps)"
all_fig['data'][1]['name']="B (fixed high precision)"
all_fig.show()


# In[ ]: