# 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

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

# 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())

chem_data.describe_reactions()

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

dynamics.describe_state()

dynamics.get_history()

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)

df = dynamics.get_history()
df

dynamics.explain_time_advance()

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

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

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

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

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

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

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

# 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()

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

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

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

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

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

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()

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

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

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

df2 = dynamics2.get_history()
df2

df.loc[1]

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

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

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

new - old

df.loc[16]

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

df.loc[70]

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

df.loc[85]

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

# 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()

# 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()

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()