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

# ## A cycle of reactions between `A`, `B` and `C` :  (1) `A <-> B <-> C `
# #### the "closing" of the above cycle (the "return" path from `C` to `A`) is coupled with an "energy donor" reaction:
# ## (2) `C + E_High <-> A + E_Low`
# #### where `E_High` and `E_Low` are, respectively, the high- and low- energy molecules that drive the cycle (for example, think of ATP/ADP).   
# Comparisons are made between results obtained with 3 different time resolutions.
# 
# All 1st-order kinetics.    
# 
# LAST REVISED: June 10, 2024 (using v. 1.0 beta33)

# 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
from src.modules.reactions.uniform_compartment import UniformCompartment
from src.modules.numerical.numerical import Numerical as num

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_2"],
                  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 = ChemData(names=["A", "B", "C", "E_high", "E_low"])

# Reaction A <-> B, mostly in forward direction (favored energetically)
# Note: all reactions in this experiment have 1st-order kinetics for all species
chem_data.add_reaction(reactants="A", products="B",
                       forward_rate=9., reverse_rate=3.)

# Reaction B <-> C, also favored energetically
chem_data.add_reaction(reactants="B", products="C",
                       forward_rate=8., reverse_rate=4.)

# Reaction C + E_High <-> A + E_Low, also favored energetically, but kinetically slow
# Note that, thanks to the energy donation from E, we can go "upstream" from C, to the higher-energy level of "A"
chem_data.add_reaction(reactants=["C" , "E_high"], products=["A", "E_low"],
                       forward_rate=1., reverse_rate=0.2)

chem_data.describe_reactions()

# Send the plot of the reaction network to the HTML log file
chem_data.plot_reaction_network("vue_cytoscape_2")


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

# In[5]:


initial_conc = {"A": 100., "E_high": 1000.}  # Note the abundant energy source "E_high"; anything not specified will default to zero
initial_conc


# In[ ]:





# In[ ]:





# # Run # 1 : FIXED time resolution, with COARSE time steps - broken up in 3 time intervals    
# (trial and error, not shown, reveals that increasing any of the time steps below, leads to "excessive time step" errors;  
# note: fixed time resolution is generelly NOT recommended, except for double-checks and error analysis)

# In[6]:


dynamics = UniformCompartment(chem_data=chem_data)
dynamics.set_conc(conc=initial_conc, snapshot=True)
dynamics.describe_state()


# In[7]:


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


# ### We'll split the simulation in three segments (apparent from the graphs of the runs, further down):  
# Time [0-0.03] fast changes   
# Time [0.03-5.] medium changes   
# Time [5.-8.] slow changes, as we approach equilibrium  

# In[8]:


dynamics.single_compartment_react(initial_step=0.0008, target_end_time=0.03, variable_steps=False)
#dynamics.get_history()


# In[9]:


dynamics.single_compartment_react(initial_step=0.001, target_end_time=5., variable_steps=False)
#dynamics.get_history()


# In[10]:


dynamics.single_compartment_react(initial_step=0.005, target_end_time=8., variable_steps=False)


# In[11]:


dynamics.get_history()


# ### Notice we created 5,609 data points, in 3 hand-picked segments at different time resolutions

# In[12]:


dynamics.explain_time_advance()


# ## Plots of changes of concentration with time

# In[13]:


dynamics.plot_history(chemicals=["E_high", "E_low"], colors=["red", "grey"])


# In[14]:


dynamics.plot_history(chemicals=["A", "B", "C"])


# ### The plots has 4 distinctive intersections; locate them and save them for later comparisons across repeated runs:

# In[15]:


run1 = []


# In[16]:


run1.append(dynamics.curve_intersect("E_high", "E_low", t_start=1., t_end=2.))


# In[17]:


run1.append(dynamics.curve_intersect("A", "B", t_start=2.31, t_end=2.33))


# In[18]:


run1.append(dynamics.curve_intersect("A", "C", t_start=3., t_end=4.))


# In[19]:


run1.append(dynamics.curve_intersect("B", "C", t_start=3., t_end=4.))


# In[20]:


run1


# ### Verify the final equilibrium state:

# In[21]:


dynamics.is_in_equilibrium()


# In[ ]:





# In[ ]:





# In[ ]:





# # _NOTE: Everything below, to the end, is JUST A REPEAT of the same experiment, with different time steps, for accuracy comparisons_

# In[ ]:





# # Run # 2. VARIABLE time resolution

# In[22]:


dynamics = UniformCompartment(chem_data=chem_data, preset="mid_inclusive")   # Note: OVER-WRITING the "dynamics" object  # heavy_brakes
dynamics.set_conc(conc=initial_conc, snapshot=True) 
dynamics.describe_state()


# In[23]:


# These will affect the dynamic automated choices of time steps, 
# and were specified by the preset used in instantiating the "dynamics" object
dynamics.show_adaptive_parameters()  


# In[24]:


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


# In[25]:


dynamics.single_compartment_react(initial_step=0.0001, target_end_time=8.0,
                                  variable_steps=True, explain_variable_steps=None)


# ### Notice we created 1,911 data points, a good deal less than in run #1

# In[26]:


dynamics.get_history()


# In[27]:


dynamics.plot_history(chemicals=["E_high", "E_low"], colors=["red", "grey"])


# In[28]:


dynamics.plot_history(chemicals=["A", "B", "C"])


# In[29]:


# Show the timestepe taken (vertical dashed lines) in a small section of the plot
dynamics.plot_history(chemicals=["A", "B", "C"], show_intervals=True, xrange=[2.5, 3.5])


# In[30]:


dynamics.explain_time_advance()  # Notice a lot of timestep adjustments, as needed


# ### Just like we saw in the earlier run, the two plots have 4 distinctive intersections among them; locate and save them:

# In[31]:


run2 = []


# In[32]:


run2.append(dynamics.curve_intersect(t_start=1., t_end=2., chem1="E_high", chem2="E_low"))  # This can be seen in the 1st plot


# In[33]:


run2.append(dynamics.curve_intersect(t_start=2.31, t_end=2.33, chem1="A", chem2="B"))       # This, and later, can be seen in the 2nd plot


# In[34]:


run2.append(dynamics.curve_intersect(t_start=3., t_end=4., chem1="A", chem2="C"))


# In[35]:


run2.append(dynamics.curve_intersect(t_start=3., t_end=4., chem1="B", chem2="C"))


# In[36]:


run2


# In[ ]:





# In[ ]:





# # Run # 3. FIXED time resolution, using FINER fixed resolution than in run #1 

# In[37]:


dynamics = UniformCompartment(chem_data=chem_data)   # Note: OVER-WRITING the "dynamics" object
dynamics.set_conc(conc=initial_conc, snapshot=True) 
dynamics.describe_state()


# In[38]:


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


# In[39]:


dynamics.single_compartment_react(initial_step=0.0004, target_end_time=0.03, variable_steps=False)


# In[40]:


dynamics.single_compartment_react(initial_step=0.0005, target_end_time=5., variable_steps=False)


# In[41]:


dynamics.single_compartment_react(initial_step=0.0025, target_end_time=8., variable_steps=False)


# In[42]:


dynamics.get_history()


# ### Notice we created 11,215 data points, A LOT MORE than in runs #1 and #2

# In[43]:


dynamics.explain_time_advance()


# ## Plots of changes of concentration with time

# In[44]:


dynamics.plot_history(chemicals=["E_high", "E_low"], colors=["red", "grey"])


# In[45]:


dynamics.plot_history(chemicals=["A", "B", "C"])


# ### Yet again, the two plots have 4 distinctive intersections among them; locate and save them:

# In[46]:


run3 = []


# In[47]:


run3.append(dynamics.curve_intersect(t_start=1., t_end=2., chem1="E_high", chem2="E_low"))


# In[48]:


run3.append(dynamics.curve_intersect(t_start=2.31, t_end=2.33, chem1="A", chem2="B"))


# In[49]:


run3.append(dynamics.curve_intersect(t_start=3., t_end=4., chem1="A", chem2="C"))


# In[50]:


run3.append(dynamics.curve_intersect(t_start=3., t_end=4., chem1="B", chem2="C"))


# In[51]:


run3


# In[ ]:





# In[ ]:





# # Finally, compare (using a Euclidean metric) the discrepancy between the runs
# Run #3, at high resolution, could be thought of as "the closest to the actual values"

# In[52]:


# Discrepancy of run1 from run3
num.compare_results(run1, run3)


# In[53]:


# Discrepancy of run2 from run3
num.compare_results(run2, run3)


# Run 2 (with a lot fewer points than run 1), not surprisingly deviates a bit more from the (presumably) highly-accurate run 3

# #### The coordinates of the 4 critical points, from the 3 different runs, are pretty similar to one another - as can be easily seen:

# In[54]:


run1


# In[55]:


run2


# In[56]:


run3


# In[ ]: