import pandas
pandas.options.display.max_rows = 12
import escher
from cameo import models, flux_variability_analysis, fba

"""
Explanation: Flux Variability Anlysis (FVA)
Load a few packages and functions.
End of explanation
"""


model = models.bigg.e_coli_core.copy()

"""
Explanation: First we load a model from the BiGG database (and make a copy of it).
End of explanation
"""


result = flux_variability_analysis(model)

"""
Explanation: Run flux variablity analysis
Calculate all flux ranges of all reactions in the model.
End of explanation
"""


result.data_frame

"""
Explanation: Inspect the result.
End of explanation
"""


result.data_frame.describe()

"""
Explanation: Get an overview of a few key statistics of the resulting flux ranges.
End of explanation
"""


result.plot(index=result.data_frame.index, height=1200)

"""
Explanation: Visualize the flux ranges.
End of explanation
"""


abs_flux_ranges = abs(result.data_frame.lower_bound - result.data_frame.upper_bound).to_dict()
escher.Builder('e_coli_core.Core metabolism', reaction_data=abs_flux_ranges).display_in_notebook()

"""
Explanation: Visualize the flux ranges on a pathway map of E. coli's central carbon metabolism.
End of explanation
"""


result.data_frame[result.data_frame.upper_bound > 500]

result_no_cyles = flux_variability_analysis(model, remove_cycles=True)

abs_flux_ranges = abs(result_no_cyles.data_frame.lower_bound - result_no_cyles.data_frame.upper_bound).to_dict()
escher.Builder('e_coli_core.Core metabolism', reaction_data=abs_flux_ranges).display_in_notebook()

"""
Explanation: Those reactions showing up in red are futile cyles.
End of explanation
"""


fba(model)

model_optimal = model.copy()

model_optimal.reactions.BIOMASS_Ecoli_core_w_GAM.lower_bound = 0.8739215069684299

result_max_obj = flux_variability_analysis(model_optimal, remove_cycles=True)

result_max_obj.plot(index=result_max_obj.data_frame.index, height=1200)

"""
Explanation: Run flux variability analysis for optimally growing E. coli
(Optimal) Flux Balance Analysis solutions are not necessariliy unique. Flux Variablity Analysis is a good tool for estimating the space of alternative optimal solutions.
End of explanation
"""


result_max_obj = flux_variability_analysis(model, fraction_of_optimum=1., remove_cycles=True)

result_max_obj.plot(index=result_max_obj.data_frame.index, height=1200)

"""
Explanation: This is actually such a common task that flux_variability_analysis provides an option for fixing the objective's flux at a certain percentage.
End of explanation
"""


sum(abs(result_max_obj.data_frame.lower_bound - result_max_obj.data_frame.upper_bound))

"""
Explanation: Turns out that in this small core metabolic model, the optimal solution is actually unique!
End of explanation
"""


percentage = (0.7 / model.solve().f) * 100
percentage

result_80perc_max_obj = flux_variability_analysis(model, fraction_of_optimum=percentage/100, remove_cycles=True)

result_80perc_max_obj.plot(index=result_80perc_max_obj.data_frame.index, height=1200)

"""
Explanation: Exercises
Exercise 1
Explore how relaxing the constraint on the growth rate affects the solution space:
1. Modify the code to explore flux ranges for $\mu \gt 0.7 \ h^{-1}$ 
1. Plot the sum of flux ranges over a range of percentages.
Exercise 2
Using FVA, determine all blocked reactions ($v = 0$) in the model.
Solutions
Solution 1
End of explanation
"""


flux_sums = []
optimum_percentages = range(50, 105, 5)
for i in optimum_percentages:
    df = flux_variability_analysis(model, fraction_of_optimum=i/100, remove_cycles=True).data_frame
    flux_sum = sum(abs(df.lower_bound - df.upper_bound))
    print("{}%: ".format(i), flux_sum)
    flux_sums.append(flux_sum)

import matplotlib.pyplot as plt

plt.plot(optimum_percentages, flux_sums)
plt.xlabel('Optimum (%)')
plt.ylabel('Flux sum [mmol gDW^-1 h^-1]')
plt.show()

"""
Explanation: Solution 2
End of explanation
"""


result = flux_variability_analysis(model, remove_cycles=True)

result.data_frame[(result.data_frame.lower_bound == 0) & (result.data_frame.upper_bound == 0)]

"""
Explanation: Solution 3
End of explanation
"""


