text:
# The simplified model of ECME-RIRR plus ROS generation by ETC ODE system
using Parameters, StaticArrays

include("common.jl")
include("cicr16.jl")
include("ck.jl")
include("force.jl")
include("ik2009.jl")
include("ina.jl")
include("inaca2009.jl")
include("ins_bg.jl")
include("katp.jl")
include("mtdna.jl")
include("oxphos-ode.jl")
include("pump.jl")
include("ros.jl")
include("tca.jl")
include("transmito.jl")
include("cabuf.jl")
include("concentrations.jl")

# Initial conditions (largely from Zhou et al. (2009))
u0_tup = (
    vm = -87.28,  # Sarcolemmal membrane potential (mV)
    dpsi = 175.318,  # Mitochondrial membrane potential (mV)
    m_na = 0.0327,  # Fast Na gating (activation)
    h_na = 0.9909,  # Fast Na gating (inactivation)
    j_na = 0.9941,  # Fast Na gating (slow inactivation)
    x_k = 1.1212e-4,  # K gating (activation)
    na_i = 8.2143,  # Cytoplasmic Na (mM)
    k_i = 150.8,  # Cytoplasmic K (mM)
    ca_i = 5.8653e-5,  # Cytoplasmic Ca (mM)
    ca_jsr = 0.1948,  # Junctional SR Ca (mM)
    ca_nsr = 0.1948,  # Network SR Ca (mM)
    ca_ss = 7.057e-5,  # Subspace calcium
    ltr_ca = 8.949E-3,  # Ca bound to low affinity troponin sites (mM)
    htr_ca = 0.1321,  # Ca bound to high affinity troponin sites (mM)
    # CICR states
    po1 = 8.309e-5,
    po2 = 7.12e-11,
    pc2 = 0.2471,
    c0 = 0.9991,
    c1 = 8.175e-5,
    c2 = 2.508e-9,
    c3 = 3.421e-14,
    c4 = 1.749e-20,
    o = 2.624e-20,
    cca0 = 1.1328e-3,
    cca1 = 1.7591e-8,
    cca2 = 6.3826e-13,
    cca3 = 1.815e-15,
    y_ca = 0.9479,
    # Tropomyosin cross bridge states
    p0 = 2.601e-5,
    p1 = 2.248e-5,
    p2 = 4.199e-5,
    p3 = 3.657e-5,
    n1 = 2.243e-5,
    adp_i = 8 - 7.8707,  # Cytosolic ATP(mM), linked to EC coupling
    adp_ic = 8 - 7.7107,  # Cytosolic ATP(mM), not linked to EC coupling
    crp_i = 5.0297,  # Cytosolic creatine phosphate(mM), linked to EC coupling
    crp_ic = 5.1291,  # Cytosolic creatine phosphate(mM), not linked to EC coupling
    # TCA cycle states
    ca_m = 2.190477831750499e-5,
    adp_m = 0.005788892127721868,
    nadh = 9.267368089681039,
    isoc = 0.4858673722382066,
    akg = 0.0001907868710833752,
    scoa = 0.28625004579156677,
    suc = 0.00011243840322880198,
    fum = 0.005275867239849815,
    mal = 0.002781903362877464,
    oaa = 2.683029248919834e-8,
    # ROS
    sox_m = 0.0004966725362784072,
    sox_i = 8.570060710246488e-9,
    h2o2_i = 2.061879987223193e-9,
    gsh_i = 0.9999899536989296,
    # ETC states
    Q_n = 1.3268400422082767,
    Qdot_n = 0.06757658200311795,
    QH2_n = 0.6002709581315435,
    QH2_p = 0.6001505595245672,
    Qdot_p = 0.07820141736951611,
    Q_p = 1.3269604407629754,
    b1 = 0.11362728970178769,
    b2 = 0.14764414671393097,
    b3 = 0.04917170786399363,
    b4 = 0.014656855720284924,
    fes_ox = 0.054765936912457895,
    cytc1_ox = 0.16423864061514593,
    cytc_ox = 0.1081945884323241)

nameLUT = Dict((k=>i) for (i, k) in enumerate(keys(u0_tup)))

# TODO: conservation of Q pool
Q_POOL = u0_tup[:Q_n] + u0_tup[:Qdot_n] + u0_tup[:QH2_n] + u0_tup[:QH2_p] + u0_tup[:Qdot_p] + u0_tup[:Q_p]
# TODO: b pool = C3 conc
B_POOL = u0_tup[:b1] + u0_tup[:b2] + u0_tup[:b3] + u0_tup[:b4]

conc = Concentrations(nadph_i=1.0,
                      pi_m=2.0,
                      mg_m=0.4,
                      mg_i=3.1,
                      AXP_M_T=1.5,
                      NAD_T=10.0,
                      AXP_I_T=8.0,
                      O2=6E-3)

@with_kw struct Params
    conc = conc
    pSODi = SODParams(K1=1.2E3, K3=24.0, K5=2.5E-4, ET=1.43E-3, KI_H2O2=0.5)
    pGPXi = GPXParams(ET=0.01, 𝚽1=5e-3, 𝚽2=0.75)
    pGRi = GRParams(K1=5E-3, ET=0.01, KM_GSSG=0.06, KM_NADPH=0.015, GSH_T=1.0)
    pCAT = CATParams(K1=17.0, ET=0.01, FR=0.05)
    pC1 = C1Params()
    pC2 = C2Params()
    pC3 = C3Params()
    pC4 = C4Params()
    pC5 = C5Params(ρF1=0.05)
    pANT = ANTParams(VMAX=5E-3)
    pMNCE = MNCEParams(VMAX=7.3E-4)
    pMCU = MCUParams(VMAX=0.028)
    pIMAC = IMACParams()
    pCS = CSParams(ACCOA=1.0, KCAT=0.5, ET=0.4, KM_ACCOA = 0.0126, KM_OAA = 6.4E-4)
    pACO = ACOParams(TCA_T=1.0, KF=1.25E-2, KEQ=2.22)
    pIDH = IDHParams(KM_ADP=0.62, KM_CA=5E-4, KI_NADH=0.19, KCAT=0.05, ET=0.109, KH1=8.1e-5, KH2=5.98E-5, KM_ISOC=1.52, NI_ISOC=2, KM_NAD=0.923)
    pKGDH = KGDHParams(KM_MG=0.0308, KM_CA=1.27E-3, ET=0.5, KCAT=0.075, KM_AKG=1.94, KM_NAD=38.7, NI_AKG=1.2)
    pSL = SLParams(KF=5E-3, KEQ=3.115, COA=0.02)
    pFH = FHParams(KF=3.32E-3, KEQ=1.0)
    pMDH = MDHParams(KH1=1.131E-5, KH2=26.7, KH3=6.68E-9, KH4=5.62E-6, K_OFFSET=3.99E-2, KCAT=0.11, ET=0.15, KM_MAL=1.49, KI_OAA=3.1E-3, KM_NAD=0.22)
    pAAT = AATParams(KF=6.44E-4, KEQ=6.6, KASP=1.5E-6, GLU=10.0)
    pDOX = DOXParams()
    pMTDNA = MTDNAParams()
    # G_H_MITO = 2E-6  # Mitochondrial membrane conductance (Gauthier 2013)
    G_H_MITO = 1E-8  # (Li 2015)
    δCA = 3E-4  # Ca buffering factor
    R_TR= inv(9.09) # Rate of Ca diffusion from NSR to JSR (1/ms)
    R_XFER = inv(0.5747)  # Rate of Ca diffusion from subspace to cytosol (1/ms)
    G_NA = 12.8  # Sodium conductance
    V_NSR = 1.4E-6  # Network SR Volume (μL)
    V_JSR = 0.16E-6  # Junctional SR Volume (μL)
    V_SS = 0.495E-9  # Total subspace Volume (μL)
    V_MYO = 25.84E-6  # Cytoplasm Volume (μL)
    V_MITO = 15.89E-6  # Mitochiondrial Volume (μL)
    A_CAP = 1.534E-4  # Cell capacitance area (cm²)
    A_CAP_V_MYO_F = A_CAP / (V_MYO * F)
    V_MITO_V_MYO = V_MITO / V_MYO
    pIK = IKParams()
    pCK = CKParams()
    pXB = XBParams()
    pLCC = LCCParams()
    pRYR = RYRParams()
    pSERCA = SERCAParams(KMF_CA=2.4e-4, KMR_CA=1.64269, NFB=1.4, NRB=1)
    pPMCA = PMCAParams()
    pNS = NSParams()
    pNKA = NKAParams()
    pNCX = NCXParams()
    pKATP = KATPParams()
    pCSQN = CABUFParams(KM_CA=0.8, ET=5.0)
    pCMDN = CABUFParams(KM_CA=2.38E-3, ET=5E-2)
    pBG = BGParams()
    BCL = 0
end

function _istim(t, bcl=0, peroid=0.5, strength=-80.0)
    if (bcl > 0) && (rem(t, bcl) <= peroid)
        return strength
    else
        return zero(strength)
    end
end

function ecme_dox(u, p, t)
    (vm, dpsi, m_na, h_na, j_na, x_k, na_i, k_i, ca_i, ca_jsr, ca_nsr, ca_ss, ltr_ca, htr_ca,
     po1, po2, pc2, c0, c1, c2, c3, c4, o, cca0, cca1, cca2, cca3, y_ca, p0, p1, p2, p3, n1,
     adp_i, adp_ic, crp_i, crp_ic, ca_m, adp_m, nadh, isoc, akg, scoa, suc, fum, mal, oaa,
     sox_m, sox_i, h2o2_i, gsh_i, Q_n, Qdot_n, QH2_n, QH2_p, Qdot_p, Q_p, b1, b2, b3, b4,
     fes_ox, cytc1_ox, cytc_ox) = u

    @unpack conc = p
	@unpack ca_o, na_o, k_o, AXP_I_T, h_i, h_m, mg_i, mg_m, ΔpH, FAC_PH, O2, nadph_i, pi_m, AXP_I_T = conc
    vfrt = vm * F_RT
    evfrtm1 = expm1(vfrt)
	ΔμH = _ΔμH(dpsi, ΔpH)  # proton motive force
    atp_i = AXP_I_T - adp_i
    atp_ic = AXP_I_T - adp_ic

    # CICR
    @unpack pLCC, pRYR, R_TR, R_XFER = p
	(d_c0, d_c1, d_c2, d_c3, d_c4, d_o, d_cca0, d_cca1, d_cca2,
     d_cca3) = lcc_system(c0, c1, c2, c3, c4, o, cca0, cca1, cca2, cca3, ca_ss, vm, pLCC)
    d_y_ca = dyca(y_ca, vm)
    (d_po1, d_po2, d_pc2) = ryr_system(po1, po2, pc2, ca_ss, pRYR)
    (iCaL, iCaK) = ical(y_ca, o, vfrt, evfrtm1, k_i, k_o, ca_o, pLCC)
    jRel = jrel(po1, po2, ca_jsr, ca_ss, pRYR.R_RYR)
    jTr = R_TR * (ca_nsr - ca_jsr)
    jXfer = R_XFER * (ca_ss - ca_i)
    # Myofibril
    @unpack pCK, pXB = p
    (d_p0, d_p1, d_p2, d_p3, d_n1, d_ltr_ca, d_htr_ca) = xb_system(p0, p1, p2, p3, n1, ltr_ca, htr_ca, ca_i, pXB)
	vAM = vam(p0, p1, p2, atp_i, adp_i, pXB)
    jTrpn = d_ltr_ca + d_htr_ca

    # CK shuttle
    (d_atp_in_ck, d_atp_ic, d_crp_i, d_crp_ic) = ck_system(atp_i, atp_ic, crp_i, crp_ic, pCK, AXP_I_T)

    # Other ion channels
    @unpack G_CAB, G_NAB = p.pBG
    @unpack G_KP, P_NA_K, G_K, G_K1 = p.pIK
    @unpack pNKA, pPMCA, pSERCA, pKATP, G_NA, pNS, pNCX = p
    d_m_na = d_mna(m_na, vm)
    d_h_na = d_hna(h_na, vm)
    d_j_na = d_jna(j_na, vm)
    eNa = _ena(na_o, na_i)
    iNa = ina(vm, eNa, m_na, h_na, j_na, G_NA)
    d_x_k = dxk(x_k, vm)
    iK = ik(vm, _enak(k_o, na_o, k_i, na_i, P_NA_K), x_k, G_K)
	eK = _ek(k_o, k_i)
    iKp = ikp(vm, eK, G_KP)
    iK1 = ik1(vm, eK, G_K1)
	iKatp = ikatp(vm, eK, na_i, atp_i, adp_i, mg_i, pKATP)
    iNsNa = insna(na_i, ca_i, vfrt, evfrtm1, na_o, pNS)
    iNaCa = inaca(na_i, ca_i, vfrt, evfrtm1, ca_o, pNCX)
    iCaB = icab(vm, _eca(ca_o, ca_i), G_CAB)
    iNaB = inab(vm, eNa, G_NAB)
    iNaK = inak(na_i, atp_i, adp_i, vm, evfrtm1, pNKA)
    iPCa = ipca(ca_i, atp_i, adp_i, pPMCA)
    jUp = jup(ca_i, ca_nsr, atp_i, adp_i, pSERCA)

    # Ion fluxes across mitochondria
    @unpack pMCU, G_H_MITO, pMNCE = p
    vUni = vuni(ca_i, dpsi, pMCU)
    vNaCa = vnaca(ca_i, ca_m, na_i, dpsi, pMNCE, pMCU)
    vHleak = vhleak(ΔμH, G_H_MITO)

    # ETC
    @unpack pC1, pC2, pC3, pC4, pC5, pDOX, pANT = p
    @unpack C2_INHIB, C3_INHIB, C4_INHIB = pDOX
    @unpack MT_PROT = p.pMTDNA
    vANT = vant(adp_m, atp_i, dpsi, conc, pANT)
    vc1, vROSC1 = c1_rates(nadh, Q_n, QH2_n, dpsi, h_m, sox_m, conc, pC1, pDOX, MT_PROT)
    vc2 = vSDH = vsdh_mm(Q_n, QH2_n, suc, fum, oaa, pC2, C2_INHIB)
    vE, vHresC4, vO2 = c4_rates(cytc_ox, dpsi, h_m, conc, pC4, C4_INHIB, MT_PROT)
    (vATPase, vHu) = c5_rates(ΔμH, pi_m, adp_m, conc, pC5, MT_PROT)
    # Complex III rates
    @unpack (KD_Q, K03, KEQ3, α, δ₁, K04, KEQ4_OX, KEQ4_RD, β, δ₂, KEQ6, K06, γ, δ₃,
        KEQ7_OX, K07_OX, KEQ7_RD, K07_RD, KEQ8_OX, K08_OX, KEQ8_RD, K08_RD, K09,
        KEQ9, K010, KEQ10, K33, KEQ33, ρC3) = pC3
    C3_CONC = MT_PROT * ρC3
	C4_CONC = MT_PROT * pC4.ρC4
    fes_rd = C3_CONC - fes_ox
    cytc1_rd = C3_CONC - cytc1_ox
	cytc_rd = C4_CONC - cytc_ox
	MYXOTHIAZOLE = ANTIMYCIN = 1

	# v1 = Q reduction
	v1 = vc1 + vc2
	# v2 = QH2 diffusion (n-side -> p-side)
    v2 = KD_Q * (QH2_n - QH2_p)
    # v3 = QH2 to FeS
	v3 = _v3(fes_ox, QH2_p, fes_rd, Qdot_p, K03, KEQ3, FAC_PH, MYXOTHIAZOLE)
	# v4 = Qdot_p to bH
	el = _el(α, δ₁, dpsi)
	er = _er(α, δ₁, dpsi)
	v4_ox = _v4(b1, Qdot_p, el, b2, Q_p, er, K04, KEQ4_OX)
	v4_rd = _v4(b3, Qdot_p, el, b4, Q_p, er, K04, KEQ4_RD)
    # v5 = Q diffusion (p-side -> n-side)
    v5 = KD_Q * (Q_p - Q_n)
    # v6 = bH to bL
	v6 = _v6(b2, b3, dpsi, K06, KEQ6, β, δ₂)
    # v7 = bL to Qn; v8: bL to Qdot_n
	el = _el(γ, δ₃, dpsi)
	er = _er(γ, δ₃, dpsi)
	v7_ox = _v7(b3, Q_n, el, b1, Qdot_n, er, K07_OX, KEQ7_OX, ANTIMYCIN) * C3_INHIB
	v7_rd = _v7(b4, Q_n, el, b2, Qdot_n, er, K07_RD, KEQ7_RD, ANTIMYCIN) * C3_INHIB
	# v8 = bL to Qdot_n and proton transfer
	v8_ox = _v8(b3, Qdot_n, el, b1, QH2_n, er, K08_OX, KEQ8_OX, FAC_PH, ANTIMYCIN) * C3_INHIB
	v8_rd = _v8(b4, Qdot_n, el, b2, QH2_n, er, K08_RD, KEQ8_RD, FAC_PH, ANTIMYCIN) * C3_INHIB
    # v9 = fes -> cytc1
    v9 = _v9(fes_rd, cytc1_ox, fes_ox, cytc1_rd, K09, KEQ9)
    # v10 = Qdot + O2 -> O2- + Q  (ROS produced by complex III)
    vROSC3 = v10 = _v10(O2, Qdot_p, sox_m, Q_p, K010, KEQ10)

	# Redox reaction between cytc1 and cytc
	v33 = _v3(cytc1_rd, cytc_ox, cytc1_ox, cytc_rd, K33, KEQ33)

    # ODE system of complex III
    d_Q_n = v5 - v7_ox - v7_rd - v1
    d_Qdot_n = v7_ox + v7_rd - v8_ox - v8_rd
    d_QH2_n = v8_ox + v8_rd - v2 + v1
    d_QH2_p = v2 - v3
    d_Qdot_p = v3 - v10 - v4_ox - v4_rd
    d_Q_p = v10 + v4_ox + v4_rd - v5
    d_b1 = v7_ox + v8_ox - v4_ox
    d_b2 = v4_ox + v7_rd + v8_rd - v6
    d_b3 = v6 - v4_rd - v7_ox - v8_ox
    d_b4 = v4_rd - v7_rd - v8_rd
    d_fes_ox = v9 - v3
    d_cytc1_ox = v33 - v9
	d_cytc_ox = vE - v33

	vROS = vROSC1 + vROSC3
    vHres = 4 * vc1 + 2 * v3 + vHresC4

    # TCA cycle rates
    @unpack pSL, pCS, pACO, pIDH, pKGDH, pFH, pMDH, pAAT = p
    nad = _nad(nadh, conc)
	vCS  = vcs(oaa, pCS)
	vACO = vaco(isoc, akg, scoa, suc, fum, mal, oaa, pACO)
	vIDH = vidh(isoc, nadh, nad, adp_m, ca_m, h_m, mg_m, pIDH)
	vKGDH = vkgdh(akg, nadh, nad, ca_m, mg_m, pKGDH)
	vSL = vsl(scoa, suc, adp_m, _atp_m(adp_m, conc), pSL)
	vFH = vfh(fum, mal, pFH)
	vMDH = vmdh(mal, oaa, nadh, nad, h_m, pMDH)
	vAAT = vaat(akg, oaa, pAAT)

    # ODE system of TCA cycle
	d_nadh = -vc1 + vIDH + vKGDH + vMDH
	d_isoc = vACO - vIDH
	d_akg = vIDH + vAAT - vKGDH
	d_scoa = vKGDH - vSL
	d_suc = vSL - vSDH
	d_fum = vSDH - vFH
	d_mal = vFH - vMDH
	d_oaa = vMDH - vCS - vAAT
    d_adp_m = vANT - vSL - vATPase

    @unpack (pSODi, pGPXi, pGRi, pCAT, pIMAC, V_MITO_V_MYO) = p
    vSOD_i = vsod(sox_i, h2o2_i, pSODi)
    vGPX_i = vgpx(h2o2_i, gsh_i, pGPXi)
    vGR_i = vgr(gsh_i, nadph_i, pGRi)
    vCAT = v_cat(h2o2_i, pCAT)
    vIMAC, vTrROS = vimac_vtrros(dpsi, sox_m, sox_i, pIMAC)

    # ODE system of ROS scavenging
    d_sox_m = vROS - vTrROS
    d_sox_i = V_MITO_V_MYO * vTrROS - vSOD_i
	d_h2o2_i = 0.5 * vSOD_i - vGPX_i - vCAT
	d_gsh_i  = vGR_i - vGPX_i

	# Remaining ODE equations
	@unpack (δCA, A_CAP_V_MYO_F, V_MYO, V_SS, V_NSR, V_JSR, V_MITO_V_MYO, BCL) = p
	iStim = _istim(t, BCL)
	d_vm = -CM_INV * (iNa + iCaL + iK + iK1 + iKp + iNaCa + iNaK + iNsNa + iPCa + iCaB + iNaB + iStim + iKatp + iCaK)
	d_na_i = -A_CAP_V_MYO_F * (iNa + iNaB + iNsNa + 3 * (iNaCa + iNaK)) # + V_MITO_V_MYO * (vNHE - 3 * vNaCa)
    d_k_i = -A_CAP_V_MYO_F * (iK + iK1 + iKp + iKatp + iStim - 2 * iNaK + iCaK)
	d_atp_i = d_atp_in_ck + V_MITO_V_MYO * vANT - vAM - 0.5 * jUp - A_CAP_V_MYO_F * (iPCa + iNaK)

	@unpack pCSQN, pCMDN, A_CAP = p
	d_ca_i = _β_ca(ca_i, pCMDN) * (jXfer - jUp - jTrpn - 0.5 * A_CAP_V_MYO_F * (iPCa + iCaB - 2 * iNaCa) + V_MITO_V_MYO * (vNaCa - vUni))
	d_ca_nsr = _β_ca(ca_nsr, pCSQN) * (V_MYO * jUp - V_JSR * jTr) / V_NSR
    d_ca_jsr = _β_ca(ca_jsr, pCSQN) * (jTr - jRel)
    d_ca_ss = _β_ca(ca_ss, pCMDN) * (V_JSR * jRel - V_MYO * jXfer - 0.5 * A_CAP / F * iCaL) / V_SS
    d_ca_m = δCA * (vUni - vNaCa)

    d_dpsi = CM_MITO_INV * (vHres - vHu - vANT - vHleak - vNaCa - 2 * vUni - vIMAC)
    # @show t, na_i, m_na, h_na, j_na, d_m_na, d_h_na, d_j_na, iNa

	# Disable rates for debugging
    #=
	d_vm = 0
	d_m_na = d_h_na = d_j_na = d_na_i = 0
	d_x_k = d_k_i = 0
	d_p0 = d_p1 = d_p2 = d_p3 = d_n1 = 0
	d_ca_i = d_ca_jsr = d_ca_nsr = d_ca_ss = d_ltr_ca = d_htr_ca = 0
	d_atp_i = d_atp_ic = d_crp_i = d_crp_ic = 0
    d_po1 = d_po2 = d_pc2 = d_c0 = d_c1 = d_c2 = d_c3 = d_c4 = d_o = d_cca0 = d_cca1 = d_cca2 = d_cca3 = d_y_ca = 0
    =#

    d_adp_i = -d_atp_i
    d_adp_ic = -d_atp_ic
    du = @SVector[d_vm, d_dpsi, d_m_na, d_h_na, d_j_na, d_x_k, d_na_i, d_k_i, d_ca_i, d_ca_jsr, d_ca_nsr, d_ca_ss, d_ltr_ca, d_htr_ca,
     d_po1, d_po2, d_pc2, d_c0, d_c1, d_c2, d_c3, d_c4, d_o, d_cca0, d_cca1, d_cca2, d_cca3, d_y_ca, d_p0, d_p1, d_p2, d_p3, d_n1,
     d_adp_i, d_adp_ic, d_crp_i, d_crp_ic, d_ca_m, d_adp_m, d_nadh, d_isoc, d_akg, d_scoa, d_suc, d_fum, d_mal, d_oaa,
     d_sox_m, d_sox_i, d_h2o2_i, d_gsh_i, d_Q_n, d_Qdot_n, d_QH2_n, d_QH2_p, d_Qdot_p, d_Q_p, d_b1, d_b2, d_b3, d_b4,
     d_fes_ox, d_cytc1_ox, d_cytc_ox,
     iCaL, iCaK, jRel, jTr, jXfer, vAM, jTrpn, iNa, iK, iKp, iK1, iKatp, iNsNa, iNaCa, iCaB,
     iNaB, iNaK, iPCa, jUp, vUni, vNaCa, vHleak, vANT, vc1, vROSC1, vSDH, vE, vHresC4, vO2,
     vATPase, vHu, vROSC3, vHres, vCS, vACO, vIDH, vKGDH, vSL, vFH, vMDH, vAAT, vSOD_i, vGPX_i,
     vGR_i, vCAT, vIMAC, vTrROS]
end

rateSymbols = [:d_vm, :d_dpsi, :d_m_na, :d_h_na, :d_j_na, :d_x_k, :d_na_i, :d_k_i, :d_ca_i, :d_ca_jsr, :d_ca_nsr, :d_ca_ss, :d_ltr_ca, :d_htr_ca,
 :d_po1, :d_po2, :d_pc2, :d_c0, :d_c1, :d_c2, :d_c3, :d_c4, :d_o, :d_cca0, :d_cca1, :d_cca2, :d_cca3, :d_y_ca, :d_p0, :d_p1, :d_p2, :d_p3, :d_n1,
 :d_adp_i, :d_adp_ic, :d_crp_i, :d_crp_ic, :d_ca_m, :d_adp_m, :d_nadh, :d_isoc, :d_akg, :d_scoa, :d_suc, :d_fum, :d_mal, :d_oaa,
 :d_sox_m, :d_sox_i, :d_h2o2_i, :d_gsh_i, :d_Q_n, :d_Qdot_n, :d_QH2_n, :d_QH2_p, :d_Qdot_p, :d_Q_p, :d_b1, :d_b2, :d_b3, :d_b4,
 :d_fes_ox, :d_cytc1_ox, :d_cytc_ox,
 :iCaL, :iCaK, :jRel, :jTr, :jXfer, :vAM, :jTrpn, :iNa, :iK, :iKp, :iK1, :iKatp, :iNsNa, :iNaCa, :iCaB,
 :iNaB, :iNaK, :iPCa, :jUp, :vUni, :vNaCa, :vHleak, :vANT, :vc1, :vROSC1, :vSDH, :vE, :vHresC4, :vO2,
 :vATPase, :vHu, :vROSC3, :vHres, :vCS, :vACO, :vIDH, :vKGDH, :vSL, :vFH, :vMDH, :vAAT, :vSOD_i, :vGPX_i,
 :vGR_i, :vCAT, :vIMAC, :vTrROS]

rateMap = Dict((k => i) for (i, k) in enumerate(rateSymbols))

# Ignore the intermediate rates
rhs(u, p, t) = ecme_dox(u, p, t)[1:length(u)]

reject_step(u, p, t) = any( x -> (x < 0.0), @view u[3:end])


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("World of machine learning")

# # Image classification with very deep convolutional networks

using Knet
cd(datapath("vgg"))
include("vgg.jl")

#-

@doc VGG

#-

VGG.main("--help")

#-

using FileIO
img = download("http://home.mweb.co.za/pa/pak04857/uniweb/animalimages/elephantthumb.jpg")
load(img)

#-

VGG.main(img)

#-

ext: jl

max_stars_repo_path: Courses/World of machine learning/vgg.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
"""
    make()
---
This is to make binaries from my C++ files.
It is supposed to run at the startup.
But thic can also be run after new C++ codes were written.
"""
function make()
    title("C++ binaries")
    isdir(rdBin) || mkdir(rdBin)
    bins = ["find-dup",
            "merge-dup",
            "split-h",
            "split-v",
            "raw2gt",
            "pedsort",
            "merge-ped"]
    width = 60
    for bin in bins
        exe = joinpath(rdBin, bin)
        cpp = joinpath(rdCpp, "$bin.cpp")
        if !isfile(exe) || stat(exe).mtime < stat(cpp).mtime
            msg = "g++ -O2 -Wall $bin.cpp -o $bin"
            print(lpad(msg, width))
            run(`g++ -O2 -Wall -o $exe $cpp`)
            done()
        else
            print(lpad(bin, width))
            done("OK")
        end
    end
end


text:
toargs(vs :: Vector{Symbol}) = Tuple(vs)
toargs(vs :: NTuple{N,Symbol} where {N}) = vs

macro model(vs::Expr,expr::Expr)
    theModule = __module__
    @assert vs.head == :tuple
    @assert expr.head == :block
    Model(theModule,Vector{Symbol}(vs.args), expr)
end

macro model(v::Symbol, expr::Expr)
    theModule = __module__
    Model(theModule,[v], expr)
end

macro model(expr :: Expr)
    theModule = __module__
    Model(theModule,Vector{Symbol}(), expr)
end


text:
module PooledArrays

export PooledArray

##############################################################################
##
## PooledArray type definition
##
##############################################################################

const DEFAULT_POOLED_REF_TYPE = UInt32

# This is used as a wrapper during PooledArray construction only, to distinguish
# arrays of pool indices from normal arrays
mutable struct RefArray{R}
    a::R
end

function _invert(d::Dict{K,V}) where {K,V}
    d1 = Dict{V,K}()
    for (k, v) in d
        d1[v] = k
    end
    d1
end

mutable struct PooledArray{T, R<:Integer, N, RA} <: AbstractArray{T, N}
    refs::RA
    pool::Dict{T,R}
    revpool::Dict{R,T}

    function PooledArray(rs::RefArray{RA}, p::Dict{T, R}, revpool=_invert(p)) where {T,R,N,RA<:AbstractArray{R, N}}
        # refs mustn't overflow pool
        if length(rs.a) > 0 && maximum(rs.a) > length(p)
            throw(ArgumentError("Reference array points beyond the end of the pool"))
        end
        new{T,R,N,RA}(rs.a,p,revpool)
    end
end
const PooledVector{T,R} = PooledArray{T,R,1}
const PooledMatrix{T,R} = PooledArray{T,R,2}

##############################################################################
##
## PooledArray constructors
##
# Algorithm:
# * Start with:
#   * A null pool
#   * A pre-allocated refs
#   * A hash from T to Int
# * Iterate over d
#   * If value of d in pool already, set the refs accordingly
#   * If value is new, add it to the pool, then set refs
##############################################################################

# Echo inner constructor as an outer constructor
function PooledArray(refs::RefArray{R}, pool::Dict{T,R}) where {T,R}
    PooledArray{T,eltype(R),ndims(R),R}(refs, pool)
end

# A no-op constructor
PooledArray(d::PooledArray) = d

function _label(xs::AbstractArray{T},
                ::Type{I}=UInt8,
                start = 1,
                labels = Array{I}(undef, size(xs)),
                pool::Dict{T,I} = Dict{T, I}(),
                nlabels = 0,
               ) where {T, I<:Integer}

    @inbounds for i in start:length(xs)
        x = xs[i]
        lbl = get(pool, x, zero(I))
        if lbl !== zero(I)
            labels[i] = lbl
        else
            if nlabels == typemax(I)
                I2 = _widen(I)
                return _label(xs, I2, i, convert(Vector{I2}, labels),
                              convert(Dict{T, I2}, pool), nlabels)
            end
            nlabels += 1
            labels[i] = convert(I, nlabels)
            pool[x] = convert(I, nlabels)
        end
    end
    labels, pool
end

_widen(::Type{UInt8}) = UInt16
_widen(::Type{UInt16}) = UInt32
_widen(::Type{UInt32}) = UInt64

# Constructor from array, pool, and ref type

"""
    PooledArray(array, [reftype])

Convert the given array to a PooledArray where each element will be referenced
as an integer of the given type. If no `reftype` is specified one is chosen
automatically based on the number of unique elements.
"""
PooledArray

function PooledArray{T}(d::AbstractArray, r::Type{R}) where {T,R<:Integer}
    refs, pool = _label(d)

    if length(pool) > typemax(R)
        throw(ArgumentError("Cannot construct a PooledArray with type $R with a pool of size $(length(pool))"))
    end

    refs1 = convert(Vector{R}, refs)
    pool1 = convert(Dict{T,R}, pool)
    PooledArray(RefArray(refs1), pool1)
end

function PooledArray{T}(d::AbstractArray) where T
    refs, pool = _label(d)
    PooledArray(RefArray(refs), pool)
end

PooledArray(d::AbstractArray{T}, r::Type{R}) where {T,R<:Integer} = PooledArray{T}(d, r)
PooledArray(d::AbstractArray{T}) where {T} = PooledArray{T}(d)

# Construct an empty PooledVector of a specific type
PooledArray(t::Type) = PooledArray(Array(t,0))
PooledArray(t::Type, r::Type{R}) where {R<:Integer} = PooledArray(Array(t,0), r)

##############################################################################
##
## Basic interface functions
##
##############################################################################

Base.size(pa::PooledArray) = size(pa.refs)
Base.length(pa::PooledArray) = length(pa.refs)
Base.lastindex(pa::PooledArray) = lastindex(pa.refs)

Base.copy(pa::PooledArray) = PooledArray(RefArray(copy(pa.refs)), copy(pa.pool))
# TODO: Implement copy_to()

function Base.resize!(pa::PooledArray{T,R,1}, n::Integer) where {T,R}
    oldn = length(pa.refs)
    resize!(pa.refs, n)
    pa.refs[oldn+1:n] .= zero(R)
    pa
end

Base.reverse(x::PooledArray) = PooledArray(RefArray(reverse(x.refs)), x.pool)

function Base.permute!!(x::PooledArray, p::AbstractVector{T}) where T<:Integer
    Base.permute!!(x.refs, p)
    x
end

function Base.invpermute!!(x::PooledArray, p::AbstractVector{T}) where T<:Integer
    Base.invpermute!!(x.refs, p)
    x
end

function Base.similar(pa::PooledArray{T,R}, S::Type, dims::Dims) where {T,R}
    PooledArray(RefArray(zeros(R, dims)), Dict{S,R}())
end

Base.findall(pdv::PooledVector{Bool}) = findall(convert(Vector{Bool}, pdv))

##############################################################################
##
## map
## Calls `f` only once per pool entry.
##
##############################################################################

function Base.map(f, x::PooledArray{T,R}) where {T,R<:Integer}
    ks = collect(keys(x.pool))
    vs = collect(values(x.pool))
    ks1 = map(f, ks)
    uks = Set(ks1)
    if length(uks) < length(ks1)
        # this means some keys have repeated
        newpool = Dict{eltype(ks1), eltype(vs)}()
        translate = Dict{eltype(vs), eltype(vs)}()
        i = 1
        for (k, k1) in zip(ks, ks1)
            if haskey(newpool, k1)
                translate[x.pool[k]] = newpool[k1]
            else
                newpool[k1] = i
                translate[x.pool[k]] = i
                i+=1
            end
        end
        refarray = map(x->translate[x], x.refs)
    else
        newpool = Dict(zip(map(f, ks), vs))
        refarray = copy(x.refs)
    end
    PooledArray(RefArray(refarray), newpool)
end

##############################################################################
##
## Sorting can use the pool to speed things up
##
##############################################################################

function groupsort_indexer(x::AbstractVector, ngroups::Integer, perm)
    # translated from Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).

    # count group sizes, location 0 for NA
    n = length(x)
    # counts = x.pool
    counts = fill(0, ngroups + 1)
    @inbounds for i = 1:n
        counts[x[i] + 1] += 1
    end
    counts[2:end] = counts[perm.+1]

    # mark the start of each contiguous group of like-indexed data
    where = fill(1, ngroups + 1)
    @inbounds for i = 2:ngroups+1
        where[i] = where[i - 1] + counts[i - 1]
    end

    # this is our indexer
    result = fill(0, n)
    iperm = invperm(perm)

    @inbounds for i = 1:n
        label = iperm[x[i]] + 1
        result[where[label]] = i
        where[label] += 1
    end
    result, where, counts
end

function Base.sortperm(pa::PooledArray; alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                       lt::Function=isless, by::Function=identity,
                       rev::Bool=false, order=Base.Sort.Forward,
                       _ord = Base.ord(lt, by, rev, order),
                       poolperm = sortperm([pa.revpool[i] for i=1:length(pa.pool)], alg=alg, order=_ord))

    groupsort_indexer(pa.refs, length(pa.pool), poolperm)[1]
end

Base.sort(pa::PooledArray; kw...) = pa[sortperm(pa; kw...)]

#type FastPerm{O<:Base.Sort.Ordering,V<:AbstractVector} <: Base.Sort.Ordering
#    ord::O
#    vec::V
#end
#Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ForwardOrdering,V}) = x[sortperm(o.vec)]
#Base.sortperm{V}(x::AbstractVector, a::Base.Sort.Algorithm, o::FastPerm{Base.Sort.ReverseOrdering,V}) = x[reverse(sortperm(o.vec))]
#Perm{O<:Base.Sort.Ordering}(o::O, v::PooledVector) = FastPerm(o, v)

##############################################################################
##
## conversions
##
##############################################################################

Base.convert(::Type{PooledArray{S,R1,N}}, pa::PooledArray{T,R2,N}) where {S,T,R1<:Integer,R2<:Integer,N} =
    PooledArray(RefArray(convert(Array{R1,N}, pa.refs)), convert(Dict{S,R1}, pa.pool))
Base.convert(::Type{PooledArray{S,R,N}}, pa::PooledArray{T,R,N}) where {S,T,R<:Integer,N} =
    PooledArray(RefArray(copy(pa.refs)), convert(Dict{S,R}, pa.pool))
Base.convert(::Type{PooledArray{T,R,N}}, pa::PooledArray{T,R,N}) where {T,R<:Integer,N} = pa
Base.convert(::Type{PooledArray{S,R1}}, pa::PooledArray{T,R2,N}) where {S,T,R1<:Integer,R2<:Integer,N} =
    convert(PooledArray{S,R1,N}, pa)
Base.convert(::Type{PooledArray{S}}, pa::PooledArray{T,R,N}) where {S,T,R<:Integer,N} =
    convert(PooledArray{S,R,N}, pa)
Base.convert(::Type{PooledArray}, pa::PooledArray{T,R,N}) where {T,R<:Integer,N} = pa

Base.convert(::Type{PooledArray{S,R,N}}, a::AbstractArray{T,N}) where {S,T,R<:Integer,N} =
    PooledArray(convert(Array{S,N}, a), R)
Base.convert(::Type{PooledArray{S,R}}, a::AbstractArray{T,N}) where {S,T,R<:Integer,N} =
    PooledArray(convert(Array{S,N}, a), R)
Base.convert(::Type{PooledArray{S}}, a::AbstractArray{T,N}) where {S,T,N} =
    PooledArray(convert(Array{S,N}, a))
Base.convert(::Type{PooledArray}, a::AbstractArray) =
    PooledArray(a)

function Base.convert(::Type{Array{S, N}}, pa::PooledArray{T, R, N}) where {S, T, R, N}
    res = Array{S}(undef, size(pa))
    for i in 1:length(pa)
        if pa.refs[i] != 0
            res[i] = pa.revpool[pa.refs[i]]
        end
    end
    return res
end

Base.convert(::Type{Vector}, pv::PooledVector{T, R}) where {T, R} = convert(Array{T, 1}, pv)

Base.convert(::Type{Matrix}, pm::PooledMatrix{T, R}) where {T, R} = convert(Array{T, 2}, pm)

Base.convert(::Type{Array}, pa::PooledArray{T, R, N}) where {T, R, N} = convert(Array{T, N}, pa)

##############################################################################
##
## indexing
##
##############################################################################

# Scalar case
function Base.getindex(pa::PooledArray, I::Integer...)
    idx = pa.refs[I...]
    iszero(idx) && throw(UndefRefError())
    return pa.revpool[idx]
end

function Base.isassigned(pa::PooledArray, I::Int...)
    !iszero(pa.refs[I...])
end

# Vector case
function Base.getindex(A::PooledArray, I::Union{Real,AbstractVector}...)
    PooledArray(RefArray(getindex(A.refs, I...)), copy(A.pool))
end

# Dispatch our implementation for these cases instead of Base
Base.getindex(A::PooledArray, I::AbstractVector) =
    PooledArray(RefArray(getindex(A.refs, I)), copy(A.pool))
Base.getindex(A::PooledArray, I::AbstractArray) =
    PooledArray(RefArray(getindex(A.refs, I)), copy(A.pool))

##############################################################################
##
## setindex!() definitions
##
##############################################################################

function getpoolidx(pa::PooledArray{T,R}, val::Any) where {T,R}
    val::T = convert(T,val)
    pool_idx = get(pa.pool, val, zero(R))
    if pool_idx == zero(R)
        pool_idx = unsafe_pool_push!(pa, val)
    end
    return pool_idx
end

function unsafe_pool_push!(pa::PooledArray{T,R}, val) where {T,R}
    _pool_idx = length(pa.pool)+1
    if _pool_idx > typemax(R)
        throw(ErrorException(string(
            "You're using a PooledArray with ref type $R, which can only hold $(Int(typemax(R))) values,\n",
            "and you just tried to add the $(typemax(R)+1)th reference.  Please change the ref type\n",
            "to a larger int type, or use the default ref type ($DEFAULT_POOLED_REF_TYPE)."
           )))
    end
    pool_idx = convert(R, _pool_idx)
    pa.pool[val] = pool_idx
    pa.revpool[pool_idx] = val
    pool_idx
end

function Base.setindex!(x::PooledArray, val, ind::Integer)
    x.refs[ind] = getpoolidx(x, val)
    return x
end

##############################################################################
##
## growing and shrinking
##
##############################################################################

function Base.push!(pv::PooledVector{S,R}, v::T) where {S,R,T}
    v = convert(S,v)
    push!(pv.refs, getpoolidx(pv, v))
    return v
end

Base.pop!(pv::PooledVector) = pv.pool[pop!(pv.refs)]

function Base.pushfirst!(pv::PooledVector{S,R}, v::T) where {S,R,T}
    v = convert(S,v)
    pushfirst!(pv.refs, getpoolidx(pv, v))
    return v
end

Base.popfirst!(pv::PooledVector) = pv.pool[popfirst!(pv.refs)]

Base.empty!(pv::PooledVector) = (empty!(pv.refs); pv)

function _vcat!(c,a,b)
    copyto!(c, 1, a, 1, length(a))
    copyto!(c, length(a)+1, b, 1, length(b))
end

function Base.vcat(a::PooledArray{<:Any, <:Integer, 1}, b::AbstractArray{<:Any, 1})
    output = similar(b, promote_type(eltype(a), eltype(b)), length(b) + length(a))
    _vcat!(output, a, b)
end

function Base.vcat(a::AbstractArray{<:Any, 1}, b::PooledArray{<:Any, <:Integer, 1})
    output = similar(a, promote_type(eltype(a), eltype(b)), length(b) + length(a))
    _vcat!(output, a, b)
end

function Base.vcat(a::PooledArray{T, <:Integer, 1}, b::PooledArray{S, <:Integer, 1}) where {T, S}
    ap = a.pool
    bp = b.pool

    U = promote_type(T,S)

    poolmap = Dict{Int, Int}()
    l = length(ap)
    newlabels = Dict{U, Int}(ap)
    for (x, i) in bp
        j = if x in keys(ap)
            poolmap[i] = ap[x]
        else
            poolmap[i] = (l+=1)
        end
        newlabels[x] = j
    end
    types = [UInt8, UInt16, UInt32, UInt64]
    tidx = findfirst(t->l < typemax(t), types)
    refT = types[tidx]
    refs2 = map(r->convert(refT, poolmap[r]), b.refs)
    newrefs = Base.typed_vcat(refT, a.refs, refs2)
    return PooledArray(RefArray(newrefs), convert(Dict{U, refT}, newlabels))
end

end


text:
@inline function available_systems(graph::PharmGraph)
    all_systems_set = Set{String}()

    function f(class)
        system = class.system::String
        push!(all_systems_set, system)
        return nothing
    end

    _for_each_node(f, graph)

    all_systems_vector = collect(all_systems_set)
    unique!(all_systems_vector)
    sort!(all_systems_vector)
    return all_systems_vector
end

@inline function system_matches(class::PharmClass, system::String)
    class_system = class.system::String
    return class_system == system
end

@inline function system_matches(class::PharmClass, system::Regex)
    class_system = class.system::String
    return occursin(system, class_system)
end

@inline function system_matches(system::String)
     function _system_matches(class::PharmClass)
         class_system = class.system::String
         return class_system == system
     end
     return _system_matches
end

@inline function system_matches(system::Regex)
    function _system_matches(class::PharmClass)
        class_system = class.system::String
        return occursin(system, class_system)
    end
    return _system_matches
end


text:
using LowLevelParticleFilters, ForwardDiff, Distributions
const LLPF = LowLevelParticleFilters
eye(n) = Matrix{Float64}(I,n,n)
n = 2 # Dinemsion of state
m = 2 # Dinemsion of input
p = 2 # Dinemsion of measurements

# Define random linenar state-space system
A = SMatrix{n,n}([0.99 0.1; 0 0.2])
B = @SMatrix randn(n,m)
C = SMatrix{p,p}(eye(p))
# C = SMatrix{p,n}([1 1])

dynamics(x,u,t) = A*x .+ B*u
measurement(x,u,t) = C*x

T = 800 # Number of time steps

d0 = MvNormal(randn(n),2.0)   # Initial state Distribution
du = MvNormal(2,1) # Control input distribution
kf = KalmanFilter(A, B, C, 0, 0.001eye(n), eye(p), d0)
ekf = LLPF.ExtendedKalmanFilter(kf, dynamics, measurement)
x,u,y = LLPF.simulate(kf,T,du)

xf,xt,R,Rt,ll = forward_trajectory(kf, u, y)
xf2,xt2,R2,Rt2,ll2 = forward_trajectory(ekf, u, y)

# When the dynamics is linear, these should be the same
@test reduce(hcat, xf) ≈ reduce(hcat, xf2)
@test reduce(hcat, xt) ≈ reduce(hcat, xt2)
@test reduce(hcat, R)  ≈ reduce(hcat, R2)
@test reduce(hcat, Rt) ≈ reduce(hcat, Rt2)
@test ll ≈ ll2

## add nonlinear dynamics
dynamics2(x,u,t) = A*x - 0.01*abs.(x) .+ B*u
ekf = LLPF.ExtendedKalmanFilter(kf, dynamics2, measurement)
x,u,y = LLPF.simulate(ekf,T,du)

xf,xt,R,Rt,ll = forward_trajectory(kf, u, y)
xf2,xt2,R2,Rt2,ll2 = forward_trajectory(ekf, u, y)

@test norm(reduce(hcat, x .- xf)) > norm(reduce(hcat, x .- xf2))
@test norm(reduce(hcat, x .- xt)) > norm(reduce(hcat, x .- xt2))
@test ll < ll2

# plot(reduce(hcat, x)', layout=2, lab="True")
# plot!(reduce(hcat, xf)', lab="Kf")
# plot!(reduce(hcat, xf2)', lab="EKF")


text:
# Propeller Aerodynamic Validation using UiUC Propeller data,

include("propwing.jl")

module SysOptInit_M
# export sysoptinit, getvar, printmaxviol

# ----- Define function for getting values within objective function ----- #
function getvar(key,x,optvars,optrange,optval,optscale)
    # Determine whether variable is a design variable, assign default if not
    if key in optvars
        value = x[optrange[key]]./optscale[key]'
    else
        value = optval[key][:,3]
    end
    # Make sure Float64 is returned instead of Array{Float64,1} where applicable
    if length(value) == 1
        value = value[1]
    end
    return value
end

# ----- Define function for setting values within objective function ----- #
function setvar(key,x,optvars,optrange,optval,optscale,value)
    # Determine whether variable is a design variable, assign default if not
    if key in optvars
        x[optrange[key]] = value.*optscale[key]'
    end

end

# ----- Define function for printing maximum constraint violations ----- #
function printmaxviol(c,name)
    if !isempty(find(c.>0.0))
        println("$name: ",maximum(c))
    end
end

function sysoptinit(optvars::Array{Symbol,1})
    # ----- Define lower, upper, and starting values ----- #
    optval = Dict{Symbol,Array{Float64,2}}()
    #---General Geometry Variables

    #---Power System Variables
    # Propeller Defaults are for the Epema N250

    #Propulsion
    optval[:velocity]  = [1.0 50.0 19.0]  #m/s
    optval[:proptwist] = [-15.0 80.0 54.9242;
    -15.0 75.0 49.6058;
    -15.0 60.0 44.6923;
    -15.0 45.0 40.1836;
    -15.0 45.0 36.0798;
    -15.0 45.0 32.3809;
    -15.0 45.0 29.0869;
    -15.0 45.0 26.1977;
    -15.0 45.0 23.7134]  #degrees # Propeller Defaults are for the Epema N250

    optval[:propchord] = [0.03 0.3 0.158788;
    0.03 0.3 0.15083;
    0.03 0.3 0.146158;
    0.03 0.3 0.145793;
    0.03 0.3 0.148978;
    0.03 0.3 0.153175;
    0.03 0.3 0.154064;
    0.03 0.3 0.145548;
    0.03 0.3 0.119747]  #chord/Rtip # Propeller Defaults are for the Epema N250

    optval[:RPMprop]   = [10.0 20000.0 4040.0] #RPM
    optval[:Rtip]      = [0.01 2.0 .203] # Radius meters, # Propeller Defaults are for the Epema N250

    # ----- Define scaling for each variable ----- #
    optscale = Dict()

    #Propulsion
    optscale[:velocity]  = [1e-3]
    optscale[:proptwist] = [1e-1]
    optscale[:propchord] = [1e3]
    optscale[:RPMprop]   = [1e-4] #RPM
    optscale[:Rtip]      = [1e0] # Radius meters
    optscale[:T_prop]      = [1e-2] # Radius meters

    # ----- Apply scaling ----- #
    for var in optvars
        optval[var] = optval[var].*optscale[var]'
    end

    # ----- Define Ranges ----- #
    optrange = Dict()
    cx = 0
    for var in optvars
        optrange[var] = cx+(1:size(optval[var],1))
        cx += size(optval[var],1)
    end
    vartotal = cx

    # ----- Assemble lower, upper, and starting optval ----- #
    lb = zeros(vartotal)
    ub = zeros(vartotal)
    x0 = zeros(vartotal)
    for var in optvars
        lb[optrange[var]] = optval[var][:,1]
        ub[optrange[var]] = optval[var][:,2]
        x0[optrange[var]] = optval[var][:,3]
    end

    # ----- Define Optimizer Options ----- #
    options = Dict{String, Any}()
    options["Derivative level"] = 0
    # options["Function precision"] = 3.00E-4
    options["Difference interval"] = 1e-5
    options["Central difference interval"] = 1e-6
    options["Iterations limit"] = 1e8
    options["Major iterations limit"] = 500
    options["Minor iterations limit"]= 1e8
    options["Major optimality tolerance"] = 1e-4
    options["Minor optimality  tolerance"] = 1e-4
    options["Major feasibility tolerance"] = 1e-4
    options["Minor feasibility tolerance"] = 1e-4
    options["Minor print level"] = 1
    options["Print frequency"] = 100
    options["Scale option"] = 2
    options["Scale tolerance"] = .95
    # options["Verify level"] = 3 #only if specifying gradients

    return x0,lb,ub,options,optrange,optval,optscale
end #function sysoptinit

end #module

# Do all calling relative to file's location
fileLoc,_ = splitdir(@__FILE__)

include("$(fileLoc)/../../hale-module-miscellany/params.jl")
push!(LOAD_PATH,"$(fileLoc)/../../hale-module-miscellany/")
using StandardAtmosphere

using CCBlade

using JLD
using PyPlot
close("all")

function prop_aero(x0,optvars,optrange,optval,optscale,T_req,af,rho,unitareas,pitch=0.0,optimize=false,numblades = 2)
    #Unpack the variables, sysoptinit takes care of which are being optimized
    theta = (SysOptInit_M.getvar(:proptwist,x0,optvars,optrange,optval,optscale)+pitch)*pi/180 #in degrees, convert to rad
    RPMp = SysOptInit_M.getvar(:RPMprop,x0,optvars,optrange,optval,optscale)
    Rtip = SysOptInit_M.getvar(:Rtip,x0,optvars,optrange,optval,optscale)
    chords = SysOptInit_M.getvar(:propchord,x0,optvars,optrange,optval,optscale)*Rtip
    velocity = SysOptInit_M.getvar(:velocity,x0,optvars,optrange,optval,optscale)

    r = linspace(.2,.911,length(chords)) * Rtip
    B = numblades
    Rhub = r[1]/2
    chord_prop = chords
    twist_prop = theta
    Vinf = velocity
    Omega = RPMp/60*2*pi
    xtowing = 0.5
    uwake, vwake, T, Q, reff = propwing.propanalysis(Rhub, Rtip, r, chord_prop, twist_prop, B, af, Vinf, Omega, rho, xtowing)
    # Tprop,Qprop,Uprop,Utprop,propvolume,Np,Tp,J_p,eff,CT,CQ = CCcall.cccall(RPMp,theta,chords,CCBlade,Rtip,r,rho,velocity,af,unitareas,numblades)
    Vhub = Vinf
    precone = 0.0
    turbine = false
    J_p = Vhub/(RPMp/60*2*Rtip)
    eff, CT, CQ = CCBlade.nondim(T, Q, Vhub, Omega, rho, Rtip, precone, turbine)
    return T,Q,uwake, vwake, 0.0, 0.0, 0.0,J_p,eff, CT, CQ
end

# push!(LOAD_PATH,"$(fileLoc)/../../Snopt.jl/src/")
# using Snopt
#--------- Initialization ---------#
optvars = [
# :proptwist;
:RPMprop;
:velocity
# :Rtip;
]

x0,lb,ub,options,optrange,optval,optscale = SysOptInit_M.sysoptinit(optvars)

#Initialize Airfoil Data
twists = SysOptInit_M.getvar(:proptwist,x0,optvars,optrange,optval,optscale)
Rtip = SysOptInit_M.getvar(:Rtip,x0,optvars,optrange,optval,optscale)
chords = SysOptInit_M.getvar(:propchord,x0,optvars,optrange,optval,optscale)*Rtip
velocity = SysOptInit_M.getvar(:velocity,x0,optvars,optrange,optval,optscale)
RPMprop = SysOptInit_M.getvar(:RPMprop,x0,optvars,optrange,optval,optscale)
r = collect(linspace(.2,.911,length(chords))) .* Rtip
J = velocity/(RPMprop/60*2*Rtip) #Vinf/(Omega/(2*pi)*Rtip*2)
alt = 222.0 #meters
numblades = 6
rho,mu,q,a = StandardAtmosphere.atmosphere(alt)
re = rho/mu.*chords.*sqrt.((RPMprop*pi/30.0*r).^2+velocity^2)

afgeomnames = ["$(fileLoc)/EpemaExperimental/R023.dat";"$(fileLoc)/EpemaExperimental/R023.dat";
"$(fileLoc)/EpemaExperimental/R0405.dat";"$(fileLoc)/EpemaExperimental/R0405.dat";
"$(fileLoc)/EpemaExperimental/R0405.dat";"$(fileLoc)/EpemaExperimental/R079.dat";
"$(fileLoc)/EpemaExperimental/R079.dat";"$(fileLoc)/EpemaExperimental/R079.dat";"$(fileLoc)/EpemaExperimental/R1.dat"]

N=100
M = 3
new_RPM=4040 #linspace(4000,20000,M)
pitch_array = [-5.0,0.0,5.0]
new_velocity = linspace(1,60,N)
J_p = zeros(N,M)
eta_p = zeros(N,M)
CT_p = zeros(N,M)
CQ_p = zeros(N,M)

SysOptInit_M.setvar(:RPMprop,x0,optvars,optrange,optval,optscale,new_RPM)

velocity = SysOptInit_M.getvar(:velocity,x0,optvars,optrange,optval,optscale)
RPMprop = SysOptInit_M.getvar(:RPMprop,x0,optvars,optrange,optval,optscale)

J = velocity/(RPMprop/60*2*Rtip)

re = rho/mu.*chords.*sqrt.((RPMprop*pi/30.0*r).^2+velocity^2)

# af,unitareas, Al, Au = CCBladeAirfoilTools.af_performance_from_geom(afgeomnames,chords,Rtip,r,J,re)
# JLD.save("$(fileLoc)/af_prop$(new_RPM).jld", "af", af, "unitareas",unitareas,"Al", Al, "Au", Au)
JLDin = JLD.load("$(fileLoc)/af_prop$(new_RPM).jld")
af = JLDin["af"]
unitareas = JLDin["unitareas"]

for j = 1:M

    T_req = 0.0
    pitch = pitch_array[j]+3 #degrees
    optimize=false
    numblades = 6
    args = (optvars,optrange,optval,optscale,T_req,af,rho,unitareas,pitch,optimize,numblades)

    for i = 1:N
        SysOptInit_M.setvar(:velocity,x0,optvars,optrange,optval,optscale,new_velocity[i])

        Tprop,Qprop,Uprop, Utprop, propvolume, Np, Tp,J_p[i,j],eta_p[i,j], CT_p[i,j], CQ_p[i,j] = prop_aero(x0,args...)

    end
end

using CSV
CT_N250_Epema_25deg = CSV.read("$(fileLoc)/EpemaExperimental/CT_N250_25deg.txt", delim = '\t',header = false)
CP_N250_Epema_25deg = CSV.read("$(fileLoc)/EpemaExperimental/CP_N250_25deg.txt", delim = '\t',header = false)
CQ_N250_Epema_25deg = copy(CP_N250_Epema_25deg)
CQ_N250_Epema_25deg[:,2] = CQ_N250_Epema_25deg[:,2]/(2*pi)
ETA_N250_Epema_25deg = CT_N250_Epema_25deg[:,2].*CT_N250_Epema_25deg[:,1]./CP_N250_Epema_25deg[:,2]

CT_N250_Epema_30deg = CSV.read("$(fileLoc)/EpemaExperimental/CT_N250_30deg.txt", delim = '\t',header = false)
CP_N250_Epema_30deg = CSV.read("$(fileLoc)/EpemaExperimental/CP_N250_30deg.txt", delim = '\t',header = false)
CQ_N250_Epema_30deg = copy(CP_N250_Epema_30deg)
CQ_N250_Epema_30deg[:,2] = CQ_N250_Epema_30deg[:,2]/(2*pi)
ETA_N250_Epema_30deg = CT_N250_Epema_30deg[:,2].*CT_N250_Epema_30deg[:,1]./CP_N250_Epema_30deg[:,2]

CT_N250_Epema_35deg = CSV.read("$(fileLoc)/EpemaExperimental/CT_N250_35deg.txt", delim = '\t',header = false)
CP_N250_Epema_35deg = CSV.read("$(fileLoc)/EpemaExperimental/CP_N250_35deg.txt", delim = '\t',header = false)
CQ_N250_Epema_35deg = copy(CP_N250_Epema_35deg)
CQ_N250_Epema_35deg[:,2] = CQ_N250_Epema_35deg[:,2]/(2*pi)
ETA_N250_Epema_35deg = CT_N250_Epema_35deg[:,2].*CT_N250_Epema_35deg[:,1]./CP_N250_Epema_35deg[:,2]

rc("figure", figsize=(6.5, 2.6))
rc("font", size=10.0)
rc("lines", linewidth=1.5)
rc("lines", markersize=3.0)
rc("legend", frameon=false)
rc("axes.spines", right=false, top=false)
rc("figure.subplot", left=0.15, bottom=0.18, top=0.97, right=0.72)
rc("axes", color_cycle=["348ABD", "A60628", "009E73", "7A68A6", "D55E00", "CC79A7"])

mygreen = "#009e73"
myblue = "#348abd"
myred = "#a60628"
plotColors = [myblue,myred,mygreen]

#Thrust Coefficient
PyPlot.figure()
for i = 1:M
    PyPlot.plot(J_p[:,i],CT_p[:,i], marker="o",color = plotColors[i],markeredgecolor = "k",label = "BEM $(30+pitch_array[i])deg RPM=$(new_RPM)")

end
PyPlot.plot(CT_N250_Epema_25deg[:,1],CT_N250_Epema_25deg[:,2],linestyle="-", color = plotColors[1],label = "Epema 25deg")
PyPlot.plot(CT_N250_Epema_30deg[:,1],CT_N250_Epema_30deg[:,2],linestyle="-", color = plotColors[2],label = "Epema 30deg")
PyPlot.plot(CT_N250_Epema_35deg[:,1],CT_N250_Epema_35deg[:,2],linestyle="-", color = plotColors[3],label = "Epema 35deg")
PyPlot.xlabel(L"J\, (V/(n*D_{prop})")
PyPlot.ylabel(L"CT\, (T/(\rho*n^2*{D_{prop}}^4))")
PyPlot.xlim([0,2.0])
PyPlot.ylim([0,.4])
PyPlot.legend(loc = 1)

#Torque Coefficient

PyPlot.figure()
for i = 1:M
    PyPlot.plot(J_p[:,i],CQ_p[:,i], marker="o",color = plotColors[i],markeredgecolor = "k",label = "BEM $(30+pitch_array[i])deg RPM=$(new_RPM)")

end
PyPlot.plot(CQ_N250_Epema_25deg[:,1],CQ_N250_Epema_25deg[:,2],linestyle="-", color = plotColors[1],label = "Epema 25deg")
PyPlot.plot(CQ_N250_Epema_30deg[:,1],CQ_N250_Epema_30deg[:,2],linestyle="-", color = plotColors[2],label = "Epema 30deg")
PyPlot.plot(CQ_N250_Epema_35deg[:,1],CQ_N250_Epema_35deg[:,2],linestyle="-", color = plotColors[3],label = "Epema 35deg")
PyPlot.xlabel(L"J\, (V/(n*D_{prop})")
PyPlot.ylabel(L"Cq\, (T/(\rho*n^2*{D_{prop}}^5)")
PyPlot.xlim([0,2.0])
PyPlot.ylim([0,.1])
PyPlot.legend(loc = 1)

#Efficiency
PyPlot.figure()
for i = 1:M
    PyPlot.plot(J_p[:,i],eta_p[:,i], linestyle="-",color = plotColors[i])#,markeredgecolor = "k",label = "BEM $(30+pitch_array[i])deg RPM=$(new_RPM)")

end
PyPlot.plot(CT_N250_Epema_25deg[:,1],ETA_N250_Epema_25deg,marker="o",linestyle="", color = plotColors[1])#,label = "Epema 25deg")
PyPlot.plot(CT_N250_Epema_30deg[:,1],ETA_N250_Epema_30deg,marker="o",linestyle="", color = plotColors[2])#,label = "Epema 30deg")
PyPlot.plot(CT_N250_Epema_35deg[:,1],ETA_N250_Epema_35deg,marker="o",linestyle="", color = plotColors[3])#,label = "Epema 35deg")

PyPlot.plot([],[],"ko",label = "Experimental")
PyPlot.plot([],[],"k-",label = "Blade Element Momentum")
PyPlot.xlabel("Advance Ratio")
PyPlot.ylabel("Efficiency")
# PyPlot.xlim([0,2.0])
# PyPlot.ylim([0,1.5])
PyPlot.legend(loc="center left", bbox_to_anchor=(.8, 0.8))


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

abstract AbstractCartesianIndex{N} # This is a hacky forward declaration for CartesianIndex
typealias NonSliceIndex Union{Colon, AbstractArray}
typealias ViewIndex Union{Real, NonSliceIndex}
typealias ScalarIndex Union{Real, AbstractCartesianIndex}

# L is true if the view itself supports fast linear indexing
immutable SubArray{T,N,P,I,L} <: AbstractArray{T,N}
    parent::P
    indexes::I
    offset1::Int       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    function SubArray(parent, indexes, offset1, stride1)
        check_parent_index_match(parent, indexes)
        new(parent, indexes, offset1, stride1)
    end
end
# Compute the linear indexability of the indices, and combine it with the linear indexing of the parent
function SubArray(parent::AbstractArray, indexes::Tuple, dims::Tuple)
    SubArray(linearindexing(viewindexing(indexes), linearindexing(parent)), parent, indexes, dims)
end
function SubArray{P, I, N}(::LinearSlow, parent::P, indexes::I, dims::NTuple{N})
    SubArray{eltype(P), N, P, I, false}(parent, indexes, 0, 0)
end
function SubArray{P, I, N}(::LinearFast, parent::P, indexes::I, dims::NTuple{N})
    # Compute the stride and offset
    stride1 = compute_stride1(parent, indexes)
    SubArray{eltype(P), N, P, I, true}(parent, indexes, compute_offset1(parent, stride1, indexes), stride1)
end

check_parent_index_match(parent, indexes) = check_parent_index_match(parent, index_ndims(indexes...))
check_parent_index_match{T,N}(parent::AbstractArray{T,N}, ::NTuple{N, Bool}) = nothing
check_parent_index_match{N}(parent, ::NTuple{N, Bool}) =
    throw(ArgumentError("number of indices ($N) must match the parent dimensionality ($(ndims(parent)))"))

# This computes the linear indexing compatability for a given tuple of indices
viewindexing() = LinearFast()
# Leading scalar indexes simply increase the stride
viewindexing(I::Tuple{ScalarIndex, Vararg{Any}}) = (@_inline_meta; viewindexing(tail(I)))
# Colons may begin a section which may be followed by any number of Colons
viewindexing(I::Tuple{Colon, Colon, Vararg{Any}}) = (@_inline_meta; viewindexing(tail(I)))
# A UnitRange can follow Colons, but only if all other indices are scalar
viewindexing(I::Tuple{Colon, UnitRange, Vararg{ScalarIndex}}) = LinearFast()
# In general, ranges are only fast if all other indices are scalar
viewindexing(I::Tuple{Union{Range, Colon}, Vararg{ScalarIndex}}) = LinearFast()
# All other index combinations are slow
viewindexing(I::Tuple{Vararg{Any}}) = LinearSlow()
# Of course, all other array types are slow
viewindexing(I::Tuple{AbstractArray, Vararg{Any}}) = LinearSlow()

# Simple utilities
size(V::SubArray) = (@_inline_meta; map(n->Int(unsafe_length(n)), indices(V)))

similar(V::SubArray, T::Type, dims::Dims) = similar(V.parent, T, dims)

parent(V::SubArray) = V.parent
parentindexes(V::SubArray) = V.indexes

parent(a::AbstractArray) = a
"""
    parentindexes(A)

From an array view `A`, returns the corresponding indexes in the parent.
"""
parentindexes(a::AbstractArray) = ntuple(i->OneTo(size(a,i)), ndims(a))

## SubArray creation
# We always assume that the dimensionality of the parent matches the number of
# indices that end up getting passed to it, so we store the parent as a
# ReshapedArray view if necessary. The trouble is that arrays of `CartesianIndex`
# can make the number of effective indices not equal to length(I).
_maybe_reshape_parent(A::AbstractArray, ::NTuple{1, Bool}) = reshape(A, Val{1})
_maybe_reshape_parent{_,N}(A::AbstractArray{_,N}, ::NTuple{N, Bool}) = A
_maybe_reshape_parent{N}(A::AbstractArray, ::NTuple{N, Bool}) = reshape(A, Val{N}) # TODO: DEPRECATE FOR #14770
"""
    view(A, inds...)

Like [`getindex`](:func:`getindex`), but returns a view into the parent array `A` with the
given indices instead of making a copy.  Calling [`getindex`](:func:`getindex`) or
[`setindex!`](:func:`setindex!`) on the returned [`SubArray`](:obj:`SubArray`) computes the
indices to the parent array on the fly without checking bounds.
"""
function view(A::AbstractArray, I::ViewIndex...)
    @_inline_meta
    @boundscheck checkbounds(A, I...)
    unsafe_view(_maybe_reshape_parent(A, index_ndims(I...)), I...)
end
# But we can simply flatten scalar `CartesianIndex`s first to make life easier
view(A::AbstractArray, I::Union{ViewIndex, AbstractCartesianIndex}...) = view(A, IteratorsMD.flatten(I)...)

function unsafe_view(A::AbstractArray, I::ViewIndex...)
    @_inline_meta
    J = to_indexes(I...)
    SubArray(A, J, map(unsafe_length, index_shape(A, J...)))
end
# When we take the view of a view, it's often possible to "reindex" the parent
# view's indices such that we can "pop" the parent view and keep just one layer
# of indirection. But we can't always do this because arrays of `CartesianIndex`
# might span multiple parent indices, making the reindex calculation very hard.
# So we use _maybe_reindex to figure out if there are any arrays of
# `CartesianIndex`, and if so, we punt and keep two layers of indirection.
unsafe_view(V::SubArray, I::ViewIndex...) = (@_inline_meta; _maybe_reindex(V, to_indexes(I...)))
_maybe_reindex(V, I) = (@_inline_meta; _maybe_reindex(V, I, I))
_maybe_reindex{C<:AbstractCartesianIndex}(V, I, ::Tuple{AbstractArray{C}, Vararg{Any}}) =
    (@_inline_meta; SubArray(V, I, map(unsafe_length, index_shape(V, I...))))
# But allow arrays of CartesianIndex{1}; they behave just like arrays of Ints
_maybe_reindex{C<:AbstractCartesianIndex{1}}(V, I, A::Tuple{AbstractArray{C}, Vararg{Any}}) =
    (@_inline_meta; _maybe_reindex(V, I, tail(A)))
_maybe_reindex(V, I, A::Tuple{Any, Vararg{Any}}) = (@_inline_meta; _maybe_reindex(V, I, tail(A)))
function _maybe_reindex(V, I, ::Tuple{})
    @_inline_meta
    idxs = reindex(V, V.indexes, to_indexes(I...))
    SubArray(V.parent, idxs, map(unsafe_length, (index_shape(V.parent, idxs...))))
end

## Re-indexing is the heart of a view, transforming A[i, j][x, y] to A[i[x], j[y]]
#
# Recursively look through the heads of the parent- and sub-indexes, considering
# the following cases:
# * Parent index is array  -> re-index that with one or more sub-indexes (one per dimension)
# * Parent index is Colon  -> just use the sub-index as provided
# * Parent index is scalar -> that dimension was dropped, so skip the sub-index and use the index as is

typealias AbstractZeroDimArray{T} AbstractArray{T, 0}

reindex(V, ::Tuple{}, ::Tuple{}) = ()

# Skip dropped scalars, so simply peel them off the parent indices and continue
reindex(V, idxs::Tuple{ScalarIndex, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1], reindex(V, tail(idxs), subidxs)...))

# Colons simply pass their subindexes straight through
reindex(V, idxs::Tuple{Colon, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (subidxs[1], reindex(V, tail(idxs), tail(subidxs))...))

# Re-index into parent vectors with one subindex
reindex(V, idxs::Tuple{AbstractVector, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1]], reindex(V, tail(idxs), tail(subidxs))...))

# Parent matrices are re-indexed with two sub-indices
reindex(V, idxs::Tuple{AbstractMatrix, Vararg{Any}}, subidxs::Tuple{Any, Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1], subidxs[2]], reindex(V, tail(idxs), tail(tail(subidxs)))...))

# In general, we index N-dimensional parent arrays with N indices
@generated function reindex{T,N}(V, idxs::Tuple{AbstractArray{T,N}, Vararg{Any}}, subidxs::Tuple{Vararg{Any}})
    if length(subidxs.parameters) >= N
        subs = [:(subidxs[$d]) for d in 1:N]
        tail = [:(subidxs[$d]) for d in N+1:length(subidxs.parameters)]
        :(@_propagate_inbounds_meta; (idxs[1][$(subs...)], reindex(V, tail(idxs), ($(tail...),))...))
    else
        :(throw(ArgumentError("cannot re-index $(ndims(V)) dimensional SubArray with fewer than $(ndims(V)) indices\nThis should not occur; please submit a bug report.")))
    end
end

# In general, we simply re-index the parent indices by the provided ones
typealias SlowSubArray{T,N,P,I} SubArray{T,N,P,I,false}
function getindex{T,N}(V::SlowSubArray{T,N}, I::Vararg{Real,N})
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[reindex(V, V.indexes, to_indexes(I...))...]
    r
end

typealias FastSubArray{T,N,P,I} SubArray{T,N,P,I,true}
function getindex(V::FastSubArray, i::Real)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + V.stride1*to_index(i)]
    r
end
# We can avoid a multiplication if the first parent index is a Colon or UnitRange
typealias FastContiguousSubArray{T,N,P,I<:Tuple{Union{Colon, UnitRange}, Vararg{Any}}} SubArray{T,N,P,I,true}
function getindex(V::FastContiguousSubArray, i::Real)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds r = V.parent[V.offset1 + to_index(i)]
    r
end

function setindex!{T,N}(V::SlowSubArray{T,N}, x, I::Vararg{Real,N})
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds V.parent[reindex(V, V.indexes, to_indexes(I...))...] = x
    V
end
function setindex!(V::FastSubArray, x, i::Real)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + V.stride1*to_index(i)] = x
    V
end
function setindex!(V::FastContiguousSubArray, x, i::Real)
    @_inline_meta
    @boundscheck checkbounds(V, i)
    @inbounds V.parent[V.offset1 + to_index(i)] = x
    V
end

linearindexing{T<:FastSubArray}(::Type{T}) = LinearFast()
linearindexing{T<:SubArray}(::Type{T}) = LinearSlow()

getindex(::Colon, i) = to_index(i)
unsafe_getindex(::Colon, i) = to_index(i)

step(::Colon) = 1
isempty(::Colon) = false
in(::Integer, ::Colon) = true

# Strides are the distance between adjacent elements in a given dimension,
# so they are well-defined even for non-linear memory layouts
strides{T,N,P,I}(V::SubArray{T,N,P,I}) = substrides(V.parent, V.indexes)

substrides(parent, I::Tuple) = substrides(1, parent, 1, I)
substrides(s, parent, dim, ::Tuple{}) = ()
substrides(s, parent, dim, I::Tuple{Real, Vararg{Any}}) = (substrides(s*size(parent, dim), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) = substrides(s, parent, dim, (I[1].I..., tail(I)...))
substrides(s, parent, dim, I::Tuple{Colon, Vararg{Any}}) = (s, substrides(s*size(parent, dim), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{Range, Vararg{Any}}) = (s*step(I[1]), substrides(s*size(parent, dim), parent, dim+1, tail(I))...)
substrides(s, parent, dim, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubArrays with indices of type $(typeof(I[1]))"))

stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

compute_stride1{N}(parent::AbstractArray, I::NTuple{N}) =
    compute_stride1(1, fill_to_length(indices(parent), OneTo(1), Val{N}), I)
compute_stride1(s, inds, I::Tuple{}) = s
compute_stride1(s, inds, I::Tuple{Real, Vararg{Any}}) =
    (@_inline_meta; compute_stride1(s*unsafe_length(inds[1]), tail(inds), tail(I)))
compute_stride1(s, inds, I::Tuple{Range, Vararg{Any}}) = s*step(I[1])
compute_stride1(s, inds, I::Tuple{Colon, Vararg{Any}}) = s
function compute_stride1{N}(s, inds, I::Tuple{AbstractCartesianIndex{N}, Vararg{Any}})
    @_inline_meta
    h, t = IteratorsMD.split(inds, Val{N})
    compute_stride1(s*prod(map(unsafe_length, h)), t, tail(I))
end
compute_stride1(s, inds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("invalid strided index type $(typeof(I[1]))"))

iscontiguous(A::SubArray) = iscontiguous(typeof(A))
iscontiguous{S<:SubArray}(::Type{S}) = false
iscontiguous{F<:FastContiguousSubArray}(::Type{F}) = true

first_index(V::FastSubArray) = V.offset1 + V.stride1 # cached for fast linear SubArrays
function first_index(V::SubArray)
    P, I = parent(V), V.indexes
    s1 = compute_stride1(P, I)
    s1 + compute_offset1(P, s1, I)
end

# Computing the first index simply steps through the indices, accumulating the
# sum of index each multiplied by the parent's stride.
# The running sum is `f`; the cumulative stride product is `s`.
# If the result is one-dimensional and it's a Colon, then linear
# indexing uses the indices along the given dimension. Otherwise
# linear indexing always starts with 1.
compute_offset1(parent, stride1::Integer, I::Tuple) = (@_inline_meta; compute_offset1(parent, stride1, find_extended_dims(I)..., I))
compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{Colon}, I::Tuple) = compute_linindex(parent, I) - stride1*first(indices(parent, dims[1]))  # index-preserving case
compute_offset1(parent, stride1::Integer, dims, inds, I::Tuple) = compute_linindex(parent, I) - stride1  # linear indexing starts with 1

function compute_linindex{N}(parent, I::NTuple{N})
    IP = fill_to_length(indices(parent), OneTo(1), Val{N})
    compute_linindex(1, 1, IP, I)
end
function compute_linindex(f, s, IP::Tuple, I::Tuple{Real, Vararg{Any}})
    @_inline_meta
    Δi = I[1]-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
# Just splat out the cartesian indices and continue
compute_linindex(f, s, IP::Tuple, I::Tuple{AbstractCartesianIndex, Vararg{Any}}) =
    (@_inline_meta; compute_linindex(f, s, IP, (I[1].I..., tail(I)...)))
compute_linindex(f, s, IP::Tuple, I::Tuple{Colon, Vararg{Any}}) =
    (@_inline_meta; compute_linindex(f, s*unsafe_length(IP[1]), tail(IP), tail(I)))
function compute_linindex(f, s, IP::Tuple, I::Tuple{Any, Vararg{Any}})
    @_inline_meta
    Δi = first(I[1])-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
compute_linindex(f, s, IP::Tuple, I::Tuple{}) = f

find_extended_dims(I) = (@_inline_meta; _find_extended_dims((), (), 1, I...))
_find_extended_dims(dims, inds, dim) = dims, inds
_find_extended_dims(dims, inds, dim, ::Real, I...) = _find_extended_dims(dims, inds, dim+1, I...)
_find_extended_dims(dims, inds, dim, i1::AbstractCartesianIndex, I...) = _find_extended_dims(dims, inds, dim, i1.I..., I...)
_find_extended_dims(dims, inds, dim, i1, I...) = _find_extended_dims((dims..., dim), (inds..., i1), dim+1, I...)

unsafe_convert{T,N,P,I<:Tuple{Vararg{RangeIndex}}}(::Type{Ptr{T}}, V::SubArray{T,N,P,I}) =
    unsafe_convert(Ptr{T}, V.parent) + (first_index(V)-1)*sizeof(T)

pointer(V::FastSubArray, i::Int) = pointer(V.parent, V.offset1 + V.stride1*i)
pointer(V::FastContiguousSubArray, i::Int) = pointer(V.parent, V.offset1 + i)
pointer(V::SubArray, i::Int) = _pointer(V, i)
_pointer{T}(V::SubArray{T,1}, i::Int) = pointer(V, (i,))
_pointer(V::SubArray, i::Int) = pointer(V, ind2sub(indices(V), i))

function pointer{T,N,P<:Array,I<:Tuple{Vararg{RangeIndex}}}(V::SubArray{T,N,P,I}, is::Tuple{Vararg{Int}})
    index = first_index(V)
    strds = strides(V)
    for d = 1:length(is)
        index += (is[d]-1)*strds[d]
    end
    return pointer(V.parent, index)
end

# indices of the parent are preserved for ::Colon indices, otherwise
# they are taken from the range/vector
# Since bounds-checking is performance-critical and uses
# indices, it's worth optimizing these implementations thoroughly
indices(S::SubArray) = (@_inline_meta; _indices_sub(S, indices(S.parent), S.indexes...))
_indices_sub(S::SubArray, pinds) = ()
function _indices_sub(S::SubArray, pinds, ::Real, I...)
    @_inline_meta
    _indices_sub(S, tail(pinds), I...)
end
function _indices_sub(S::SubArray, pinds, ::Colon, I...)
    @_inline_meta
    (pinds[1], _indices_sub(S, tail(pinds), I...)...)
end
function _indices_sub(S::SubArray, pinds, i1::AbstractArray, I...)
    @_inline_meta
    (unsafe_indices(i1)..., _indices_sub(S, tail(pinds), I...)...)
end
function _indices_sub{N}(S::SubArray, pinds, ::AbstractCartesianIndex{N}, I...)
    @_inline_meta
    _indices_sub(S, IteratorsMD.split(pinds, Val{N})[2], I...)
end
# _indices_sub for arrays of CartesianIndex is defined in multidimensional.jl

## Compatability
# deprecate?
function parentdims(s::SubArray)
    nd = ndims(s)
    dimindex = Array{Int}(nd)
    sp = strides(s.parent)
    sv = strides(s)
    j = 1
    for i = 1:ndims(s.parent)
        r = s.indexes[i]
        if j <= nd && (isa(r,Union{Colon,Range}) ? sp[i]*step(r) : sp[i]) == sv[j]
            dimindex[j] = i
            j += 1
        end
    end
    dimindex
end

"""
    replace_ref_end!(ex)

Recursively replace occurrences of the symbol :end in a "ref" expression (i.e. A[...]) `ex`
with the appropriate function calls (`endof`, `size` or `trailingsize`). Replacement uses
the closest enclosing ref, so

    A[B[end]]

should transform to

    A[B[endof(B)]]

"""
function replace_ref_end!(ex,withex=nothing)
    if isa(ex,Symbol) && ex == :end
        withex === nothing && error("Invalid use of end")
        return withex
    elseif isa(ex,Expr)
        if ex.head == :ref
            S = ex.args[1] = replace_ref_end!(ex.args[1],withex)
            # new :ref, so redefine withex
            nargs = length(ex.args)-1
            if nargs == 0
                return ex
            elseif nargs == 1
                # replace with endof(S)
                ex.args[2] = replace_ref_end!(ex.args[2],:(Base.endof($S)))
            else
                n = 1
                J = endof(ex.args)
                for j = 2:J-1
                    exj = ex.args[j] = replace_ref_end!(ex.args[j],:(Base.size($S,$n)))
                    if isa(exj,Expr) && exj.head == :...
                        # splatted object
                        exjs = exj.args[1]
                        n = :($n + length($exjs))
                    elseif isa(n, Expr)
                        # previous expression splatted
                        n = :($n + 1)
                    else
                        # an integer
                        n += 1
                    end
                end
                ex.args[J] = replace_ref_end!(ex.args[J],:(Base.trailingsize($S,$n)))
            end
        else
            # recursive search
            for i = eachindex(ex.args)
                ex.args[i] = replace_ref_end!(ex.args[i],withex)
            end
        end
    end
    ex
end

"""
    @view A[inds...]

Creates a `SubArray` from an indexing expression. This can only be applied directly to a
reference expression (e.g. `@view A[1,2:end]`), and should *not* be used as the target of
an assignment (e.g. `@view(A[1,2:end]) = ...`).
"""
macro view(ex)
    if isa(ex, Expr) && ex.head == :ref
        ex = replace_ref_end!(ex)
        Expr(:&&, true, esc(Expr(:call,:(Base.view),ex.args...)))
    else
        throw(ArgumentError("Invalid use of @view macro: argument must be a reference expression A[...]."))
    end
end

max_stars_repo_name: MichaelHatherly/julia

max_issues_repo_name: MichaelHatherly/julia

max_forks_repo_name: MichaelHatherly/julia


text:
@testset "Stats" begin
    autocor(DataArray([1, 2, 3, 4, 5]))
end

@testset "describe" begin
    io = IOBuffer()
    describe(io, @data(collect(1:10)))
    @test String(take!(io)) == """
                               Summary Stats:
                               Mean:           5.500000
                               Minimum:        1.000000
                               1st Quartile:   3.250000
                               Median:         5.500000
                               3rd Quartile:   7.750000
                               Maximum:        10.000000
                               Length:         10
                               Type:           $Int
                               Number Missing: 0
                               % Missing:      0.000000
                               """
    describe(io, @data([1, missing]))
    @test String(take!(io)) == """
                               Summary Stats:
                               Mean:           1.000000
                               Minimum:        1.000000
                               1st Quartile:   1.000000
                               Median:         1.000000
                               3rd Quartile:   1.000000
                               Maximum:        1.000000
                               Length:         1
                               Type:           $Int
                               Number Missing: 1
                               % Missing:      50.000000
                               """
    describe(io, @data(["s"]))
    @test String(take!(io)) == """
                               Summary Stats:
                               Length:         1
                               Type:           String
                               Number Unique:  1
                               Number Missing: 0
                               % Missing:      0.000000
                               """
    describe(io, @data(["s", missing]))
    @test String(take!(io)) == """
                               Summary Stats:
                               Length:         2
                               Type:           String
                               Number Unique:  2
                               Number Missing: 1
                               % Missing:      50.000000
                               """
    describe(io, DataArray(Any, 5))
    @test String(take!(io)) == """
                               Summary Stats:
                               Length:         5
                               Type:           Any
                               Number Unique:  1
                               Number Missing: 5
                               % Missing:      100.000000
                               """
    describe(io, DataArray(Float64, 5))
    @test String(take!(io)) == """
                               Summary Stats:
                               Type:           Float64
                               Number Missing: 5
                               % Missing:      100.000000
                               """
end


text:
using LabConnections.BeagleBone
import LabConnections.BeagleBone: initdev, listdev, closedev, printdev, write!, read, gpio_channels

using Base.Test

@testset "GPIO tests" begin

    @testset "Inialization/Termination" begin
        # Initialize three devices
        initdev("gpio", Int32(1))
        @test sum(listdev()) == 1

        initdev("gpio", Int32(3))
        @test sum(listdev()) == 2

        initdev("gpio", Int32(5))
        @test sum(listdev()) == 3

        #printdev("gpio", 3)

        # Attempt to initialize a device which has already been initialized
        @test_throws ErrorException initdev("gpio", Int32(1))
        @test_throws ErrorException initdev("gpio", Int32(3))
        @test_throws ErrorException initdev("gpio", Int32(5))

        # Attempt to initialize a device with a very high index (no matching channel)
        @test_throws ErrorException initdev("gpio", Int32(1000))

        # Attempt to remove devices which have not been initialized
        @test_throws ErrorException closedev("gpio", Int32(2))
        @test_throws ErrorException closedev("gpio", Int32(4))
        @test_throws ErrorException closedev("gpio", Int32(6))

        #printdev("gpio", 3)

        # Remove devices from TOC
        closedev("gpio", Int32(1))
        @test sum(listdev()) == 2

        closedev("gpio", Int32(3))
        @test sum(listdev()) == 1

        closedev("gpio", Int32(5))
        @test sum(listdev()) == 0

    end

    @testset "Exceptions" begin

        # Fixture
        device = initdev("gpio", Int32(1))

        # Test that an exception is thrown when an invalid operation is given
        # supported operations are 1,2,3
        @test_throws ErrorException write!(device, (Int32(0), "something"))
        @test_throws ErrorException write!(device, (Int32(4), "something"))

        # Test that exceptions are thrown for each individual operation
        @test_throws ErrorException write!(device, (Int32(1), "bad_entry"))
        @test_throws ErrorException write!(device, (Int32(2), "bad_entry"))
        @test_throws ErrorException write!(device, (Int32(3), "bad_entry"))

        # Test operation 3
        @test_throws ErrorException write!(device, (Int32(3), "none"))
        @test_throws ErrorException write!(device, (Int32(3), "rising"))
        @test_throws ErrorException write!(device, (Int32(3), "falling"))
        @test_throws ErrorException write!(device, (Int32(3), "both"))

        # Close Gpio
        closedev("gpio", Int32(1))
    end

    @testset "All channels" begin
        # Configure the GPIO for output usage
        devices = []
        for ii = 1:length(gpio_channels)
            device = initdev("gpio", Int32(ii))
            # Operation 2 -> in/out, set out
            write!(device, (Int32(2), "out"))
            # Append to list
            append!(devices, [device])
        end

        # Sets all available GPIO pins high/low in an alternating pattern
        for i = 1:20
            for ii = 1:length(gpio_channels)
                state = "$(i%2)"
                write!(devices[ii], (Int32(1), state))
                @test read(devices[ii], Int32(1)) == state
            end
            sleep(0.10)
        end

        # Closes all devices
        for ii = 1:length(gpio_channels)
            closedev("gpio", Int32(ii))
        end
    end

    @testset "Read/Write" begin

        # Fixture
        device = initdev("gpio", Int32(1))

        # Test operation 1
        write!(device, (Int32(1), "1"))
        @test read(device, Int32(1)) == "1"
        write!(device, (Int32(1), "0"))
        @test read(device, Int32(1)) == "0"
        write!(device, (Int32(1), "1"))
        @test read(device, Int32(1)) == "1"
        write!(device, (Int32(1), "0"))
        @test read(device, Int32(1)) == "0"

        write!(device, (Int32(2), "in"))
        @test read(device, Int32(2)) == "in"
        write!(device, (Int32(2), "out"))
        @test read(device, Int32(2)) == "out"
        write!(device, (Int32(2), "in"))
        @test read(device, Int32(2)) == "in"
        write!(device, (Int32(2), "out"))
        @test read(device, Int32(2)) == "out"

        # Close Gpio
        closedev("gpio", Int32(1))
    end
end


text:
# This file was generated, do not modify it. # hide
rc = range(svc_mdl, :cost, lower=0.1, upper=5)
tm = TunedModel(model=svc_mdl, ranges=[rc], tuning=Grid(resolution=10),
                resampling=CV(nfolds=3, rng=33), measure=misclassification_rate)
mtm = machine(tm, X, y)

fit!(mtm)

ypred = predict(mtm, X)
misclassification_rate(ypred, y)


text:
#Script for saving atomic energy levels
#24/07/2017

using JLD

atom = "87Rb";

D1trans = 1.56051536; #D1 line transition energy in eV
D2trans = 1.58999085; #D2 line transition energy in eV
n4l2j5_2en = 2.40116159;
n4l2j3_2en = 2.40121692;
n6l0j1_2en = 2.49759249;
n6l1j1_2en = 2.94203794;
n6l1j3_2en = 2.95165365;
n5l0j1_2en = 0;
ionlim = 4.1771270; #Rb groundstate energy in eV

save(string(atom,"levels.jld"),"D1trans",D1trans,"D2trans",D2trans,"n4l2j5_2en",n4l2j5_2en,"n4l2j3_2en",n4l2j3_2en,"n6l0j1_2en",n6l0j1_2en,"n6l1j1_2en",n6l1j1_2en,"n6l1j3_2en",n6l1j3_2en,"n5l0j1_2en",n5l0j1_2en,"ionlim",ionlim);

max_stars_repo_name: CSChisholm/Rydberg

max_issues_repo_name: CSChisholm/Rydberg

max_forks_repo_name: CSChisholm/Rydberg


text:
# Sample postprocessing script for customized figure.
#
# Reading 2D simulation data, zooming-in to a region of interest, and
# adding streamlines on top of colored mesh.
#
# Hongyang Zhou, hyzhou@umich.edu

using Vlasiator, Plots

file = "bulk.0000501.vlsv"
nameρ = "rho"
nameV = "rho_v"

boxcoords = [0, 20, -15, 15]

meta = load(file)

heatmap(meta, nameρ,
   xlim=(boxcoords[1], boxcoords[2]),
   ylim=(boxcoords[3], boxcoords[4]),
   aspect_ratio=:equal,
   c=:turbo)

#=
# Attributes can be modified afterwards, but it's slower.
heatmap(meta, nameρ, c=:turbo)
xlims!(boxcoords[1], boxcoords[2])
ylims!(boxcoords[3], boxcoords[4])
=#

#streamplot(meta, nameV, comp="xy", color="w", density=2.0)


text:
function extract_cells(file, ::Type{TI}=Int) where {TI}
    line = readline(file)

    cell_idx_pattern = r"^(\d+)\s*,"
    m = match(cell_idx_pattern, line)
    first_cell_idx = parse(TI, m[1])

    node_idx_pattern = r",\s*(\d+)"
    local cells
    nodes = TI[]
    for m in eachmatch(node_idx_pattern, line)
        push!(nodes, parse(TI, m[1]))
    end
    cells = [Tuple(nodes)]

    nextline = _extract_cells!(cells, file, first_cell_idx)
    return cells, first_cell_idx - TI(1), nextline
end

function _extract_cells!(
    cells::AbstractVector{NTuple{nnodes,TI}}, file, prev_cell_idx::TI
) where {nnodes,TI}
    cell_idx_pattern = r"^(\d+)\s*,"
    node_idx_pattern = r",\s*(\d+)"
    nodes = zeros(TI, nnodes)

    line = readline(file)
    m = match(stopping_pattern, line)
    while m isa Nothing
        m = match(cell_idx_pattern, line)
        if m != nothing
            cell_idx = parse(TI, m[1])
            cell_idx == prev_cell_idx + TI(1) || throw("Cell indices are not consecutive.")

            nodes .= zero(TI)
            for (i, m) in enumerate(eachmatch(node_idx_pattern, line))
                nodes[i] = parse(TI, m[1])
            end
            all(nodes .!= zero(TI)) ||
                throw("Cell $cell_idx has fewer nodes than it should.")
            push!(cells, Tuple(nodes))

            prev_cell_idx = cell_idx
        end
        line = readline(file)
        m = match(stopping_pattern, line)
    end

    return line
end


text:
#=
Given a list of numbers and a number k, return whether any two numbers from the list add up to k.

For example, given [10, 15, 3, 7] and k of 17, return true since 10 + 7 is 17.
=#

function addupto(X, k)
    for i in 1:length(X)-1
        for j in i+1:length(X)
             if X[i]+X[j] == k
                 return true
             end
        end
    end
    return false
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
# Julia program to toggle the k-th bit of a number.
function toggle(num, k)
    return (num ⊻ (1 << (k-1)))
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
print("Enter the value of k(where you need to toggle the k'th bit): ")
k = readline()
k = parse(Int, k)
new_num = toggle(num, k)
println("The given number, after toggling the k-th bit is $new_num.")

"""
Time Complexity: O(1)
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 24
Enter the value of k(where you need to toggle the k'th bit): 3
The given number, after toggling the k-th bit is 28.

SAMPLE 2
Enter the number: 33
Enter the value of k(where you need to toggle the k'th bit): 12
The given number, after toggling the k-th bit is 2081.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
@testset "Spectral functions" begin
    function testKernelT(isFermi, symmetry, sign, β, τ, ε)
        println("Testing $symmetry  for β=$β, τ=$τ, ε=$ε")
        n(ε, β) = 1.0 / (exp(ε * β) - sign * 1.0)
        @test Spectral.density(isFermi, ε, β) ≈ n(ε, β)
        @test Spectral.kernelT(Val(isFermi), Val(symmetry), eps(0.0), ε, β) ≈ 1.0 + sign * n(ε, β)

        # @test Spectral.kernelT(type, 0.0, ε, β) ≈ sign * n(ε, β) # τ=0.0 should be treated as the 0⁻
        @test Spectral.kernelT(Val(isFermi), Val(symmetry), -eps(0.0), ε, β) ≈ sign * n(ε, β)

        @test Spectral.kernelT(Val(isFermi), Val(symmetry), -τ, ε, β) ≈ sign * Spectral.kernelT(Val(isFermi), Val(symmetry), β - τ, ε, β)
        @test Spectral.kernelT(Val(isFermi), Val(symmetry), -eps(0.0), 1000.0, 1.0) ≈ 0.0
        @test Spectral.kernelT(Val(isFermi), Val(symmetry), -eps(0.0), -1000.0, 1.0) ≈ -1.0
        if isFermi
            # ω can not be zero for boson
            @test Spectral.kernelT(Val(true), Val(:none), τ, 0.0, β) ≈ n(0.0, β)
        end
    end

    testKernelT(true, :none, -1.0, 10.0, 1.0, 1.0)
    testKernelT(false, :none, 1.0, 10.0, 1.0, 1.0)

    testKernelT(true, :none, -1.0, 10.0, 1e-10, 1.0)
    testKernelT(false, :none, 1.0, 10.0, 1e-10, 1.0)

    testKernelT(true, :none, -1.0, 10.0, 1.0, 1.0e-6)
    testKernelT(false, :none, 1.0, 10.0, 1.0, 1.0e-6) # small ϵ for bosonic case is particularly dangerous because the kernal diverges ~1/ϵ

    function testAccuracy(isFermi, symmetry, τGrid, ωGrid, β)
        maxErr = BigFloat(0.0)
        τ0, ω0, macheps = 0.0, 0.0, 0.0
        for (τi, τ) in enumerate(τGrid)
            for (ωi, ω) in enumerate(ωGrid)
                ker1 = Spectral.kernelT(Val(isFermi), Val(symmetry), τ, ω, β)
                ker2 = Spectral.kernelT(Val(isFermi), Val(symmetry), BigFloat(τ), BigFloat(ω), BigFloat(β))
                if abs(ker1 - ker2) > maxErr
                    maxErr = abs(ker1 - ker2)
                    τ0, ω0, macheps = τ, ω, eps(ker1) # get the machine accuracy for the float number ker1
                end
            end
        end
        @test maxErr < 2macheps
        return maxErr, τ0, ω0, macheps
    end

    β, Euv = 10000.0, 100.0
    τGrid = [t for t in LinRange(-β + 1e-10, β, 100)]
    ωGrid = [w for w in LinRange(-Euv, Euv, 100)]
    maxErr, τ0, ω0, macheps = testAccuracy(true, :none, τGrid, ωGrid, β)
    maxErr, τ0, ω0, macheps = testAccuracy(false, :none, τGrid, ωGrid, β)

    τGrid = [t for t in LinRange(-β + 1e-6, β, 100)]
    ωGrid = [-1e-6, -1e-8, -1e-10, 1e-10, 1e-8, 1e-6]
    maxErr, τ0, ω0, macheps = testAccuracy(true, :none, τGrid, ωGrid, β)
    maxErr, τ0, ω0, macheps = testAccuracy(false, :none, τGrid, ωGrid, β)

    τGrid = [-1e-6, -1e-8, -1e-10, 1e-10, 1e-8, 1e-6]
    ωGrid = [w for w in LinRange(-Euv, Euv, 100)]
    maxErr, τ0, ω0, macheps = testAccuracy(true, :none, τGrid, ωGrid, β)
    maxErr, τ0, ω0, macheps = testAccuracy(false, :none, τGrid, ωGrid, β)

    τGrid = [β - 1e-6, β - 1e-8, β - 1e-10, -β + 1e-10, -β + 1e-8, -β + 1e-6]
    ωGrid = [w for w in LinRange(-Euv, Euv, 100)]
    maxErr, τ0, ω0, macheps = testAccuracy(true, :none, τGrid, ωGrid, β)
    maxErr, τ0, ω0, macheps = testAccuracy(false, :none, τGrid, ωGrid, β)
end


text:
base_path = "/Volumes/davidkunz/Documents/bakalarka/mri_head_homography"
target_images_path = joinpath(base_path, "images", "target")
target_images_fnames = filter(x -> !isfile(x), readdir(target_images_path))
target_images_paths = map(fname -> joinpath(target_images_path, fname), target_images_fnames)

for target_img_path in target_images_paths
    k = parse(Int, target_img_path[findfirst(r"[0-9]+", target_img_path)])
    #path
    source_img_path = joinpath(base_path, "images", "source", "img$(k).png")
    target = load(target_img_path)
    source = load(source_img_path)

    source_jpg_path = source_img_path[1:end-3]*"jpg"
    target_jpg_path = target_img_path[1:end-3]*"jpg"
    @info k source_jpg_path target_jpg_path
    # resave
    save(target_jpg_path, target)
    save(source_jpg_path, source)
end

source_images_path = joinpath(base_path, "images", "source")

for img in filter(x -> occursin(".png", x), readdir(source_images_path))
    run(`rm $(joinpath(source_images_path, img))`)
end

## remane to jpg

tmp_df = CSV.read("/Volumes/davidkunz/Documents/bakalarka/mri_head_homography/loc_table_png.csv") |> DataFrame

names(tmp_df)

collect(tmp_df["Source image"])
new_names_source = map(x -> replace(x, "png" => "jpg", count=1), collect(tmp_df["Source image"]))
new_names_target = map(x -> replace(x, "png" => "jpg", count=1), collect(tmp_df["Target image"]))

df = deepcopy(tmp_df)
df["Source image"] = new_names_source
df["Target image"] = new_names_target

CSV.write("/Volumes/davidkunz/Documents/bakalarka/mri_head_homography/loc_table.csv", df)


text:
using Private
using Base.Test
using Logging
using Lint
using EEG
using TimeModels

Logging.configure(level=DEBUG)
Logging.configure(output=open("logfile.log", "a"))

include("source_analysis/beamformer.jl")
include("types/Leadfield/Leadfield.jl")
include("statistics/hotelling.jl")
include("statistics/kalman.jl")
include("reshaping/subsample.jl")
include("reshaping/detrend.jl")
include("plotting/srf.jl")

@printf("\n\n\n All tests passed. Now checking lint\n\n\n")

lintpkg( "Private" )


text:
@doc raw"""
    get_reason(c)

return the current reason stored within a [`StoppingCriterion`](@ref) `c`.
This reason is empty if the criterion has never been met.
"""
get_reason(c::sC) where {sC<:StoppingCriterion} = c.reason

@doc raw"""
    get_stopping_criteria(c)
return the array of internally stored [`StoppingCriterion`](@ref)s for a
[`StoppingCriterionSet`](@ref) `c`.
"""
function get_stopping_criteria(c::S) where {S<:StoppingCriterionSet}
    return error("get_stopping_criteria() not defined for a $(typeof(c)).")
end

@doc raw"""
    StopAfterIteration <: StoppingCriterion

A functor for an easy stopping criterion, i.e. to stop after a maximal number
of iterations.

# Fields
* `maxIter` – stores the maximal iteration number where to stop at
* `reason` – stores a reason of stopping if the stopping criterion has one be
  reached, see [`get_reason`](@ref).

# Constructor

    StopAfterIteration(maxIter)

initialize the stopafterIteration functor to indicate to stop after `maxIter`
iterations.
"""
mutable struct StopAfterIteration <: StoppingCriterion
    maxIter::Int
    reason::String
    StopAfterIteration(mIter::Int) = new(mIter, "")
end
function (c::StopAfterIteration)(::P, ::O, i::Int) where {P<:Problem,O<:Options}
    if i >= c.maxIter
        c.reason = "The algorithm reached its maximal number of iterations ($(c.maxIter)).\n"
        return true
    end
    return false
end
"""
    StopWhenTimeElapsed <: StoppingCriterion

stores a threshold indicating when to stop after some elapsed time  [`GradientProblem`](@ref).
"""
mutable struct StopWhenTimeElapsed <: StoppingCriterion
    threshold::Float64
    reason::String
    StopWhenTimeElapsed(T::Float64) = new(T, "")
end
function (c::StopWhenTimeElapsed)(p::P, o::O, iter::Int) where {P<:Problem,O<:Options}
    if o.timer > c.threshold # let's see if o.timer works, perhaps will have to do another one for RTR too
        c.reason = "The algorithm reached maximum time allowed after $iter iterations; time elapsed ($(o.timer)) is more than $(c.threshold).\n"
        return true
    end
    return false
end

"""
    StopWhenGradientNormLess <: StoppingCriterion

stores a threshold when to stop looking at the norm of the gradient from within
a [`GradientProblem`](@ref).
"""
mutable struct StopWhenGradientNormLess <: StoppingCriterion
    threshold::Float64
    reason::String
    StopWhenGradientNormLess(ε::Float64) = new(ε, "")
end
function (c::StopWhenGradientNormLess)(p::Problem, o::AbstractGradientOptions, iter::Int)
    if norm(p.M, o.x, o.gradient) < c.threshold
        c.reason = "The algorithm reached approximately critical point after $iter iterations; the gradient norm ($(norm(p.M,o.x,o.gradient))) is less than $(c.threshold).\n"
        return true
    end
    return false
end
"""
    StopWhenChangeLess <: StoppingCriterion

stores a threshold when to stop looking at the norm of the change of the
optimization variable from within a [`Options`](@ref), i.e `o.x`.
For the storage a [`StoreOptionsAction`](@ref) is used

# Constructor

    StopWhenChangeLess(ε[, a])

initialize the stopping criterion to a threshold `ε` using the
[`StoreOptionsAction`](@ref) `a`, which is initialized to just store `:x` by
default.
"""
mutable struct StopWhenChangeLess <: StoppingCriterion
    threshold::Float64
    reason::String
    storage::StoreOptionsAction
    function StopWhenChangeLess(ε::Float64, a::StoreOptionsAction=StoreOptionsAction((:x,)))
        return new(ε, "", a)
    end
end
function (c::StopWhenChangeLess)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if has_storage(c.storage, :x)
        xOld = get_storage(c.storage, :x)
        if distance(p.M, o.x, xOld) < c.threshold && i > 0
            c.reason = "The algorithm performed a step with a change ($(distance(p.M, o.x, xOld))) less than $(c.threshold).\n"
            c.storage(p, o, i)
            return true
        end
    end
    c.storage(p, o, i)
    return false
end

"""
    StopWhenStepsizeLess <: StoppingCriterion

stores a threshold when to stop looking at the last step size determined or found
during the last iteration from within a [`Options`](@ref).

# Constructor

    StopWhenStepsizeLess(ε)

initialize the stopping criterion to a threshold `ε`.
"""
mutable struct StopWhenStepSizeLess <: StoppingCriterion
    threshold::Float64
    reason::String
    function StopWhenStepSizeLess(ε::Float64)
        return new(ε, "")
    end
end
function (c::StopWhenStepSizeLess)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    s = get_last_stepsize(p, o, i)
    if s < c.threshold && i > 0
        c.reason = "The algorithm computed a step size ($s) less than $(c.threshold).\n"
        return true
    end
    return false
end

@doc raw"""
    StopWhenAll <: StoppingCriterion

store an array of [`StoppingCriterion`](@ref) elements and indicates to stop,
when _all_ indicate to stop. The `reason` is given by the concatenation of all
reasons.

# Constructor
    StopWhenAll(c::NTuple{N,StoppingCriterion} where N)
    StopWhenAll(c::StoppingCriterion,...)
"""
mutable struct StopWhenAll{TCriteria<:Tuple} <: StoppingCriterionSet
    criteria::TCriteria
    reason::String
    StopWhenAll(c::Vector{StoppingCriterion}) = new{typeof(tuple(c...))}(tuple(c...), "")
    StopWhenAll(c...) = new{typeof(c)}(c, "")
end
function (c::StopWhenAll)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if all(subC -> subC(p, o, i), c.criteria)
        c.reason = string([get_reason(subC) for subC in c.criteria]...)
        return true
    end
    return false
end
get_stopping_criteria(c::StopWhenAll) = c.criteria

"""
    &(s1,s2)
    s1 & s2

Combine two [`StoppingCriterion`](@ref) within an [`StopWhenAll`](@ref).
If either `s1` (or `s2`) is already an [`StopWhenAll`](@ref), then `s2` (or `s1`) is
appended to the list of [`StoppingCriterion`](@ref) within `s1` (or `s2`).

# Example
    a = StopAfterIteration(200) & StopWhenChangeLess(1e-6)
    b = a & StopWhenGradientNormLess(1e-6)

Is the same as

    a = StopWhenAll(StopAfterIteration(200), StopWhenChangeLess(1e-6))
    b = StopWhenAll(StopAfterIteration(200), StopWhenChangeLess(1e-6), StopWhenGradientNormLess(1e-6))
"""
function Base.:&(s1::S, s2::T) where {S<:StoppingCriterion,T<:StoppingCriterion}
    return StopWhenAll(s1, s2)
end
function Base.:&(s1::S, s2::StopWhenAll) where {S<:StoppingCriterion}
    return StopWhenAll(s1, s2.criteria...)
end
function Base.:&(s1::StopWhenAll, s2::T) where {T<:StoppingCriterion}
    return StopWhenAll(s1.criteria..., s2)
end

@doc raw"""
    StopWhenAny <: StoppingCriterion

store an array of [`StoppingCriterion`](@ref) elements and indicates to stop,
when _any_ single one indicates to stop. The `reason` is given by the
concatenation of all reasons (assuming that all non-indicating return `""`).

# Constructor
    StopWhenAny(c::NTuple{N,StoppingCriterion} where N)
    StopWhenAny(c::StoppingCriterion...)
"""
mutable struct StopWhenAny{TCriteria<:Tuple} <: StoppingCriterionSet
    criteria::TCriteria
    reason::String
    StopWhenAny(c::Vector{StoppingCriterion}) = new{typeof(tuple(c...))}(tuple(c...), "")
    StopWhenAny(c::StoppingCriterion...) = new{typeof(c)}(c)
end
function (c::StopWhenAny)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if any(subC -> subC(p, o, i), c.criteria)
        c.reason = string([get_reason(subC) for subC in c.criteria]...)
        return true
    end
    return false
end
get_stopping_criteria(c::StopWhenAny) = c.criteria

"""
    |(s1,s2)
    s1 | s2

Combine two [`StoppingCriterion`](@ref) within an [`StopWhenAny`](@ref).
If either `s1` (or `s2`) is already an [`StopWhenAny`](@ref), then `s2` (or `s1`) is
appended to the list of [`StoppingCriterion`](@ref) within `s1` (or `s2`)

# Example
    a = StopAfterIteration(200) | StopWhenChangeLess(1e-6)
    b = a | StopWhenGradientNormLess(1e-6)

Is the same as

    a = StopWhenAny(StopAfterIteration(200), StopWhenChangeLess(1e-6))
    b = StopWhenAny(StopAfterIteration(200), StopWhenChangeLess(1e-6), StopWhenGradientNormLess(1e-6))
"""
function Base.:|(s1::S, s2::T) where {S<:StoppingCriterion,T<:StoppingCriterion}
    return StopWhenAny(s1, s2)
end
function Base.:|(s1::S, s2::StopWhenAny) where {S<:StoppingCriterion}
    return StopWhenAny(s1, s2.criteria...)
end
function Base.:|(s1::StopWhenAny, s2::T) where {T<:StoppingCriterion}
    return StopWhenAny(s1.criteria..., s2)
end

"""
    StopWhenCostLess <: StoppingCriterion

store a threshold when to stop looking at the cost function of the
optimization problem from within a [`Problem`](@ref), i.e `get_cost(p,o.x)`.

# Constructor

    StopWhenCostLess(ε)

initialize the stopping criterion to a threshold `ε`.
"""
mutable struct StopWhenCostLess <: StoppingCriterion
    threshold::Float64
    reason::String
    StopWhenCostLess(ε::Float64) = new(ε, "")
end
function (c::StopWhenCostLess)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if i > 0 && get_cost(p, o.x) < c.threshold
        c.reason = "The algorithm reached a cost function value ($(get_cost(p,o.x))) less than the threshold ($(c.threshold)).\n"
        return true
    end
    return false
end

"""
    StopAfter <: StoppingCriterion

store a threshold when to stop looking at the complete runtime. It uses
`time_ns()` to measure the time and you provide a `Period` as a time limit,
i.e. `Minute(15)`

# Constructor

    StopAfter(t)

initialize the stopping criterion to a `Period t` to stop after.
"""
mutable struct StopAfter <: StoppingCriterion
    threshold::Period
    reason::String
    start::Nanosecond
    function StopAfter(t::Period)
        return if value(t) < 0
            error("You must provide a positive time period")
        else
            new(t, "", Nanosecond(0))
        end
    end
end
function (c::StopAfter)(p::P, o::O, i::Int) where {P<:Problem,O<:Options}
    if value(c.start) == 0 || i <= 0 # (re)start timer
        c.start = Nanosecond(time_ns())
    else
        cTime = Nanosecond(time_ns()) - c.start
        if i > 0 && (cTime > Nanosecond(c.threshold))
            c.reason = "The algorithm ran for about $(floor(cTime, typeof(c.threshold))) and has hence reached the threshold of $(c.threshold).\n"
            return true
        end
    end
    return false
end

"""
    are_these_stopping_critera_active(c::StoppingCriterion, cond)

Return `true` if any criterion from the given set is both active and fulfils the given
condition `cond` (`cond(c)` returns `true`).
"""
function are_these_stopping_critera_active(cond::Function, c::StoppingCriterionSet)
    return any(cc -> are_these_stopping_critera_active(cond, cc), get_stopping_criteria(c))
end
function are_these_stopping_critera_active(cond::Function, c::StoppingCriterion)
    return !isempty(c.reason) && cond(c)
end

@doc raw"""
    get_active_stopping_criteria(c)

returns all active stopping criteria, if any, that are within a
[`StoppingCriterion`](@ref) `c`, and indicated a stop, i.e. their reason is
nonempty.
To be precise for a simple stopping criterion, this returns either an empty
array if no stop is indicated or the stopping criterion as the only element of
an array. For a [`StoppingCriterionSet`](@ref) all internal (even nested)
criteria that indicate to stop are returned.
"""
function get_active_stopping_criteria(c::sCS) where {sCS<:StoppingCriterionSet}
    c = get_active_stopping_criteria.(get_stopping_criteria(c))
    return vcat(c...)
end
# for non-array containing stopping criteria, the recursion ends in either
# returning nothing or an 1-element array containing itself
function get_active_stopping_criteria(c::sC) where {sC<:StoppingCriterion}
    if c.reason != ""
        return [c] # recursion top
    else
        return []
    end
end


text:
function gcd(a, b)
    if b == 0
        return a
    end
    return gcd(b, a % b)
end

max_stars_repo_path: Utility Functions/Algorithms/GCD.jl

max_issues_repo_path: Utility Functions/Algorithms/GCD.jl

max_forks_repo_path: Utility Functions/Algorithms/GCD.jl


text:
"""
Defines common methods whose implementation belongs to a number of different
packages.

Should aim to reduce the number of methods in this package.
"""
module common
  importall Base
  export  compose!, rep, prepare, prepare!, get, set, scale!, breed!,
          compare, execute

  compare(::Any) = error("Common stub method: compare")

  execute(::Any) = error("Common stub method: execute")

  """
  Stub method for scale! functions.
  """
  scale!(::Any) = error("Common stub method: scale!")

  """
  Stub method for prepare functions.
  """
  prepare(::Any) = error("Common stub method: prepare")

  """
  Stub method for prepare! functions.
  """
  prepare!(::Any) = error("Common stub method: prepare!")

  """
  Stub method for breed! functions.
  """
  breed!(::Any) = error("Common stub method: breed!")

  """
  Should compose an object from its provided definition.
  """
  compose!(x::Any) = x

  """
  Should return the representation used by some given object.
  """
  rep(::Any) = error("Common stub method: rep")

  get(::Any) = error("Common stub method: get")

  set(::Any) = error("Common stub method: set")
end


text:
#find out the test statistic when m0=1
import LatentGaussianMixtureModel
import Distributions, StatsBase
import Yeppp

@everywhere using LatentGaussianMixtureModel
@everywhere using Distributions, Yeppp, StatsBase, Random

#Brun calculate the statistic for one data set;
#b is the the random number seed, from 1 to 100
@everywhere function Brun(b::Integer, X::Matrix{Float64},
    groupindex::Vector{Int}; debuginfo::Bool=false, ntrials::Int=5)

    mu_true = [-2.528,-0.969]
    wi_true =  [.05, .95]
    sigmas_true = [.234, .244]
    betas_true=[0.019498134778826972,0.007110782868926634,0.030919753747641616,-0.26957245719168493,-0.5257588591653842,-0.6314135703754931,-0.799509687993618,0.07786105411635978,0.12059632526169425,0.22504173963135948]

    #Randomly data generation based on the setting on datagen.jl
    Random.seed!(b * 100)
    n = maximum(groupindex)
    m = MixtureModel(map((u, v) -> Normal(u, v), mu_true, sigmas_true), wi_true)
    gamma_true = rand(m, n)

    prob = exp(gamma_true[groupindex] .+ X*betas_true)
    prob= prob ./ (1 .+ prob)
    Y = Bool[rand(Binomial(1, prob[i])) == 1 for i in 1:length(groupindex)];

    lr=LatentGaussianMixtureModel.EMtest(X, Y, groupindex, 1, ntrials=ntrials, debuginfo=debuginfo, ctauparallel=false, ngh=100)
    println("Mission $b completed with lr=$lr")
    lr
end