text:
function diffeq(du,u,p,t)
    #=Integrated model
    - Birtwistle et al., Mol. Syst. Biol. (2007)
    - Nakakuki et al., Cell (2010)
    =#

    # extracellular volume to cytoplasmic volume ratio
    VeVc = 33.3

    # fraction definitions
    if u[V.sigmaS] + u[V.sigmaSP] + u[V.sigmaSP_G] > 0.0
        fsigmaS = u[V.sigmaS]/(u[V.sigmaS] + u[V.sigmaSP] + u[V.sigmaSP_G])
    else
        fsigmaS = 0.0
    end

    if u[V.sigmaG] + u[V.sigmaG_A] + u[V.sigmaG_O] + u[V.A_sigmaG_O] > 0.0
        fsigmaG = u[V.sigmaG]/(u[V.sigmaG] + u[V.sigmaG_A] + u[V.sigmaG_O] + u[V.A_sigmaG_O])
    else
        fsigmaG = 0.0
    end

    if u[V.sigmaA] + u[V.sigmaAP] + u[V.sigmaAP_S] + u[V.sigmaAP_R] + u[V.sigmaAP_I] + u[V.sigmaAP_T] > 0.0
        fsigmaA = u[V.sigmaA]/(u[V.sigmaA] + u[V.sigmaAP] + u[V.sigmaAP_S] + u[V.sigmaAP_R] + u[V.sigmaAP_I] + u[V.sigmaAP_T])
    else
        fsigmaA = 0.0
    end

    if u[V.sigmaR] + u[V.sigmaRP] > 0.0
        fsigmaR = u[V.sigmaR]/(u[V.sigmaR] + u[V.sigmaRP])
    else
        fsigmaR = 0.0
    end

    sigmaEP = u[V.E11P] + u[V.E12P] + u[V.E23P] + u[V.E24P] + u[V.E34P] + u[V.E44P] + u[V.E13P] + u[V.E14P]
    if sigmaEP > 0.0
        f11 = u[V.E11P]/sigmaEP
    else
        f11 = 0.0
    end

    v = Dict{Int64,Float64}()

    # Birtwistle et al., MSB (2007)
    v[1] = (p[C.kon1]*u[V.E]*u[V.E1] - p[C.EGF_off]*u[V.E_E1])
    v[2] = (p[C.kon2]*u[V.H]*u[V.E3] - p[C.HRGoff_3]*u[V.H_E3])
    v[3] = (p[C.kon3]*u[V.H]*u[V.E4] - p[C.HRGoff_4]*u[V.H_E4])
    v[4] = (p[C.kon4]*u[V.E_E1]*u[V.E_E1] - p[C.koff4]*u[V.E11])
    v[5] = (p[C.kon5]*u[V.E_E1]*u[V.E2] - p[C.koff5]*u[V.E12])
    v[6] = (p[C.kon6]*u[V.H_E3]*u[V.E2] - p[C.koff6]*u[V.E23])
    v[7] = (p[C.kon7]*u[V.H_E3]*u[V.H_E4] - p[C.koff7]*u[V.E34])
    v[8] = (p[C.kon8]*u[V.H_E4]*u[V.E2] - p[C.koff8]*u[V.E24])
    v[9] = (p[C.kon9]*u[V.H_E4]*u[V.H_E4] - p[C.koff9]*u[V.E44])
    v[10] = (p[C.kf10]*u[V.E11] - p[C.VmaxPY]*u[V.E11P]/(p[C.KmPY] + u[V.E11P]) - p[C.kPTP10]*u[V.sigmaT]*u[V.E11P])
    v[11] = (p[C.kf11]*u[V.E12] - p[C.VmaxPY]*u[V.E12P]/(p[C.KmPY] + u[V.E12P]) - p[C.kPTP11]*u[V.sigmaT]*u[V.E12P])
    v[12] = (p[C.kf12]*u[V.E23] - p[C.VmaxPY]*u[V.E23P]/(p[C.KmPY] + u[V.E23P]) - p[C.kPTP12]*u[V.sigmaT]*u[V.E23P])
    v[13] = (p[C.kf13]*u[V.E34] - p[C.VmaxPY]*u[V.E34P]/(p[C.KmPY] + u[V.E34P]) - p[C.kPTP13]*u[V.sigmaT]*u[V.E34P])
    v[14] = (p[C.kf14]*u[V.E24] - p[C.VmaxPY]*u[V.E24P]/(p[C.KmPY] + u[V.E24P]) - p[C.kPTP14]*u[V.sigmaT]*u[V.E24P])
    v[15] = (p[C.kf15]*u[V.E44] - p[C.VmaxPY]*u[V.E44P]/(p[C.KmPY] + u[V.E44P]) - p[C.kPTP15]*u[V.sigmaT]*u[V.E44P])
    v[16] = (4*p[C.kon16]*u[V.E11P]*u[V.G] - p[C.koff16]*fsigmaG*u[V.E11G])
    v[17] = (8*p[C.kon17]*u[V.E11P]*u[V.S] - p[C.koff17]*fsigmaS*u[V.E11S])
    v[18] = (2*p[C.kon18]*u[V.E11P]*u[V.R] - p[C.koff18]*fsigmaR*u[V.E11R])
    v[19] = (3*p[C.kon19]*u[V.E12P]*u[V.G] - p[C.koff19]*fsigmaG*u[V.E12G])
    v[20] = (6*p[C.kon20]*u[V.E12P]*u[V.S] - p[C.koff20]*fsigmaS*u[V.E12S])
    v[21] = (2*p[C.kon21]*u[V.E12P]*u[V.R] - p[C.koff21]*fsigmaR*u[V.E12R])
    v[22] = (3*p[C.kon22]*u[V.E23P]*u[V.G] - p[C.koff22]*fsigmaG*u[V.E23G])
    v[23] = (3*p[C.kon23]*u[V.E23P]*u[V.S] - p[C.koff23]*fsigmaS*u[V.E23S])
    v[24] = (3*p[C.kon24]*u[V.E23P]*u[V.I] - p[C.koff24]*u[V.E23I])
    v[25] = (2*p[C.kon25]*u[V.E23P]*u[V.R] - p[C.koff25]*fsigmaR*u[V.E23R])
    v[26] = (4*p[C.kon26]*u[V.E34P]*u[V.G] - p[C.koff26]*fsigmaG*u[V.E34G])
    v[27] = (3*p[C.kon27]*u[V.E34P]*u[V.S] - p[C.koff27]*fsigmaS*u[V.E34S])
    v[28] = (4*p[C.kon28]*u[V.E34P]*u[V.I] - p[C.koff28]*u[V.E34I])
    v[29] = (2*p[C.kon29]*u[V.E34P]*u[V.R] - p[C.koff29]*fsigmaR*u[V.E34R])
    v[30] = (3*p[C.kon30]*u[V.E24P]*u[V.G] - p[C.koff30]*fsigmaG*u[V.E24G])
    v[31] = (4*p[C.kon31]*u[V.E24P]*u[V.S] - p[C.koff31]*fsigmaS*u[V.E24S])
    v[32] = (1*p[C.kon32]*u[V.E24P]*u[V.I] - p[C.koff32]*u[V.E24I])
    v[33] = (2*p[C.kon33]*u[V.E24P]*u[V.R] - p[C.koff33]*fsigmaR*u[V.E24R])
    v[34] = (4*p[C.kon34]*u[V.E44P]*u[V.G] - p[C.koff34]*fsigmaG*u[V.E44G])
    v[35] = (4*p[C.kon35]*u[V.E44P]*u[V.S] - p[C.koff35]*fsigmaS*u[V.E44S])
    v[36] = (2*p[C.kon36]*u[V.E44P]*u[V.I] - p[C.koff36]*u[V.E44I])
    v[37] = (2*p[C.kon37]*u[V.E44P]*u[V.R] - p[C.koff37]*fsigmaR*u[V.E44R])
    v[38] = (p[C.kf38]*u[V.sigmaS]*sigmaEP - p[C.VmaxPY]*u[V.sigmaSP]/(p[C.KmPY] + u[V.sigmaSP]) - p[C.kPTP38]*u[V.sigmaT]*u[V.sigmaSP])
    v[39] = (p[C.kf39]*u[V.sigmaA]*sigmaEP - p[C.VmaxPY]*u[V.sigmaAP]/(p[C.KmPY] + u[V.sigmaAP]) - p[C.kPTP39]*u[V.sigmaT]*u[V.sigmaAP])
    v[40] = (p[C.kon40]*u[V.sigmaG]*u[V.O] - p[C.koff40]*u[V.sigmaG_O])
    v[41] = (p[C.kon41]*u[V.sigmaG]*u[V.A] - p[C.koff41]*u[V.sigmaG_A]*fsigmaA)
    v[42] = (p[C.kon42]*u[V.sigmaSP]*u[V.G] - p[C.koff42]*u[V.sigmaSP_G]*fsigmaG)
    v[43] = (3*p[C.kon43]*u[V.sigmaAP]*u[V.S] - p[C.koff43]*u[V.sigmaAP_S]*fsigmaS)
    v[44] = (3*p[C.kon44]*u[V.sigmaAP]*u[V.I] - p[C.koff44]*u[V.sigmaAP_I])
    v[45] = (2*p[C.kon45]*u[V.sigmaAP]*u[V.R] - p[C.koff45]*u[V.sigmaAP_R]*fsigmaR)
    v[46] = (p[C.kon46]*u[V.P3]*u[V.A] - p[C.koff46]*u[V.P3_A]*fsigmaA)
    v[47] = (p[C.kf47]*u[V.P3]*u[V.Akt]/(p[C.Kmf47] + u[V.Akt]) - p[C.Vmaxr47]*u[V.Aktstar]/(p[C.Kmr47] + u[V.Aktstar]))
    v[48] = (p[C.kf48]*(1 - u[V.fint]*f11)*u[V.sigmaI]*u[V.P2]/(p[C.Kmf48] + u[V.P2]) - 3*p[C.PTEN]*u[V.P3]/(p[C.Kmr48] + u[V.P3]))
    v[49] = (p[C.kf49]*u[V.sigmaO]*u[V.RsD]/(p[C.Kmf49] + u[V.RsD]) - p[C.kr49]*u[V.sigmaR]*u[V.RsT]/(p[C.Kmr49] + u[V.RsT]) - p[C.kr49b]*u[V.sigmaRP]*u[V.RsT]/(p[C.Kmr49b] + u[V.RsT]) - p[C.kcon49]*u[V.RsT])
    v[50] = (p[C.kf50]*u[V.sigmaR]*sigmaEP - p[C.VmaxPY]*u[V.sigmaRP]/(p[C.KmPY] + u[V.sigmaRP]) - p[C.kPTP50]*u[V.sigmaT]*u[V.sigmaRP])
    v[51] = (p[C.kf51]*u[V.RsT]*u[V.Raf]/(p[C.Kmf51] + u[V.Raf]) - p[C.Vmaxr51]*u[V.Rafstar]/(p[C.Kmrb51] + u[V.Rafstar]))
    v[52] = (p[C.kf52]*u[V.Rafstar]*u[V.MEK]/(p[C.Kmf52] + u[V.MEK]) - p[C.Vmaxr52]*u[V.ppMEKc]/(p[C.Kmr52] + u[V.ppMEKc]))
    # v[53] No reaction
    v[54] = (p[C.kf54]*u[V.O]*u[V.ppERKc]/(p[C.Kmf54] + u[V.O]) - p[C.Vmaxr54]*u[V.OP]/(p[C.Kmr54] + u[V.OP]))
    v[55] = (p[C.kf55]*u[V.A]*u[V.ppERKc]/(p[C.Kmf55] + u[V.A]) - p[C.Vmaxr55]*u[V.AP]/(p[C.Kmr55] + u[V.AP]))
    # v[56] No reaction
    v[57] = (p[C.kon57]*u[V.P3_A]*u[V.G] - p[C.koff57]*u[V.sigmaA_G])
    v[58] = (p[C.kon58]*u[V.sigmaA_G]*u[V.O] - p[C.koff58]*u[V.sigmaA_G_O])
    v[59] = (p[C.kon59]*u[V.sigmaG_O]*u[V.A] - p[C.koff59]*u[V.A_sigmaG_O]*fsigmaA)
    v[60] = (p[C.kon60]*u[V.sigmaG_A]*u[V.O] - p[C.koff60]*u[V.A_sigmaG_O])
    v[61] = (p[C.kon61]*u[V.H_E3]*u[V.E_E1] - p[C.koff61]*u[V.E13])
    v[62] = (p[C.kon62]*u[V.H_E4]*u[V.E_E1] - p[C.koff62]*u[V.E14])
    v[63] = (p[C.kf63]*u[V.E13] - p[C.VmaxPY]*u[V.E13P]/(p[C.KmPY] + u[V.E13P]) - p[C.kPTP63]*u[V.sigmaT]*u[V.E13P])
    v[64] = (p[C.kf64]*u[V.E14] - p[C.VmaxPY]*u[V.E14P]/(p[C.KmPY] + u[V.E14P]) - p[C.kPTP64]*u[V.sigmaT]*u[V.E14P])
    v[65] = (4*p[C.kon65]*u[V.E13P]*u[V.G] - p[C.koff65]*fsigmaG*u[V.E13G])
    v[66] = (5*p[C.kon66]*u[V.E13P]*u[V.S] - p[C.koff66]*fsigmaS*u[V.E13S])
    v[67] = (3*p[C.kon67]*u[V.E13P]*u[V.I] - p[C.koff67]*u[V.E13I])
    v[68] = (2*p[C.kon68]*u[V.E13P]*u[V.R] - p[C.koff68]*fsigmaR*u[V.E13R])
    v[69] = (4*p[C.kon69]*u[V.E14P]*u[V.G] - p[C.koff69]*fsigmaG*u[V.E14G])
    v[70] = (6*p[C.kon70]*u[V.E14P]*u[V.S] - p[C.koff70]*fsigmaS*u[V.E14S])
    v[71] = (1*p[C.kon71]*u[V.E14P]*u[V.I] - p[C.koff71]*u[V.E14I])
    v[72] = (2*p[C.kon72]*u[V.E14P]*u[V.R] - p[C.koff72]*fsigmaR*u[V.E14R])
    v[73] = (4*p[C.kon73]*u[V.E11P]*u[V.T] - p[C.koff73]*u[V.E11T])
    v[74] = (3*p[C.kon74]*u[V.E12P]*u[V.T] - p[C.koff74]*u[V.E12T])
    v[75] = (2*p[C.kon75]*u[V.E23P]*u[V.T] - p[C.koff75]*u[V.E23T])
    v[76] = (2*p[C.kon76]*u[V.E34P]*u[V.T] - p[C.koff76]*u[V.E34T])
    v[77] = (2*p[C.kon77]*u[V.E24P]*u[V.T] - p[C.koff77]*u[V.E24T])
    v[78] = (2*p[C.kon78]*u[V.E44P]*u[V.T] - p[C.koff78]*u[V.E44T])
    v[79] = (3*p[C.kon79]*u[V.E13P]*u[V.T] - p[C.koff79]*u[V.E13T])
    v[80] = (3*p[C.kon80]*u[V.E14P]*u[V.T] - p[C.koff80]*u[V.E14T])
    v[81] = (p[C.kf81]*u[V.E1]*u[V.ppERKc]/(p[C.Kmf81] + u[V.E1]) - p[C.Vmaxr81]*u[V.E1_PT]/(p[C.Kmr81] + u[V.E1_PT]))
    v[82] = (p[C.kf82]*u[V.E2]*u[V.ppERKc]/(p[C.Kmf82] + u[V.E2]) - p[C.Vmaxr82]*u[V.E2_PT]/(p[C.Kmr82] + u[V.E2_PT]))
    v[83] = (p[C.kf83]*u[V.E4]*u[V.ppERKc]/(p[C.Kmf83] + u[V.E4]) - p[C.Vmaxr83]*u[V.E4_PT]/(p[C.Kmr83] + u[V.E4_PT]))
    v[84] = (p[C.kf84]*u[V.E_E1]*u[V.ppERKc]/(p[C.Kmf84] + u[V.E_E1]) - p[C.Vmaxr84]*u[V.E_E1_PT]/(p[C.Kmr84] + u[V.E_E1_PT]))
    v[85] = (p[C.kf85]*u[V.H_E4]*u[V.ppERKc]/(p[C.Kmf85] + u[V.H_E4]) - p[C.Vmaxr85]*u[V.H_E4_PT]/(p[C.Kmr85] + u[V.H_E4_PT]))
    v[86] = (p[C.kon86]*u[V.E]*u[V.E1_PT] - p[C.EGF_off]*u[V.E_E1_PT])
    v[87] = (p[C.kon87]*u[V.H]*u[V.E4_PT] - p[C.HRGoff_4]*u[V.H_E4_PT])
    v[88] = (2*p[C.kon88]*u[V.sigmaAP]*u[V.T] - p[C.koff88]*u[V.sigmaAP_T])
    # E11 degradation
    v[89] = p[C.kdeg]*u[V.E11P]
    v[90] = p[C.kdeg]*u[V.E11G]
    v[91] = p[C.kdeg]*u[V.E11S]
    v[92] = p[C.kdeg]*u[V.E11R]
    v[93] = p[C.kdeg]*u[V.E11T]

    # Nakakuki et al., Cell (2010)
    v[94] = p[C.V1] * u[V.ppMEKc] * u[V.ERKc] /  ( p[C.Km1] * (1 + u[V.pERKc] / p[C.Km2]) + u[V.ERKc] )
    v[95] = p[C.V2] * u[V.ppMEKc] * u[V.pERKc] /  ( p[C.Km2] * (1 + u[V.ERKc] / p[C.Km1]) + u[V.pERKc] )
    v[96] = p[C.V3] * u[V.pERKc] /  ( p[C.Km3] * (1 + u[V.ppERKc] / p[C.Km4]) + u[V.pERKc] )
    v[97] = p[C.V4] * u[V.ppERKc] /  ( p[C.Km4]* (1 + u[V.pERKc] / p[C.Km3]) + u[V.ppERKc] )
    v[98] = p[C.V5] * u[V.pERKn] /  ( p[C.Km5] * (1 + u[V.ppERKn] / p[C.Km6]) + u[V.pERKn] )
    v[99] = p[C.V6] * u[V.ppERKn] /  ( p[C.Km6] * (1 + u[V.pERKn] / p[C.Km5]) + u[V.ppERKn] )
    v[100] = p[C.KimERK] * u[V.ERKc] - p[C.KexERK] * (p[C.Vn]/p[C.Vc]) * u[V.ERKn]
    v[101] = p[C.KimpERK] * u[V.pERKc] - p[C.KexpERK] * (p[C.Vn]/p[C.Vc]) * u[V.pERKn]
    v[102] = p[C.KimppERK] * u[V.ppERKc] - p[C.KexppERK] * (p[C.Vn]/p[C.Vc]) * u[V.ppERKn]
    v[103] = p[C.V10] * u[V.ppERKn]^p[C.n10] / ( p[C.Km10]^p[C.n10] + u[V.ppERKn]^p[C.n10] )
    v[104] = p[C.p11] * u[V.PreduspmRNAn]
    v[105] = p[C.p12] * u[V.duspmRNAc]
    v[106] = p[C.p13] * u[V.duspmRNAc]
    v[107] = p[C.V14] * u[V.ppERKc] * u[V.DUSPc] / ( p[C.Km14] + u[V.DUSPc] )
    v[108] = p[C.V15] * u[V.pDUSPc] / ( p[C.Km15] + u[V.pDUSPc] )
    v[109] = p[C.p16] * u[V.DUSPc]
    v[110] = p[C.p17] * u[V.pDUSPc]
    v[111] = p[C.KimDUSP] * u[V.DUSPc] - p[C.KexDUSP] * (p[C.Vn]/p[C.Vc]) * u[V.DUSPn]
    v[112] = p[C.KimpDUSP] * u[V.pDUSPc] - p[C.KexpDUSP] * (p[C.Vn]/p[C.Vc]) * u[V.pDUSPn]
    v[113] = p[C.V20] * u[V.ppERKn] * u[V.DUSPn] / ( p[C.Km20] + u[V.DUSPn] )
    v[114] = p[C.V21] * u[V.pDUSPn] / ( p[C.Km21] + u[V.pDUSPn] )
    v[115] = p[C.p22] * u[V.DUSPn]
    v[116] = p[C.p23] * u[V.pDUSPn]
    v[117] = p[C.V24] * u[V.ppERKc] * u[V.RSKc] / ( p[C.Km24] + u[V.RSKc] )
    v[118] = p[C.V25] * u[V.pRSKc] / ( p[C.Km25] + u[V.pRSKc] )
    v[119] = p[C.KimRSK] * u[V.pRSKc] - p[C.KexRSK] * (p[C.Vn]/p[C.Vc]) * u[V.pRSKn]
    v[120] = p[C.V27] * u[V.pRSKn] * u[V.CREBn] / ( p[C.Km27] + u[V.CREBn] )
    v[121] = p[C.V28] * u[V.pCREBn] / ( p[C.Km28] + u[V.pCREBn] )
    v[122] = p[C.V29] * u[V.ppERKn] * u[V.Elk1n] / ( p[C.Km29] + u[V.Elk1n] )
    v[123] = p[C.V30] * u[V.pElk1n] / ( p[C.Km30] + u[V.pElk1n] )
    v[124] = p[C.V31] * (u[V.pCREBn] * u[V.pElk1n])^p[C.n31] / ( p[C.Km31]^p[C.n31] + (u[V.pCREBn] * u[V.pElk1n])^p[C.n31] + (u[V.Fn] / p[C.KF31])^p[C.nF31] )
    v[125] = p[C.p32] * u[V.PrecfosmRNAn]
    v[126] = p[C.p33] * u[V.cfosmRNAc]
    v[127] = p[C.p34] * u[V.cfosmRNAc]
    v[128] = p[C.V35] * u[V.ppERKc] * u[V.cFOSc] / ( p[C.Km35] + u[V.cFOSc] )
    v[129] = p[C.V36] * u[V.pRSKc] * u[V.cFOSc] / ( p[C.Km36] + u[V.cFOSc] )
    v[130] = p[C.V37] * u[V.pcFOSc] / ( p[C.Km37] + u[V.pcFOSc] )
    v[131] = p[C.p38] * u[V.cFOSc]
    v[132] = p[C.p39] * u[V.pcFOSc]
    v[133] = p[C.KimFOS] * u[V.cFOSc] - p[C.KexFOS] * (p[C.Vn]/p[C.Vc]) * u[V.cFOSn]
    v[134] = p[C.KimpcFOS] * u[V.pcFOSc] - p[C.KexpcFOS] * (p[C.Vn]/p[C.Vc]) * u[V.pcFOSn]
    v[135] = p[C.V42] * u[V.ppERKn] * u[V.cFOSn] / ( p[C.Km42] + u[V.cFOSn] )
    v[136] = p[C.V43] * u[V.pRSKn] * u[V.cFOSn] / ( p[C.Km43] + u[V.cFOSn] )
    v[137] = p[C.V44] * u[V.pcFOSn] / ( p[C.Km44] + u[V.pcFOSn] )
    v[138] = p[C.p45] * u[V.cFOSn]
    v[139] = p[C.p46] * u[V.pcFOSn]
    v[140] = p[C.p47] * u[V.DUSPn] * u[V.ppERKn] - p[C.m47] * u[V.DUSPn_ppERKn]
    v[141] = p[C.p48] * u[V.DUSPn_ppERKn]
    v[142] = p[C.p49] * u[V.DUSPn] * u[V.pERKn] - p[C.m49] * u[V.DUSPn_pERKn]
    v[143] = p[C.p50] * u[V.DUSPn_pERKn]
    v[144] = p[C.p51] * u[V.DUSPn] * u[V.ERKn] - p[C.m51] * u[V.DUSPn_ERKn]
    v[145] = p[C.p52] * u[V.pDUSPn] * u[V.ppERKn] - p[C.m52] * u[V.pDUSPn_ppERKn]
    v[146] = p[C.p53] * u[V.pDUSPn_ppERKn]
    v[147] = p[C.p54] * u[V.pDUSPn] * u[V.pERKn] - p[C.m54] * u[V.pDUSPn_pERKn]
    v[148] = p[C.p55] * u[V.pDUSPn_pERKn]
    v[149] = p[C.p56] * u[V.pDUSPn] * u[V.ERKn] - p[C.m56] * u[V.pDUSPn_ERKn]
    v[150] = p[C.V57] * u[V.pcFOSn]^p[C.n57] / ( p[C.Km57]^p[C.n57] + u[V.pcFOSn]^p[C.n57] )
    v[151] = p[C.p58] * u[V.PreFmRNAn]
    v[152] = p[C.p59] * u[V.FmRNAc]
    v[153] = p[C.p60] * u[V.FmRNAc]
    v[154] = p[C.p61] * u[V.Fc]
    v[155] = p[C.KimF] * u[V.Fc] - p[C.KexF] * (p[C.Vn]/p[C.Vc]) * u[V.Fn]
    v[156] = p[C.p63] * u[V.Fn]

    du[V.E1] = -v[1] - v[81]
    du[V.E2] = -v[5] - v[6] - v[8] - v[82]
    du[V.E3] = -v[2]
    du[V.E4] = -v[3] - v[83]
    du[V.E_E1] = v[1] - v[4] - v[4] - v[5] - v[61] - v[62] - v[84]
    du[V.H_E3] = v[2] - v[6] - v[7] - v[61]
    du[V.H_E4] = v[3] - v[7] - v[8] - v[9] - v[9] - v[62] - v[85]
    du[V.E11] = v[4] - v[10]
    du[V.E12] = v[5] - v[11]
    du[V.E23] = v[6] - v[12]
    du[V.E34] = v[7] - v[13]
    du[V.E24] = v[8] - v[14]
    du[V.E44] = v[9] - v[15]
    du[V.E11P] = v[10] - v[16] - v[17] - v[18] - v[73] - v[89]
    du[V.E12P] = v[11] - v[19] - v[20] - v[21] - v[74]
    du[V.E23P] = v[12] - v[22] - v[23] - v[24] - v[25] - v[75]
    du[V.E34P] = v[13] - v[26] - v[27] - v[28] - v[29] - v[76]
    du[V.E24P] = v[14] - v[30] - v[31] - v[32] - v[33] - v[77]
    du[V.E44P] = v[15] - v[34] - v[35] - v[36] - v[37] - v[78]
    du[V.G] = -v[16] - v[19] - v[22] - v[26] - v[30] - v[34] - v[42] - v[57] - v[65] - v[69] + v[90]
    du[V.S] = -v[17] - v[20] - v[23] - v[27] - v[31] - v[35] - v[43] - v[66] - v[70] + v[91]
    du[V.I] = -v[24] - v[28] - v[32] - v[36] - v[44] - v[67] - v[71]
    du[V.R] = -v[18] - v[21] - v[25] - v[29] - v[33] - v[37] - v[45] - v[68] - v[72] + v[92]
    du[V.O] = -v[40] - v[54] - v[58] - v[60]
    du[V.A] = -v[41] - v[46] - v[55] - v[59]
    du[V.E11G] = v[16] - v[90]
    du[V.E11S] = v[17] - v[91]
    du[V.E11R] = v[18] - v[92]
    du[V.E12G] = v[19]
    du[V.E12S] = v[20]
    du[V.E12R] = v[21]
    du[V.E23G] = v[22]
    du[V.E23S] = v[23]
    du[V.E23I] = v[24]
    du[V.E23R] = v[25]
    du[V.E34G] = v[26]
    du[V.E34S] = v[27]
    du[V.E34I] = v[28]
    du[V.E34R] = v[29]
    du[V.E24G] = v[30]
    du[V.E24S] = v[31]
    du[V.E24I] = v[32]
    du[V.E24R] = v[33]
    du[V.E44G] = v[34]
    du[V.E44S] = v[35]
    du[V.E44I] = v[36]
    du[V.E44R] = v[37]
    du[V.sigmaG] = v[16] + v[19] + v[22] + v[26] + v[30] + v[34] - v[40] - v[41] + v[42] + v[65] + v[69] - v[90]
    du[V.sigmaS] = v[17] + v[20] + v[23] + v[27] + v[31] + v[35] - v[38] + v[43] + v[66] + v[70] - v[91]
    du[V.sigmaI] = v[24] + v[28] + v[32] + v[36] + v[44] + v[67] + v[71]
    du[V.sigmaR] = v[18] + v[21] + v[25] + v[29] + v[33] + v[37] + v[45] - v[50] + v[68] + v[72] - v[92]
    du[V.sigmaA] = -v[39] + v[41] + v[46] + v[59]
    du[V.sigmaSP] = v[38] - v[42]
    du[V.sigmaAP] = v[39] - v[43] - v[44] - v[45] - v[88]
    du[V.sigmaG_O] = v[40] - v[59]
    du[V.sigmaG_A] = v[41] - v[60]
    du[V.sigmaSP_G] = v[42]
    du[V.sigmaAP_S] = v[43]
    du[V.sigmaAP_I] = v[44]
    du[V.sigmaAP_R] = v[45]
    du[V.P3_A] = v[46] - v[57]
    du[V.P2] = -v[48]
    du[V.P3] = -v[46] + v[48]
    du[V.Akt] = -v[47]
    du[V.RsD] = -v[49]
    du[V.RsT] = v[49]
    du[V.sigmaRP] = v[50]
    du[V.Raf] = -v[51]
    du[V.Rafstar] = v[51]
    du[V.MEK] = -v[52]
    du[V.ppMEKc] = v[52]
    du[V.OP] = v[54]
    du[V.AP] = v[55]
    du[V.A_sigmaG_O] = v[59] + v[60]
    du[V.sigmaA_G] = v[57] - v[58]
    du[V.sigmaA_G_O] = v[58]
    du[V.sigmaO] = v[40] + v[58] + v[60]
    du[V.E13] = v[61] - v[63]
    du[V.E14] = v[62] - v[64]
    du[V.E13P] = v[63] - v[65] - v[66] - v[67] - v[68] - v[79]
    du[V.E14P] = v[64] - v[69] - v[70] - v[71] - v[72] - v[80]
    du[V.E13G] = v[65]
    du[V.E13S] = v[66]
    du[V.E13I] = v[67]
    du[V.E13R] = v[68]
    du[V.E14G] = v[69]
    du[V.E14S] = v[70]
    du[V.E14I] = v[71]
    du[V.E14R] = v[72]
    du[V.T] = -v[73] - v[74] - v[75] - v[76] - v[77] - v[78] - v[79] - v[80] - v[88] + v[93]
    du[V.E11T] = v[73] - v[93]
    du[V.E12T] = v[74]
    du[V.E23T] = v[75]
    du[V.E34T] = v[76]
    du[V.E24T] = v[77]
    du[V.E44T] = v[78]
    du[V.E13T] = v[79]
    du[V.E14T] = v[80]
    du[V.sigmaT] = v[73] + v[74] + v[75] + v[76] + v[77] + v[78] + v[79] + v[80] + v[88] - v[93]
    du[V.E1_PT] = v[81] - v[86]
    du[V.E2_PT] = v[82]
    du[V.E4_PT] = v[83] - v[87]
    du[V.E_E1_PT] = v[84] + v[86]
    du[V.H_E4_PT] = v[85] + v[87]
    du[V.Aktstar] = v[47]
    du[V.sigmaAP_T] = v[88]
    du[V.E] = (-v[1] - v[86])/VeVc
    du[V.H] = (-v[2] - v[3] - v[87])/VeVc
    du[V.fint] = p[C.a98]*(-u[V.fint] + p[C.b98])

    du[V.CREBn] = -v[120] + v[121]
    du[V.pCREBn] = v[120] - v[121]
    du[V.ERKc] = -v[94] + v[96] - v[100]
    du[V.ERKn] = v[98] + v[100]*(p[C.Vc]/p[C.Vn]) + v[143] -v[144] + v[148] -v[149]
    du[V.pERKc] = v[94] - v[95] -v[96] +v[97]-v[101]
    du[V.pERKn] = -v[98] + v[99] + v[101]*(p[C.Vc]/p[C.Vn]) + v[141] - v[142] + v[146] - v[147]
    du[V.ppERKc] =  v[95] - v[97] - v[102]
    du[V.ppERKn] = -v[99] + v[102]*(p[C.Vc]/p[C.Vn]) - v[140] - v[145]
    du[V.Elk1n] = -v[122] + v[123]
    du[V.pElk1n] = v[122] - v[123]
    du[V.cFOSc] = v[127] - v[128] - v[129] + v[130] - v[131] - v[133]
    du[V.cFOSn] = v[133]*(p[C.Vc]/p[C.Vn]) - v[135] - v[136] + v[137] - v[138]
    du[V.pcFOSc] = v[128] + v[129] - v[130] - v[132] - v[134]
    du[V.pcFOSn] = v[134]*(p[C.Vc]/p[C.Vn]) + v[135] + v[136] - v[137] - v[139]
    du[V.DUSPc] = v[106] - v[107] + v[108] - v[109] - v[111]
    du[V.DUSPn] = v[111]*(p[C.Vc]/p[C.Vn]) - v[113] + v[114] - v[115] - v[140] + v[141] - v[142] + v[143] - v[144]
    du[V.pDUSPc] = v[107] - v[108] - v[110] - v[112]
    du[V.pDUSPn] = v[112]*(p[C.Vc]/p[C.Vn]) + v[113] - v[114] - v[116] - v[145] + v[146] - v[147] + v[148] - v[149]
    du[V.DUSPn_ERKn] = v[144]
    du[V.DUSPn_pERKn] = v[142] - v[143]
    du[V.DUSPn_ppERKn] = v[140] - v[141]
    du[V.pDUSPn_ERKn] = v[149]
    du[V.pDUSPn_pERKn] = v[147] - v[148]
    du[V.pDUSPn_ppERKn] = v[145] - v[146]
    du[V.RSKc] = -v[117] + v[118]
    du[V.pRSKc] = v[117] - v[118] - v[119]
    du[V.pRSKn] = v[119]*(p[C.Vc]/p[C.Vn])
    du[V.PrecfosmRNAn] = v[124] - v[125]
    du[V.PreduspmRNAn] = v[103] - v[104]
    du[V.cfosmRNAc] = v[125]*(p[C.Vn]/p[C.Vc]) - v[126]
    du[V.duspmRNAc] = v[104]*(p[C.Vn]/p[C.Vc]) - v[105]
    du[V.Fc] = v[153] - v[154] - v[155]
    du[V.Fn] = v[155]*(p[C.Vc]/p[C.Vn]) - v[156]
    du[V.FmRNAc] = v[151]*(p[C.Vn]/p[C.Vc]) - v[152]
    du[V.PreFmRNAn] = v[150] - v[151]

end

function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.n_parameters)

    # Birtwistle et al. MSB (2007)
    p[C.VmaxPY] = 223.8776
    p[C.KmPY] = 486.1398
    p[C.kdeg] = 0.0259
    p[C.kf47] = 24.6048
    p[C.Vmaxr47] = 590.5058
    p[C.Kmf47] = 698.6022
    p[C.Kmr47] = 483.8622
    p[C.kf48] = 16.833
    p[C.Kmf48] = 715.5688
    p[C.Kmr48] = 324.9294
    p[C.PTEN] = 1 #  693.5786
    p[C.kf49] = 44.3501
    p[C.kr49] = 552.6746
    p[C.Kmf49] = 343.2483
    p[C.Kmr49] = 753.1667
    p[C.Kmr49b] = 381.2208
    p[C.kr49b] = 640.8212
    p[C.kf51] = 3.6515
    p[C.Vmaxr51] = 16.737
    p[C.Kmf51] = 599.7076
    p[C.Kmr51] = 346.4779
    p[C.Kmrb51] = 988.4496
    p[C.kf52] = 0.7745
    p[C.Vmaxr52] = 199.2773
    p[C.Kmf52] = 545.4408
    p[C.Kmr52] = 675.2994
    p[C.kf54] = 0.0538
    p[C.Vmaxr54] = 588.2671
    p[C.Kmf54] = 457.9645
    p[C.Kmr54] = 336.183
    p[C.kf55] = 0.2256
    p[C.Vmaxr55] = 646.9003
    p[C.Kmf55] = 460.9446
    p[C.Kmr55] = 643.9247
    p[C.kf38] = 279.9929
    p[C.kf39] = 385.7428
    p[C.kf50] = 389.1061
    p[C.a98] = 0.0849
    p[C.b98] = 0.1833
    p[C.koff46] = 0.5194
    p[C.EGF_off] = 0.0175
    p[C.HRGoff_3] = 9.0E-4
    p[C.HRGoff_4] = 0.0973
    p[C.koff4] = 0.1717
    p[C.koff5] = 4.3985
    p[C.koff6] = 2.6619
    p[C.koff7] = 8.0557
    p[C.koff8] = 9.1034
    p[C.koff9] = 5.5425
    p[C.koff57] = 0.4526
    p[C.koff16] = 0.5737
    p[C.koff17] = 4.6874
    p[C.koff18] = 2.2768
    p[C.koff19] = 2.3361
    p[C.koff20] = 0.6761
    p[C.koff21] = 4.7291
    p[C.koff22] = 3.6962
    p[C.koff23] = 2.3619
    p[C.koff24] = 4.4226
    p[C.koff25] = 2.225
    p[C.koff26] = 0.0103
    p[C.koff27] = 1.8922
    p[C.koff28] = 4.6432
    p[C.koff29] = 2.0148
    p[C.koff30] = 4.9936
    p[C.koff31] = 1.2204
    p[C.koff32] = 3.8752
    p[C.koff33] = 1.2817
    p[C.koff34] = 3.2036
    p[C.koff35] = 1.8696
    p[C.koff36] = 1.2567
    p[C.koff37] = 0.4059
    p[C.koff65] = 0.1185
    p[C.koff66] = 2.2988
    p[C.koff67] = 1.6142
    p[C.koff40] = 3.1051
    p[C.koff41] = 7.0487
    p[C.koff42] = 3.5195
    p[C.koff43] = 0.5441
    p[C.koff44] = 0.4265
    p[C.koff45] = 3.9967
    p[C.koff58] = 6.3059
    p[C.koff59] = 9.172
    p[C.koff68] = 2.8871
    p[C.kPTP10] = 29.8531
    p[C.kPTP11] = 78.204
    p[C.kPTP12] = 11.4211
    p[C.kPTP13] = 55.2104
    p[C.kPTP14] = 57.7506
    p[C.kPTP15] = 60.2628
    p[C.kPTP63] = 7.4766
    p[C.kPTP64] = 48.6335
    p[C.koff73] = 2.2768
    p[C.koff19] = 2.3361
    p[C.koff20] = 0.6761
    p[C.koff21]  = 3.0048
    p[C.koff74] = 1.2496
    p[C.koff75] = 1.4323
    p[C.koff76] = 2.1542
    p[C.koff77] = 1.2237
    p[C.koff78] = 0.2007
    p[C.koff79] = 1.1852
    p[C.koff80] = 2.9373
    p[C.kPTP38] = 83.4465
    p[C.kPTP39] = 79.6132
    p[C.koff88] = 3.9255
    p[C.kPTP50] = 96.5716
    p[C.kf81] = 1.361
    p[C.Vmaxr81] = 242.6034
    p[C.Kmf81] = 485.2626
    p[C.Kmr81] = 323.4012
    p[C.kf82] = 6.9987
    p[C.Vmaxr82] = 398.1931
    p[C.Kmf82] = 781.4374
    p[C.Kmr82] = 595.8395
    p[C.kf83] = 1.7628
    p[C.Vmaxr83] = 534.0531
    p[C.Kmf83] = 609.4766
    p[C.Kmr83] = 653.5184
    p[C.kf84] = 4.6894
    p[C.Vmaxr84] = 634.1626
    p[C.Kmf84] = 622.3847
    p[C.Kmr84] = 258.4637
    p[C.kf85] = 6.7591
    p[C.Vmaxr85] = 369.2261
    p[C.Kmf85] = 179.6486
    p[C.Kmr85] = 290.7667
    p[C.kcon49] = 9.9783
    p[C.kon1] = 1.0086E-4
    p[C.kon86] = 0.0038
    p[C.kon2] = 0.0059
    p[C.kon3] = 0.0334
    p[C.kon87] = 8.0E-4
    p[C.kon4] = 0.5005
    p[C.kon5] = 2.5427
    p[C.kon6] = 0.2283
    p[C.kon7] = 1.0606
    p[C.kon8] = 1.0259
    p[C.kon9] = 2.2868
    p[C.kon57] = 0.0039
    p[C.kf10] = 0.6496
    p[C.kf11] = 0.3721
    p[C.kf12] = 1.8012
    p[C.kf13] = 0.8875
    p[C.kf14] = 6.1726
    p[C.kf15] = 1.3565
    p[C.kf63] = 0.9297
    p[C.kf64] = 1.2083
    p[C.kon16] = 0.0097
    p[C.kon17] = 0.0049
    p[C.kon18] = 0.0117
    p[C.kon73] = 0.0116
    p[C.kon19] = 0.0896
    p[C.kon20] = 0.0478
    p[C.kon21] = 0.0114
    p[C.kon74] = 0.0133
    p[C.kon22] = 7.0E-4
    p[C.kon23] = 0.0138
    p[C.kon24] = 0.005
    p[C.kon25] = 0.0995
    p[C.kon75] = 0.0137
    p[C.kon26] = 0.0355
    p[C.kon27] = 0.0201
    p[C.kon28] = 0.0074
    p[C.kon29] = 0.0346
    p[C.kon76] = 0.0053
    p[C.kon30] = 0.002
    p[C.kon31] = 0.0032
    p[C.kon32] = 9.0E-4
    p[C.kon33] = 0.0335
    p[C.kon77] = 0.0101
    p[C.kon34] = 1.0E-4
    p[C.kon35] = 0.0602
    p[C.kon36] = 0.0043
    p[C.kon37] = 0.0791
    p[C.kon78] = 0.0076
    p[C.kon79] = 0.0078
    p[C.kon65] = 0.0123
    p[C.kon66] = 1.9264E-4
    p[C.kon67] = 6.6667E-5
    p[C.kon80] = 2.0E-4
    p[C.kon40] = 0.0191
    p[C.kon41] = 0.0051
    p[C.kon42] = 0.0023
    p[C.kon43] = 0.0127
    p[C.kon44] = 0.0122
    p[C.kon45] = 0.0028
    p[C.kon88] = 0.0108
    p[C.kon46] = 0.0148
    p[C.kon58] = 0.0215
    p[C.kon59] = 0.0077
    p[C.kon60] = 1.1994E-4
    p[C.koff60] = 4.9981
    p[C.koff61] = 5.229
    p[C.kon61] = 0.8048
    p[C.kon62] = 1.782
    p[C.koff62] = 5.5142
    p[C.kon68] = 0.0045
    p[C.kon69] = 0.0084
    p[C.koff69] = 3.97
    p[C.kon70] = 0.0116
    p[C.koff70] = 2.6069
    p[C.kon71] = 0.0078
    p[C.koff71] = 2.2988
    p[C.kon72] = 0.0355
    p[C.koff72] = 0.907

    # Nakakuki et al., Cell (2010)
    p[C.V1] = 0.34284837
    p[C.Km1] = 307.0415253
    p[C.V2] = 2.20e-01
    p[C.Km2] = 3.50e+02
    p[C.V3] = 7.20e-01
    p[C.Km3] = 1.60e+02
    p[C.V4] = 6.48e-01
    p[C.Km4] = 6.00e+01
    p[C.V5] = 19.49872346
    p[C.Km5]  =29.94073716
    p[C.V6] = p[C.V5]
    p[C.Km6] = p[C.Km5]
    p[C.KimERK] = 1.20e-02
    p[C.KexERK] = 1.80e-02
    p[C.KimpERK] = 1.20e-02
    p[C.KexpERK] = 1.80e-02
    p[C.KimppERK] = 1.10e-02
    p[C.KexppERK] = 1.30e-02
    p[C.V10] = 29.24109258
    p[C.Km10] = 169.0473748
    p[C.n10] = 1#3.970849295
    p[C.p11] = 0.000126129
    p[C.p12] = 0.007875765
    p[C.p13] = 0.001245747
    p[C.V14] = 5.636949216
    p[C.Km14] = 34180.48
    p[C.V15] = 2.992346912
    p[C.Km15] = 0.001172165
    p[C.p16] = 2.57e-04
    p[C.p17] = 9.63e-05
    p[C.KimDUSP] = 0.024269764
    p[C.KexDUSP] = 0.070467899
    p[C.KimpDUSP] = p[C.KimDUSP]
    p[C.KexpDUSP] = p[C.KexDUSP]
    p[C.V20] = 0.157678678
    p[C.Km20] = 735598.6967
    p[C.V21] = 0.005648117
    p[C.Km21] = 387.8377182
    p[C.p22] = 2.57e-04
    p[C.p23] = 9.63e-05
    p[C.V24] = 0.550346114
    p[C.Km24] = 29516.06587
    p[C.V25] = 10.09063736
    p[C.Km25] = 0.913939859
    p[C.KimRSK] = 0.025925065
    p[C.KexRSK] = 0.129803956
    p[C.V27] = 19.23118154
    p[C.Km27] = 441.5834425
    p[C.V28] = 6.574759504
    p[C.Km28] = 14.99180922
    p[C.V29] = 0.518529841
    p[C.Km29] = 21312.69109
    p[C.V30] = 13.79479021
    p[C.Km30] = 15.04396629
    p[C.V31] = 0.655214248
    p[C.Km31] = 185.9760682
    p[C.n31] = 1#1.988003164
    p[C.KF31] = 0.013844393
    p[C.nF31] = 1#2.800340453
    p[C.p32] = 0.003284434
    p[C.p33] = 0.000601234
    p[C.p34] = 7.65E-05
    p[C.V35] = 8.907637012
    p[C.Km35] = 8562.744184
    p[C.V36] = 0.000597315
    p[C.Km36] = 528.552427
    p[C.V37] = 1.745848179
    p[C.Km37] = 0.070379236
    p[C.p38] = 2.57e-04
    p[C.p39] = 9.63e-05
    p[C.KimFOS] = 0.54528521
    p[C.KexFOS] = 0.133249762
    p[C.KimpcFOS] = p[C.KimFOS]
    p[C.KexpcFOS] = p[C.KexFOS]
    p[C.V42] = 0.909968714
    p[C.Km42] = 3992.061328
    p[C.V43] = 0.076717457
    p[C.Km43] = 1157.116021
    p[C.V44] = 0.078344305
    p[C.Km44] = 0.051168202
    p[C.p45] = 2.57e-04
    p[C.p46] = 9.63e-05
    p[C.p47] = 0.001670815
    p[C.m47] = 15.80783969
    p[C.p48] = 0.686020478
    p[C.p49] = 0.314470502
    p[C.m49] = 2.335459127
    p[C.p50] = 26.59483436
    p[C.p51] = 0.01646825
    p[C.m51] = 9.544308421
    p[C.p52] = p[C.p47]
    p[C.m52] = p[C.m47]
    p[C.p53] = p[C.p48]
    p[C.p54] = p[C.p49]
    p[C.m54] = p[C.m49]
    p[C.p55] = p[C.p50]
    p[C.p56] = p[C.p51]
    p[C.m56] = p[C.m51]
    p[C.V57] = 1.026834758
    p[C.Km57] = 0.637490056
    p[C.n57] = 1#3.584464176
    p[C.p58] = 0.000270488
    p[C.p59] = 0.001443889
    p[C.p60] = 0.002448164
    p[C.p61] = 3.50E-05
    p[C.KimF] = 0.019898797
    p[C.KexF] = 0.396950616
    p[C.p63] = 4.13E-05

    p[C.Vn] = 0.22
    p[C.Vc] = 0.94

    p[C.scale_E1] = 10.0
    p[C.scale_E2] = 10.0
    p[C.scale_E3] = 10.0
    p[C.scale_E4] = 10.0
    p[C.scale_G] = 10.0
    #p[C.scale_S] = 10.0
    p[C.a_SHC1] = 10.0
    p[C.a_SHC2] = 10.0
    p[C.a_SHC3] = 10.0
    p[C.a_SHC4] = 10.0
    #p[C.scale_I] = 10.0
    p[C.a_PIK3CA] = 10.0
    p[C.a_PIK3CB] = 10.0
    p[C.a_PIK3CD] = 10.0
    p[C.a_PIK3CG] = 10.0
    # ver.5
    p[C.scale_PTEN] = 10.0
    #p[C.scale_R] = 10.0
    p[C.a_RASA1] = 10.0
    p[C.a_RASA2] = 10.0
    p[C.a_RASA3] = 10.0
    #p[C.scale_O] = 10.0
    p[C.a_SOS1] = 10.0
    p[C.a_SOS2] = 10.0
    p[C.scale_A] = 10.0
    #p[C.scale_Akt] = 10.0
    p[C.a_AKT1] = 10.0
    p[C.a_AKT2] = 10.0
    #p[C.scale_RsD] = 10.0
    p[C.a_HRAS] = 10.0
    p[C.a_KRAS] = 10.0
    p[C.a_NRAS] = 10.0
    #p[C.scale_Raf] = 10.0
    p[C.a_ARAF] = 10.0
    p[C.a_BRAF] = 10.0
    p[C.a_RAF1] = 10.0
    #p[C.scale_MEK] = 10.0
    p[C.a_MAP2K1] = 10.0
    p[C.a_MAP2K2] = 10.0
    p[C.scale_T] = 10.0
    p[C.scale_CREBn] = 10.0
    #p[C.scale_ERKc] = 10.0
    p[C.a_MAPK1] = 10.0
    p[C.a_MAPK3] = 10.0
    p[C.scale_Elk1n] = 10.0
    #p[C.scale_RSKc] = 10.0
    p[C.a_RPS6KA1] = 10.0
    p[C.a_RPS6KA2] = 10.0
    p[C.a_RPS6KA3] = 10.0
    p[C.a_RPS6KA6] = 10.0

    return p
end

function initial_values()::Vector{Float64}
    u0::Vector{Float64} = zeros(V.n_species)

    u0[V.E1] = 1
    u0[V.E2] = 1
    u0[V.E3] = 1
    u0[V.E4] = 1
    u0[V.G] = 1
    u0[V.S] = 1
    u0[V.I] = 1
    u0[V.R] = 1
    u0[V.O] = 1
    u0[V.A] = 1
    u0[V.P2] = 1
    u0[V.Akt] = 1
    u0[V.RsD] = 1
    u0[V.Raf] = 1
    u0[V.MEK] = 1
    u0[V.T] = 1
    u0[V.CREBn] = 1
    u0[V.ERKc] = 1
    u0[V.Elk1n] = 1
    u0[V.RSKc] = 1

    return u0
end


text:
using Example
hello("it looks like you were dreaming.")

ext: jl

max_stars_repo_name: luca-mm/Bachelor-Thesis

max_issues_repo_name: luca-mm/Bachelor-Thesis

max_forks_repo_name: luca-mm/Bachelor-Thesis


text:
function parsehtml(input::String; strict=false)
    result_ptr = ccall((:gumbo_parse,libgumbo),
                       Ptr{CGumbo.Output},
                       (Ptr{Uint8},),
                       input)
    goutput::CGumbo.Output = unsafe_load(result_ptr)

    if strict && goutput.errors.length > 0
        throw(InvalidHTMLException("input html was invalid"))
    end
    doc = document_from_gumbo(goutput)
    gumbo_dl = dlopen(libgumbo)
    default_options = dlsym(gumbo_dl, :kGumboDefaultOptions)
    ccall((:gumbo_destroy_output,libgumbo),
          Void,
          (Ptr{Void}, Ptr{CGumbo.Output}),
          default_options, result_ptr)
    return doc
end

# turn a gumbo vector into a Julia vector
# of Ptr{T} where T is the struct contained
# by the gumbo vector
gvector_to_jl(T,gv::CGumbo.Vector) = pointer_to_array(convert(Ptr{Ptr{T}},gv.data),
                                                     (int(gv.length),))

# convert a vector of pointers to GumboAttributes to
# a Dict String => String
function attributes(av::Vector{Ptr{CGumbo.Attribute}})
    result = Dict{String,String}()
    for ptr in av
        ga::CGumbo.Attribute = unsafe_load(ptr)
        result[bytestring(ga.name)] = bytestring(ga.value)
    end
    return result
end

function element(parent::HTMLNode, ge::CGumbo.Element)
    tag = CGumbo.TAGS[ge.tag+1]  # +1 is for 1-based julia indexing
    # TODO when tag is unknown get out actual tag ??
    attrs = attributes(gvector_to_jl(CGumbo.Attribute,ge.attributes))
    children = HTMLNode[]
    res = HTMLElement{tag}(children, parent, attrs)
    for nodeptr in gvector_to_jl(CGumbo.Node,ge.children)
        node::CGumbo.Node = unsafe_load(nodeptr)
        if node.gntype == CGumbo.ELEMENT
            push!(children, element(res, node.v))  # already a GumboElement
        elseif node.gntype == CGumbo.TEXT
            push!(children,HTMLText(res,
                                    bytestring(reinterpret(CGumbo.Text,node.v).text)))
        end
        # TODO handle CDATA, comments, etc.
    end
    res
end

element(ge::CGumbo.Element) = element(NullNode(), ge)

# transform gumbo output into Julia data
function document_from_gumbo(goutput::CGumbo.Output)
    # TODO convert some of these typeasserts to better error messages?
    gnode::CGumbo.Node = unsafe_load(goutput.document)
    if gnode.gntype != CGumbo.DOCUMENT
        error("Document appears to have incorrent GumboNodeType")
    end
    gdoc = reinterpret(CGumbo.Document,gnode.v)
    doctype = bytestring(gdoc.name)
    groot::CGumbo.Node = unsafe_load(goutput.root)
    if groot.gntype != CGumbo.ELEMENT  # should always be true
        error("root appears to have wrong GumboNodeType")
    end
    grootnode::CGumbo.Node = unsafe_load(goutput.root)
    root = element(grootnode.v)  # already an element
    HTMLDocument(doctype, root)
end


text:
using Random

# TODO: use Macros and Symbols

function divide_on_feature(x, idx, threshold)
    """
    Divide dataset based on if sample value on feature index is larger than
    the given threshold
    """

    x1, x2 = [], []

    for i in 1:size(x, 1)
        if x[i, idx] >= threshold
            push!(x1, x[i,:]')
        else
            push!(x2, x[i,:]')
        end
    end

    return vcat(x1...), vcat(x2...)

end

function sample(data::Array{T,2}, n::Int64; seed=123::Int64) where T <: Any
    # TODO: accept seed!
    n_samples = size(data, 1)

    # Shuffle data
    idx = shuffle(MersenneTwister(seed), 1:n_samples)
    data = data[idx, :]

    data[1:n, :]
end

function split_data(data, fraction)
    n_samples = size(data, 1)

    # Suffle data
    shuffle!(data)

    s = round(Integer, n_samples * fraction)

    return data[1:s, :], data[s:end, :]
end

# From https://bit.ly/2QnEacc

function entropy(y::AbstractArray)
    """ Calculate the entropy of label array y """
    unique_labels = unique(y)
    e = 0
    n = length(y)

    for label in unique_labels
        c = length(filter(x-> x==label, y))
        p = c / n
        e += -p * log2(p)
    end

    return e
end

function gini_index(y::AbstractArray)
    """ Calculate the Gini index of label array y """
    unique_labels = unique(y)
    e = 0
    n = length(y)

    for label in unique_labels
        c = length(filter(x-> x==label, y))
        p = c / n
        e += p^2
    end

    return 1 - e
end

function classification_error(y::AbstractArray)
    """ Calculate the Classification Error of label array y """
    unique_labels = unique(y)
    pp = []
    n = length(y)

    for label in unique_labels
        c = length(filter(x-> x==label, y))
        p = c / n
        push!(pp, p)
    end

    return 1 - max(pp)
end

function modes(values)
    # From https://rosettacode.org/wiki/Averages/Mode#Julia
    dict = Dict() # Values => Number of repetitions
    modesArray = typeof(values[1])[] # Array of the modes so far
    max = 0 # Max of repetitions so far

    for v in values
        # Add one to the dict[v] entry (create one if none)
        if v in keys(dict)
            dict[v] += 1
        else
            dict[v] = 1
        end

        # Update modesArray if the number of repetitions
        # of v reaches or surpasses the max value
        if dict[v] >= max
            if dict[v] > max
                empty!(modesArray)
                max += 1
            end
            append!(modesArray, [v])
        end
    end

    return modesArray
end

function count_distinct(seq)
    dict = Dict{Any, Int64}();
    for item in seq
        if !haskey(dict, item)
            dict[item] = 0;
        end
        dict[item] += 1
    end

    dict
end


text:
using StanSample, Test

TestDir = @__DIR__

tmpdir = mktempdir()

original_stanrun_tests = [
  "original_stanrun_test/test_stanrun.jl",
]

for test in original_stanrun_tests
  println("\nTesting: $test.")
  include(test)
end
println()

basic_run_tests = [
  "test_basic_runs/test_bernoulli_dict.jl",
  "test_basic_runs/test_bernoulli_nt.jl",
  "test_basic_runs/test_bernoulli_array_dict_1.jl",
  "test_basic_runs/test_bernoulli_array_dict_2.jl",
  "test_basic_runs/test_parse_interpolate.jl"
]

@testset "Bernoulli basic run tests" begin
  for test in basic_run_tests
    println("\nTesting: $test.")
    include(joinpath(TestDir, test))
    @test sdf[:theta, :mean][1] ≈ 0.34 atol=0.1
  end
  println()
end

sample_settings_tests = [
  "test_sample_settings/test_bernoulli.jl"
]

@testset "Bernoulli Sample() settings tests" begin
  for test in sample_settings_tests
    println("\nTesting: $test.")
    include(joinpath(TestDir, test))
    @test stanmodel.method.adapt.delta ≈ 0.85 atol=0.01
  end
  println()
end

generate_quantities_tests = [
  "test_generate_quantities/test_generate_quantities.jl"
]

@testset "Generate_quantities tests" begin
  for test in generate_quantities_tests
    println("\nTesting: $test.")
    include(joinpath(TestDir, test))
  end
  println()
end


text:
module ReactiveMPProdFinalTest

using Test
using ReactiveMP
using Random
using LinearAlgebra

import ReactiveMP: getdist

@testset "ProdFinal" begin

    uni_distributions = [
        Gamma(1.0, 2.0),
        NormalMeanVariance(-10.0, 3.0)
    ]

    mv_distributions = [
        MvNormalMeanCovariance(ones(3)),
        MvNormalMeanPrecision(3ones(3))
    ]

    mxv_distributions = [
        Wishart(3, Matrix(Diagonal(ones(3)))),
        Wishart(4, Matrix(Diagonal(ones(4))))
    ]

    sets = (uni_distributions, mv_distributions, mxv_distributions)

    for dset in sets
        for i in 1:length(dset), j in 1:length(dset)
            left = dset[i]
            right = dset[j]

            @test (@inferred getdist(prod(ProdAnalytical(), ProdFinal(left), right))) === left
            @test (@inferred getdist(prod(ProdAnalytical(), left, ProdFinal(right)))) === right
            @test_throws ErrorException prod(ProdAnalytical(), ProdFinal(left), ProdFinal(right))

        end
    end

    # Check errors for different variate_forms
    for i in 1:length(sets), j in 1:length(sets)
        if i !== j
            dset1 = sets[i]
            dset2 = sets[j]

            for dist1 in dset1, dist2 in dset2
                @test_throws ErrorException prod(ProdAnalytical(), ProdFinal(dist1), dist2)
                @test_throws ErrorException prod(ProdAnalytical(), dist1, ProdFinal(dist2))
            end
        end
    end

end

end


text:
# Algorithm 2.6
# Newton's method for solving f(x) = x given an approximation

function bisect_approx(f, p0, tolerance, maxIter)
	for i ∈ 1:maxIter
		p1 = f(p0)
		p2 = f(p1)
		p = p0 - ((p1 - p0)^2)/(p2 - 2p1 + p0)
		if abs(p - p0) < tolerance
			return p
		end
		p0 = p
	end
	throw(ArgumentError("Method failed"))
end

println(bisect_approx(x->((10/(4 + x))^(1/2)), 1.5, 0.0001, 5))


text:
# Try to open at least all downloaded datasets
@testset "Reading rasters" begin
    for f in remotefiles
        @testset "Reading $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f))
        end
    end
end

@testset "Reading rasters first band" begin
    for f in remotefiles
        @testset "Reading $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), band=1)
            @test last(size(ga)) == 1
        end
    end
end

@testset "Reading rasters streaming" begin
    for f in remotefiles
        @testset "Reading $f streaming" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), masked=false)
            ga[1,1,1]
        end
            end
end

@testset "Reading rasters streaming first band" begin
    for f in remotefiles
        @testset "Reading $f streaming" begin
            ga = GeoArrays.read(joinpath(testdatadir, f), masked=false, band=1)
            @test last(size(ga)) == 1
            ga[1,1,1]
        end
    end
end

@testset "Reading rasters and writing" begin
    for (i, f) in enumerate(remotefiles)
        @testset "Reading and writing $f" begin
            ga = GeoArrays.read(joinpath(testdatadir, f))
            GeoArrays.write!(joinpath(tempdir(), "test_$i.tif"), ga)
            ga_copy = GeoArrays.read(joinpath(tempdir(), "test_$i.tif"))
            @test ga[1,1,1] === ga_copy[1,1,1]
            @test ga[end,end,end] === ga_copy[end,end,end]
        end
    end
end

@testset "Read second band" begin
    fn = joinpath(testdatadir, remotefiles[end - 1])
    GeoArrays.read(fn, band=2)
end

@testset "Writing rasters" begin

    @testset "Simplest version" begin
        ga = GeoArray(rand(100, 200, 3))
        fn = GeoArrays.write!(joinpath(tempdir(), "test.tif"), ga)
        GeoArrays.read(fn)
    end
    @testset "Simplest version" begin
        ga = GeoArray(rand(100, 200, 3))
        fn = GeoArrays.write!(joinpath(tempdir(), "test.img"), ga)
        GeoArrays.read(fn)
    end
    @testset "Nodata" begin
        ga = GeoArray(Array{Union{Missing,Int32}}(rand(1:10, 100, 200, 3)))
        fn = GeoArrays.write!(joinpath(tempdir(), "test_nodata.tif"), ga, 1)
        GeoArrays.read(fn)
    end
    @testset "Nodata" begin
        ga = GeoArray(Array{Union{Missing,Int32}}(rand(1:10, 100, 200, 3)))
        fn = GeoArrays.write!(joinpath(tempdir(), "test_nodata.img"), ga, 1)
        GeoArrays.read(fn)
    end

end


text:
using BinaryBuilder

name = "MUMPS"
version = v"5.2.1"

# Collection of sources required to build PARMETIS.
# The patch prevents building the source of METIS that ships with PARMETIS;
# we rely on METIS_jll instead.
sources = [
  ArchiveSource("http://mumps.enseeiht.fr/MUMPS_5.2.1.tar.gz",
                "d988fc34dfc8f5eee0533e361052a972aa69cc39ab193e7f987178d24981744a"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir -p ${libdir}
cd $WORKSPACE/srcdir/MUMPS_5.2.1

OPENBLAS=(-lopenblas)
FFLAGS=()
CFLAGS=()
if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
  FFLAGS+=(-ffixed-line-length-none)  # replacing symbols below sometimes makes lines > 72 chars
  OPENBLAS=(-lopenblas64_)
  if [[ "${target}" == powerpc64le-linux-gnu ]]; then
    OPENBLAS+=(-lgomp)
  fi

  syms=(DGEMM IDAMAX ISAMAX SNRM2 XERBLA ccopy cgemm cgemv cgeru cher clarfg cscal cswap ctrsm ctrsv cungqr cunmqr dcopy dgemm dgemv dger dlamch dlarfg dorgqr dormqr dnrm2 dscal dswap dtrsm dtrsv dznrm2 idamax isamax ilaenv scnrm2 scopy sgemm sgemv sger slamch slarfg sorgqr sormqr snrm2 sscal sswap strsm strsv xerbla zcopy zgemm zgemv zgeru zlarfg zscal zswap ztrsm ztrsv zungqr zunmqr)
  for sym in ${syms[@]}
  do
    FFLAGS+=("-D${sym}=${sym}_64")
    CFLAGS+=("-D${sym}=${sym}_64")
  done
fi

makefile="Makefile.G95.PAR"
cp Make.inc/${makefile} Makefile.inc

make_args+=(OPTF=-O
            CDEFS=-DAdd_
            LMETISDIR=${libdir}
            IMETIS=-I${prefix}/include
            LMETIS='-L$(LMETISDIR) -lparmetis -lmetis'
            ORDERINGSF="-Dpord -Dparmetis"
            CC="mpicc -fPIC"
            FC="mpif90 -fPIC ${FFLAGS[@]}"
            FL="mpif90 -fPIC ${CFLAGS[@]}"
            SCALAP=-lscalapack
            INCPAR=  # Let MPI compilers fill in the blanks
            LIBPAR=-lscalapack
            LIBBLAS=${OPENBLAS})

if [[ "${target}" == *-apple* ]]; then
  make_args+=(RANLIB=echo)
fi

# NB: parallel build fails
make alllib "${make_args[@]}"

# build shared libs
all_load="--whole-archive"
noall_load="--no-whole-archive"
extra=""
if [[ "${target}" == *-apple-* ]]; then
    all_load="-all_load"
    noall_load="-noall_load"
    extra="-Wl,-undefined -Wl,dynamic_lookup -headerpad_max_install_names"
fi

cd lib
libs=(-lparmetis -lmetis -lscalapack ${OPENBLAS})
mpif90 -fPIC -shared -Wl,${all_load} libpord.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o libpord.${dlext}
cp libpord.${dlext} ${libdir}

libs+=(-L${libdir} -lpord)
mpif90 -fPIC -shared -Wl,${all_load} libmumps_common.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o libmumps_common.${dlext}
cp libmumps_common.${dlext} ${libdir}

libs+=(-lmumps_common)
for libname in cmumps dmumps smumps zmumps
do
  mpif90 -fPIC -shared -Wl,${all_load} lib${libname}.a ${libs[@]} -Wl,${noall_load} ${extra[@]} -o lib${libname}.${dlext}
done
cp *.${dlext} ${libdir}
cd ..

cp include/* ${prefix}/include
"""

# OpenMPI and MPICH are not precompiled for Windows
# SCALAPACK doesn't build on PowerPC
platforms = expand_gfortran_versions(filter!(p -> !Sys.iswindows(p) && arch(p) != "powerpc64le", supported_platforms()))

# The products that we will ensure are always built
products = [
    LibraryProduct("libsmumps", :libsmumps),
    LibraryProduct("libdmumps", :libdmumps),
    LibraryProduct("libcmumps", :libcmumps),
    LibraryProduct("libzmumps", :libzmumps),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("MPICH_jll"),
    Dependency("METIS_jll"),
    Dependency("PARMETIS_jll"),
    Dependency("SCALAPACK_jll"),
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
__precompile__()

module SpatialGrids

using StaticArrays

export Raster, rasterize_points
export SparseVoxelGrid, in_cuboid, voxel_center, make_voxel_id

include("raster.jl")
include("sparse_voxels.jl")

end # module


text:
platform = load_platform("$HetaSimulatorDir/test/examples/single_comp", rm_out = false);
model = platform.models[:nameless];
@test test_show(platform)
@test test_show(model)

@test isfile("$HetaSimulatorDir/test/examples/single_comp/_julia/model.jl")

# Wrong input tests
## no saveat or tspan
@test_throws ErrorException("Please, provide either `saveat` or `tspan` value.") Scenario(model)
## wrong observables
@test_throws ErrorException("The following observables have not been found in the model: [:GG]") Scenario(model; observables=[:r1, :GG], tspan=(0,120))

# Simulate default scenario
s = Scenario(model, tspan = (0., 200.), saveat = [0.0, 150., 250.], observables=[:r1]) |> sim
@test test_show(s)
@test typeof(s) <: HetaSimulator.SimResult
@test status(s) == :Success
@test times(s)[end] == 250.
@test vals(s)[1][:r1] == 0.1

# tests related to saveat, tspan behavior
s1 = Scenario(model, tspan = (0., 200.)) |> sim
s2 = Scenario(model, saveat = [150., 200.]) |> sim
@test times(s1)[1] == 0.0
@test times(s1)[end] == 200.0
@test times(s2) == [150., 200.]

# Simulate scenario
scn1 = Scenario(model; parameters = [:k1=>0.02], tspan = (0., 200.), saveat = [0.0, 150., 250.], observables=[:r1])
scn2 = Scenario(model; parameters = [:k1=>0.015], tspan = (0., 200.), observables=[:r1])
@test typeof(scn1) <: Scenario
@test test_show(scn1)
@test saveat(scn1) == [0.0, 150., 250.]
@test saveat(scn2) == Float64[]
@test tspan(scn1) == (0., 250.)
@test tspan(scn2) == (0., 200.)
@test parameters(scn1)[:k1] == 0.02
@test parameters(scn2)[:k1] == 0.015

cs1 = sim(scn1)
cs2 = sim([:one=>scn1,:two=>scn2], parameters_upd=[:k1=>0.03])
@test typeof(cs1) <: HetaSimulator.SimResult
@test typeof(cs2) <: Vector{Pair{Symbol, HetaSimulator.SimResult}}
@test length(parameters(cs1))==0
@test parameters(last(cs2[:one]))[:k1] == 0.03

# Monte-Carlo simulation tests
mciter = 100
mc1 =  mc(model, [:k1=>Normal(0.02,1e-3)], mciter; tspan = (0., 200.), observables=[:r1])
mc2 = mc([:one=>scn1,:two=>scn2], [:k1=>Normal(0.02,1e-3)], mciter)
gsar = gsa(mc1,200)
@test typeof(mc1) <: HetaSimulator.MCResult
@test test_show(mc1)
@test typeof(mc2) <: Vector{Pair{Symbol, HetaSimulator.MCResult}}
@test typeof(mc1[1]) <: HetaSimulator.Simulation
@test typeof(mc2[1]) <: Pair{Symbol, HetaSimulator.MCResult}
@test length(mc1) == mciter
@test typeof(parameters(mc1)) <: Vector
@test times(mc1[1])[end] == 200.
@test keys((parameters(mc1[1]))) == (:k1,)
@test size(pearson(gsar)) == (1,1)
@test size(partial(gsar)) == (1,1)
@test size(standard(gsar)) == (1,1)

# Fitting tests
fscn1 = Scenario(model; parameters = [:k1=>0.02], tspan = (0., 200.), observables=[:A, :B, :r1])
fscn2 = Scenario(model; parameters = [:k1=>0.015], tspan = (0., 200.), observables=[:A, :B, :r1])
data = read_measurements("$HetaSimulatorDir/test/examples/single_comp/single_comp_data.csv")
add_measurements!(fscn1, data; subset = [:scenario => :one])
add_measurements!(fscn2, data; subset = [:scenario => :two])
fres = fit([:one=>fscn1, :two=>fscn2], [:k1=>0.01])

@test typeof(fres) <: HetaSimulator.FitResult
@test test_show(fres)
@test length(measurements(fscn1)) == 24
@test length(measurements(fscn2)) == 24
@test status(fres) == :FTOL_REACHED
@test obj(fres) ≈ 146.056244
@test typeof(optim(fres)) == Vector{Pair{Symbol, Float64}}
@test length(optim(fres)) == 1


text:
@testset "Message dispatcher" begin

    if Sys.iswindows()
        global_socket_name1 = "\\\\.\\pipe\\jsonrpc-testrun1"
    elseif Sys.isunix()
        global_socket_name1 = joinpath(tempdir(), "jsonrpc-testrun1")
    else
        error("Unknown operating system.")
    end

    request1_type = Stratum.RequestType("request1", Foo, String)
    request2_type = Stratum.RequestType("request2", Nothing, String)
    notify1_type = Stratum.NotificationType("notify1", String)

    global g_var = ""

    server_is_up = Base.Condition()

    server_task = @async try
        server = listen(global_socket_name1)
        notify(server_is_up)
        sock = accept(server)
        global conn = Stratum.StratumEndpoint(sock, sock)
        global msg_dispatcher = Stratum.MsgDispatcher()

        msg_dispatcher[request1_type] = (conn, params) -> begin
            @test Stratum.is_currently_handling_msg(msg_dispatcher)
            params.fieldA == 1 ? "YES" : "NO"
        end
        msg_dispatcher[request2_type] = (conn, params) -> Stratum.StratumError(-32600, "Our message", nothing)
        msg_dispatcher[notify1_type] = (conn, params) -> g_var = params

        run(conn)

        for msg in conn
            Stratum.dispatch_msg(conn, msg_dispatcher, msg)
        end
    catch err
        Base.display_error(stderr, err, catch_backtrace())
    end

    wait(server_is_up)

    sock2 = connect(global_socket_name1)
    conn2 = StratumEndpoint(sock2, sock2)

    run(conn2)

    Stratum.send(conn2, notify1_type, "TEST")

    res = Stratum.send(conn2, request1_type, Foo(fieldA=1, fieldB="FOO"))

    @test res == "YES"
    @test g_var == "TEST"

    @test_throws Stratum.StratumError(-32600, "Our message", nothing) Stratum.send(conn2, request2_type, nothing)

    close(conn2)
    close(sock2)
    close(conn)

    fetch(server_task)

    # Now we test a faulty server

    if Sys.iswindows()
        global_socket_name2 = "\\\\.\\pipe\\jsonrpc-testrun2"
    elseif Sys.isunix()
        global_socket_name2 = joinpath(tempdir(), "jsonrpc-testrun2")
    else
        error("Unknown operating system.")
    end

    server_is_up = Base.Condition()

    server_task2 = @async try
        server = listen(global_socket_name2)
        notify(server_is_up)
        sock = accept(server)
        global conn = Stratum.StratumEndpoint(sock, sock)
        global msg_dispatcher = Stratum.MsgDispatcher()

        msg_dispatcher[request2_type] = (conn, params)->34 # The request type requires a `String` return, so this tests whether we get an error.

        run(conn)

        for msg in conn
            @test_throws ErrorException("The handler for the 'request2' request returned a value of type $Int, which is not a valid return type according to the request definition.") Stratum.dispatch_msg(conn, msg_dispatcher, msg)
        end
    catch err
        Base.display_error(stderr, err, catch_backtrace())
    end

    wait(server_is_up)

    sock2 = connect(global_socket_name2)
    conn2 = StratumEndpoint(sock2, sock2)

    run(conn2)

    @test_throws Stratum.StratumError(-32603, "The handler for the 'request2' request returned a value of type $Int, which is not a valid return type according to the request definition.", nothing) Stratum.send(conn2, request2_type, nothing)

    close(conn2)
    close(sock2)
    close(conn)

    fetch(server_task)

end


text:
InputFolder = './Images/NucleiNoisy/';
OutputFolder = './Results/Images/NucleiNoisy/';
Rescale = 1;

@iA = '*.tif';

@fxg_sLoGLocMax [iA] > [S];
params.GRad = 7;
params.LocalMaxBox = [9 9];
params.ThrLocalMax = 0.3;
params.MinLocThr = 0.2;
/endf

/show iA > S;
/keep S > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
using ArgParse

export append_module_info

function append_module_info(input_jl_file::String, output_jl_file::String, module_name::String )

    open(output_jl_file, "w") do out_file

        open(input_jl_file) do in_file

            module_inserted = false
            prev_line_comment = true

            for ln in eachline(in_file)

                if isempty(ln) && !(module_inserted)
                    write(out_file, "\n")
                    write(out_file, "module $module_name \n")
                    module_inserted = true

                else
                    write(out_file, "$(ln) \n")

                end

            end

            write(out_file, "\nend \n")

        end

    end

end

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "input_jl_file"
            help = "Path to the input Julia file."
            required = true
        "output_jl_file"
            help = "Path to write the output (julia file) to."
            required = true
        "module_name"
            help = "The module name to insert into the Julia code."
            required = true
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    append_module_info(parsed_args["input_jl_file"], parsed_args["output_jl_file"], parsed_args["module_name"])
end

main()

# usage:
# julia insert_module_name.jl ../src/DiscreteStochErlangEpidModel.jl ../src/DiscreteStochErlangEpidModel_v2.jl DiscreteStochErlangEpidModel


text:
@testset "Plots" begin

    ## RM Plotting
    game = MatrixGame([
        (0,0) (-1,1) (1,-1);
        (1,-1) (0,0) (-1,1);
        (-1,1) (1,-1) (0,0)
    ])

    p1 = player(game, 1, [0.2,0.3,0.5])
    p2 = player(game, 2, [0.2,0.3,0.5])
    train_both!(p1, p2, 100)

    l1 = length(p1.hist)
    l2 = length(first(p1.hist))
    mat = CFR.cumulative_strategies(p1)
    @test size(mat) == (l1, l2)

    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), p1,p2) ≠ nothing
    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), p1) ≠ nothing

    ## Matrix CFR Plotting

    game = SimpleIIGame([
        (0,0) (-1,1) (1,-1);
        (1,-1) (0,0) (-1,1);
        (-1,1) (1,-1) (0,0)
    ])

    p1 = player(game, 1, [0.2,0.3,0.5])
    p2 = player(game, 2, [0.2,0.3,0.5])
    train_both!(p1, p2, 100)

    l1 = length(p1.hist)
    l2 = length(first(p1.hist))
    mat = CFR.cumulative_strategies(p1)
    @test size(mat) == (l1, l2)

    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), p1,p2) ≠ nothing
    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), p1) ≠ nothing

    ## CFR plotting

    game = IIEMatrixGame()

    sol = CFRSolver(game; debug=true)
    train!(sol, 100)

    hist = sol.I[0].hist
    l1 = length(hist)
    l2 = length(first(hist))
    mat = CFR.cumulative_strategies(hist)
    @test size(mat) == (l1, l2)

    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), sol) ≠ nothing
    @test RecipesBase.apply_recipe(Dict{Symbol,Any}(), sol.I[0]) ≠ nothing
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

## interpolate from a set of elements

element1 = Element(Seg2, (1, 2))
element2 = Element(Seg2, (2, 3))
X = Dict(1 => [0.0], 2 => [1.0], 3 => [2.0])
T = Dict(1 => [0.0], 2 => [1.0], 3 => [0.0])
problem = Problem(Heat, "foo", 1)
add_elements!(problem, element1, element2)
problem_elements = get_elements(problem)
update!(problem_elements, "geometry", X)
update!(problem_elements, "temperature", T)
@test isnan(problem("temperature", [-0.1], 0.0))
@test isapprox(problem("temperature", [0.0], 0.0), [0.0])
@test isapprox(problem("temperature", [0.5], 0.0), [0.5])
@test isapprox(problem("temperature", [1.0], 0.0), [1.0])
@test isapprox(problem("temperature", [1.5], 0.0), [0.5])
@test isapprox(problem("temperature", [2.0], 0.0), [0.0])
@test isnan(problem("temperature", [ 2.1], 0.0))


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using DifferentialEquations, DiffEqCallbacks

using FMIImport: fmi3Instance

############ Model-Exchange ############

# Read next time event from FMU and provide it to the integrator
function time_choice(c::fmi3Instance, integrator)
    @debug "time_choice(_, _): Time event @ t=$(integrator.t)"

    discreteStatesNeedUpdate, terminateSimulation, nominalsOfContinuousStatesChanged, valuesOfContinuousStatesChanged, nextEventTimeDefined, nextEventTime = fmi3UpdateDiscreteStates(c)
    fmi3EnterContinuousTimeMode(c)
    if Bool(nextEventTimeDefined)
        c.timeEvent = integrator.t >= nextEventTime
        nextEventTime
    else
        Inf
    end
end

# Handles events and returns the values and nominals of the changed continuous states.
function handleEvents(c::fmi3Instance, enterEventMode::Bool, exitInContinuousMode::Bool)
    nominalsChanged = fmi3False
    valuesChanged = fmi3False
    if enterEventMode

        fmi3EnterEventMode(c, c.stepEvent, c.stateEvent, c.rootsFound, Csize_t(c.fmu.modelDescription.numberOfEventIndicators), c.timeEvent)
        # TODO inputEvent handling

        discreteStatesNeedUpdate = fmi3True
        while discreteStatesNeedUpdate == fmi3True

            # update discrete states
            discreteStatesNeedUpdate, terminateSimulation, nominalsOfContinuousStatesChanged, valuesOfContinuousStatesChanged, nextEventTimeDefined, nextEventTime = fmi3UpdateDiscreteStates(c)

            if valuesOfContinuousStatesChanged == fmi3True
                valuesChanged = true
            end

            if nominalsOfContinuousStatesChanged == fmi3True
                nominalsChanged = true
            end

            if terminateSimulation == fmi3True
                @error "fmi3UpdateDiscreteStates returned error!"
            end
        end
    end

    if exitInContinuousMode
        fmi3EnterContinuousTimeMode(c)
    end
    @debug "handleEvents(_, $(enterEventMode), $(exitInContinuousMode)): rootsFound: $(c.rootsFound)   valuesChanged: $(valuesChanged)   continuousStates: $(fmi3GetContinuousStates(c))",
    return valuesChanged, nominalsChanged

end

# Returns the event indicators for an FMU.
function condition(c::fmi3Instance, out, x, t, integrator, inputFunction, inputValues::Array{fmi3ValueReference}) # Event when event_f(u,t) == 0
    if inputFunction !== nothing
        fmi3SetFloat64(c, inputValues, inputFunction(integrator.t))
    end

    c.stateEvent = fmi3False
    fmi3SetTime(c, t)
    fmi3SetContinuousStates(c, x)
    indicators = fmi3GetEventIndicators(c)
    if length(indicators) > 0
        for i in 1:length(indicators)
            if c.z_prev[i] < 0 && indicators[i] >= 0
                c.rootsFound[i] = 1
            elseif c.z_prev[i] > 0 && indicators[i] <= 0
                c.rootsFound[i] = -1
            else
                c.rootsFound[i] = 0
            end
            c.stateEvent |= (c.rootsFound[i] != 0)
            # c.z_prev[i] = indicators[i]
        end
    end
    @debug "condition(_, _, $(x), $(t), _, _, _): eventIndicators $indicators   rootsFound $(c.rootsFound)   stateEvent $(c.stateEvent)"
    copy!(out, indicators)
end

# Handles the upcoming events.
function affectFMU!(c::fmi3Instance, integrator, idx, inputFunction, inputValues::Array{fmi3ValueReference}, force=false)
    # Event found - handle it

    @debug "affectFMU!(_, _, $(idx), _, _): x:$(integrator.u)   [before handle events]"

    fmi3SetContinuousStates(c, integrator.u)

    continuousStatesChanged, nominalsChanged = handleEvents(c, true, Bool(sign(idx)))

    @debug "affectFMU!(_, _, $(idx), _, _): continuousStatesChanged=$(continuousStatesChanged)   x_int:$(integrator.u)   x_fmu:$(fmi3GetContinuousStates(c))   [after handle events]"

    if inputFunction !== nothing
        fmi3SetFloat64(c, inputValues, inputFunction(integrator.t))
    end

    if continuousStatesChanged == fmi3True
        integrator.u = fmi3GetContinuousStates(c)
        @debug "affectFMU!(_, _, $(idx), _, _): Set new state $(integrator.u)"
    end

    if nominalsChanged == fmi3True
        x_nom = fmi3GetNominalsOfContinuousStates(c)
    end
    # timeEventCb = PresetTimeCallback(2.0, (integrator) -> affect!(c, integrator, 0))
end

# Does one step in the simulation.
function stepCompleted(c::fmi3Instance, x, t, integrator, inputFunction, inputValues::Array{fmi3ValueReference})

    fmi3SetContinuousStates(c, x)

    indicators = fmi3GetEventIndicators(c)
        if length(indicators) > 0
        c.stateEvent = fmi3False

        for i in 1:length(indicators)
            if c.z_prev[i] < 0 && indicators[i] >= 0
                c.rootsFound[i] = 1
            elseif c.z_prev[i] > 0 && indicators[i] <= 0
                c.rootsFound[i] = -1
            else
                c.rootsFound[i] = 0
            end
            c.stateEvent |= (c.rootsFound[i] != 0)
            c.z_prev[i] = indicators[i]
        end
    end

    (status, c.stepEvent, terminateSimulation) = fmi3CompletedIntegratorStep(c, fmi3True)
    @debug "stepCompleted(_, $(x), $(t), _, _, _): stepEvent $(c.stepEvent)"
    @assert terminateSimulation == fmi3False "completed Integratorstep failed!"

end

# Returns the state derivatives of the FMU.
function fx(c::fmi3Instance, x, p, t)
    @debug "fx($(x), _, $(t))"
    fmi3SetTime(c, t)
    fmi3SetContinuousStates(c, x)
    dx = fmi3GetContinuousStateDerivatives(c)
end

# save FMU values
function saveValues(c::fmi3Instance, recordValues, u, t, integrator)
    fmi3SetTime(c, t)
    x = integrator.sol(t)
    fmi3SetContinuousStates(c, x)

    (fmi3GetFloat64(c, recordValues)...,)
end

"""
Source: FMISpec3.0, Version D5ef1c1: 3.3. Code Example

Simulates a FMU instance for the given simulation time interval.
State- and Time-Events are handled correctly.

Returns a tuple of type (ODESolution, DiffEqCallbacks.SavedValues).
If keyword `recordValues` is not set, a tuple of type (ODESolution, nothing) is returned for consitency.
"""
function fmi3SimulateME(c::fmi3Instance, t_start::Union{Real, Nothing} = nothing, t_stop::Union{Real, Nothing} = nothing;
    solver = nothing,
    customFx = nothing,
    recordValues::fmi3ValueReferenceFormat = nothing,
    saveat = [],
    saveAtEvent::Bool = false,
    setup::Bool = true,
    reset = nothing, # nothing = auto
    inputValues::fmi3ValueReferenceFormat = nothing,
    inputFunction = nothing,
    rootSearchInterpolationPoints = 100,
    kwargs...)

    if t_start == nothing
        t_start = c.fmu.modelDescription.defaultStartTime
    end

    if t_stop == nothing
        t_stop = c.fmu.modelDescription.defaultStopTime
    end

    recordValues = prepareValueReference(c, recordValues)
    inputValues = prepareValueReference(c, inputValues)
    solution = nothing
    callbacks = []
    savedValues = nothing

    savingValues = (length(recordValues) > 0)
    hasInputs = (length(inputValues) > 0)

    if savingValues
        savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})

        savingCB = SavingCallback((u,t,integrator) -> saveValues(c, recordValues, u, t, integrator),
                                  savedValues,
                                  saveat=saveat)
        push!(callbacks, savingCB)
    end

    # auto correct reset if only setup is given
    if reset === nothing
        reset = setup
    end
    @assert !(setup==false && reset==true) "fmi3SimulateME(...): keyword argument `setup=false`, but `reset=true`. This may cause a FMU crash."

    if reset
        fmi3Reset(c)
    end

    if setup
        fmi3EnterInitializationMode(c, t_start, t_stop)
        fmi3ExitInitializationMode(c)
    end

    eventHandling = c.fmu.modelDescription.numberOfEventIndicators > 0
    timeEventHandling = false

    if eventHandling
        discreteStatesNeedUpdate = fmi3True
        while discreteStatesNeedUpdate == fmi3True
            discreteStatesNeedUpdate, terminateSimulation, nominalsOfContinuousStatesChanged, valuesOfContinuousStatesChanged, nextEventTimeDefined, nextEventTime = fmi3UpdateDiscreteStates(c)
            @assert terminateSimulation == fmi3False ["Initial Event handling failed!"]
            timeEventHandling |= (nextEventTimeDefined == fmi3True)
        end
    end

    if customFx === nothing
        customFx = (x, p, t) -> fx(c, x, p, t)
    end

    fmi3EnterContinuousTimeMode(c)

    if eventHandling
        c.z_prev = fmi3GetEventIndicators(c)
    end

    # First evaluation of the FMU
    x0 = fmi3GetContinuousStates(c)
    x0_nom = fmi3GetNominalsOfContinuousStates(c)

    p = []
    problem = ODEProblem(customFx, x0, (t_start, t_stop), p,)

    # callback functions

    # use step callback always if we have inputs or need evenet handling
    if hasInputs || eventHandling
        stepCb = FunctionCallingCallback((x, t, integrator) -> stepCompleted(c, x, t, integrator, inputFunction, inputValues);
                                            func_everystep = true,
                                            func_start = true)
        push!(callbacks, stepCb)
    end

    if eventHandling

        eventCb = VectorContinuousCallback((out, x, t, integrator) -> condition(c, out, x, t, integrator, inputFunction, inputValues),
                                           (integrator, idx) -> affectFMU!(c, integrator, idx, inputFunction, inputValues, true),
                                           Int64(c.fmu.modelDescription.numberOfEventIndicators);
                                           rootfind = RightRootFind,
                                           save_positions=(saveAtEvent,saveAtEvent),
                                           interp_points=rootSearchInterpolationPoints)#,abstol=1e-16, reltol=1e-12, repeat_nudge=1//100)
        push!(callbacks, eventCb)

        if timeEventHandling
            timeEventCb = IterativeCallback((integrator) -> time_choice(c, integrator),
                                            (integrator) -> affectFMU!(c, integrator, 0, inputFunction, inputValues), Float64;
                                            initial_affect = true,
                                            save_positions=(saveAtEvent,saveAtEvent))
            push!(callbacks, timeEventCb)
        end
    end

    if solver === nothing
        solution = solve(problem; callback = CallbackSet(callbacks...), saveat = saveat, kwargs...)
    else
        solution = solve(problem, solver; callback = CallbackSet(callbacks...), saveat = saveat, kwargs...)
    end

    if savingValues
        return solution, savedValues
    else
        return solution
    end
end

"""
Source: FMISpec3.0, Version D5ef1c1: 4.3. Code Examples

Starts a simulation of the Co-Simulation FMU instance.

Returns a tuple of (success::Bool, DiffEqCallbacks.SavedValues) with success = `true` or `false`.
If keyword `recordValues` is not set, a tuple of type (success::Bool, nothing) is returned for consitency.

ToDo: Improve Documentation.
"""
function fmi3SimulateCS(c::fmi3Instance, t_start::Real, t_stop::Real;
                        recordValues::fmi3ValueReferenceFormat = nothing,
                        saveat = [],
                        setup::Bool = true,
                        reset = nothing,
                        inputValues::fmi3ValueReferenceFormat = nothing,
                        inputFunction = nothing)

    recordValues = prepareValueReference(c, recordValues)
    inputValues = prepareValueReference(c, inputValues)
    variableSteps = c.fmu.modelDescription.CScanHandleVariableCommunicationStepSize

    success = false
    savedValues = nothing

    # default setup
    if length(saveat) == 0
        saveat = LinRange(t_start, t_stop, 100)
    end
    dt = (t_stop - t_start) / length(saveat)

    # setup if no variable steps
    if variableSteps == false
        if length(saveat) >= 2
            dt = saveat[2] - saveat[1]
        end
    end

    # auto correct reset if only setup is given
    if reset === nothing
        reset = setup
    end
    @assert !(setup==false && reset==true) "fmi3SimulateME(...): keyword argument `setup=false`, but `reset=true`. This may cause a FMU crash."

    if reset
        fmi3Reset(c)
    end

    if setup
        fmi3EnterInitializationMode(c, t_start, t_stop)
        fmi3ExitInitializationMode(c)
    end

    t = t_start

    record = length(recordValues) > 0

    #numDigits = length(string(round(Integer, 1/dt)))
    noSetFMUStatePriorToCurrentPoint = fmi3False
    eventEncountered = fmi3False
    terminateSimulation = fmi3False
    earlyReturn = fmi3False
    lastSuccessfulTime = fmi3Float64(0.0)

    if record
        savedValues = SavedValues(Float64, Tuple{collect(Float64 for i in 1:length(recordValues))...})

        i = 1

        values = (fmi3GetFloat64(c, recordValues)...,)
        DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
        DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, values, Val{false})

        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            if inputFunction !== nothing
                fmi3SetFloat64(c, inputValues, inputFunction(t))
            end

            fmi3DoStep(c, t, dt, fmi3True, eventEncountered, terminateSimulation, earlyReturn, lastSuccessfulTime)
            if eventEncountered == fmi3True
                @warn "Event handling"
            end
            if terminateSimulation == fmi3True
                @error "fmi3DoStep returned error!"
            end
            if earlyReturn == fmi3True
                @warn "early Return"
            end
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1

            values = (fmi3GetFloat64(c, recordValues)...,)
            DiffEqCallbacks.copyat_or_push!(savedValues.t, i, t)
            DiffEqCallbacks.copyat_or_push!(savedValues.saveval, i, values, Val{false})
        end

        success = true
        return success, savedValues
    else
        i = 1
        while t < t_stop
            if variableSteps
                if length(saveat) > i
                    dt = saveat[i+1] - saveat[i]
                else
                    dt = t_stop - saveat[i]
                end
            end

            if inputFunction !== nothing
                fmi3SetFloat64(c, inputValues, inputFunction(t))
            end

            fmi3DoStep(c, t, dt, fmi3True, eventEncountered, terminateSimulation, earlyReturn, lastSuccessfulTime)
            if eventEncountered == fmi3True
                @warn "Event handling"
            end
            if terminateSimulation == fmi3True
                @error "fmi3DoStep returned error!"
            end
            if earlyReturn == fmi3True
                @warn "early Return"
            end
            t = t + dt #round(t + dt, digits=numDigits)
            i += 1
        end

        success = true
        return success
    end
end

# TODO simulate ScheduledExecution
function fmi3SimulateSE(c::fmi3Instance, t_start::Real, t_stop::Real;
    recordValues::fmi3ValueReferenceFormat = nothing,
    saveat = [],
    setup::Bool = true,
    reset = nothing,
    inputValues::fmi3ValueReferenceFormat = nothing,
    inputFunction = nothing)
    @assert false "Not implemented"
end

"""
Starts a simulation of the fmu instance for the matching fmu type, if both types are available, CS is preferred.

Returns:
    - a tuple of (success::Bool, DiffEqCallbacks.SavedValues) with success = `true` or `false` for CS-FMUs
    - a tuple of (ODESolution, DiffEqCallbacks.SavedValues) for ME-FMUs
    - if keyword `recordValues` is not set, a tuple of type (..., nothing)

ToDo: Improve Documentation.
"""
function fmi3Simulate(c::fmi3Instance, t_start::Real = 0.0, t_stop::Real = 1.0;kwargs...)

    if fmi3IsCoSimulation(c.fmu)
        return fmi3SimulateCS(c, t_start, t_stop; kwargs...)
    elseif fmi3IsModelExchange(c.fmu)
        return fmi3SimulateME(c, t_start, t_stop; kwargs...)
    else
        error(unknownFMUType)
    end
end


text:
tree = Any[1, Any[2, Any[4, Any[7, Any[],
                                          Any[]],
                                   Any[]],
                            Any[5, Any[],
                                   Any[]]],
                     Any[3, Any[6, Any[8, Any[],
                                          Any[]],
                                   Any[9, Any[],
                                          Any[]]],
                            Any[]]]

preorder(t, f) = if !isempty(t)
                     f(t[1]); preorder(t[2], f); preorder(t[3], f)
                 end

inorder(t, f) = if !isempty(t)
                    inorder(t[2], f); f(t[1]); inorder(t[3], f)
                end

postorder(t, f) = if !isempty(t)
                      postorder(t[2], f); postorder(t[3], f); f(t[1])
                  end

levelorder(t, f) = while !isempty(t)
                       t = mapreduce(x -> isa(x, Number) ? (f(x); []) : x, vcat, t)
                   end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function get_dzhanibekov(;
    timestep=0.01,
    gravity=-9.81,
    color=RGBA(0.9,0.9,0.9,1),
    T=Float64)

    radius = 0.1
    body_length = 1.0
    body_mass = 1.0
    origin = Origin{Float64}()
    main_body = Capsule(radius, body_length, body_mass,
        color=color, name=:main)
    side_body = Capsule(0.5 * radius, 0.35 * body_length, 0.5 * body_mass,
        orientation_offset=RotY(0.5 * π),
        color=color, name=:side)
    links = [main_body, side_body]

    # Joint Constraints
    joint_float = JointConstraint(Floating(origin, links[1]), name=:floating)
    joint_fixed = JointConstraint(Fixed(links[1], links[2];
        parent_vertex=szeros(3),
        child_vertex=[-0.25 * body_length; 0.0; 0.0]), name=:fixed)
    joints = [joint_float, joint_fixed]

    links[1].inertia = Diagonal([3e-2, 1e-3, 1e-1])
    # Mechanism
    return Mechanism(origin, links, joints;
        gravity,
        timestep)
end

function initialize_dzhanibekov!(mechanism::Mechanism{T,Nn,Ne,Nb};
    linear_velocity=zeros(3),
    angular_velocity=zeros(3)) where {T,Nn,Ne,Nb}

    zero_velocity!(mechanism)
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :floating);
        xmin=[0;0;1;0;0;0; linear_velocity; angular_velocity])
    set_minimal_coordinates_velocities!(mechanism, get_joint(mechanism, :fixed))
end


text:
const obj = Dict([
    "children"=> 3,
    "cats"=> 7,
    "samoyeds"=> 2,
    "pomeranians"=> 3,
    "akitas"=> 0,
    "vizslas"=> 0,
    "goldfish"=> 5,
    "trees"=> 3,
    "cars"=> 2,
    "perfumes"=> 1
])

function parse_line(l)
    n, o1, v1, o2, v2, o3, v3 = match(r"Sue (\d+): (\w+): (\d+), (\w+): (\d+), (\w+): (\d+)", l).captures
    (parse(Int, n), Dict([
        o1=>parse(Int, v1),
        o2=>parse(Int, v2),
        o3=>parse(Int, v3),
    ]))
end

const data = parse_line.(readlines())

function correspond(d)
    i, infos = d
    all(v == obj[k] for (k, v) in infos)
end

function correspond2(d)
    i, infos = d
    b = true
    for (k, v) in infos
        if k in ("cats", "trees")
            b = b && v > obj[k]
        elseif k in ("pomeranians", "goldfish")
            b = b && v < obj[k]
        else
            b = b && v == obj[k]
        end
    end
    b
end

filter(correspond, data)[1] |> println
filter(correspond2, data)[1] |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
using DataFrames
# Add all custom functions here
include("1.4.medline.jl")
include("fetch_full_article.jl")

# diagnosis for the year 2000 to  2018 (limit is 1000 for now)
df = @time  searchBreastCancerArticles("breast neoplasms", "diagnosis",2015,2018,100)
println(head(df))

# Then fetch full artcile from PMC and store it as dataframe
df_full_text = fetchFullArticleFromPmc(df)
println(head(df_full_text))

# save it as excel
output_file = open("output/diagnosis_full_article.csv", "w")
close(output_file)
writetable("output/diagnosis_full_article.csv",df_full_text)

max_stars_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_issues_repo_name: jmulloka/Breast-Cancer-Text-Mining

max_forks_repo_name: jmulloka/Breast-Cancer-Text-Mining


text:
using DashBootstrapComponents

row = html_div(
    [
        dbc_row(
            [
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
            ],
            align = "start",
        ),
        dbc_row(
            [
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
            ],
            align = "center",
        ),
        dbc_row(
            [
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
                dbc_col(html_div("One of three columns")),
            ],
            align = "end",
        ),
        dbc_row([
            dbc_col(html_div("One of three columns"), align = "start"),
            dbc_col(html_div("One of three columns"), align = "center"),
            dbc_col(html_div("One of three columns"), align = "end"),
        ]),
    ],
    className = "pad-row",
);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
# indicator of an affine set

using LinearAlgebra
using SparseArrays
using SuiteSparse

export IndAffine

### ABSTRACT TYPE

abstract type IndAffine <: ProximableFunction end

is_affine(f::IndAffine) = true
is_generalized_quadratic(f::IndAffine) = true

fun_name(f::IndAffine) = "Indicator of an affine subspace"

### CONSTRUCTORS

"""
**Indicator of an affine subspace**

    IndAffine(A, b; iterative=false)

If `A` is a matrix (dense or sparse) and `b` is a vector, returns the indicator function of the set
```math
S = \\{x : Ax = b\\}.
```
If `A` is a vector and `b` is a scalar, returns the indicator function of the set
```math
S = \\{x : \\langle A, x \\rangle = b\\}.
```
By default, a direct method (QR factorization of matrix `A'`) is used to evaluate `prox!`.
If `iterative=true`, then `prox!` is evaluated approximately using an iterative method instead.
"""
function IndAffine(A::M, b::V; iterative=false) where {M, V}
    if iterative == false
        IndAffineDirect(A, b)
    else
        IndAffineIterative(A, b)
    end
end

### INCLUDE CONCRETE TYPES

using LinearAlgebra
using SparseArrays
using SuiteSparse

include("indAffineDirect.jl")
include("indAffineIterative.jl")


text:
"""
    struct ComplexSpace <: EuclideanSpace{ℂ}

A standard complex vector space ℂ^d with Euclidean inner product and no additional
structure. It is completely characterised by its dimension and whether its the normal space
or its dual (which is canonically isomorphic to the conjugate space).
"""
struct ComplexSpace <: EuclideanSpace{ℂ}
  d::Int
  dual::Bool
end
ComplexSpace(d::Integer = 0; dual = false) = ComplexSpace(Int(d), dual)
function ComplexSpace(dim::Pair; dual = false)
    if dim.first === Trivial()
        return ComplexSpace(dim.second; dual = dual)
    else
        msg = "$(dim) is not a valid dimension for ComplexSpace"
        throw(SectorMismatch(msg))
    end
end
function ComplexSpace(dims::AbstractDict; kwargs...)
    if length(dims) == 0
        return ComplexSpace(0; kwargs...)
    elseif length(dims) == 1
        return ComplexSpace(first(dims); kwargs...)
    else
        msg = "$(dims) is not a valid dimension dictionary for ComplexSpace"
        throw(SectorMismatch(msg))
    end
end

# convenience constructor
Base.getindex(::ComplexNumbers) = ComplexSpace
Base.:^(::ComplexNumbers, d::Int) = ComplexSpace(d)

# Corresponding methods:
#------------------------
dim(V::ComplexSpace) = V.d
isdual(V::ComplexSpace) = V.dual
Base.axes(V::ComplexSpace) = Base.OneTo(dim(V))

Base.conj(V::ComplexSpace) = ComplexSpace(dim(V), !isdual(V))

Base.oneunit(::Type{ComplexSpace}) = ComplexSpace(1)
⊕(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(dim(V1)+dim(V2), isdual(V1)) :
    throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist"))
fuse(V1::ComplexSpace, V2::ComplexSpace) = ComplexSpace(V1.d*V2.d)
flip(V::ComplexSpace) = dual(V)

infimum(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(min(dim(V1), dim(V2)), isdual(V1)) :
    throw(SpaceMismatch("Infimum of space and dual space does not exist"))
supremum(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(max(dim(V1), dim(V2)), isdual(V1)) :
    throw(SpaceMismatch("Supremum of space and dual space does not exist"))

Base.show(io::IO, V::ComplexSpace) = print(io, isdual(V) ? "(ℂ^$(V.d))'" : "ℂ^$(V.d)")


text:
using NNLS
using BenchmarkTools
import NonNegLeastSquares
using PyCall
using Random

const pyopt_nnls = pyimport_conda("scipy.optimize", "scipy")[:nnls]

function runbenchmarks()
    Random.seed!(1)
    println("M\tN\tt_nnls\tt_nnls_with_workspace_reuse")
    for m in 50:50:250
        for n_over_m in [0.5, 1, 2]
            n = round(Int, m * n_over_m)
            inputs = [(randn(m, n), randn(m)) for i in 1:50]

            for (A, b) in inputs
                # @assert nnls(A, b) ≈ NonNegLeastSquares.nnls(A, b)
                @assert nnls(A, b) == pyopt_nnls(A, b)[1]
            end

            t0 = @elapsed begin
                for (A, b) in inputs
                    pyopt_nnls(A, b)
                end
            end

            x1, t1, m1, gc1, mem1 = @timed begin
                for (A, b) in inputs
                    nnls(A, b)
                end
            end
            x2, t2, m2, gc2, mem2 = @timed begin
                work = NNLSWorkspace(m, n)
                for (A, b) in inputs
                    solve!(work, A, b)
                end
            end
            # @show (t0, t1, t2)
            @printf("%d\t%d\t%.2f\t%.2f\n", m, n, t1 / t0, t2 / t0)
            # b1 = median(@benchmark nnls(A, b) setup=(A = randn($m, $n); b = randn($m)))
            # b2 = median(@benchmark NonNegLeastSquares.nnls(A, b) setup=(A = randn($m, $n); b = randn($m)))
            # r = ratio(b2, b1)
            # @printf "%d\t%d\t%.1f\t%.1f\n" m n r.time r.memory
        end
    end
end

runbenchmarks()


text:
type Path <: AbstractPath
    walks::Vector
end

Path() = Path(Any[])

length(p::Path) = length(p.walks)

has_node(p::Path, v) = v in p.walks

==(p1::Path, p2::Path) = (p1.walks == p2.walks)

# find the shortest path via depth first search
function _DFS_shortest_path{V}(g::DiGraph{V},
                               v1::V,
                               v2::V,
                               path = Path(),
                               shortest =  Union{};
                               verbose = false)

    path = deepcopy(path)
    path.walks = [path.walks; v1]

    if verbose
        println("current:", path)
    end

    if v1 == v2
        return path
    end

    for node in forward_neighbors(g, v1)
        if !has_node(path, node)
            if (shortest ==  Union{}) || length(path) < length(shortest)
                newPath = _DFS_shortest_path(g, node, v2, path, shortest, verbose = verbose)
                if !(newPath ==  Union{})
                    shortest = newPath
                end
            end
        end
    end
    return shortest
end

@doc doc"""
`shortest_path(g, v1, v2 [, verbose = false])` finds the shortest path from `v1` to `v2` on
the time graph `g`. if `verbose = true`, prints the current path at each search step.
"""->
shortest_path{V}(g::DiGraph{Node{V}}, v1::V, v2::V; verbose = false) =
_DFS_shortest_path(g, make_node(g,v1), make_node(g,v2), verbose = verbose)

#shortest_path{V}(g::DiGraph{AttributeNode{V}}, v1::V, v2::V; verbose = false) =
#_DFS_shortest_path(g, make_node(g,v1), make_node(g,v2), verbose = verbose)

@doc doc"""
`shortest_distance(g, v1, v2)` finds the shortest distance from `v1` to `v2`
on the time graph `g`.
"""->
shortest_distance{V}(g::DiGraph{Node{V}}, v1::V, v2::V) =
          shortest_path(g, v1, v2) == Union{} ? Inf : length(shortest_path(g, v1, v2)) - 1

#shortest_distance{V}(g::DiGraph{AttributeNode{V}}, v1::V, v2::V) =
#            shortest_path(g, v1, v2) ==  Union{} ? Inf : length(shortest_path(g, v1, v2)) - 1


text:
# test whether prox satisfies some conditions

using LinearAlgebra
using SparseArrays
using ProximalOperators
using Test

stuff = [
  Dict( "constr" => LeastSquares,
        "params" => ( (randn(20, 10), randn(20)), (randn(15, 40), randn(15), rand()), (rand(Complex{Float64}, 15, 40), rand(Complex{Float64}, 15), rand()), (sprandn(100,1000,0.05), randn(100), rand()) ),
        "args"   => ( randn(10), randn(40), rand(Complex{Float64}, 40), randn(1000) ),
        "gammas" => ( rand(), rand(), rand(), rand() ),
        # test first order optimality condition of prox
        "test"   => (f, x, gamma, y) -> norm(y + gamma*f.lambda*(f.A'*(f.A*y-f.b)) - x) <= 1e-10
      ),

  Dict( "constr" => IndSimplex,
        "params" => ( (), (2,), (5.0,), (rand(),) ),
        "args"   => ( randn(10), randn(10), randn(10), randn(10) ),
        "gammas" => ( rand(), rand(), rand(), rand() ),
        # test y belonging to the simplex
        "test"   => (f, x, gamma, y) -> all(y .>= 0.0) && abs(sum(y)-f.a) <= (1+f.a)*1e-12
      ),

  Dict( "constr" => IndBallL1,
        "params" => ( (), (1.7,), (5.0,), (rand(), ) ),
        "args"   => ( randn(10), randn(10), randn(10), randn(10) ),
        "gammas" => ( rand(), rand(), rand(), rand() ),
        # test y belonging to the L1 ball
        "test"   => (f, x, gamma, y) -> norm(y, 1) <= (1+1e-12)*f.r
      ),

  Dict( "constr" => HuberLoss,
        "params" => ( (), (rand(),), (rand(), rand()), (rand(), rand()) ),
        "args"   => ( randn(10), randn(8, 10), randn(20), rand(Complex{Float64}, 12, 15) ),
        "gammas" => ( rand(), rand(), rand(), rand() ),
        # test optimality condition of prox
        "test"   => (f, x, gamma, y) -> isapprox((x-y)/gamma, (norm(y) <= f.rho ? f.mu*y : f.rho*f.mu*y/norm(y)))
      ),

  Dict( "constr" => SqrHingeLoss,
        "params" => ( (randn(5),), (randn(5), 0.1+rand()), (randn(3, 5), 0.1+rand()) ),
        "args"   => ( randn(5), randn(5), randn(3, 5) ),
        "gammas" => ( 0.1+rand(), 0.1+rand(), 0.1+rand() ),
        # test optimality condition of prox
        "test"   => (f, x, gamma, z) -> isapprox((x-z)/gamma, -2 .* f.mu.*f.y.*max.(0, 1 .- f.y.*z))
      ),
]

@testset "$(d["constr"])" for d in stuff
      constr = d["constr"]
      params = d["params"]
      args = d["args"]
      gammas = d["gammas"]
      test = d["test"]
      for i = 1:length(params)
            f = constr(params[i]...)
            y, fy = prox(f, args[i], gammas[i])
            @test test(f, args[i], gammas[i], y)
      end
end


text:
using MAT
using jInv.Mesh
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.Utils
using OED

model = 1;
IPmodes = [0;1;2;3];
lambda = .1;

r       = 48
l       = 2
pad     = 4
nex     = 20
matfile = "data/OED_shepp.mat"
resfile = "results/plotObjectiveB/OED_shepp"
matfile = matread(matfile)
n       = Int(matfile["n"])
L       = .1*speye(n*n)
Xtrue   = matfile["x_true"]

Ce = matfile["Ce"]
ce = Ce*Xtrue
Ci = matfile["Ci"]
ci = matfile["ci"]

if nex!==-1
	idx = randperm(size(Xtrue,2))[1:nex]
	Xtrue = Xtrue[:,idx]
	ci   = ci[:,idx]
	ce   = ce[:,idx]
end

Mimg   = getRegularMesh([-1 1 -1 1.],[n,n])
domain = copy(Mimg.domain);
domain[1:2:end] -= pad*Mimg.h
domain[2:2:end] += pad*Mimg.h
Mgrid = getRegularMesh(domain,[n+2*pad,n+2*pad])
wGrid = getRegularMesh([0 180 0 180],[32 32]);

# divide problems among available workers
ids = mod.(0:size(Xtrue,2)-1,nworkers())+1
OEDparam = Array{OEDTypeBParam}(nworkers())
for k=1:nworkers()
	OEDparam[k] = getOEDTypeBParam(r,l,Xtrue[:,ids.==k],Mimg,Mgrid)
	OEDparam[k].noiseLevel = 1e-5
end

ipDict  = Dict(0=>:unconstr, 1=>:eqconstr, 2=>:nonnegconstr, 3=>:boxconstr)

for k=1:length(IPmodes)
	IPmode  = ipDict[IPmodes[k]]
	if IPmode==:unconstr
		pFor = Array{UnconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = UnconstrainedOEDParam(OEDparam[j],L,zeros(0,0),[])
		end
	elseif IPmode==:eqconstr
		if isempty(Ce) || isempty(ce)
			error("equality constraints must be provided. To run unconstrained OED, use IPmode=:unconstr.")
		end
		pFor = Array{EqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
			pFor[j] = EqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],zeros(0,0),[])
		end
	elseif IPmode==:nonnegconstr
		I = speye(size(Xtrue,1))
		pFor = Array{IneqconstrainedOEDParam}(nworkers())
		for j=1:nworkers()
				pFor[j] = getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],I,0*OEDparam[j].Ftrue,
				          ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end

	elseif IPmode==:boxconstr
		if isempty(Ci) || isempty(ci)
			error("inequality constraints must be provided. To run equality constrained OED, use IPmode=:eqconstr.")
		end
		pFor = Array{IneqconstrainedOEDParam}(nworkers())

		for j=1:nworkers()
		pFor[j] = 	getIneqconstrainedOEDParam(OEDparam[j],L,Ce,ce[:,ids.==j],Ci,ci[:,ids.==j],
		ftolIPQP=1e-3,mutolIPQP=1e-3,outIPQP=0)
		end
	else
		error("Invalid IPmode = $IPmode")
	end
	display(pFor[1])

	if nworkers()>1
		pMis  = Array{RemoteChannel}(nworkers())
		for j=1:nworkers()
			pMis[j] = initRemoteChannel(getMisfitParam,workers()[j], pFor[j],ones(size(OEDparam[j].Ftrue)),Xtrue[:,ids.==j],SSDFun)
		end
	else
		pMis  = getMisfitParam(pFor[1],ones(size(Xtrue)),Xtrue,SSDFun);
	end
	w1,w2 = getCellCenteredAxes(wGrid)
	F     = zeros(length(w1),length(w2))
	for k1=1:length(w1)
	    println("k1=$k1")
	    for k2=1:k1
	        Dc,F[k1,k2],dF,d2F,pMis,tMis = computeMisfit([w1[k1],w2[k2]],pMis);
	    end
	end
	outfile = @sprintf "%s-%s.mat" resfile IPmode
	matwrite(outfile,Dict("F"=>F,"domain"=>wGrid.domain,"n"=>wGrid.n))
end


text:
# This file was generated, do not modify it. # hide
test, train = partition(eachindex(y), 0.8);

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
import DataFrames
import Dates
import TimeSeries

function verify_forecasts(sys::System, num_initial_times, num_forecasts, horizon)
    initial_times = get_forecast_initial_times(sys)
    if length(initial_times) != num_initial_times
        @error "count of initial_times doesn't match" num_initial_times initial_times
        return false
    end

    total_forecasts = 0
    forecasts = iterate_forecasts(sys)
    for forecast in forecasts
        if IS.get_horizon(forecast) != horizon
            @error "horizon doesn't match" IS.get_horizon(forecast) horizon
            return false
        end
        total_forecasts += 1
    end

    if num_forecasts != total_forecasts
        @error "num_forecasts doesn't match" num_forecasts total_forecasts
        return false
    end

    return true
end

@testset "Test read_time_series_metadata" begin
    for file in ("timeseries_pointers.json", "timeseries_pointers.csv")
        filename = joinpath(RTS_GMLC_DIR, file)
        forecasts = IS.read_time_series_metadata(filename)
        @test length(forecasts) == 282

        for forecast in forecasts
            @test isfile(forecast.data_file)
        end
    end
end

@testset "Test forecast normalization" begin
    component_name = "122_HYDRO_1"
    timeseries_file = joinpath(
        DATA_DIR,
        "forecasts",
        "RTS_GMLC_forecasts",
        "gen",
        "Hydro",
        "DAY_AHEAD_hydro.csv",
    )
    timeseries = IS.read_time_series(timeseries_file)[Symbol(component_name)]
    max_value = maximum(TimeSeries.values(timeseries))

    metadata = IS.TimeseriesFileMetadata(
        "DAY_AHEAD",
        "Generator",
        "122_HYDRO_1",
        "activepower",
        1.0,
        timeseries_file,
        [],
        "DeterministicInternal",
    )

    # Test code path where no normalization occurs.
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries)

    # Test code path where timeseries is normalized by dividing by the max value.
    metadata.scaling_factor = "Max"
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries ./ max_value)

    # Test code path where timeseries is normalized by dividing by a custom value.
    sf = 95.0
    metadata.scaling_factor = sf
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries ./ sf)
end

@testset "Test single forecast addition" begin
    component_name = "122_HYDRO_1"
    label = "activepower"
    timeseries_file = joinpath(
        DATA_DIR,
        "forecasts",
        "RTS_GMLC_forecasts",
        "gen",
        "Hydro",
        "DAY_AHEAD_hydro.csv",
    )
    timeseries = IS.read_time_series(timeseries_file)[Symbol(component_name)]

    # Test with a filename.
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    component = get_component(HydroEnergyReservoir, sys, component_name)
    add_forecast!(sys, timeseries_file, component, label, 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.timestamp(get_data(forecast)) == TimeSeries.timestamp(timeseries)
    @test TimeSeries.values(get_data(forecast)) == TimeSeries.values(timeseries)

    # Test with TimeSeries.TimeArray.
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    component = get_component(HydroEnergyReservoir, sys, component_name)
    add_forecast!(sys, timeseries, component, label, 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(get_data(forecast)) == TimeSeries.values(timeseries)

    # Test with DataFrames.DataFrame.
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "RTS_GMLC.m")))
    component = get_component(HydroEnergyReservoir, sys, component_name)
    df = DataFrames.DataFrame(timeseries)
    add_forecast!(sys, df, component, label, 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
end

@testset "Forecast data matpower" begin
    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))
    forecasts_metadata = joinpath(FORECASTS_DIR, "5bus_ts", "timeseries_pointers_da.json")
    add_forecasts!(sys, forecasts_metadata)
    @test verify_forecasts(sys, 1, 5, 24)

    # Add the same files.
    # This will fail because the component-label pairs will be duplicated.
    @test_throws ArgumentError add_forecasts!(sys, forecasts_metadata)

    forecasts_metadata = joinpath(FORECASTS_DIR, "5bus_ts", "timeseries_pointers_rt.json")

    ## This will fail because the resolutions are different.
    @test_throws PowerSystems.DataFormatError add_forecasts!(sys, forecasts_metadata)

    ## TODO: need a dataset with same resolution but different horizon.

    sys = System(PowerSystems.PowerModelsData(joinpath(MATPOWER_DIR, "case5_re.m")))
    add_forecasts!(sys, forecasts_metadata)
    @test verify_forecasts(sys, 1, 5, 288)
end


text:
# inplace
function constrain_all(
    m::AbstractMetricParams{T},
    us::AbstractArray{T,2},
    vs::AbstractArray{T,2},
    μ
) where {T<:Number}
    us_cols = eachcol(us)
    vs_cols = eachcol(vs)
    curried(u, v) = constrain(m, u, v, μ = μ)
    @. vs[1, :] = curried(us_cols, vs_cols)
    vs
end

# SMatrix -- do we actually want to support this?
function constrain_all(
    m::AbstractMetricParams{T},
    us::SMatrix{4,M,T,L},
    vs::SMatrix{4,M,T,L},
    μ
) where {M,T<:Number,L}
    alloc_vs = constrain_all(m, us, similar(vs), μ)
    SMatrix{4,M}(alloc_vs)
end

# vectors of vectors
function constrain_all(m::AbstractMetricParams{T}, us, vs, μ) where {T<:Number}
    curried(u, v) = constrain_all(m, u, v, μ)
    curried.(us, vs)
end

function constrain_all(
    m::AbstractMetricParams{T},
    u::AbstractVector{T},
    v::AbstractVector{T},
    μ
) where {T<:Number}
    v[1] = constrain(m, u, v, μ = μ)
    v
end

function constrain_all(
    m::AbstractMetricParams{T},
    u::StaticVector{S,T},
    v::StaticVector{S,T},
    μ
) where {S,T<:Number}
    mut = MVector{S,T}(v)
    mut[1] = constrain(m, u, v, μ = μ)
    SVector{S,T}(mut)
end


text:
# Matern 3/2 isotropic covariance function

"""
    Mat32Iso <: MaternIso

Isotropic Matern 3/2 kernel (covariance)
```math
k(x,x') = σ²(1 + √3|x-x'|/ℓ)\\exp(-√3|x-x'|/ℓ)
```
with length scale ``ℓ`` and signal standard deviation ``σ``.
"""
mutable struct Mat32Iso{T<:Real} <: MaternIso
    "Length scale"
    ℓ::T
    "Signal variance"
    σ2::T
    "Priors for kernel parameters"
    priors::Array
end

"""
Matern 3/2 isotropic covariance function

    Mat32Iso(ll::T, lσ::T)

# Arguments
  - `ll::Real`: length scale (given on log scale)
  - `lσ::Real`: signal standard deviation (given on log scale)
"""
Mat32Iso(ll::T, lσ::T) where T = Mat32Iso{T}(exp(ll), exp(2 * lσ), [])

function set_params!(mat::Mat32Iso, hyp::AbstractVector)
    length(hyp) == 2 || throw(ArgumentError("Matern 3/2 has two parameters, received $(length(hyp))."))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2 * hyp[2])
end

get_params(mat::Mat32Iso{T}) where T = T[log(mat.ℓ), log(mat.σ2) / 2 ]
get_param_names(mat::Mat32Iso) = [:ll, :lσ]
num_params(mat::Mat32Iso) = 2

cov(mat::Mat32Iso, r::Number) =
    (s = √3 * r / mat.ℓ; mat.σ2 * (1 + s) * exp(-s))

@inline dk_dll(mat::Mat32Iso, r::Real) =
    (s = √3 * r / mat.ℓ; mat.σ2 * s^2 * exp(-s))