# Size of variable arrays:
sizeAlgebraic = 24
sizeStates = 24
sizeConstants = 58
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_voi = "time in component environment (minute)"
    legend_states[0] = "Glc_x_0 in component Glc_x_0 (millimolar)"
    legend_algebraic[0] = "inGlc in component inGlc (flux)"
    legend_states[1] = "Glc_x in component Glc_x (millimolar)"
    legend_algebraic[1] = "GlcTrans in component GlcTrans (flux)"
    legend_states[2] = "Glc in component Glc (millimolar)"
    legend_algebraic[2] = "HK in component HK (flux)"
    legend_states[3] = "G6P in component G6P (millimolar)"
    legend_algebraic[3] = "PGI in component PGI (flux)"
    legend_algebraic[10] = "storage in component storage (flux)"
    legend_states[4] = "F6P in component F6P (millimolar)"
    legend_algebraic[4] = "PFK in component PFK (flux)"
    legend_states[5] = "FBP in component FBP (millimolar)"
    legend_algebraic[5] = "ALD in component ALD (flux)"
    legend_states[6] = "GAP in component GAP (millimolar)"
    legend_algebraic[7] = "GAPDH in component GAPDH (flux)"
    legend_algebraic[6] = "TIM in component TIM (flux)"
    legend_states[7] = "DHAP in component DHAP (millimolar)"
    legend_algebraic[14] = "lpGlyc in component lpGlyc (flux)"
    legend_states[8] = "BPG in component BPG (millimolar)"
    legend_algebraic[8] = "lpPEP in component lpPEP (flux)"
    legend_states[9] = "PEP in component PEP (millimolar)"
    legend_algebraic[9] = "PK in component PK (flux)"
    legend_states[10] = "Pyr in component Pyr (millimolar)"
    legend_algebraic[11] = "PDC in component PDC (flux)"
    legend_states[11] = "ATP in component ATP (millimolar)"
    legend_algebraic[13] = "consum in component consum (flux)"
    legend_algebraic[15] = "AK in component AK (flux)"
    legend_states[12] = "ADP in component ADP (millimolar)"
    legend_states[13] = "AMP in component AMP (millimolar)"
    legend_states[14] = "CN_x_0 in component CN_x_0 (millimolar)"
    legend_algebraic[23] = "inCN in component inCN (flux)"
    legend_states[15] = "CN_x in component CN_x (millimolar)"
    legend_algebraic[21] = "lacto in component lacto (flux)"
    legend_states[16] = "ACA in component ACA (millimolar)"
    legend_algebraic[18] = "difACA in component difACA (flux)"
    legend_algebraic[12] = "ADH in component ADH (flux)"
    legend_states[17] = "ACA_x in component ACA_x (millimolar)"
    legend_algebraic[16] = "outACA in component outACA (flux)"
    legend_states[18] = "EtOH in component EtOH (millimolar)"
    legend_algebraic[19] = "difEtOH in component difEtOH (flux)"
    legend_states[19] = "EtOH_x in component EtOH_x (millimolar)"
    legend_algebraic[17] = "outEtOH in component outEtOH (flux)"
    legend_states[20] = "Glyc in component Glyc (millimolar)"
    legend_algebraic[22] = "difGlyc in component difGlyc (flux)"
    legend_states[21] = "Glyc_x in component Glyc_x (millimolar)"
    legend_algebraic[20] = "outGlyc in component outGlyc (flux)"
    legend_states[22] = "NADH in component NADH (millimolar)"
    legend_states[23] = "NAD in component NAD (millimolar)"
    legend_constants[0] = "k0 in component model_parameters (first_order_rate_constant)"
    legend_constants[1] = "y_vol in component model_parameters (dimensionless)"
    legend_constants[2] = "K2_Glc in component GlcTrans (millimolar)"
    legend_constants[3] = "K2_IG6P in component GlcTrans (millimolar)"
    legend_constants[4] = "K2_IIG6P in component GlcTrans (millimolar)"
    legend_constants[5] = "V_2m in component GlcTrans (flux)"
    legend_constants[6] = "P2 in component GlcTrans (dimensionless)"
    legend_constants[7] = "K3_DGlc in component HK (millimolar)"
    legend_constants[8] = "K3_Glc in component HK (millimolar)"
    legend_constants[9] = "K3_ATP in component HK (millimolar)"
    legend_constants[10] = "V_3m in component HK (flux)"
    legend_constants[11] = "V_4m in component PGI (flux)"
    legend_constants[12] = "K4_G6P in component PGI (millimolar)"
    legend_constants[13] = "K4_F6P in component PGI (millimolar)"
    legend_constants[14] = "K4_eq in component PGI (dimensionless)"
    legend_constants[15] = "K5 in component PFK (millimolar2)"
    legend_constants[16] = "kappa5 in component PFK (dimensionless)"
    legend_constants[17] = "V_5m in component PFK (flux)"
    legend_constants[18] = "K6_eq in component ALD (millimolar)"
    legend_constants[19] = "K6_FBP in component ALD (millimolar)"
    legend_constants[20] = "K6_DHAP in component ALD (millimolar)"
    legend_constants[21] = "K6_GAP in component ALD (millimolar)"
    legend_constants[22] = "K6_IGAP in component ALD (millimolar)"
    legend_constants[23] = "V_6r in component ALD (flux)"
    legend_constants[24] = "V_6f in component ALD (flux)"
    legend_constants[25] = "K7_eq in component TIM (dimensionless)"
    legend_constants[26] = "K7_DHAP in component TIM (millimolar)"
    legend_constants[27] = "K7_GAP in component TIM (millimolar)"
    legend_constants[28] = "V_7m in component TIM (flux)"
    legend_constants[29] = "K8_NAD in component GAPDH (millimolar)"
    legend_constants[30] = "K8_NADH in component GAPDH (millimolar)"
    legend_constants[31] = "K8_GAP in component GAPDH (millimolar)"
    legend_constants[32] = "K8_BPG in component GAPDH (millimolar)"
    legend_constants[33] = "K8_eq in component GAPDH (dimensionless)"
    legend_constants[34] = "V_8m in component GAPDH (flux)"
    legend_constants[35] = "k9f in component lpPEP (second_order_rate_constant)"
    legend_constants[36] = "k9r in component lpPEP (second_order_rate_constant)"
    legend_constants[37] = "K10_PEP in component PK (millimolar)"
    legend_constants[38] = "K10_ADP in component PK (millimolar)"
    legend_constants[39] = "V_10m in component PK (flux)"
    legend_constants[40] = "K11 in component PDC (millimolar)"
    legend_constants[41] = "V_11m in component PDC (flux)"
    legend_constants[42] = "K12_NADH in component ADH (millimolar)"
    legend_constants[43] = "K12_ACA in component ADH (millimolar)"
    legend_constants[44] = "V_12m in component ADH (flux)"
    legend_constants[45] = "k13 in component difEtOH (first_order_rate_constant)"
    legend_constants[46] = "K15_NADH in component lpGlyc (millimolar)"
    legend_constants[47] = "K15_INADH in component lpGlyc (millimolar)"
    legend_constants[48] = "K15_INAD in component lpGlyc (millimolar)"
    legend_constants[49] = "K15_DHAP in component lpGlyc (millimolar)"
    legend_constants[50] = "V_15m in component lpGlyc (flux)"
    legend_constants[51] = "k16 in component difGlyc (first_order_rate_constant)"
    legend_constants[52] = "k18 in component difACA (first_order_rate_constant)"
    legend_constants[53] = "k20 in component lacto (second_order_rate_constant)"
    legend_constants[54] = "k22 in component storage (second_order_rate_constant)"
    legend_constants[55] = "k23 in component consum (first_order_rate_constant)"
    legend_constants[56] = "k24f in component AK (second_order_rate_constant)"
    legend_constants[57] = "k24r in component AK (second_order_rate_constant)"
    legend_rates[0] = "d/dt Glc_x_0 in component Glc_x_0 (millimolar)"
    legend_rates[1] = "d/dt Glc_x in component Glc_x (millimolar)"
    legend_rates[2] = "d/dt Glc in component Glc (millimolar)"
    legend_rates[3] = "d/dt G6P in component G6P (millimolar)"
    legend_rates[4] = "d/dt F6P in component F6P (millimolar)"
    legend_rates[5] = "d/dt FBP in component FBP (millimolar)"
    legend_rates[6] = "d/dt GAP in component GAP (millimolar)"
    legend_rates[7] = "d/dt DHAP in component DHAP (millimolar)"
    legend_rates[8] = "d/dt BPG in component BPG (millimolar)"
    legend_rates[9] = "d/dt PEP in component PEP (millimolar)"
    legend_rates[10] = "d/dt Pyr in component Pyr (millimolar)"
    legend_rates[11] = "d/dt ATP in component ATP (millimolar)"
    legend_rates[12] = "d/dt ADP in component ADP (millimolar)"
    legend_rates[13] = "d/dt AMP in component AMP (millimolar)"
    legend_rates[14] = "d/dt CN_x_0 in component CN_x_0 (millimolar)"
    legend_rates[15] = "d/dt CN_x in component CN_x (millimolar)"
    legend_rates[16] = "d/dt ACA in component ACA (millimolar)"
    legend_rates[17] = "d/dt ACA_x in component ACA_x (millimolar)"
    legend_rates[18] = "d/dt EtOH in component EtOH (millimolar)"
    legend_rates[19] = "d/dt EtOH_x in component EtOH_x (millimolar)"
    legend_rates[20] = "d/dt Glyc in component Glyc (millimolar)"
    legend_rates[21] = "d/dt Glyc_x in component Glyc_x (millimolar)"
    legend_rates[22] = "d/dt NADH in component NADH (millimolar)"
    legend_rates[23] = "d/dt NAD in component NAD (millimolar)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = 24.0
    states[1] = 6.7
    states[2] = 0.573074
    states[3] = 4.2
    states[4] = 0.49
    states[5] = 4.64
    states[6] = 0.115
    states[7] = 2.95
    states[8] = 0.00027
    states[9] = 0.04
    states[10] = 8.7
    states[11] = 2.1
    states[12] = 1.5
    states[13] = 0.33
    states[14] = 5.60
    states[15] = 5.20358
    states[16] = 1.48153
    states[17] = 1.28836
    states[18] = 19.2379
    states[19] = 16.4514
    states[20] = 4.196
    states[21] = 1.68478
    states[22] = 0.33
    states[23] = 0.65
    constants[0] = 0.048
    constants[1] = 59.0
    constants[2] = 1.7
    constants[3] = 1.2
    constants[4] = 7.2
    constants[5] = 1014.96
    constants[6] = 1.0
    constants[7] = 0.37
    constants[8] = 0.0
    constants[9] = 0.1
    constants[10] = 5.7547
    constants[11] = 496.042
    constants[12] = 0.8
    constants[13] = 0.15
    constants[14] = 0.13
    constants[15] = 0.021
    constants[16] = 0.15
    constants[17] = 45.4327
    constants[18] = 0.081
    constants[19] = 0.3
    constants[20] = 2.0
    constants[21] = 4.0
    constants[22] = 10.0
    constants[23] = 1.10391E4
    constants[24] = 2.20782E3
    constants[25] = 0.055
    constants[26] = 1.23
    constants[27] = 1.27
    constants[28] = 1.16365E2
    constants[29] = 0.1
    constants[30] = 0.06
    constants[31] = 0.6
    constants[32] = 0.01
    constants[33] = 0.0055
    constants[34] = 8.33858E2
    constants[35] = 4.43866E5
    constants[36] = 1.52862E3
    constants[37] = 0.2
    constants[38] = 0.17
    constants[39] = 3.43096E2
    constants[40] = 0.3
    constants[41] = 5.31328E1
    constants[42] = 0.1
    constants[43] = 0.71
    constants[44] = 8.98023E1
    constants[45] = 16.72
    constants[46] = 0.13
    constants[47] = 0.034
    constants[48] = 0.13
    constants[49] = 25.0
    constants[50] = 8.14797E1
    constants[51] = 1.9
    constants[52] = 24.7
    constants[53] = 2.83828E-3
    constants[54] = 2.25932
    constants[55] = 3.20760
    constants[56] = 4.32900E2
    constants[57] = 1.33333E2
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[0] = constants[0]*(states[0]-states[1])
    rates[0] = algebraic[0]
    algebraic[1] = (constants[5]/constants[1])*((states[1]/constants[2])/(1.00000+states[1]/constants[2]+((constants[6]*(states[1]/constants[2])+1.00000)/(constants[6]*(states[2]/constants[2])+1.00000))*(1.00000+states[2]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))-(constants[5]/constants[1])*((states[2]/constants[2])/(1.00000+states[2]/constants[2]+((constants[6]*(states[2]/constants[2])+1.00000)/(constants[6]*(states[1]/constants[2])+1.00000))*(1.00000+states[1]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))
    rates[1] = algebraic[0]+algebraic[1]
    algebraic[2] = (constants[10]*states[11]*states[2])/(constants[7]*constants[9]+constants[8]*states[11]+constants[9]*states[2]+states[2]*states[11])
    rates[2] = algebraic[1]+algebraic[2]
    algebraic[3] = (constants[11]*states[3])/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])-(constants[11]*(states[4]/constants[14]))/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])
    algebraic[4] = (constants[17]*(power(states[4], 2.00000)))/(constants[15]*(1.00000+constants[16]*(states[11]/states[13]))+power(states[4], 2.00000))
    rates[4] = algebraic[3]+algebraic[4]
    algebraic[5] = (constants[24]*states[5])/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))-(constants[24]*((states[6]*states[7])/constants[18]))/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))
    rates[5] = algebraic[4]+algebraic[5]
    algebraic[7] = (constants[34]*states[6]*states[23])/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))-(constants[34]*((states[8]*states[22])/constants[33]))/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))
    algebraic[6] = (constants[28]*states[7])/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])-(constants[28]*(states[6]/constants[25]))/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])
    rates[6] = algebraic[5]+algebraic[6]+algebraic[7]
    algebraic[8] = constants[35]*states[8]*states[12]-constants[36]*states[9]*states[11]
    rates[8] = algebraic[7]+algebraic[8]
    algebraic[9] = (constants[39]*states[12]*states[9])/((constants[37]+states[9])*(constants[38]+states[12]))
    rates[9] = algebraic[8]+algebraic[9]
    algebraic[10] = constants[54]*states[11]*states[3]
    rates[3] = algebraic[2]+algebraic[3]+algebraic[10]
    algebraic[11] = (constants[41]*states[10])/(constants[40]+states[10])
    rates[10] = algebraic[9]+algebraic[11]
    algebraic[14] = (constants[50]*states[7])/(constants[49]*(1.00000+(constants[47]/states[22])*(1.00000+states[23]/constants[48]))+states[7]*(1.00000+(constants[46]/states[22])*(1.00000+states[23]/constants[48])))
    rates[7] = algebraic[5]+algebraic[6]+algebraic[14]
    algebraic[13] = constants[55]*states[11]
    algebraic[15] = constants[56]*states[13]*states[11]-constants[57]*(power(states[12], 2.00000))
    rates[11] = algebraic[4]+algebraic[9]+algebraic[8]+algebraic[10]+algebraic[2]+algebraic[13]+algebraic[15]
    rates[12] = algebraic[4]+algebraic[9]+algebraic[8]+algebraic[10]+algebraic[2]+algebraic[13]+algebraic[15]
    rates[13] = algebraic[15]
    algebraic[12] = (constants[44]*states[16]*states[22])/((constants[42]+states[22])*(constants[43]+states[16]))
    rates[22] = algebraic[7]+algebraic[12]+algebraic[14]
    rates[23] = algebraic[7]+algebraic[12]+algebraic[14]
    algebraic[18] = (constants[52]/constants[1])*(states[16]-states[17])
    rates[16] = algebraic[18]+algebraic[11]+algebraic[12]
    algebraic[19] = (constants[45]/constants[1])*(states[18]-states[19])
    rates[18] = algebraic[19]+algebraic[12]
    algebraic[17] = constants[0]*states[19]
    rates[19] = algebraic[19]+algebraic[17]
    algebraic[21] = constants[53]*states[17]*states[15]
    algebraic[16] = constants[0]*states[17]
    rates[17] = algebraic[21]+algebraic[18]+algebraic[16]
    algebraic[22] = (constants[51]/constants[1])*(states[20]-states[21])
    rates[20] = algebraic[22]+algebraic[14]
    algebraic[20] = constants[0]*states[21]
    rates[21] = algebraic[22]+algebraic[20]
    algebraic[23] = constants[0]*(states[14]-states[15])
    rates[14] = algebraic[23]
    rates[15] = algebraic[21]+algebraic[23]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[0] = constants[0]*(states[0]-states[1])
    algebraic[1] = (constants[5]/constants[1])*((states[1]/constants[2])/(1.00000+states[1]/constants[2]+((constants[6]*(states[1]/constants[2])+1.00000)/(constants[6]*(states[2]/constants[2])+1.00000))*(1.00000+states[2]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))-(constants[5]/constants[1])*((states[2]/constants[2])/(1.00000+states[2]/constants[2]+((constants[6]*(states[2]/constants[2])+1.00000)/(constants[6]*(states[1]/constants[2])+1.00000))*(1.00000+states[1]/constants[2])+states[3]/constants[3]+(states[2]*states[3])/(constants[2]*constants[4])))
    algebraic[2] = (constants[10]*states[11]*states[2])/(constants[7]*constants[9]+constants[8]*states[11]+constants[9]*states[2]+states[2]*states[11])
    algebraic[3] = (constants[11]*states[3])/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])-(constants[11]*(states[4]/constants[14]))/(constants[12]+states[3]+(constants[12]/constants[13])*states[4])
    algebraic[4] = (constants[17]*(power(states[4], 2.00000)))/(constants[15]*(1.00000+constants[16]*(states[11]/states[13]))+power(states[4], 2.00000))
    algebraic[5] = (constants[24]*states[5])/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))-(constants[24]*((states[6]*states[7])/constants[18]))/(constants[19]+states[5]+(states[6]*constants[20]*constants[24])/(constants[18]*constants[23])+(states[7]*constants[21]*constants[24])/(constants[18]*constants[23])+(states[5]*states[6])/constants[22]+(states[6]*states[7]*constants[24])/(constants[18]*constants[23]))
    algebraic[7] = (constants[34]*states[6]*states[23])/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))-(constants[34]*((states[8]*states[22])/constants[33]))/(constants[31]*constants[29]*(1.00000+states[6]/constants[31]+states[8]/constants[32])*(1.00000+states[23]/constants[29]+states[22]/constants[30]))
    algebraic[6] = (constants[28]*states[7])/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])-(constants[28]*(states[6]/constants[25]))/(constants[26]+states[7]+(constants[26]/constants[27])*states[6])
    algebraic[8] = constants[35]*states[8]*states[12]-constants[36]*states[9]*states[11]
    algebraic[9] = (constants[39]*states[12]*states[9])/((constants[37]+states[9])*(constants[38]+states[12]))
    algebraic[10] = constants[54]*states[11]*states[3]
    algebraic[11] = (constants[41]*states[10])/(constants[40]+states[10])
    algebraic[14] = (constants[50]*states[7])/(constants[49]*(1.00000+(constants[47]/states[22])*(1.00000+states[23]/constants[48]))+states[7]*(1.00000+(constants[46]/states[22])*(1.00000+states[23]/constants[48])))
    algebraic[13] = constants[55]*states[11]
    algebraic[15] = constants[56]*states[13]*states[11]-constants[57]*(power(states[12], 2.00000))
    algebraic[12] = (constants[44]*states[16]*states[22])/((constants[42]+states[22])*(constants[43]+states[16]))
    algebraic[18] = (constants[52]/constants[1])*(states[16]-states[17])
    algebraic[19] = (constants[45]/constants[1])*(states[18]-states[19])
    algebraic[17] = constants[0]*states[19]
    algebraic[21] = constants[53]*states[17]*states[15]
    algebraic[16] = constants[0]*states[17]
    algebraic[22] = (constants[51]/constants[1])*(states[20]-states[21])
    algebraic[20] = constants[0]*states[21]
    algebraic[23] = constants[0]*(states[14]-states[15])
    return algebraic

def solve_model():
    """Solve model with ODE solver"""
    from scipy.integrate import ode
    # Initialise constants and state variables
    (init_states, constants) = initConsts()

    # Set timespan to solve over
    voi = linspace(0, 10, 500)

    # Construct ODE object to solve
    r = ode(computeRates)
    r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1)
    r.set_initial_value(init_states, voi[0])
    r.set_f_params(constants)

    # Solve model
    states = array([[0.0] * len(voi)] * sizeStates)
    states[:,0] = init_states
    for (i,t) in enumerate(voi[1:]):
        if r.successful():
            r.integrate(t)
            states[:,i+1] = r.y
        else:
            break

    # Compute algebraic variables
    algebraic = computeAlgebraic(constants, states, voi)
    return (voi, states, algebraic)

def plot_model(voi, states, algebraic):
    """Plot variables against variable of integration"""
    import pylab
    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
    pylab.figure(1)
    pylab.plot(voi,vstack((states,algebraic)).T)
    pylab.xlabel(legend_voi)
    pylab.legend(legend_states + legend_algebraic, loc='best')
    pylab.show()

if __name__ == "__main__":
    (voi, states, algebraic) = solve_model()
    plot_model(voi, states, algebraic)