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

# ## Association/Dissociation reaction `A + B <-> C`
# #### with 1st-order kinetics for each species, taken to equilibrium.
# #### Exploration of debugging and diagnostics options
# (Adaptive variable time steps are used)
# 
# _See also the experiment "1D/reactions/reaction_4"_  
# 
# LAST REVISED: May 22, 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.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 <-> C , with 1st-order kinetics for each species
chem_data.add_reaction(reactants=["A" , "B"], products=["C"],
                       forward_rate=5., reverse_rate=2.)

chem_data.describe_reactions()


# In[5]:


# 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[6]:


dynamics = ReactionDynamics(reaction_data=chem_data)


# In[7]:


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


# In[8]:


dynamics.describe_state()


# In[9]:


dynamics.get_history()


# ## Run the reaction

# In[10]:


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

# All of these settings are currently close to the default values... but subject to change; set for repeatability
dynamics.set_thresholds(norm="norm_A", low=0.5, high=0.8, abort=1.44)
dynamics.set_thresholds(norm="norm_B", low=0.08, high=0.5, abort=1.5)
dynamics.set_step_factors(upshift=1.5, downshift=0.5, abort=0.5)
dynamics.set_error_step_factor(0.5)

dynamics.single_compartment_react(initial_step=0.004, reaction_duration=0.06,
                                  variable_steps=True, explain_variable_steps=False,
                                  snapshots={"initial_caption": "1st reaction step",
                                             "final_caption": "last reaction step"})


# In[11]:


dynamics.get_history()


# ## Note: "A" (now largely depleted) is the limiting reagent

# In[12]:


dynamics.explain_time_advance()


# ### Notice how the reaction proceeds in smaller steps in the early times, when the concentrations are changing much more rapidly.
# #### The argument argument _variable_steps=True_ dynamically adjusts the initial_step (which is initially found to be too large, leading to some backtracking

# In[13]:


dynamics.plot_step_sizes(show_intervals=True)


# ## Plots of changes of concentration with time

# In[14]:


dynamics.plot_curves(colors=['red', 'violet', 'green'], show_intervals=True)


# ### Check the final equilibrium

# In[15]:


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


# In[ ]:





# # Everthing below is just for diagnostic insight 
# ### into the adaptive variable time steps

# In[16]:


dynamics.get_diagnostic_decisions_data()


# In[17]:


dynamics.get_diagnostic_rxn_data(rxn_index=0)


# In[18]:


dynamics.get_diagnostic_conc_data()


# In[ ]: