# To reproduce Figure 1 in the associated Physiome paper,
# execute this script from the command line:
#
#   cd [PathToThisFile]
#   [PathToOpenCOR]/pythonshell Figure2.py

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# matplotlib.use('agg')

import opencor as opencor

# load the reference model
simulation = opencor.open_simulation("model.sedml")
data = simulation.data()


results = dict()

def run_sim(Q_flow):

    # reset everything in case we are running interactively and have existing results
    simulation.reset(True)
    simulation.clear_results()
    data.constants()["Apical_concentrations/glucose_m"] = 0.1
    data.constants()["Apical_concentrations/Na_m"] = 0.09
    data.constants()["Apical_concentrations/Cl_m"] = 0.09
    data.constants()["Apical_concentrations/K_m"] = 0.006

    data.constants()["A_GLUT2/n_GLUT"] = 2e9
    data.constants()["GLUT2/n_GLUT"] = 2e9
    data.constants()["phenomonological_constants/n_SGLT"] = 10e7
    data.constants()["ENaC/G_ENaC"] = 0.5
    data.constants()["Cell_concentration/L_A"] = 6e-5
    data.constants()["Cell_concentration/L_B"] = 6e-5
    # data.constants()["Blood_concentrations/Q_in"] = 9e-18
    data.constants()["Blood_concentrations/Q_in"] = Q_flow
    data.set_starting_point(0)
    data.set_ending_point(3000)
    data.set_point_interval(10)
    simulation.run()
    ds = simulation.results().data_store()
    J_GLUT = ds.voi_and_variables()["GLUT2/J_GLUT"].values()[-1]
    J_A_Na = ds.voi_and_variables()["Cell_concentration/J_A_Na"].values()[-1]
    J_W_A = ds.voi_and_variables()["Cell_concentration/J_w_A"].values()[-1]


    return (J_GLUT, J_A_Na, J_W_A)

if  __name__ == '__main__':
    # different values for Q_flow
    Q_flow = [6e-18, 1e-17, 3e-17, 1e-15]
    results1 = run_sim(Q_flow[0])
    results2 = run_sim(Q_flow[1])
    results3 = run_sim(Q_flow[2])
    results4 = run_sim(Q_flow[3])


    plt.figure(figsize=(12, 8))
    plt.subplot(131)

    X1 = 0
    X2 = 0.2
    Y1 = results1[0]
    Y2 = results2[0]
    Y3 = results3[0]
    Y4 = results4[0]
    Y_A = ((Y1 + Y2 + Y3 + Y4) / 4) * 0.001 * 60
    ymin, ymax = 0, 6e-11
    Y_exp1 = [4.5e-11]
    A1 = [5e-12]
    A2 = [4.5e-12]
    plt.bar(X1, Y_A, yerr=A1, align='center', color='b', alpha=0.5, ecolor='black', capsize=10, width=0.2,
            label='Model')
    plt.bar(X2, Y_exp1, yerr=A2, align='center', color='r', alpha=0.5, ecolor='black', capsize=10, width=0.2,
            label='Experiment')
    plt.ylim(ymin, ymax)
    plt.tick_params(axis='x', which='major', labelsize=15)
    plt.tick_params(axis='y', which='major', labelsize=15)
    # plt.xticks(y_pos, objects)
    plt.ylabel('Glucose Absorption Rate (mmol/min)', fontsize='14')
    plt.xticks([])
    plt.legend(loc='best', fontsize='14')
    plt.title('A')

    plt.subplot(132)

    X1 = 0
    X2 = 0.2
    Y5 = results1[1]
    Y6 = results2[1]
    Y7 = results3[1]
    Y8 = results4[1]
    Y_A = ((Y1 + Y2 + Y3 + Y4) / 4) * 0.001 * 60
    A3 = [2e-12]
    A4 = [7e-12]
    Y_B = ((Y5 + Y6 + Y7 + Y8) / 4) * 0.001 * 60
    Y_exp2 = [4.2e-11]
    ymin, ymax = 0, 8e-11
    plt.bar(X1, Y_B, align='center', yerr=A3, color='b', ecolor='black', capsize=10, alpha=0.5, width=0.2,
            label='Model')

    plt.bar(X2, Y_exp2, align='center', yerr=A4, color='r', ecolor='black', capsize=10, alpha=0.5, width=0.2,
            label='Experiment')

    plt.ylim(ymin, ymax)
    plt.tick_params(axis='x', which='major', labelsize=15)
    plt.tick_params(axis='y', which='major', labelsize=15)
    # ~ plt.xticks(y_pos, objects)
    plt.ylabel('Sodium Absorption Rate (mmol/min)', fontsize='14')
    plt.xticks([])
    plt.legend(loc='best', fontsize='14')
    plt.title('B')

    plt.subplot(133)

    ymin, ymax = 0, 4e-8

    Y9 = -results1[2] * 0.001 * 60
    Y10 = -results2[2] * 0.001 * 60
    Y11 = -results3[2] * 0.001 * 60
    Y12 = -results4[2] * 0.001 * 60
    A5 = [3e-9]
    A6 = [2e-9]
    Y_C = ((Y9 + Y10 + Y11 + Y12) / 4)
    print(Y_C)
    Y_exp3 = [2e-8]
    plt.bar(X1, Y_C, yerr=A5, align='center', color='b', alpha=0.5, ecolor='black', capsize=10, width=0.2,
            label='Model')
    plt.bar(X2, Y_exp3, yerr=A6, align='center', color='r', alpha=0.5, ecolor='black', capsize=10, width=0.2,
            label='Experiment')
    plt.ylim(ymin, ymax)
    # plt.xticks(y_pos, objects)
    plt.ylabel('Water Absorption Rate (mmol/min)', fontsize='14')
    plt.tick_params(axis='x', which='major', labelsize=15)
    plt.tick_params(axis='y', which='major', labelsize=15)
    plt.xticks([])
    plt.legend(loc='best', fontsize='14')
    plt.title('C')
    plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25,
                        wspace=0.35)

    plt.savefig('Figure02.png')
    plt.show()