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

# ## Exploration of errors in the implementation of substeps, for the simulation of the coupled reactions `2 S <-> U` and `S <-> X`   
# Both mostly forward.  1st-order kinetics throughout.   
# 
# Same as `variable_steps_1`, but with substeps.
# 
# LAST REVISED: Mar. 6, 2023    **THIS IS AN ARCHIVED EXPERIMENT**

# # IMPORTANT: DO NOT ATTEMPT TO RUN THIS NOTEBOOK!   
# ## This is a **"frozen run"** that depends on an old version of Life123, for demonstration purposes  
# (newer versions don't contain this implementation of substeps.)  
# If you bypass the execution exit in the first cell, and run the other cells, you WILL NOT REPLICATE the results below!

# In[1]:


# To stop the current and subsequent cells: USED TO PREVENT ACCIDENTAL RUNS OF THIS NOTEBOOK!

class StopExecution(Exception):
    def _render_traceback_(self):
        return []

raise StopExecution     # See: https://stackoverflow.com/a/56953105/5478830


# In[ ]:





# 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
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 system
chem_data = chem(names=["U", "X", "S"])

# Reaction 2 S <-> U , with 1st-order kinetics for all species (mostly forward)
chem_data.add_reaction(reactants=[(2, "S")], products="U",
                       forward_rate=8., reverse_rate=2.)

# Reaction S <-> X , with 1st-order kinetics for all species (mostly forward)
chem_data.add_reaction(reactants="S", products="X",
                       forward_rate=6., reverse_rate=3.)

chem_data.describe_reactions()

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


# ### Set the initial concentrations of all the chemicals

# In[5]:


dynamics = ReactionDynamics(reaction_data=chem_data)
dynamics.set_conc(conc={"U": 50., "X": 100., "S": 0.})
dynamics.describe_state()


# In[6]:


dynamics.set_diagnostics()       # To save diagnostic information about the call to single_compartment_react()
dynamics.verbose_list = ["substeps", "variable_steps"]

dynamics.single_compartment_react(time_step=0.01, stop_time=0.3, 
                                  variable_steps=True, thresholds={"low": 0.25, "high": 0.64},
                                  dynamic_substeps=2, abs_fast_threshold=100.)

df = dynamics.get_history()
df


# In[7]:


(transition_times, step_sizes) = dynamics.explain_time_advance(return_times=True)


# In[8]:


np.array(step_sizes)


# In[9]:


np.array(transition_times)    # Note: there will be one more transition time (the end time) than step sizes


# ## Plots of changes of concentration with time

# In[10]:


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


# # Notice the irregular "bumpiness" on the plot of S.  See comments at very end.

# In[11]:


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


# In[12]:


# Show the "critical values", i.e. times when the step size changes
dynamics.plot_curves(colors=['green', 'orange', 'blue'], vertical_lines=transition_times, 
                     title="Critical values of time-step changes for reactions `2 S <-> U` and `S <-> X`")


# ## Note: the dashed lines in the plots immediatly above and below are NOT the steps; they are the "critical values", i.e. times when the step size changes.   
# The time steps were shown in an earlier plots

# In[13]:


dynamics.plot_step_sizes(show_intervals=True)


# In[14]:


dynamics.get_diagnostic_rxn_data(rxn_index=0)


# In[15]:


dynamics.get_diagnostic_rxn_data(rxn_index=1)


# In[16]:


dynamics.get_diagnostic_conc_data()


# In[17]:


dynamics.get_diagnostic_decisions_data()


# #### Notice how the first step got aborted, and re-run, because of the large adjusted L2 value in the concentrations 

# # At time t=0.185 :   

# ## Data from this simulation (with substeps):   
# [it took 39 computation steps]

# In[19]:


df.iloc[39]


# In[20]:


substep_data = df.iloc[39][['U', 'X', 'S']].to_numpy()
substep_data


# ## Data from experiment `variable_steps_1` (no substeps, and same main steps as in here):   
# [it took 21 computation steps]

# In[21]:


no_substep_data = np.array([54.631805, 71.023480, 19.712910])
no_substep_data


# ## Data from experiment `substeps_2` (with tiny fixed steps, used as a proxy for exact solution):

# In[22]:


exact_data = np.array([54.529071, 71.351769, 19.590088])
exact_data


# ## Comparing the results:

# In[23]:


substep_error = abs(substep_data - exact_data)
substep_error


# In[24]:


no_substep_error = abs(no_substep_data - exact_data)
no_substep_error


# ## Conclusions:  
# In this scenario, the substeps option made the errors WORSE - in spite of a substantial number of extra steps.  
# Also, substeps led to irregularies in the plot of the concentrations of S.   
# At any rates, substeps are far less relevant ever since the introduction of variable (main) step.

# # => This implementation of substeps was deprecated - and no longer part of current releases

# In[ ]: