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

# ### Adaptive time substeps (variable time resolution) for reaction `A <-> B`,
# with 1st-order kinetics in both directions, taken to equilibrium
# 
# Same as the experiment _"react_1"_ , but with adaptive time substeps
# 
# LAST REVISED: Feb. 5, 2023

# 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"])

# 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(reaction_data=chem_data)


# In[8]:


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

dynamics.describe_state()


# In[9]:


dynamics.history.get()


# ## Run the reaction

# In[10]:


dynamics.set_diagnostics()       # To save diagnostic information about the call to single_compartment_react()

dynamics.single_compartment_react(time_step=0.1, n_steps=11,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"},
                                  dynamic_substeps=2, rel_fast_threshold=100)


# ## The argument _dynamic_step=2_ splits the time steps in 2 whenever the reaction is "fast" (as determined using the specified value of _fast_threshold_ )

# In[11]:


df = dynamics.get_history()   # The system's history, saved during the run of single_compartment_react()
df


# In[12]:


dynamics.explain_time_advance()


# ## Notice how the reaction proceeds in smaller steps (half steps) in the early times, when [A] and [B] are changing much more rapidly
# ### That resulted from passing the argument _dynamic_steps=2_ to single_compartment_react()
# 
# * For example, upon completing the half step to t=0.30, i.e. **going from 0.25 to 0.30**, the last change in [A] was (21.508301 - 20.677734) = 0.830567  
# The threshold we specified for a reaction to be considered fast is 100% per full step, for any of the involved chemicals.  For a half step, that corresponds to 50%, i.e. 0.5   
# Since abs(0.830567) > 0.5 , the reaction is therefore marked "FAST" (as it has been so far), and the simulation then proceeds in a half step, to t=0.35

# * (Note: at t=0, in the absence of any simulation data, ALL reactions are _assumed_ to be fast)

# * By contrast, upon completing the half step to t=0.40, i.e. **going from 0.35 to 0.40**, the following changes occur in [A] and [B]:  

# In[13]:


df.iloc[7]


# In[14]:


df.iloc[8]


# In[15]:


s_0_35 = df.iloc[7][['A', 'B']].to_numpy()
s_0_35     # Concentrations of A and B at t=0.35


# In[16]:


s_0_40 = df.iloc[8][['A', 'B']].to_numpy()
s_0_40     # Concentrations of A and B at t=0.40


# In[17]:


(s_0_40 - s_0_35)


# Again, the threshold we specified for a reaction to be considered fast is 100% per full step, for any of the involved chemicals.  
# For a half step, that corresponds to 50%, i.e. 0.5    
# BOTH A's change of abs(0.46) AND B's change of abs(-0.46) are SMALLER than that.   
# The reaction is therefore marked "SLOW", and the simulation then proceeds in a _full time step_ of 0.1 (i.e. a more relaxed time resolution), to t=0.50

# ### Check the final equilibrium

# In[18]:


dynamics.get_system_conc()


# NOTE: Consistent with the 3/2 ratio of forward/reverse rates (and the 1st order reactions),
#  the systems settles in the following equilibrium:
# 
# [A] = 23.99316406
#  
# [B] = 36.00683594
# 
# The presence of equilibrium may be automatically checked withe the following handy function:
# 

# In[19]:


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


# ## Plots of changes of concentration with time

# In[20]:


fig = px.line(data_frame=dynamics.get_history(), x="SYSTEM TIME", y=["A", "B"], 
              title="Reaction A <-> B .  Changes in concentrations with time",
              color_discrete_sequence = ['navy', 'darkorange'],
              labels={"value":"concentration", "variable":"Chemical"})
fig.show()


# ## Note how the left-hand side of this plot is much smoother than it was in experiment "react_1", where no adaptive time substeps were used!

# In[ ]:





# # Diagnostics of the run may be investigated as follows:

# In[21]:


dynamics.diagnostic_data_baselines.get()


# In[22]:


dynamics.get_diagnostic_data(rxn_index=0)      # For the 0-th reaction (the only reaction in our case)


# In[ ]: