text:
#=
Detailed Oxidative phosphorylation model by Gauthier et al. (2013)
With ROS generation
Default parameter values from Kembro et al. and Gauthier et al.
Some are adjusted by An-chi Wei to prevent complex number
=#
using Parameters
include("common.jl")
include("concentrations.jl")

const E_FMN = -0.375  # FMN redox potential (mV)
const E_SOX = -0.15  # Superoxide redox potential (mV)
const E_ROS = _ra(E_SOX - E_FMN)  # Voltage factor between Superoxide and FMN

# Doxorubicin effects
@with_kw struct DOXParams
    DOX = 0.0  # Doxorubicin (DOX) concentration (mM)
    NI = 3  # Hills coefficient of DOX inhibition on OXPHOS complex
    KI_COX_C1 = 400e-3  # DOX inhibition concentration (IC50) on complex I
    KI_COX_C2 = 2000e-3 # DOX inhibition concentration (IC50) on complex II
    KI_COX_C3 = 185e-3  # DOX inhibition concentration (IC50) on complex III
    KI_COX_C4 = 165e-3  # DOX inhibition concentration (IC50) on complex IV
    K_RC = 1E3 / 15  # DOX redox cycling constant
    C1_INHIB = _hills(KI_COX_C1, DOX, NI)  # complex I inhibitory scale
    C2_INHIB = _hills(KI_COX_C2, DOX, NI) # complex II inhibitory scale
    C3_INHIB = _hills(KI_COX_C3, DOX, NI)  # complex III inhibitory scale
    C4_INHIB = _hills(KI_COX_C4, DOX, NI)  # complex IV inhibitory scale
    E_LEAK_C1 = 1 + K_RC * DOX  # Electron leak scaling factor from complex I
end

# OXPHOS (complex I)
@with_kw struct C1Params
    SCALE = 1E3  # Convert 1/s to 1/ms
	ρC1 = 24.9302  # Concentration of complex I (mM), from Gauthier et al. (2013)
    dpsi_B = 50.0  # Phase boundary potential (mV)
    K12 = 6.3396E11 / SCALE
    K21 = 5 / SCALE
    K56 = 100 / SCALE
    K65 = 2.5119E13 / SCALE
    K61 = 3.329E7 / SCALE
    K16 = 432.7642 / SCALE
    K23 = 3886.7 / SCALE
    K32 = 9.1295E6 / SCALE
    K34 = 639.1364 / SCALE
    K43 = 3.2882 / SCALE
    K47 = 1.5962E7 / SCALE
    K74 = 65.2227 / SCALE
    K75 = 2.4615E4 / SCALE
    K57 = 1.1667E3 / SCALE
    K42 = 6.0318 / SCALE
end

function c1_rates(nadh, nad, Q_n, QH2_n, dpsi, h_m, sox_m,
                  dpsi_B, K12, K21, K65, K56, K61, K16, K23, K32,
                  K34, K43, K47, K74, K57, K75, K42, ρC1,
                  MT_PROT=1, E_LEAK_C1=1, C1_INHIB=1, h_i=1E-4, O2=6E-3)
    C1_CONC = ρC1 * MT_PROT
	# A hack to wrokaround negative concentration
	nadh = max(0.0, nadh)
	nad = max(0.0, nad)
	Q_n = max(0.0, Q_n)
	QH2_n = max(0.0, QH2_n)
    a = _ra(dpsi - dpsi_B)
    a12 = K12 * h_m^2
    a21 = K21
    a65 = K65 * h_i^2
    a56 = K56
    a61 = K61 / a
    a16 = K16 * a
    a23 = K23 * sqrt(nadh)
    a32 = K32
    a34 = K34
    a43 = K43 * sqrt(nad)
    a47 = C1_INHIB * K47 * sqrt(Q_n * h_m)
    a74 = K74
    a57 = C1_INHIB * K57 * sqrt(QH2_n)
    a75 = K75
    a42 = K42 * E_LEAK_C1 * O2
    a24 = K42 * O2 * E_ROS * sox_m
	a52 = a25 = 0

    # Fraction of each state in Complex I from the matlab code, derived from KA pattern
	e1 = a21*a32*a42*a56*a61*a74 + a21*a32*a42*a56*a61*a75 + a21*a32*a42*a57*a61*a74 + a21*a32*a42*a57*a65*a74 + a21*a32*a43*a56*a61*a74 + a21*a32*a43*a56*a61*a75 + a21*a32*a43*a57*a61*a74 + a21*a32*a43*a57*a65*a74 + a21*a32*a47*a56*a61*a75 + a21*a34*a42*a56*a61*a74 + a21*a34*a42*a56*a61*a75 + a21*a34*a42*a57*a61*a74 + a21*a34*a42*a57*a65*a74 + a21*a34*a47*a56*a61*a75 + a23*a34*a47*a56*a61*a75 + a24*a32*a47*a56*a61*a75 + a24*a34*a47*a56*a61*a75
	e2 = a12*a32*a42*a56*a61*a74 + a12*a32*a42*a56*a61*a75 + a12*a32*a42*a57*a61*a74 + a12*a32*a42*a57*a65*a74 + a12*a32*a43*a56*a61*a74 + a12*a32*a43*a56*a61*a75 + a12*a32*a43*a57*a61*a74 + a12*a32*a43*a57*a65*a74 + a12*a32*a47*a56*a61*a75 + a12*a34*a42*a56*a61*a74 + a12*a34*a42*a56*a61*a75 + a12*a34*a42*a57*a61*a74 + a12*a34*a42*a57*a65*a74 + a12*a34*a47*a56*a61*a75 + a16*a32*a42*a57*a65*a74 + a16*a32*a43*a57*a65*a74 + a16*a34*a42*a57*a65*a74
	e3 = a12*a23*a42*a56*a61*a74 + a12*a23*a42*a56*a61*a75 + a12*a23*a42*a57*a61*a74 + a12*a23*a42*a57*a65*a74 + a12*a23*a43*a56*a61*a74 + a12*a23*a43*a56*a61*a75 + a12*a23*a43*a57*a61*a74 + a12*a23*a43*a57*a65*a74 + a12*a23*a47*a56*a61*a75 + a12*a24*a43*a56*a61*a74 + a12*a24*a43*a56*a61*a75 + a12*a24*a43*a57*a61*a74 + a12*a24*a43*a57*a65*a74 + a16*a21*a43*a57*a65*a74 + a16*a23*a42*a57*a65*a74 + a16*a23*a43*a57*a65*a74 + a16*a24*a43*a57*a65*a74
	e4 = a12*a23*a34*a56*a61*a74 + a12*a23*a34*a56*a61*a75 + a12*a23*a34*a57*a61*a74 + a12*a23*a34*a57*a65*a74 + a12*a24*a32*a56*a61*a74 + a12*a24*a32*a56*a61*a75 + a12*a24*a32*a57*a61*a74 + a12*a24*a32*a57*a65*a74 + a12*a24*a34*a56*a61*a74 + a12*a24*a34*a56*a61*a75 + a12*a24*a34*a57*a61*a74 + a12*a24*a34*a57*a65*a74 + a16*a21*a32*a57*a65*a74 + a16*a21*a34*a57*a65*a74 + a16*a23*a34*a57*a65*a74 + a16*a24*a32*a57*a65*a74 + a16*a24*a34*a57*a65*a74
	e5 = a12*a23*a34*a47*a61*a75 + a12*a23*a34*a47*a65*a75 + a12*a24*a32*a47*a61*a75 + a12*a24*a32*a47*a65*a75 + a12*a24*a34*a47*a61*a75 + a12*a24*a34*a47*a65*a75 + a16*a21*a32*a42*a65*a74 + a16*a21*a32*a42*a65*a75 + a16*a21*a32*a43*a65*a74 + a16*a21*a32*a43*a65*a75 + a16*a21*a32*a47*a65*a75 + a16*a21*a34*a42*a65*a74 + a16*a21*a34*a42*a65*a75 + a16*a21*a34*a47*a65*a75 + a16*a23*a34*a47*a65*a75 + a16*a24*a32*a47*a65*a75 + a16*a24*a34*a47*a65*a75
	e6 = a12*a23*a34*a47*a56*a75 + a12*a24*a32*a47*a56*a75 + a12*a24*a34*a47*a56*a75 + a16*a21*a32*a42*a56*a74 + a16*a21*a32*a42*a56*a75 + a16*a21*a32*a42*a57*a74 + a16*a21*a32*a43*a56*a74 + a16*a21*a32*a43*a56*a75 + a16*a21*a32*a43*a57*a74 + a16*a21*a32*a47*a56*a75 + a16*a21*a34*a42*a56*a74 + a16*a21*a34*a42*a56*a75 + a16*a21*a34*a42*a57*a74 + a16*a21*a34*a47*a56*a75 + a16*a23*a34*a47*a56*a75 + a16*a24*a32*a47*a56*a75 + a16*a24*a34*a47*a56*a75
	e7 = a12*a23*a34*a47*a56*a61 + a12*a23*a34*a47*a57*a61 + a12*a23*a34*a47*a57*a65 + a12*a24*a32*a47*a56*a61 + a12*a24*a32*a47*a57*a61 + a12*a24*a32*a47*a57*a65 + a12*a24*a34*a47*a56*a61 + a12*a24*a34*a47*a57*a61 + a12*a24*a34*a47*a57*a65 + a16*a21*a32*a42*a57*a65 + a16*a21*a32*a43*a57*a65 + a16*a21*a32*a47*a57*a65 + a16*a21*a34*a42*a57*a65 + a16*a21*a34*a47*a57*a65 + a16*a23*a34*a47*a57*a65 + a16*a24*a32*a47*a57*a65 + a16*a24*a34*a47*a57*a65
    denom = e1 + e2 + e3 + e4 + e5 + e6 + e7
    f2 = e2 / denom
    f4 = e4 / denom
    f7 = e7 / denom
	# Electrons flow through C1
    jResC1 = C1_CONC * (f4 * a47 - f7 * a74)
    # ROS produced in C1
    vROSC1 = C1_CONC * (f4 * a42 - f2 * a24)
    # Q reduction rate from C1
    vc1 = 0.5 * jResC1
    return vc1, vROSC1
end

function c1_rates(nadh, Q_n, QH2_n, dpsi, h_m, sox_m, conc::Concentrations, pC1::C1Params, pDOX::DOXParams, MT_PROT=1)
	@unpack  (dpsi_B, K12, K21, K65, K56, K61, K16, K23, K32,
	K34, K43, K47, K74, K57, K75, K42, ρC1) = pC1
	@unpack NAD_T, h_i, O2 = conc
	@unpack E_LEAK_C1, C1_INHIB = pDOX
	nad = _nad(nadh, conc)
	return c1_rates(nadh, nad, Q_n, QH2_n, dpsi, h_m, sox_m,
	         dpsi_B, K12, K21, K65, K56, K61, K16, K23, K32,
	         K34, K43, K47, K74, K57, K75, K42, ρC1,
	         MT_PROT, E_LEAK_C1, C1_INHIB, h_i, O2)
end

# OXPHOS (complex II)
@with_kw struct C2Params
    VMAX = 250 / 60E3  # Maximal rate of SDH = complex II (mM/ms)
    KM_Q = 0.6  # MM constant for CoQ
    KI_OAA = 0.15  # MM constant for OAA
	FAC_SUC = 0.085
	KI_FUM = 1.3
	KM_SUC = 0.03
	ET = 0.5
	KCAT = 3E-3
end

#=
Rate of complex II, aka succinate dehydrogenase (SDH)
SUC + Q -> (via FAD/FADH2) -> FUM + QH2
with a magical SUC_SCALE term from the matlab code
adapted from Demin et al. (2000)
=#
function vsdh(Q_n, QH2_n, suc, fum, oaa, KI_OAA, VMAX, KM_Q, FAC_SUC, C2_INHIB=1)
    f_q = _mm(_mm(Q_n, QH2_n), KM_Q)
    f_oaa = _mm(KI_OAA, oaa)
    f_suc = FAC_SUC * sqrt(suc / fum)
    vSDH = VMAX * C2_INHIB * f_suc * f_oaa * f_q
end

function vsdh(Q_n, QH2_n, suc, fum, oaa, pC2::C2Params, C2_INHIB=1)
	@unpack KI_OAA, VMAX, KM_Q, FAC_SUC = pC2
	vSDH = vsdh(Q_n, QH2_n, suc, fum, oaa, KI_OAA, VMAX, KM_Q, FAC_SUC, C2_INHIB)
end

# SDH rate combines Cortassa et al. (2003) and Demin et al. (2000)
function vsdh_mm(Q_n, QH2_n, suc, fum, oaa, pC2::C2Params, C2_INHIB=1)
	@unpack KI_OAA, VMAX, KM_Q, KI_FUM, KM_SUC = pC2
	f_q = _mm(_mm(Q_n, QH2_n), KM_Q)
	f_oaa = _mm(KI_OAA, oaa)
	f_fum = _mm(KI_FUM, fum)
	f_suc = _mm(suc * f_oaa * f_fum, KM_SUC)
	vSDH = C2_INHIB * VMAX * f_suc * f_q
end

# OXPHOS (complex IV)
@with_kw struct C4Params
	# Concentration of complex IV (mM)
    ρC4 = 0.325
    SCALE = 60E3  # 1/min to 1/ms
    δ₅ = 0.5
    K34 = 1.7667E28 / SCALE
    K43 = 1.7402 / SCALE
    K35 = 45000 / SCALE
    K36 = 2.8955E25 / SCALE
    K63 = 2.8955E10 / SCALE
    K37 = 1.7542E12 / SCALE
    K73 = 1.7542E4 / SCALE
end

# C4 rates (scalar version)
function c4_rates(cytc_ox, dpsi, h_m, conc::Concentrations, pC4::C4Params, C4_INHIB=1, MT_PROT=1)
	@unpack K34, K73, K43, K35, δ₅, K36, K37, K63, ρC4 = pC4
	@unpack h_i, O2 = conc
    C4_CONC = ρC4 * MT_PROT

    aδ = _ra(-δ₅ * dpsi)
    a1mδ = _ra((1 - δ₅) * dpsi)
	f_hm = h_m * aδ
	f_hi = h_i * a1mδ
	f_cr = cytc_rd = C4_CONC - cytc_ox
	f_co = cytc_ox * a1mδ
    # Transition rates
    a12 = K34 * f_cr^3 * f_hm^4  # a12 = K34 * exp(-δ₅ * 4 * vfrt) * cytc_rd^3 * h_m^4
    a14 = K73 * f_hi  # K73 * exp((1 - δ₅) * F_RT * dpsi) * h_i
    a21 = K43 * f_co^3 * f_hi  # K43 * exp((1 - δ₅) * 4 * vfrt) * cytc_ox^3 * h_i
    a23 = K35 * O2 * C4_INHIB
    a34 = K36 * f_cr * f_hm^3  # K36 * exp(-δ₅ * 3 * vfrt) * cytc_rd * h_m^3
    a41 = K37 * f_hm  # K37 * exp(-δ₅ * vfrt) * h_m
    a43 = K63 * f_co * f_hi^2  # K63 * exp((1 - δ₅) * 3 * vfrt) * cytc_ox * h_i^2

    # Weight of each state (from KA pattern)
    e1 = a21 * a41 * a34 + a41 * a34 * a23
    e2 = a12 * a41 * a34
    e3 = a23 * a12 * a41 + a43 * a14 * a21 + a23 * a43 * a12 + a23 * a43 * a14
    e4 = a14 * a34 * a21 + a34 * a23 * a12 + a34 * a23 * a14
    den = e1 + e2 + e3 + e4

    # Fraction of each state
    y = e1 / den
    yr = e2 / den
    yo = e3 / den
    yoh = e4 / den
    # Reaction rates
    v34 = C4_CONC * (y * a12 - yr * a21)
    v35 = C4_CONC * yr * a23
    v36 = C4_CONC * (a34 * yo - a43 * yoh)
    v37 = C4_CONC * (a41 * yoh - a14 * y)

    # cytc -> complex IV electron flux
    vE = 3 * v34 + v35
    # Proton flux by C4
	# proton pumped by C1 + C4
	vH = 4 * v34 + 3 * v36 + v37
	# proton pumped by C4 alone
    vHresC4 = v34 + 2 * v36 + v37
    # O2 consumption
    vO2 = v35
    return vE, vHresC4, vO2
end

# Complex V (ATP synthase) parameters
@with_kw struct C5Params
    ρF1 = 5.0  # Concentration of complex V (mM)
    P1 = 1.346E-8
    P2 = 7.739E-7
    P3 = 6.65E-15
    PA = 1.656E-8
    PB = 3.373E-10
    PC1 = 9.651E-17
    PC2 = 4.585E-17
    KEQ = 1.71E6  # Equilibrium constant of ATP synthase
	E150 = _ra(3 * 50)
end

# Without Binding corrections
function _vaf1(pi_m, adp_m, atp_m, KEQ)
	vaf1 = KEQ * atp_m / (pi_m * adp_m)
end

# Relative electrochemical activit of ATP <-> ADP + Pi, corrected by dissociation
function _vaf1(h_m, pi_m, adp_m, atp_m, mg_m, KEQ)
    # Dissociation
    (adp3_m, hadp_m, _, poly_adp_m) = breakdown_adp(adp_m, h_m, mg_m)
    (_, _, mgatp_m, poly_atp_m) = breakdown_atp(atp_m, h_m, mg_m)
	ke_app = KEQ * h_m * poly_atp_m * _h2o_poly(h_m) / (poly_adp_m * _pi_poly(h_m))
	vaf1 = _vaf1(pi_m, adp3_m + hadp_m, mgatp_m, ke_app)
end

function _c5_rates(ΔμH, vaf1, pC5::C5Params, MT_PROT=1, C5_INHIB=1)
	@unpack ρF1, P1, P2, P3, PA, PB, PC1, PC2, KEQ, E150 = pC5
	vh = _ra(3 * ΔμH)
    common = -ρF1 * MT_PROT * C5_INHIB / ((1 + P1 * vaf1) * E150 + (P2 + P3 * vaf1) * vh)
    vATPase = common * ((100 * PA + PC1 * E150) * vaf1 - (PA + PC2 * vaf1) * vh)
    vHu = 3 * common * (100 * PA * (1 + vaf1) - vh * (PA + PB))
	return (vATPase, vHu)
end

# ATP synthase (cpmplex V) and related rates from Kembro et. al (2013)
# vATPase: ATP production rate
# vHu : protons flux through ATP synthase
function c5_rates(dpsi, ΔμH, h_m, pi_m, adp_m, conc::Concentrations, pC5::C5Params, MT_PROT=1, C5_INHIB=1)
	@unpack KEQ = pC5
	@unpack mg_m = conc
    # Potential difference of ATP <-> ADP + Pi
	vaf1 = _vaf1(h_m, pi_m, adp_m, _atp_m(adp_m, conc), mg_m, KEQ)
	return _c5_rates(ΔμH, vaf1, pC5, MT_PROT, C5_INHIB)
end

# From Li et al. 2015
function c5_rates(ΔμH, pi_m, adp_m, conc::Concentrations, pC5::C5Params, MT_PROT=1, C5_INHIB=1)
	@unpack KEQ = pC5
	vaf1 = _vaf1(pi_m, adp_m, _atp_m(adp_m, conc), KEQ)
	return _c5_rates(ΔμH, vaf1, pC5, MT_PROT, C5_INHIB)
end

@with_kw struct ANTParams
	VMAX = 0.435  # Max rate of ANT (mM/ms), (from AC Wei)
    H = 0.5  # Voltage steepness
end

# Rate of ATP/ADP translocator (ANT)
function _vant(atp4_m, adp3_m, atp4_i, adp3_i, dpsi, H, VMAX)
	f_i = atp4_i / adp3_i
	f_m = adp3_m / atp4_m
	# Magnus_Keizer_1997
	# vANT = VMAX * (1 - f_i * f_m * _ra(-dpsi)) / ((1 + f_i * _ra(-H * dpsi)) * (1 + f_m))
	# Cortassa_2003
	vANT = VMAX * (1 - f_i * f_m) / ((1 + f_i * _ra(-H * dpsi)) * (1 + f_m))
	# Li_2015
	# vANT = VMAX * _ra(-dpsi) * (1 - f_i * f_m) / ((1 + f_i * _ra(-H * dpsi)) * (1 + f_m))
end

# Rate of ATP/ADP translocator (ANT), with auto correction from dissociation
function vant(atp_m, adp_m, atp_i, adp_i, mg_m, mg_i, h_m, h_i, dpsi, H, VMAX)
	(atp4_i, _, _, _) = breakdown_atp(atp_i, h_i, mg_i)
	(atp4_m, _, _, _) = breakdown_atp(atp_m, h_m, mg_m)
	(adp3_i, _, _, _) = breakdown_adp(adp_i, h_i, mg_i)
	(adp3_m, _, _, _) = breakdown_adp(adp_m, h_m, mg_m)
	vANT = _vant(atp4_m, adp3_m, atp4_i, adp3_i, dpsi, H, VMAX)
end

function vant(adp_m, atp_i, h_m, dpsi, conc::Concentrations, pANT::ANTParams)
	@unpack H, VMAX = pANT
	@unpack h_i, mg_m, mg_i = conc
	atp_m = _atp_m(adp_m, conc)
	adp_i = _adp_i(atp_i, conc)
	vANT = vant(atp_m, adp_m, atp_i, adp_i, mg_m, mg_i, h_m, h_i, dpsi, H, VMAX)
end

# Without correction
function vant(adp_m, atp_i, dpsi, conc::Concentrations, pANT::ANTParams)
	@unpack H, VMAX = pANT
	atp_m = _atp_m(adp_m, conc)
	adp_i = _adp_i(atp_i, conc)
	atp4_i = 0.25 * atp_i
	atp4_m = 0.025 * atp_m
	adp3_i = 0.45 * adp_i
	adp3_m = 0.17 * adp_m
	return _vant(atp4_m, adp3_m, atp4_i, adp3_i, dpsi, H, VMAX)
end

# Proton leak across inner mitochondrial membrane
vhleak(ΔμH, G_H_MITO) = G_H_MITO * ΔμH

# OXPHOS (complex III)
@with_kw struct C3Params
	SCALE = 60e3  # 1/min to 1/ms
	K03 = 9.9998E4 / SCALE
    KEQ3 = 0.6877
    K04 = 3.6402E3 / SCALE
    KEQ4_OX = 129.9853
    KEQ4_RD = 13.7484
    δ₁ = 0.5
    δ₂ = 0.5
    δ₃ = 0.5
    β = 0.5006
	α = 0.5 * (1 - β)
    γ = α
    KD_Q = 1.32E6 / SCALE
    K06 = 10000 / SCALE
    KEQ6 = 9.4546
    K07_OX = 800 / SCALE
    K07_RD = 100 / SCALE
    KEQ7_OX = 3.0748
    KEQ7_RD = 29.0714
    K08_OX = 5000 / SCALE
    K08_RD = 500 / SCALE
    KEQ8_OX = 129.9853
    KEQ8_RD = 9.4546
    K09 = 4.9949E4 / SCALE
    KEQ9 = 0.2697
    K010 = 1700 / SCALE
    KEQ10 = 1.4541
    K33 = 148148 / SCALE
    KEQ33 = 2.1145
    ρC3 = 0.325
    Q_T = 4.0  # Total CoQ pool
end

_v3(A1, Q1, A2, Q2, K0, KEQ) = K0 * (KEQ * A1 * Q1 - A2 * Q2)
_v3(fes_ox, QH2_p, fes_rd, Qdot_p, K03, KEQ3, FAC_PH, MYXOTHIAZOLE) = _v3(fes_ox * FAC_PH, QH2_p, fes_rd * MYXOTHIAZOLE, Qdot_p, K03, KEQ3)
_el(α, δ, dpsi) = _ra(-α * δ * dpsi)
_er(α, δ, dpsi) = _ra(α * (1 - δ) * dpsi)
_v4(bleft, Qdot_p, el, bright, Q_p, er, K04, KEQ4) = _v3(bleft, Qdot_p * el, bright, Q_p * er, K04, KEQ4)
_v6(b2, b3, dpsi, K06, KEQ6, β, δ₂) = _v3(b2, _el(β, δ₂, dpsi), b3, _er(β, δ₂, dpsi), K06, KEQ6)
_v7(bleft, Q_n, el, bright, Qdot_n, er, K07, KEQ7, ANTIMYCIN) = _v4(bleft, Q_n, el, bright, Qdot_n * ANTIMYCIN, er, K07, KEQ7)
_v8(bleft, Qdot_n, el, bright, QH2_n, er, K08, KEQ8, FAC_PH, ANTIMYCIN) = _v7(bleft, Qdot_n * FAC_PH^2, el, bright, QH2_n, er, K08, KEQ8, ANTIMYCIN)
_v9(fes_rd, cytc1_ox, fes_ox, cytc1_rd, K09, KEQ9) = _v3(fes_rd, cytc1_ox, fes_ox, cytc1_rd, K09, KEQ9)
_v10(O2, Qdot_p, sox_m, Q_p, K010, KEQ10) = _v3(O2, Qdot_p, sox_m, Q_p, K010, KEQ10)

# OXPHOS ODE system
# TODO:
# 1. conservation of Q pool (4.0 mM)
# 1.5 Fast equilibrium of (Q_n, Q_p) and (QH2_n, QH2_p) (v2 and v5 are much faster)
# 2. conservation of BH/BL pool (same as C3_CONC)
function oxphos_system(nadh, Q_n, Qdot_n, QH2_n, QH2_p, Qdot_p, Q_p, b1, b2, b3, b4, fes_ox,
	cytc1_ox, cytc_ox, dpsi, ΔμH, h_m, sox_m, suc, fum, oaa, pi_m, adp_m, atp_i, conc::Concentrations,
	pC1::C1Params, pC2::C2Params, pC3::C3Params, pC4::C4Params, pC5::C5Params,
	pDOX::DOXParams, pANT::ANTParams, MT_PROT=1)
	@unpack O2, FAC_PH = conc
	@unpack C2_INHIB, C3_INHIB, C4_INHIB = pDOX
	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(dpsi, ΔμH, h_m, pi_m, adp_m, conc, pC5, MT_PROT)
	vANT = vant(adp_m, atp_i, h_m, dpsi, conc, pANT)
	# Complex III ODE system
	@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
	d_nadh_ox = -vc1
    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
	vHres = 4 * vc1 + 2 * v3 + vHresC4
	return (d_nadh_ox, 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, vHres, vROSC1, vSDH, vROSC3, vO2,
			vATPase, vHu, vANT)
end


text:
# SPDX-License-Identifier: MIT

struct GSM{ModelType} <: AbstractModel
    transform::AffineTransform
    f::Function
    vec2params::Function
end

function newmodel(gsm::GSM{ModelType}, embeddings) where ModelType
    E = [statistics(e) for e in embeddings]
    H = gsm.f.(E |> gsm.transform)
    [ModelType(gsm.vec2params(η)...) for η in H]
end

#######################################################################
# Model interface

function loglikelihood(gsm::GSM, embeddings, Xs)
    models = newmodel(gsm, embeddings)
    [sum(loglikelihood(models[k], Xs[k])) for k in 1:length(models)]
end

ext: jl

max_stars_repo_name: lucasondel/PPCA

max_issues_repo_name: lucasondel/PPCA

max_forks_repo_name: lucasondel/PPCA


text:
"""
    RRTPFlash{T}(;K0 = nothing,rtol= 1e-10,atol = 1e-10,max_iters = 100)

Method to solve non-reactive multicomponent flash problem by using successive substitution in the Rachford-Rice equation.

Only two phases are supported. if K0 is `nothing`, it will be calculated via the Wilson correlation.

The optimizer will stop at `max_iters` evaluations, when the absolute tolerance is less than `atol` or when the relative tolerance is less than `rtol`
"""
Base.@kwdef struct RRTPFlash{T} <: TPFlashMethod
    K0::T = nothing
    rtol::Float64 = 1e-10
    atol::Float64 = 1e-10
    max_iters::Int = 100
    spec::Symbol = :unknown
end

#z is the original feed composition, x is a matrix with molar fractions, n is a matrix with molar amounts

function tp_flash_impl(model::EoSModel, p, T, n, method::RRTPFlash)
    if method.K0===nothing
        pure = split_model.(model)
        crit = crit_pure.(pure)
        Tc = [crit[i][1] for i ∈ @comps]
        pc = [crit[i][2] for i ∈ @comps]
        ω = acentric_factor.(pure)
        K0 = @. exp(log(pc/p)+5.373*(1+ω)*(1-Tc/T))
    else
        K0 = method.K0
    end

    atol = method.atol
    rtol = method.rtol
    max_iters = method.max_iters

    x = zero(n)
    y = zero(n)
    α = zeros(eltype(K0),1)
    φ_α = zero(K0)
    φ_β = zero(K0)
    f(Kout,Kin) = RR_fixpoint(Kout,Kin,model,p,T,n,x,y,φ_α,φ_β,α,method.spec)

    Solvers.fixpoint(f,K0,Solvers.SSFixPoint();atol,rtol,max_iters)
    α₀ = α[1]
    G = (gibbs_free_energy(model,p,T,x)*(1-α₀)+gibbs_free_energy(model,p,T,y)*α₀)/R̄/T

    X = hcat(x,y)'
    nvals = X.*[1-α₀
                α₀]
    return (X, nvals, G)
end

function RR_fixpoint(K1,K0,model,p,T,n,x,y,φ_α,φ_β,α,spec)
    α₀ = rr_vle_vapor_fraction(K0,n)
    α[1] = α₀
    x = rr_flash_liquid!(x,K0,n,α₀)
    y .= K0 .* x
    if spec == :lle
        phaseα = :l
        phaseβ = :l
    elseif spec ==:vle
        phaseα = :l
        phaseβ = :v
    else
        phaseα = :unknown
        phaseβ = :unknown
    end

    φ_α .= fugacity_coefficient(model,p,T,x,phase = phaseα)
    φ_β .= fugacity_coefficient(model,p,T,y,phase = phaseβ)
    K1 .= φ_α./φ_β
    K1
end

export RRTPFlash


text:
function plmdca(filename::String;
                decimation::Bool=false,
                boolmask::Union(Array{Bool,2},Nothing)=nothing,
                fracmax::Float64 = 0.3,
                fracdec::Float64 = 0.1,
                remove_dups::Bool = true,
                max_gap_fraction::Real = 0.9,
                theta = :auto,
                lambdaJ::Real=0.005,
                lambdaH::Real=0.01,
                gaugecol::Int=-1,
                epsconv::Real=1.0e-5,
                maxit::Int=1000,
                verbose::Bool=true,
                method::Symbol=:LD_LBFGS)

    W,Z,N,M,q = ReadFasta(filename,max_gap_fraction, theta, remove_dups)

    boolmask != nothing && size(boolmask) != (N,N) && error("size boolmask different from ( $N, $N )")

    plmalg = PlmAlg(method,verbose, epsconv ,maxit, boolmask)
    plmvar = PlmVar(N,M,q,q*q,gaugecol,lambdaJ,lambdaH,Z,W)

    if decimation  == false
        Jmat, pslike = MinimizePLAsym(plmalg,plmvar)
    else
        plmalg.boolmask != nothing && println("Warning: decimation with boolmask not implemented. Continuing ...")
        decvar = DecVar{2}(fracdec, fracmax, ones(Bool, (N-1)*q*q, N))
        Jmat, pslike = DecimateAsym!(plmvar, plmalg, decvar)
    end

    score, FNAPC, Jtensor = ComputeScore(Jmat, plmvar)
    return output = PlmOut{4}(sdata(pslike), Jtensor, score)
end

function MinimizePLAsym(alg::PlmAlg, var::PlmVar)

    LL = (var.N - 1) * var.q2 + var.q
    x0 = zeros(Float64, LL)
    vecps = SharedArray(Float64,var.N)
    Jmat = @parallel hcat for site=1:var.N #1:12
        function f(x::Vector, g::Vector)
            g === nothing && (g = zeros(Float64, length(x)))
            return PLsiteAndGrad!(x, g, site,  var)
        end

        opt = Opt(alg.method, length(x0))
        ftol_abs!(opt, alg.epsconv)
        maxeval!(opt, alg.maxit)
        if alg.boolmask != nothing # constrain to zero boolmask variables
            lb,ub=ComputeUL(alg,var,site,LL)
            lower_bounds!(opt, lb)
            upper_bounds!(opt, ub)
        end
        min_objective!(opt, f)
        elapstime = @elapsed  (minf, minx, ret) = optimize(opt, x0)
        alg.verbose && @printf("site = %d\t pl = %.4f\t time = %.4f\t", site, minf, elapstime)
        alg.verbose && println("exit status = $ret")
        vecps[site] = minf
        minx
    end
    return Jmat, vecps
end

function ComputeUL(alg::PlmAlg, var::PlmVar, site::Int, LL::Int)

    boolmask = alg.boolmask
    N  = var.N
    q2 = var.q2
    lb = -Inf * ones(Float64, LL)
    ub =  Inf * ones(Float64, LL)
    tiny::Float64 = 1.0e-6
    offset::Int = 0

    for i=1:site-1
        if boolmask[i,site] == false
            for s = 1:q2
                lb[offset + s] = -tiny
                ub[offset + s] =  tiny
            end
        end
        offset += q2
    end
    for i=site+1:N
        if boolmask[i, site] == false
            for s = 1:q2
                lb[offset + s] = -tiny
                ub[offset + s] =  tiny
            end
        end
        offset += q2
    end
    return lb,ub
end

function PLsiteAndGrad!(vecJ::Array{Float64,1},  grad::Array{Float64,1}, site::Int, plmvar::PlmVar)

    LL = length(vecJ)
    q2 = plmvar.q2
    q = plmvar.q
    gaugecol = plmvar.gaugecol
    N = plmvar.N
    M = plmvar.M
    Z = sdata(plmvar.Z)
    W = sdata(plmvar.W)

    for i=1:LL-q
        grad[i] = 2.0 * plmvar.lambdaJ * vecJ[i]
    end
    for i=(LL-q+1):LL
       grad[i] = 2.0 * plmvar.lambdaH * vecJ[i]
    end

    vecene = zeros(Float64,q)
    expvecenesunorm = zeros(Float64,q)
    pseudolike = 0.0

    @inbounds begin
        for a = 1:M
            fillvecene!(vecene, vecJ,site,a, q, Z,N)
            norm = sumexp(vecene)
            expvecenesunorm = exp(vecene .- log(norm))
            pseudolike -= W[a] * ( vecene[Z[site,a]] - log(norm) )
            offset = 0
            for i = 1:site-1
                @simd for s = 1:q
                    grad[ offset + s + q * ( Z[i,a] - 1 ) ] += W[a] *  expvecenesunorm[s]
                end
                grad[ offset + Z[site,a] + q * ( Z[i,a] - 1 ) ] -= W[a]
                offset += q2
            end
	    for i = site+1:N
                @simd for s = 1:q
                    grad[ offset + s + q * ( Z[i,a] - 1 ) ] += W[a] *  expvecenesunorm[s]
                end
                grad[ offset + Z[site,a] + q * ( Z[i,a] - 1 ) ] -= W[a]
                offset += q2
            end

            @simd for s = 1:q
                grad[ offset + s ] += W[a] *  expvecenesunorm[s]
            end
	    grad[ offset + Z[site,a] ] -= W[a]
        end
    end

    if 1 <= gaugecol <= q
        offset = 0;
        @inbounds begin
            for i=1:N-1
                for s=1:q
                    grad[offset + gaugecol + q * (s - 1)  ] = 0.0; # Gauge!!! set gradJ[a,q] = 0
                    grad[offset + s + q * (gaugecol - 1) ] = 0.0; # Gauge!!! set gradJ[q,a] = 0
                end
                offset += q2
            end
            grad[offset + gaugecol] = 0.0 # Gauge!!! set gradH[q] = 0
        end
    end

    pseudolike += L2norm_asym(vecJ, plmvar)

    return pseudolike
end

function fillvecene!(vecene::Array{Float64,1}, vecJ::Array{Float64,1},site::Int, a::Int, q::Int, sZ::DenseArray{Int,2},N::Int)

    q2 = q*q
    Z = sdata(sZ)
    @inbounds begin
        for l = 1:q
            offset::Int = 0
            scra::Float64 = 0.0
            for i = 1:site-1 # Begin sum_i \neq site J
                scra += vecJ[offset + l + q * (Z[i,a]-1)]
                offset += q2
            end
            # skipping sum over residue site
    	    for i = site+1:N
                scra += vecJ[offset + l + q * (Z[i,a]-1)]
                offset += q2
            end # End sum_i \neq site J
            scra += vecJ[offset + l] # sum H
            vecene[l] = scra
        end
    end
end

function L2norm_asym(vec::Array{Float64,1}, plmvar::PlmVar)
    q = plmvar.q
    N = plmvar.N
    lambdaJ = plmvar.lambdaJ
    lambdaH = plmvar.lambdaH

    LL = length(vec)

    mysum1 = 0.0
    @inbounds @simd for i=1:(LL-q)
        mysum1 += vec[i] * vec[i]
    end
    mysum1 *= lambdaJ

    mysum2 = 0.0
    @inbounds @simd for i=(LL-q+1):LL
        mysum2 += vec[i] * vec[i]
    end
    mysum2 *= lambdaH

    return mysum1+mysum2
end


text:
using Documenter, Flux, NNlib

makedocs(modules=[Flux, NNlib],
         sitename = "Flux",
         pages = ["Home" => "index.md",
                  "Building Models" =>
                    ["Basics" => "models/basics.md",
                     "Recurrence" => "models/recurrence.md",
                     "Regularisation" => "models/regularisation.md",
                     "Model Reference" => "models/layers.md",
                     "Advanced Model Building" => "models/advanced.md",
                     "NNlib" => "models/nnlib.md"],
                  "Handling Data" =>
                    ["One-Hot Encoding" => "data/onehot.md",
                     "DataLoader" => "data/dataloader.md"],
                  "Training Models" =>
                    ["Optimisers" => "training/optimisers.md",
                     "Training" => "training/training.md"],
                  "GPU Support" => "gpu.md",
                  "Saving & Loading" => "saving.md",
                  "The Julia Ecosystem" => "ecosystem.md",
                  "Performance Tips" => "performance.md",
                  "Community" => "community.md"],
         format = Documenter.HTML(assets = ["assets/flux.css"],
                                  analytics = "UA-36890222-9",
                                  prettyurls = haskey(ENV, "CI")))

deploydocs(repo = "github.com/FluxML/Flux.jl.git",
           target = "build",
           push_preview = true)


text:
srcfile(src::LineNumberNode) = abspath(string(src.file))
srcline(src::LineNumberNode) = src.line

function _check_file(obj, srcnode)

    # book
    book = _check_currbook()
    bdir = bookdir(book)

    src = srcfile(srcnode)
    if !isfile(src)
        startswith(basename(src), "REPL") &&
            error("A new object can't be created from the repl. Use a `.br.jl` file at the books folder!!!")
        error("A new object must be defined in a `.br.jl` file under the book directory!!")
    end
    abspath(srcfile(obj)) != abspath(src)

    return src
end

function _extract_expr_pairs(ex)

    ex = (ex isa Tuple) ? ex : tuple(ex)

    pairs = Dict()
    for subex in ex
        (subex isa LineNumberNode) && continue
        if Meta.isexpr(subex, :(=))
            key, val = subex.args

            if key isa Symbol
                pairs[key] = val
                continue
            end
        end
        error("Unsupported syntax at:\n", subex,
            "\n\nUsage:\n@macro(A = \"Bla\", B = [1,2,3]) or",
            "\n@macro A = \"Bla\"  B = [1,2,3]\n"
        )
    end
    return pairs
end

function _capture_output(dofunc::Function)
    tfile = tempname()
    try
        touch(tfile)
        open(tfile, "w") do out
            redirect_stdout(out) do
                redirect_stderr(out) do
                    dofunc()
                end
            end
        end
        return read(tfile, String)
    finally
        rm(tfile; force = true)
    end
end

function _spaceless_unsensitive_regex(hint)
    hint = strip(hint)
    hint = replace(hint, isspace => "\\s*")
    return Regex(string("^(?:", "(?i)", "\\s*", hint, "\\s*", ")", "\$"))
end

# ---------------------------------------------------
function _has_tag(tags::Set, hint)
    for tag in tags
        has_match(tag, hint) && return true
    end
    return false
end
_has_tag(hint) = (tags) ->  _has_tag(tags, hint)


text:
__precompile__()

module PlotRecipes

using Reexport
@reexport using Plots
@reexport using StatPlots

include("graphs.jl")
include("finance.jl")
include("misc.jl")

end # module


text:
# Object Definition
mutable struct Account
    account_number::String
    balance::Float64
    date_opened::Date
end

# Accessors
account_number(a::Account) = a.account_number
balance(a::Account) = a.balance
date_opened(a::Account) = a.date_opened

# Functions

function deposit!(a::Account, amount::Real)
    a.balance += amount
    return a.balance
end

function withdraw!(a::Account, amount::Real)
    a.balance -= amount
    return a.balance
end

function transfer!(from::Account, to::Account, amount::Real)
    # println("Transferring ", amount, " from account ",
    #     account_number(from), " to account ", account_number(to))
    withdraw!(from, amount)
    deposit!(to, amount)
    return amount
end


text:
#
#Simple example of a Poisson problem:
#
#           y
#             |
#             |           u_y = 0
#          Ly ----------------------------------
#             |                                |
#             |                                |
#             |                                |
#       u = 0 |    u - u_xx - u_yy = f         | u = 0
#             |                                |
#             |                                |
#             |                                |
#           -----------------------------------|-----  x
#                         u_y = 0               Lx
#
#Here,
#
#      u = sin(omegax x) * cos(omegay y)
#      f = ( 1 + omegax^2 + omegay^2 ) * sin(omegax x) * cos(omegay y)
#      omegax = nx pi / Lx
#      omegay = ny pi / Ly
#

using FinElt
using FinElt.PlanarPoisson
using Printf

include("params.jl")

const omegax = 3 * pi / Lx
const omegay = pi / Ly

function exact_u(x)
    return sin(omegax*x[1]) * cos(omegay*x[2])
end

function f(x)
    c = 1 + omegax^2 + omegay^2
    return c * exact_u(x)
end

ell_f(v, z) = P1SourceTimesFunc!(v, z, f)
lin_functionals = [ ("Omega", ell_f) ]

maxerr = zeros(refinements+1)
@printf("%10s  %12s  %8s  %8s\n\n",
        "N", "max error", "rate", "seconds")
for k = 0:refinements
    start = time()
    mesh = read_msh_file("rect$k.msh")
    vp = VariationalProblem(mesh, ["Left", "Right"])
    add_bilin_form!(vp, "Omega", grad_dot_grad!)
    add_bilin_form!(vp, "Omega", func_times_func!)
    add_lin_functnl!(vp, "Omega", source_times_func!, f)
    # Implicit zero boundary conditions.
    A, b = assembled_linear_system(vp)
    ufree = A \ b
    uh = complete_soln(ufree, vp)
    u = get_nodal_vals(exact_u, mesh)
    finish = time()
    maxerr[k+1] = maximum(abs.(uh-u))
    N = length(ufree)
    if k == 0
        @printf("%10d  %12.4e\n", N, maxerr[k+1])
    else
        rate = log2(maxerr[k]/maxerr[k+1])
        elapsed = finish - start
        @printf("%10d  %12.4e  %8.4f  %8.4f\n",
                N, maxerr[k+1], rate, elapsed)
    end
end


text:
module DMUStudent

using Nettle: hexdigest

export HW6
#     HW1,
#     HW2,
#     HW3,
#     HW4,
#     HW5,

is_identikey_colorado_email(email) = !isnothing(match(r"^[a-z]{4}\d{4}@colorado.edu$", email))

function hash_score(hw, email, score, key)
    return hexdigest("sha256", string(hw)*email*string(score)*key)
end

# include("HW1.jl")
# include("HW2.jl")
# include("HW3.jl")
# include("HW4.jl")
# include("HW5.jl")
include("HW6.jl")

end # module


text:
############################################################################
# main call function
############################################################################
"""
    calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT[; grad_list=[undef], rfphase_increment=[π], m0=:antiperiodic, output=:complexsignal])

Calculate the signal or magnetization evolution with the linear approximation of the generalized Bloch model assuming a super-Loretzian lineshape.

The simulation assumes a sequence of rectangluar RF-pulses with varying flip angles α and RF-pulse durations TRF, but a fixed repetition time TR. Further, it assumes balanced gradient moments.

# Arguments
- `α::Vector{<:Number}`: Array of flip angles in radians
- `TRF::Vector{<:Number}`: Array of the RF-pulse durations in seconds
- `TR::Number`: Repetition time in seconds
- `ω0::Number`: Off-resonance frequency in rad/s
- `B1::Number`: Normalized transmit B1 field, i.e. B1 = 1 corresponds to a well-calibrated B1 field
- `m0s::Number`: Fractional size of the semi-solid pool; should be in range of 0 to 1
- `R1f::Number`: Longitudinal relaxation rate of the free pool in 1/seconds
- `R2f::Number`: Transversal relaxation rate of the free pool in 1/seconds
- `Rx::Number`: Exchange rate between the two spin pools in 1/seconds
- `R1f::Number`: Longitudinal relaxation rate of the semi-solid pool in 1/seconds
- `T2s::Number`: Transversal relaxationt time of the semi-solid pool in seconds
- `R2slT::NTuple{3, Function}`: Tuple of three functions: R2sl(TRF, ω1, B1, T2s), dR2sldB1(TRF, ω1, B1, T2s), and R2sldT2s(TRF, ω1, B1, T2s). Can be generated with [`precompute_R2sl`](@ref)

Optional:
- `grad_list=[undef]`: Vector to indicate which gradients should be calculated; the vector elements can either be `undef` for no gradient, or any subset/order of `grad_list=[grad_m0s(), grad_R1f(), grad_R2f(), grad_Rx(), grad_R1s(), grad_T2s(), grad_ω0(), grad_B1()]`; the derivative wrt. to apparent `R1a = R1f = R1s` can be calculated with `grad_R1a()`
- `rfphase_increment=[π]::Vector{<:Number}`: Increment of the RF phase between consequtive pulses. The default value `π`, together with ``ω0=0`` corresponds to the on-resonance condition. When more than one value is supplied, their resulting signal is stored along the second dimension of the output array.
- `m0=:antiperiodic`: With the default keyword `:antiperiodic`, the signal and their derivatives are calcualted assuming ``m(0) = -m(T)``, where `T` is the duration of the RF-train. With the keyword :thermal, the magnetization ``m(0)`` is initialized with thermal equilibrium `[xf, yf, zf, xs, zs] = [0, 0, 1-m0s, 0, m0s]`, followed by a α[1]/2 - TR/2 prep pulse; and with the keyword `:IR`, this initalization is followed an inversion pulse of duration `TRF[1]`, (set `α[1]=π`) and a α[2]/2 - TR/2 prep pulse.
- `output=:complexsignal`: The defaul keywords triggers the function to output a complex-valued signal (`xf + 1im yf`); the keyword `output=:realmagnetization` triggers an output of the entire (real valued) vector `[xf, yf, zf, xs, zs]`

# Examples
```jldoctest
julia> R2slT = precompute_R2sl();

julia> calculatesignal_linearapprox(ones(100)*π/2, ones(100)*5e-4, 4e-3, 0, 1, 0.1, 1, 15, 30, 6.5, 10e-6, R2slT)
100×1×1 Array{ComplexF64, 3}:
[:, :, 1] =
 -0.029305987774458298 - 2.609023307487357e-18im
  0.004329424678273616 - 9.22858709360654e-19im
  -0.02276140937384383 + 1.3684406715212983e-18im
  0.008280330224850318 - 3.913888063486245e-18im
  -0.01675130572786822 + 4.6409553994856164e-18im
  0.011921649143708501 - 6.302663734240743e-18im
 -0.011190579890418322 + 7.310583447081588e-18im
   0.01530560844095237 - 8.175490153983983e-18im
 -0.006026791818066622 + 9.464626759459701e-18im
  0.018463027499697637 - 9.608537136498097e-18im
                       ⋮
   0.06153147350966063 + 1.496274325112204e-18im
  0.060814007683572674 + 5.6874647068595874e-18im
  0.061725751593187314 + 1.6576028430814674e-18im
   0.06106421622262947 + 5.557937883636948e-18im
   0.06190745705315383 + 1.808443885793045e-18im
   0.06129750535994443 + 5.436837755621598e-18im
  0.062077399758002784 + 1.949397602286564e-18im
   0.06151502191944254 + 5.3236925282029116e-18im
  0.062236337703062715 + 2.0810405785135138e-18im
```
"""
function calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, T2s, R2slT; grad_list=[undef], rfphase_increment=[π], m0=:antiperiodic, output=:complexsignal)

    R2s_vec = evaluate_R2sl_vector(abs.(α), TRF, B1, T2s, R2slT, grad_list)

    return calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, R2s_vec; grad_list=grad_list, rfphase_increment=rfphase_increment, m0=m0, output=output)
end

function calculatesignal_linearapprox(α, TRF, TR, ω0, B1, m0s, R1f, R2f, Rx, R1s, R2s_vec; grad_list=[undef], rfphase_increment=[π], m0=:antiperiodic, output=:complexsignal)

    ω1 = α ./ TRF

    if output == :complexsignal && grad_list == [undef]
        S = Array{ComplexF64}(undef, length(ω1), length(rfphase_increment), 1)
        jj = 1
    elseif output == :complexsignal
        S = Array{ComplexF64}(undef, length(ω1), length(rfphase_increment), 1 + length(grad_list))
        jj = [[1,j + 1] for j = 1:length(grad_list)]
    elseif output == :realmagnetization && grad_list == [undef]
        S = Array{SVector{6,Float64}}(undef, length(ω1), length(rfphase_increment), length(grad_list))
        jj = 1
    elseif output == :realmagnetization
        S = Array{SVector{11,Float64}}(undef, length(ω1), length(rfphase_increment), length(grad_list))
        jj = 1:length(grad_list)
    end

    if m0 == :thermal || m0 == :IR
        if grad_list == [undef]
            _m0 = @SVector [0, 0, 1 - m0s, 0, m0s, 1]
        else
            _m0 = @SVector [0, 0, 1 - m0s, 0, m0s, 0, 0, 0, 0, 0, 1]
        end
    elseif isa(m0, AbstractVector)
        if grad_list == [undef] && length(m0) < 6
            _m0 = SVector{6}([m0; zeros(5 - length(m0)); 1])
        elseif length(m0) < 11
            _m0 = SVector{11}([m0; zeros(10 - length(m0)); 1])
        else
            error("the m0 vector of length <= 5 w/o gradients or <= 10 with gradients")
        end
    elseif m0 != :antiperiodic
        error("m0 must either be :antiperiodic, :IR, :thermal, or a vector")
    end

    for j in eachindex(grad_list), i in eachindex(rfphase_increment)
        if m0 == :antiperiodic
            m = antiperiodic_boundary_conditions_linear(ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, R2s_vec[1], R2s_vec[2], R2s_vec[3], grad_list[j], rfphase_increment[i])
        elseif m0 == :thermal || isa(m0, AbstractVector)
            # this implements the α[1]/2 - TR/2 preparation (TR/2 is implemented in propagate_magnetization_linear!)
            u_pr = exp(hamiltonian_linear(ω1[1] / 2, B1, ω0, TRF[1], m0s, R1f, R2f, Rx, R1s, R2s_vec[1][1], R2s_vec[2][1], R2s_vec[3][1], grad)) # R2sl is actually wrong for the prep pulse
            m = u_pr * _m0
        elseif m0 == :IR
            # this implements the π - spoiler - TR - α[2]/2 - TR/2 preparation (TR/2 is implemented in propagate_magnetization_linear!)
            u_ip = propagator_linear_inversion_pulse(ω1[1], TRF[1], B1, R2s_vec[1][1], R2s_vec[2][1], R2s_vec[3][1], grad)
            m = u_ip * _m0
            u_fp = xs_destructor(grad) * exp(hamiltonian_linear(0, B1, ω0, TR, m0s, R1f, R2f, Rx, R1s, 0, 0, 0, grad))
            m = u_fp * m
            u_pr = exp(hamiltonian_linear(ω1[2] / 2, B1, ω0, TRF[2], m0s, R1f, R2f, Rx, R1s, R2s_vec[1][2], R2s_vec[2][2], R2s_vec[3][2], grad)) # R2sl is actually wrong for the prep pulse
            m = u_pr * m
        end

        Mij = @view S[:, i, jj[j]]
        propagate_magnetization_linear!(Mij, ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, R2s_vec[1], R2s_vec[2], R2s_vec[3], rfphase_increment[i], m, grad_list[j])
    end
    return S
end

############################################################################
# helper functions
############################################################################
function evolution_matrix_linear(ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, _R2s, _dR2sdT2s, _dR2sdB1, grad=undef, rfphase_increment=π)
    # put inversion pulse at the end (this defines m as the magnetization at the first TE after the inversion pulse)
    u_rot = z_rotation_propagator(rfphase_increment, grad)
    u_fp = xs_destructor(grad) * exp(hamiltonian_linear(0, B1, ω0, TR / 2, m0s, R1f, R2f, Rx, R1s, 0, 0, 0, grad))
    u_pl = propagator_linear_inversion_pulse(ω1[1], TRF[1], B1, _R2s[1], _dR2sdT2s[1], _dR2sdB1[1], grad)
    A = u_fp * u_pl * u_rot * u_fp

    for i = length(ω1):-1:2
        u_fp = xs_destructor(grad) * exp(hamiltonian_linear(0, B1, ω0, (TR - TRF[i]) / 2, m0s, R1f, R2f, Rx, R1s, 0, 0, 0, grad))
        u_pl = exp(hamiltonian_linear(ω1[i], B1, ω0, TRF[i], m0s, R1f, R2f, Rx, R1s, _R2s[i], _dR2sdT2s[i], _dR2sdB1[i], grad))
        A = A * u_fp * u_pl * u_rot * u_fp
    end
    return A
end

function antiperiodic_boundary_conditions_linear(ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, _R2s, _dR2sdT2s, _dR2sdB1, grad=undef, rfphase_increment=π)
    A = evolution_matrix_linear(ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, _R2s, _dR2sdT2s, _dR2sdB1, grad, rfphase_increment)

    Q = A - A0(grad)
    m = Q \ C(grad)
    return m
end

function propagate_magnetization_linear!(S, ω1, B1, ω0, TRF, TR, m0s, R1f, R2f, Rx, R1s, _R2s, _dR2sdT2s, _dR2sdB1, rfphase_increment, m, grad=undef)

    u_rot = z_rotation_propagator(rfphase_increment, grad)
    ms_setindex!(S, m, 1, grad)
    for i = 2:length(ω1)
        u_fp = xs_destructor(grad) * exp(hamiltonian_linear(0, B1, ω0, (TR - TRF[i]) / 2, m0s, R1f, R2f, Rx, R1s, 0, 0, 0, grad))
        u_pl = exp(hamiltonian_linear(ω1[i], B1, ω0, TRF[i], m0s, R1f, R2f, Rx, R1s, _R2s[i], _dR2sdT2s[i], _dR2sdB1[i], grad))

        m = u_fp * (u_pl * (u_rot * (u_fp * m)))

        ms_setindex!(S, m, i, grad)
    end
    return S
end

function ms_setindex!(S::AbstractArray{<:AbstractVector}, y, i, grad)
    S[i] = y
    return S
end
function ms_setindex!(S::AbstractArray{<:Complex}, y, i, grad)
    S[i] = y[1] + 1im * y[2]
    return S
end
function ms_setindex!(S::AbstractArray{<:Complex}, y, i, grad::grad_param)
    S[i,1] = y[1] + 1im * y[2]
    S[i,2] = y[6] + 1im * y[7]
    return S
end


text:
@testset "Biconnect" begin
    gint = SimpleGraph(12)
    add_edge!(gint, 1, 2)
    add_edge!(gint, 2, 3)
    add_edge!(gint, 2, 4)
    add_edge!(gint, 3, 4)
    add_edge!(gint, 3, 5)
    add_edge!(gint, 4, 5)
    add_edge!(gint, 2, 6)
    add_edge!(gint, 1, 7)
    add_edge!(gint, 6, 7)
    add_edge!(gint, 6, 8)
    add_edge!(gint, 6, 9)
    add_edge!(gint, 8, 9)
    add_edge!(gint, 9, 10)
    add_edge!(gint, 11, 12)

    a = [[Edge(3, 5), Edge(4, 5), Edge(2, 4), Edge(3, 4), Edge(2, 3)],
         [Edge(9, 10)],
         [Edge(6, 9), Edge(8, 9), Edge(6, 8)],
         [Edge(1, 7), Edge(6, 7), Edge(2, 6), Edge(1, 2)],
         [Edge(11, 12)]]

    for g in testgraphs(gint)
        bcc = @inferred(biconnected_components(g))
        @test bcc == a
        @test typeof(bcc) === Vector{Vector{Edge{eltype(g)}}}
    end

    g = SimpleGraph(4)
    add_edge!(g, 1, 2)
    add_edge!(g, 2, 3)
    add_edge!(g, 3, 4)
    add_edge!(g, 1, 4)

    h = SimpleGraph(4)
    add_edge!(h, 1, 2)
    add_edge!(h, 2, 3)
    add_edge!(h, 3, 4)
    add_edge!(h, 1, 4)

    gint = blockdiag(g, h)
    add_edge!(gint, 4, 5)

    a = [[Edge(5, 8), Edge(7, 8), Edge(6, 7), Edge(5, 6)], [Edge(4, 5)], [Edge(1, 4), Edge(3, 4), Edge(2, 3), Edge(1, 2)]]

    for g in testgraphs(gint)
        bcc = @inferred(biconnected_components(g))
        @test bcc == a
        @test typeof(bcc) === Vector{Vector{Edge{eltype(g)}}}
    end
end


text:
JuMP.owner_model(cref::RobustConstraintRef)::RobustModel = cref.model

JuMP.index(cref::RobustConstraintRef)::RobustConstraintIndex = cref.index
function Base.:(==)(v::RobustConstraintRef, w::RobustConstraintRef)::Bool
    return v.model === w.model && v.index == w.index
end
Base.broadcastable(cref::RobustConstraintRef) = Ref(cref)

function JuMP.build_constraint(
    _error::Function,
    fun::GenericAffExpr,
    set::MOI.AbstractScalarSet,
    con_type::RobustConstraintType;
    extra_kwargs...
)

    print("buildConst","\n")
    if typeof(con_type) == UncertaintySetConstraint
        return UncertaintySetConst(con_type.name,fun,set)
    elseif typeof(con_type) == UncertainConstraint
        return UncertainConst(con_type.name,fun,set)
    else
        _error("Weird Constraint Type!")
    end
end

function JuMP.add_constraint(
    model::RobustModel,
    con::RobustConstraint,
    name::String,
)

    print("addConst","\n")
    model.nextconidx+=1
    if  occursin("UncertaintySetConst",string(typeof(con)))
        cref = UncertaintySetConstraintRef(model,model.nextconidx,length(model.uncertaintySetConstraints)+1,'S')
        model.uncertaintySetConstraints[cref.type_idx] = con
        model.con_index_to_type_index[model.nextconidx] = length(model.uncertaintySetConstraints)
        model.type_index_to_con_index[length(model.uncertaintySetConstraints)]= model.nextconidx
    elseif occursin("UncertainConst", string(typeof(con)))
        cref = UncertainConstraintRef(model,model.nextconidx,length(model.uncertainConstraints)+1,'U')
        model.uncertainConstraints[cref.type_idx] = con
        model.con_index_to_type_index[model.nextconidx] = length(model.uncertainConstraints)
        model.type_index_to_con_index[length(model.uncertainConstraints)]= model.nextconidx
    end
    JuMP.set_name(cref,name)
    return cref
end


text:
#=
P15 (**) Replicate the elements of a list a given number of times.
Example:
* (repli '(a b c) 3)
(A A A B B B C C C)
=#

# From P12
function my_repeat(s, num)
    function _my_repeat(s, num, repeated)
        num == 0 ? repeated : _my_repeat(s, num - 1, vcat(repeated, s))
    end
    _my_repeat(s, num, [])
end

function repli(coll, num)
    reduce(vcat, [], map(item -> my_repeat(item, num), coll))
end

@assert repli([:a, :b, :c], 3) == [:a, :a, :a, :b, :b, :b, :c, :c, :c]
@assert repli([:a], 5) == [:a, :a, :a, :a, :a]
@assert repli([], 1) == []

println("Tests passed: JL-15.jl")


text:
#
# ADDITION
#
"""
    double_add(a::Float64, b::Float64)

Add two `Float64`s with `DoubleFloat64` precision.
"""
@inline function double_add(a::Float64, b::Float64)
    hi, lo =  two_sum(a, b)

    FastDouble(hi, lo)
end

@inline function +(a::DoubleFloat64{T}, b::Float64) where T
    hi, lo = two_sum(a.hi, b)
    lo += a.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleFloat64{T}(hi, lo)
end
+(a::Float64, b::DoubleFloat64) = b + a
+(a::Integer, b::DoubleFloat64) = b + float(a)
+(a::DoubleFloat64, b::Integer) = a + float(b)

@inline function accurate_add(a::DoubleFloat64, b::DoubleFloat64)
    s1, s2 = two_sum(a.hi, b.hi)
    t1, t2 = two_sum(a.lo, b.lo)
    s2 += t1
    s1, s2 = quick_two_sum(s1, s2);
    s2 += t2
    s1, s2 = quick_two_sum(s1, s2);

    AccurateDouble(s1, s2);
end

@inline function fast_add(a::DoubleFloat64, b::DoubleFloat64)
    hi, lo = two_sum(a.hi, b.hi)
    lo += (a.lo + b.lo)
    hi, lo = quick_two_sum(hi, lo)

    FastDouble(hi, lo);
end

+(a::FastDouble, b::FastDouble) = fast_add(a, b)
+(a::FastDouble, b::AccurateDouble) = accurate_add(a, b)
+(a::AccurateDouble, b::FastDouble) = accurate_add(a, b)
+(a::AccurateDouble, b::AccurateDouble) = accurate_add(a, b)

#
# SUBSTRACTION
#
"""
    double_sub(a::Float64, b::Float64)

Subtract two `Float64`s with `DoubleFloat64` precision.
"""
@inline function double_sub(a::Float64, b::Float64)
    hi, lo = two_diff(a, b)

    FastDouble(hi, lo)
end

@inline function -(a::DoubleFloat64{T}, b::Float64) where T
    hi, lo = two_diff(a.hi, b)
    lo += a.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleFloat64{T}(hi, lo)
end

@inline function -(a::Float64, b::DoubleFloat64{T}) where T
    hi, lo = two_diff(a, b.hi)
    lo -= b.lo
    hi, lo = quick_two_sum(hi, lo)

    DoubleFloat64{T}(hi, lo)
end

@inline function fast_sub(a::DoubleFloat64, b::DoubleFloat64)
    hi, lo = two_diff(a.hi, b.hi)
    lo += a.lo
    lo -= b.lo
    hi, lo = quick_two_sum(hi, lo)

    FastDouble(hi, lo)
end

@inline function accurate_sub(a::DoubleFloat64, b::DoubleFloat64)
    s1, s2 = two_diff(a.hi, b.hi)
    t1, t2 = two_diff(a.lo, b.lo)
    s2 += t1
    s1, s2 = quick_two_sum(s1, s2);
    s2 += t2
    s1, s2 = quick_two_sum(s1, s2);

    AccurateDouble(s1, s2);
end

-(a::FastDouble, b::FastDouble) = fast_sub(a, b)
-(a::FastDouble, b::AccurateDouble) = accurate_sub(a, b)
-(a::AccurateDouble, b::FastDouble) = accurate_sub(a, b)
-(a::AccurateDouble, b::AccurateDouble) = accurate_sub(a, b)

-(a::DoubleFloat64{T}) where T = DoubleFloat64{T}(-a.hi, -a.lo)

#
# MULTIPLICATION
#
"""
    double_mutiply(a::Float64, b::Float64)

Multiply two `Float64`s with `DoubleFloat64` precision.
"""
@inline function double_mul(a::Float64, b::Float64)
    hi, lo = two_prod(a, b)
    FastDouble(hi, lo)
end

@inline function *(a::DoubleFloat64{T}, b::Float64) where T
    p1, p2 = two_prod(a.hi, b)
    p2 += a.lo * b
    p1, p2 = quick_two_sum(p1, p2)

    DoubleFloat64{T}(p1, p2)
end
*(a::Float64, b::DoubleFloat64) = b * a
*(a::Integer, b::DoubleFloat64) = b * float(a)
*(a::Bool, b::DoubleFloat64) = b * convert(Float64, a)
*(a::DoubleFloat64, b::Integer) = a * float(b)

@inline function *(a::DoubleFloat64{T}, b::DoubleFloat64{S}) where {T, S}
    p1, p2 = two_prod(a.hi, b.hi)
    p2 += a.hi * b.lo + a.lo * b.hi
    p1, p2 = quick_two_sum(p1, p2)

    DoubleFloat64{promote_type(T, S)}(p1, p2)
end

#
# DIVISION
#
"""
    double_div(a::Float64, b::Float64)

Divide two `Float64`s with `DoubleFloat64` precision.
"""
@inline function double_div(a::Float64, b::Float64)
    q1 = a / b
    # Compute a - q1 * b
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a, p1)
    e -= p2

    # get next approximation
    q2 = (s + e) / b

    s, e = quick_two_sum(q1, q2)

    FastDouble(s, e)
end

@inline function /(a::DoubleFloat64{T}, b::Float64) where T
    q1 = a.hi / b

    # Compute  this - q1 * d
    p1, p2 = two_prod(q1, b)
    s, e = two_diff(a.hi, p1)
    e += a.lo
    e -= p2

    # get next approximation.
    q2 = (s + e) / b

    # renormalize
    hi, lo = quick_two_sum(q1, q2)

    DoubleFloat64{T}(hi, lo)
end

@inline function fast_div(a::DoubleFloat64, b::DoubleFloat64)
    q1 = a.hi / b.hi  # approximate quotient

    # compute  this - q1 * dd
    r = b * q1
    s1, s2 = two_diff(a.hi, r.hi)
    s2 -= r.lo
    s2 += a.lo

    # get next approximation
    q2 = (s1 + s2) / b.hi

    # renormalize
    hi, lo = quick_two_sum(q1, q2)

    FastDouble(hi, lo)
end

@inline function accurate_div(a::DoubleFloat64, b::DoubleFloat64)
    q1 = a.hi / b.hi  # approximate quotient

    r = a -  b * q1

    q2 = r.hi / b.hi
    r -= q2 * b

    q3 = r.hi / b.hi

    q1, q2 = quick_two_sum(q1, q2)

    AccurateDouble(q1, q2) + q3
end

/(a::FastDouble, b::FastDouble) = fast_div(a, b)
/(a::FastDouble, b::AccurateDouble) = accurate_div(a, b)
/(a::AccurateDouble, b::FastDouble) = accurate_div(a, b)
/(a::AccurateDouble, b::AccurateDouble) = accurate_div(a, b)

/(a::Float64, b::DoubleFloat64) = DoubleFloat64(a) / b
/(a::Integer, b::DoubleFloat64) = float(a) / b
/(a::DoubleFloat64, b::Integer) = a / float(b)

Base.inv(a::DoubleFloat64) = 1.0 / a

Base.rem(a::DoubleFloat64, b::DoubleFloat64) = a - round(a / b) * b
@inline function Base.divrem(a::DoubleFloat64, b::DoubleFloat64)
    n = round(a / b)
    n, a - n * b
end

function Base.mod(x::DoubleFloat64, y::DoubleFloat64)
    n = round(a / b)
    return (a - b * n)
end

#
# POWERS
#

"""
    square(x::DoubleFloat64)

Compute `x * x` in a more efficient way.
"""
@inline function square(a::DoubleFloat64{T}) where T
    p1, p2 = two_square(a.hi)
    p2 += 2.0 * a.hi * a.lo
    p2 += a.lo * a.lo
    hi, lo = quick_two_sum(p1, p2)
    return DoubleFloat64{T}(hi, lo)
end
square(a::Float64) = double_square(a)
"""
    double_square(x::Float64)

Convert `x` to a DoubleFloat64 and then compute `x*x`.
"""
@inline function double_square(a::Float64)
    hi, lo = two_square(a)
    return FastDouble(hi, lo)
end
# Implementation adapted from Base
@inline function power_by_squaring(x::DoubleFloat64, p::Integer)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return square(x)
    end
    P = abs(p)
    t = trailing_zeros(P) + 1
    P >>= t
    while (t -= 1) > 0
        x = square(x)
    end
    y = x
    while P > 0
        t = trailing_zeros(P) + 1
        P >>= t
        while (t -= 1) >= 0
            x = square(x)
        end
        y *= x
    end
    if (p < 0)
        return 1.0 / y
    end
    y
end

^(a::DoubleFloat64, p::Integer) = power_by_squaring(a, p)
^(a::DoubleFloat64, b::Real) = iszero(a) ? one(a) : exp(b * log(a))

@inline function Base.sqrt(a::DoubleFloat64{T}) where T
    #= Strategy:  Use Karp's trick:  if x is an approximation
   to sqrt(a), then

      sqrt(a) = a*x + [a - (a*x)^2] * x / 2   (approx)

   The approximation is accurate to twice the accuracy of x.
   Also, the multiplication (a*x) and [-]*x can be done with
   only half the precision.
   =#
    if a.hi < 0
        throw(DomainError("sqrt will only return a complex result if called with a complex argument."))
    end
    if iszero(a)
        return zero(a)
    end

    x = inv(sqrt(a.hi))
    ax = a.hi * x

    convert(DoubleFloat64{T}, double_add(ax, (a - square(ax)).hi * (x * 0.5)))
end

"""
    double_sqrt(x::Float64)

Convert `x` to a DoubleFloat64 and then compute `sqrt(x)`.
"""
double_sqrt(a::Float64) = sqrt(FastDouble(a))

# TODO: Add n-th root?


text:
module MEstimation # general estimating equations with or without bias-reducting adjustments

using NLsolve
using Optim
using FiniteDiff
using ForwardDiff
using LinearAlgebra
using Distributions
# using InvertedIndices

import Base: show, print
import StatsBase: fit, aic, vcov, coef, coeftable, stderror, CoefTable

export objective_function
export objective_function_template

export estimating_function
export get_estimating_function
export estimating_function_template

export aic
export tic
export vcov
export coef
export coeftable
export fit
export stderror

export slice
# export profile

include("estimating_functions.jl")
include("objective_functions.jl")
include("fit.jl")
include("result_methods.jl")
include("slice.jl")
# include("profile.jl")

end # module


text:
using JuMP
using Gurobi

include("../MMKP/MMKP.jl")
include("../src/JOH.jl")
include("SSIT.jl")

ssit_method(m=60, n_threads=1) = JOH.Matheur.SSIT.make_SSIT_method(
	[.001, .003, .005],
	[m*2,	  m, 	m],
	"ssit", n_threads)

base_method(m=60, n_threads=1) = JOH.Matheur.SSIT.make_SSIT_method(
	[.0001],
	[m*5],
	"base", n_threads)

_make_result_table(res_dir) = try mkpath(res_dir) catch IOError println("dir already exists") end
function run_trial(res_dir, problems)
	bd = joinpath(res_dir, "base")
	sd = joinpath(res_dir, "ssit")
	_make_result_table(bd)
	_make_result_table(sd)

	base_problems = deepcopy(problems)

	for (method, probs, dir) in [
			(ssit_method(), problems, sd),
			(base_method(), base_problems, bd)]
		map(p->set_optimizer(p.model, Gurobi.Optimizer), probs)
		SE.generate_comparison_data2(method, probs, results_dir=dir)
	end
end

problems =  MMKP.load_problems()
resdir="./results/long_mmkp"
mkpath(resdir)

run_trial(resdir, problems)

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
abstract type SymStore end

mutable struct ModuleStore <: SymStore
    name::String
    vals::Dict{String,Any}
    exported::Set{String}
    doc::String
end
ModuleStore(name::String) = ModuleStore(name, Dict{String,Any}(), Set{String}(), "")

struct MethodStore <: SymStore
    file::String
    line::Int
    args::Vector{Tuple{String,String}}
end

struct FunctionStore <: SymStore
    methods::Vector{MethodStore}
    doc::String
end

struct abstractStore <: SymStore
    params::Vector{String}
    doc::String
end

struct primitiveStore <: SymStore
    params::Vector{String}
    doc::String
end

struct structStore <: SymStore
    params::Vector{String}
    fields::Vector{String}
    ts::Vector{String}
    methods::Vector{MethodStore}
    doc::String
end

struct genericStore <: SymStore
    t::String
    params::Vector{String}
    doc::String
end

function _getdoc(x)
    string(Docs.doc(x))
end

function read_methods(x)
    map(methods(x)) do m
        path = isabspath(String(m.file)) ? String(m.file) : Base.find_source_file(String(m.file))
        if path == nothing
            path = ""
        end
        MethodStore(path,
                    m.line,
                    Base.arg_decl_parts(m)[2][2:end])
    end
end

function collect_params(t, params = [])
    if t isa UnionAll
        push!(params, t.var)
        return collect_params(t.body, params)
    else
        return t, params
    end
end

function load_module(m, pkg, depot, out)
    out.doc = string(Docs.doc(m))
    out.exported = Set{String}(string.(names(m)))
    if haskey(depot["manifest"], first(pkg))
        for pkg1 in depot["manifest"][first(pkg)]
            if pkg1["uuid"] == last(pkg)
                for dep in get(pkg1, "deps", [])
                    try
                        depm = getfield(m, Symbol(first(dep)))
                        if !haskey(depot["packages"], last(dep))
                            depot["packages"][last(dep)] = ModuleStore(first(dep))
                            load_module(depm, dep, depot, depot["packages"][last(dep)])
                            out.vals[first(dep)] = first(dep)
                        else
                            out.vals[first(dep)] = first(dep)
                        end
                        # the above make reference to the name of the module, may have to change to uuid
                    catch err
                    end
                end
            end
        end
    end
    for n in names(m, all = true)
        !isdefined(m, n) && continue
        startswith(string(n), "#") && continue
        if Base.isdeprecated(m, n)
        else
            x = getfield(m, n)
            if x isa Function
                out.vals[String(n)] = FunctionStore(read_methods(x), _getdoc(x))
            elseif x isa DataType
                t, p = collect_params(x)
                if t.abstract
                    out.vals[String(n)] = abstractStore(string.(p), _getdoc(x))
                elseif t.isbitstype
                        out.vals[String(n)] = primitiveStore(string.(p), _getdoc(x))
                elseif !(isempty(t.types) || Base.isvatuple(t))
                        out.vals[String(n)] = structStore(string.(p),
                                                     collect(string.(fieldnames(t))),
                                                     string.(collect(t.types)),
                                                     read_methods(x),
                                                     _getdoc(x))
                else
                    out.vals[String(n)] = genericStore("DataType", string.(p), _getdoc(x))
                end
            elseif x isa Module && x != m # include reference to current module
                if parentmodule(x) == m # load non-imported submodules

                    out.vals[String(n)] = ModuleStore(String(n))
                    load_module(x, pkg, depot, out.vals[String(n)])
                end
            else
                out.vals[String(n)] = genericStore(string(typeof(x)), [], _getdoc(x))
            end
        end
    end
    out
end

function import_package(pkg, depot)
    depot["packages"][last(pkg)] = ModuleStore(first(pkg))
    try
        Main.eval(:(import $(Symbol(first(pkg)))))
        m = getfield(Main, Symbol(first(pkg)))
        load_module(m, pkg, depot, depot["packages"][last(pkg)])
    catch err
    end
    return depot["packages"][last(pkg)]
end

function load_core()
    c = Pkg.Types.Context()
    depot = Dict("manifest" => c.env.manifest,
                 "installed" => (VERSION < v"1.1.0-DEV.857" ? c.env.project["deps"] : c.env.project.deps),
                 "packages" => Dict{String,Any}("Base" => ModuleStore("Base"), "Core" => ModuleStore("Core")))

    load_module(Base, "Base"=>"Base", depot, depot["packages"]["Base"])
    load_module(Core, "Core"=>"Core", depot, depot["packages"]["Core"])
    push!(depot["packages"]["Base"].exported, "include")
    # Add special case macros
    depot["packages"]["Base"].vals["@."] = depot["packages"]["Base"].vals["@__dot__"]
    push!(depot["packages"]["Base"].exported, "@.")

    return depot
end

function save_store_to_disc(store, file)
    io = open(file, "w")
    serialize(io, store)
    close(io)
end


text:
"""
    ModeOccupation(state)

    fields:
         - n::Int
         - m::Int
         - state::Vector{Int}
"""
struct ModeOccupation
    n::Int
    m::Int
    state::Vector{Int}
    ModeOccupation(state) = all(state[:] .>= 0) ? new(sum(state), length(state), state) : error("negative photon counts")
end

Base.show(io::IO, i::ModeOccupation) = print(io, "state = ", i.state)

at_most_one_photon_per_bin(state) = all(state[:] .<= 1)
at_most_one_photon_per_bin(r::ModeOccupation) = at_most_one_photon_per_bin(r.state)

isa_subset(subset_modes::Vector{Int}) = (at_most_one_photon_per_bin(subset_modes) && sum(subset_modes) != 0)
isa_subset(subset_modes::ModeOccupation) = isa_subset(subset_modes.state)

"""
    first_modes(n::Int, m::Int)

Create a [`ModeOccupation`](@ref) with `n` photons in the first sites of `m` modes.
"""
first_modes(n::Int,m::Int) = n<=m ? ModeOccupation([i <= n ? 1 : 0 for i in 1:m]) : error("n>m")

"""
    Subset(state::Vector{Int})

Create a mode occupation list with at most one count per mode.

    Fields:
         - n::Int
         - m::Int
         - subset::Vector{Int}
"""
struct Subset
        # basically a mode occupation list with at most one count per mode
        n::Int
        m::Int
        subset::Vector{Int}
        function Subset(state)

                isa_subset(state) ? new(sum(state), length(state), state) : error("invalid subset")
        end
        function Subset(modeocc::ModeOccupation)

                state = modeocc.state
                isa_subset(state) ? new(sum(state), length(state), state) : error("invalid subset")
        end
end

Base.show(io::IO, s::Subset) = print(io, "subset = ", convert(Vector{Int},occupancy_vector_to_partition(s.subset)))

function check_disjoint_subsets(s1::Subset, s2::Subset)
        @argcheck s1.m == s2.m "subsets do not have the same dimension"
        @argcheck all(s1.subset .* s2.subset .== 0) "subsets overlap"
end

function check_subset_overlap(subsets::Vector{Subset})

        for (i,subset_1) in enumerate(subsets)
                for (j,subset_2) in enumerate(subsets)
                        if i>j
                                check_disjoint_subsets(subset_1, subset_2)
                        end
                end
        end

end

"""
    Partition(subsets::Vector{Subset})

Create a partition from multiple [`Subset`](@ref).
"""
struct Partition
        subsets::Vector{Subset}
        n_subset::Int
        m::Int
        function Partition(subsets)
                check_subset_overlap(subsets)
                new(subsets, length(subsets), subsets[1].m)
        end
end

Base.show(io::IO, part::Partition) = begin

    println(io, "partition =")
    for s in part.subsets
        println(io, s)
    end
end

"""
    partition_from_subset_lengths(subset_lengths)

Return a partition from a vector of subset lengths.
"""
function partition_from_subset_lengths(subset_lengths)

        """returns a partition from a vector of subset lengths, such as [2,1] gives Partition([[1,1],[1]])"""

    m = sum(subset_lengths)
    subsets = []

    mode = 1

    for subset_length in subset_lengths
        subset_vector = zeros(Int, m)

        for j in 0:subset_length-1
            if mode+j <= m
                subset_vector[mode+j ] = 1
            end
        end

        Subset(subset_vector)
        mode += subset_length
        push!(subsets, Subset(subset_vector))
    end

    Partition(convert(Vector{Subset}, subsets))

end

"""
    equilibrated_partition_vector(m,n_subsets)

Returns a (most) equilibrated partition possible by euclidian division.

(a problem is that euclidian distribution may give n_subsets or n_subsets+1 if not done like below - here it is the most obvious thing I could think of to get a somewhat equilibrated partition with a constant number of subsets)
"""
function equilibrated_partition_vector(m,n_subsets)

    q = div(m,n_subsets)
    y = n_subsets
    r = rem(m,n_subsets)

    first_part = [q for i in 1:y]
    first_part[1] += r

    first_part

end

equilibrated_mode_occupation(m,n_subsets) = ModeOccupation(equilibrated_partition_vector(m,n_subsets))

"""
    equilibrated_partition(m,n_subsets)

Returns a most equilibrated_partition according to the principles of [`equilibrated_partition_vector`](@ref).
"""
function equilibrated_partition(m,n_subsets)

        partition_from_subset_lengths(equilibrated_partition_vector(m,n_subsets))

end

"""
    occupies_all_modes(part::Partition)

Check wether a partition occupies all modes or not.
"""
function occupies_all_modes(part::Partition)

        """checks if a partition occupies all m modes"""

        occupied_modes = zeros(Int, part.m)
        for s in part.subsets
                occupied_modes .+= s.subset
        end

        all(occupied_modes .== 1)

end

"""
    PartitionOccupancy(counts::ModeOccupation, n::Int, partition::Partition)

   Fields:
        - counts::ModeOccupation
        - partition::Partition
        - n::Int
        - m::Int
"""
struct PartitionOccupancy
        counts::ModeOccupation
        partition::Partition
        n::Int
        m::Int
        function PartitionOccupancy(counts::ModeOccupation, n::Int, partition::Partition)

                @argcheck counts.m == partition.n_subset "counts do not have as many modes as parition has subsets"
                @argcheck sum(counts.state) <= n "more photons at the output than input"
                if occupies_all_modes(partition)
                        @argcheck sum(counts.state) == n "photons lost in a partition that occupies all modes"
                end

                new(counts, partition, n, partition.subsets[1].m)

        end
end

Base.show(io::IO, part_occ::PartitionOccupancy) = begin

        for (i, count) in enumerate(part_occ.counts.state)
                println(io, count," in ", part_occ.partition.subsets[i])
        end

end

function partition_occupancy_to_partition_size_vector_and_counts(part_occ::PartitionOccupancy)

    part = part_occ.partition

    @argcheck occupies_all_modes(part) "need to have a partition that occupies all modes if using Shsnovitch's formulas in partition_expectation_values.jl"

    partition_size_vector = [part.subsets[i].n for i in 1:length(part.subsets)]
    partition_counts = part_occ.counts.state

    partition_size_vector, partition_counts

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CyrusSASL_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CyrusSASL")
JLLWrappers.@generate_main_file("CyrusSASL", UUID("6422fedd-75a7-50c2-a7c3-a11dad25a896"))
end  # module CyrusSASL_jll


text:
function Legendre(n,x)
    if n == 0
        1
    elseif n == 1
        x
    elseif n == 2
        (1/2)*(3x^2-1)
    elseif n == 3
        (1/2)*(5x^3-3x)
    elseif n == 4
        (1/8)*(35x^4-30x^2 + 3)
    elseif n == 5
        (1/8)*(63x^5 - 70x^3 + 15)
    elseif n == 6
        (1/16)*(231x^6 - 315x^4 + 105x^2 - 5)
    elseif n == 7
        (1/16)*(-35x + 315x^3 - 693x^5 + 429x^7)
    elseif n == 8
        (1/128)*(35 - 1260x^2 + 6930x^4 - 12012x^6 + 6435x^8)
    end
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
"""
Contains some functions for calculating cross-correlations, autocorrelations,
and Pearson correlations among rate units in a rate network
"""

using Statistics
import DSP.xcorr

"""
	xcorr_unbiased(x::Vector, y::Vector, maxlag::Int64)

Returns the cross-correlation of two signals normalized by the length of the
available convolution window.

Written by Rainer Engelken (re2365@columbia.edu).
"""
function xcorr_unbiased(x, y, maxlag)
	l = length(x) # we assume same length for x and y
	lags = -maxlag:maxlag
	scale = l .- abs.(lags)
	scale[scale .<= 0] .= 1 # avoid zero division
    xcorr(x, y)[l .+ lags] ./ scale
end

"""
	avg_autocorrs(A::Matrix{T}, maxlag::Int64, wait::Int64=0) where T <: AbstractFloat

Returns the average autocorrelation for the activity of each unit described in
the `N x NT` matrix `A`.

Note the autocorrelations are not mean-subtracted. The `maxlag` parameter
gives the maximum shift of the signal, and the `wait` parameter gives the number
of timesteps in the beginning to disregard.
"""
function avg_autocorrs(A::Matrix{T}, maxlag::Int64, wait::Int64=0) where T <: AbstractFloat

	(N, NT) = size(A)

	unit_autocorrs = mapslices(view(A, :, (wait + 1):NT), dims=2) do x
		return xcorr_unbiased(x, x, maxlag)
	end

	return vec(mean(unit_autocorrs, dims=1))
end


text:
layers(::Type{CHELSA{BioClim}}) = layers(BioClim)
layerkeys(::Type{CHELSA{BioClim}}, args...) = layerkeys(BioClim, args...)

"""
    getraster(source::Type{CHELSA{BioClim}}, [layer]) => Union{Tuple,String}

Download [`CHELSA`](@ref) [`BioClim`](@ref) data from [chelsa-climate.org](https://chelsa-climate.org/).

# Arguments
- `layer`: `Integer` or tuple/range of `Integer` from `$(layers(BioClim))`,
    or `Symbol`s form `$(layerkeys(BioClim))`. Without a `layer` argument, all layers
    will be downloaded, and a `NamedTuple` of paths returned.

Returns the filepath/s of the downloaded or pre-existing files.
"""
getraster(T::Type{CHELSA{BioClim}}, layer::Union{Tuple,Int,Symbol}) = _getraster(T, layer)

_getraster(T::Type{CHELSA{BioClim}}, layers::Tuple) = _map_layers(T, layers)
_getraster(T::Type{CHELSA{BioClim}}, layer::Symbol) = _getraster(T, bioclim_int(layer))
function _getraster(T::Type{CHELSA{BioClim}}, layer::Integer)
    _check_layer(T, layer)
    path = rasterpath(T, layer)
    url = rasterurl(T, layer)
    return _maybe_download(url, path)
end

rastername(::Type{CHELSA{BioClim}}, layer::Integer) = "CHELSA_bio10_$(lpad(layer, 2, "0")).tif"

rasterpath(::Type{CHELSA{BioClim}}) = joinpath(rasterpath(CHELSA), "BioClim")
rasterpath(T::Type{CHELSA{BioClim}}, layer::Integer) = joinpath(rasterpath(T), rastername(T, layer))

rasterurl(::Type{CHELSA{BioClim}}) = joinpath(rasterurl(CHELSA), "chelsa_V1/climatologies/bio/")
rasterurl(T::Type{CHELSA{BioClim}}, layer::Integer) = joinpath(rasterurl(T), rastername(T, layer))

max_forks_repo_name: cesaraustralia/RasterDataSources


text:
```Implement binary search.```

function binary_search(arr,l,r,x)
    if r ≥ l
        mid = l + (r - l) ÷ 2

            if arr[mid] == x
                return mid
            elseif arr[mid] > x
                return binary_search(arr, l, mid-1, x)
            else
                return binary_search(arr,  mid+1, r,  x)
            end
    else
        return -1
    end
end

arr = [1,2,3,6,8,10,12,15]
x=10

# Function call
result = binary_search(arr, 0, length(arr)-1, x)

arr[result]

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
include("ToeplitzOperator.jl")
include("ConstantOperator.jl")
include("TridiagonalOperator.jl")
include("PermutationOperator.jl")

## Operators overrided for spaces

include("Conversion.jl")
include("Multiplication.jl")
include("CalculusOperator.jl")
include("Reverse.jl")


text:
function testCreateComputePipeline(device, outPipelineRef::Ref)
    @scope begin
        shaderFolder = String(@__DIR__) * "/shaders/"
        shader = @autodestroy device.createShaderFromFile(shaderFolder * "test_comp.spv")
        stage = lava.createComputeStage(shader)
        plLayout = @autodestroy lava.createPipelineLayout(device)
        createInfo = lava.ComputePipelineCreateInfo(stage, plLayout)
        pi = lava.ComputePipeline(createInfo)
        @test pi.getLayout() == plLayout
        outPipelineRef[] = pi
    end
    return true
end

function testComputePipeline(device)
    # TODO Add support of compute pipeline to Queue
    # queue = device.namedQueue("compute")
    # cmd = queue.beginCommandBuffer()
    pipRef = Ref{lava.ComputePipeline}()
    @test testCreateComputePipeline(device, pipRef)
    # cmd.bindPipeline(pipRef[])
    # cmd.endCommandBuffer()

    return true
end


text:
include("../src/point_system.jl")

# Score an offensive player
# Missing (not in stats, so can't be considered):
#  - 2pt Conversions
function score_off(; rush_yds = 0, rush_td = 0, rec_yds = 0, rec_td = 0, pass_yds = 0, pass_td = 0, pass_int = 0, sack = 0, fum_l = 0, pr_td = 0, kr_td = 0
)
  (
    div(rush_yds, RUSH_YDS_BATCH) * RUSH_YDS_PTS +
    rush_td * RUSH_TD_PTS +
    div(rec_yds, REC_YDS_BATCH) * REC_YDS_PTS +
    rec_td * REC_TD_PTS +
    div(pass_yds, PASS_YDS_BATCH) * PASS_YDS_PTS +
    pass_td * PASS_TD_PTS +
    pass_int * PASS_INT_PTS +
    sack * SACK_PTS +
    fum_l * FUM_L_PTS +
    pr_td * PR_TD_PTS +
    kr_td * KR_TD_PTS
  )
end

function get_points_allowed_pts(papg::Real)
  cur = PTS_ALLOWED_PTS[1]
  for bin in PTS_ALLOWED_PTS
    if papg < bin[1]
      return cur[2]
    else
      cur = bin
    end
  end
  return cur[2]
end

# Score a team's defense
#  - No blocks
#  - Points allowed is just an average
function score_def_st(; games = 0, papg = 0, kr_td = 0, pr_td = 0, int_td = 0, fum_td = 0, sack = 0, sfty = 0)
  (
    games * get_points_allowed_pts(papg) +
    kr_td * DEF_KR_TD_PTS +
    pr_td * DEF_PR_TD_PTS +
    int_td * INT_TD_PTS +
    fum_td * FUM_TD_PTS +
    sack * DEF_SACK_PTS +
    sfty * DEF_SFTY_PTS
  )
end

function score_k(; xp_made = 0, xp_att = 0, fg_made = 0, fg_att = 0, f0019 = 0, f2029 = 0, f3039 = 0, f4049 = 0, f50 = 0)
  (
    xp_made * XP_MADE_PTS +
    (max(xp_att - xp_made, 0) + max(fg_att - fg_made, 0)) * FG_MISS_PTS +
    (f0019 + f2029 + f3039) * FG_0039_PTS +
    f4049 * FG_4049_PTS +
    f50 * FG_50_PTS
  )
end

max_stars_repo_name: jnu/FantasyFootball

max_issues_repo_name: jnu/FantasyFootball

max_forks_repo_name: jnu/FantasyFootball


text:
"""
    cnse_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{FreeSlip}, ::CNSE2D)

apply free slip boundary condition for momentum
sets reflective ghost point
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{FreeSlip},
    ::CNSE2D,
    ::TurbulenceClosure,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = @SVector [state⁻.ρu[1], state⁻.ρu[2], -0]
    ρu⁺ = ρu⁻ - 2 * n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    state⁺.ρu = @SVector [ρu⁺[1], ρu⁺[2]]

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxGradient, NumericalFluxSecondOrder}, ::Impenetrable{FreeSlip}, ::CNSE2D)

no second order flux computed for linear drag
"""
cnse_boundary_state!(
    ::Union{NumericalFluxGradient, NumericalFluxSecondOrder},
    ::MomentumBC,
    ::CNSE2D,
    ::LinearDrag,
    _...,
) = nothing

"""
    cnse_boundary_state!(::NumericalFluxGradient, ::Impenetrable{FreeSlip}, ::CNSE2D)

apply free slip boundary condition for momentum
sets non-reflective ghost point
"""
function cnse_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{FreeSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ

    ρu⁻ = @SVector [state⁻.ρu[1], state⁻.ρu[2], -0]
    ρu⁺ = ρu⁻ - n⁻ ⋅ ρu⁻ .* SVector(n⁻)

    state⁺.ρu = @SVector [ρu⁺[1], ρu⁺[2]]

    return nothing
end

"""
    shallow_normal_boundary_flux_second_order!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::CNSE2D)

apply free slip boundary condition for momentum
apply zero numerical flux in the normal direction
"""
function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxFirstOrder, ::Impenetrable{NoSlip}, ::CNSE2D)

apply no slip boundary condition for momentum
sets reflective ghost point
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxFirstOrder,
    ::Impenetrable{NoSlip},
    ::CNSE2D,
    ::TurbulenceClosure,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρ = state⁻.ρ
    state⁺.ρu = -state⁻.ρu

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxGradient, ::Impenetrable{NoSlip}, ::CNSE2D)

apply no slip boundary condition for momentum
set numerical flux to zero for U
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxGradient,
    ::Impenetrable{NoSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    state⁺,
    aux⁺,
    n⁻,
    state⁻,
    aux⁻,
    t,
    args...,
)
    FT = eltype(state⁺)
    state⁺.ρu = @SVector zeros(FT, 2)

    return nothing
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::CNSE2D)

apply no slip boundary condition for momentum
sets ghost point to have no numerical flux on the boundary for U
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = -state⁻.ρu
    gradflux⁺.ν∇u = gradflux⁻.ν∇u

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{FreeSlip}, ::CNSE2D)

no mass boundary condition for penetrable
"""
cnse_boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{FreeSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    _...,
) = nothing

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{FreeSlip}, ::CNSE2D)

apply free slip boundary condition for momentum
apply zero numerical flux in the normal direction
"""
function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Penetrable{FreeSlip},
    ::CNSE2D,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
    args...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0])'

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Impenetrable{MomentumFlux}, ::HBModel)

apply kinematic stress boundary condition for momentum
applies free slip conditions for first-order and gradient fluxes
"""
function cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Impenetrable{<:MomentumFlux},
    model::CNSE2D,
    turb::TurbulenceClosure,
    args...,
)
    return cnse_boundary_state!(
        nf,
        Impenetrable(FreeSlip()),
        model,
        turb,
        args...,
    )
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{MomentumFlux}, ::HBModel)

apply kinematic stress boundary condition for momentum
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{<:MomentumFlux},
    model::CNSE2D,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * bc.drag(state⁻, aux⁻, t)'

    return nothing
end

"""
    cnse_boundary_state!(::Union{NumericalFluxFirstOrder, NumericalFluxGradient}, ::Penetrable{MomentumFlux}, ::HBModel)

apply kinematic stress boundary condition for momentum
applies free slip conditions for first-order and gradient fluxes
"""
function cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    ::Penetrable{<:MomentumFlux},
    model::CNSE2D,
    turb::TurbulenceClosure,
    args...,
)
    return cnse_boundary_state!(
        nf,
        Penetrable(FreeSlip()),
        model,
        turb,
        args...,
    )
end

"""
    cnse_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{MomentumFlux}, ::HBModel)

apply kinematic stress boundary condition for momentum
sets ghost point to have specified flux on the boundary for ν∇u
"""
@inline function cnse_boundary_state!(
    ::NumericalFluxSecondOrder,
    bc::Penetrable{<:MomentumFlux},
    shallow::CNSE2D,
    ::TurbulenceClosure,
    state⁺,
    gradflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    aux⁻,
    t,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * bc.drag(state⁻, aux⁻, t)'

    return nothing
end


text:
using ModelingToolkit, OrdinaryDiffEq, Test, LinearAlgebra

@parameters t σ ρ β
@variables x(t) y(t) z(t) k(t)
D = Differential(t)

eqs = [D(D(x)) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

sys′ = ODESystem(eqs)
sys = ode_order_lowering(sys′)

eqs2 = [0 ~ x*y - k,
        D(D(x)) ~ σ*(y-x),
        D(y) ~ x*(ρ-z)-y,
        D(z) ~ x*y - β*z]
sys2 = ODESystem(eqs2, t, [x, y, z, k], parameters(sys′))
sys2 = ode_order_lowering(sys2)
# test equation/varible ordering
ModelingToolkit.calculate_massmatrix(sys2) == Diagonal([1, 1, 1, 1, 0])

u0 = [D(x) => 2.0,
      x => 1.0,
      y => 0.0,
      z => 0.0]

p  = [σ => 28.0,
      ρ => 10.0,
      β => 8/3]

tspan = (0.0,100.0)
prob = ODEProblem(sys,u0,tspan,p,jac=true)
probexpr = ODEProblemExpr(sys,u0,tspan,p,jac=true)
sol = solve(prob,Tsit5())
solexpr = solve(eval(prob),Tsit5())
@test all(x->x==0,Array(sol - solexpr))
#using Plots; plot(sol,vars=(:x,:y))

@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

lorenz1 = ODESystem(eqs,name=:lorenz1)
lorenz2 = ODESystem(eqs,name=:lorenz2)

@variables α(t)
@parameters γ
connections = [0 ~ lorenz1.x + lorenz2.y + α*γ]
connected = ODESystem(connections,t,[α],[γ],systems=[lorenz1,lorenz2])

u0 = [lorenz1.x => 1.0,
      lorenz1.y => 0.0,
      lorenz1.z => 0.0,
      lorenz2.x => 0.0,
      lorenz2.y => 1.0,
      lorenz2.z => 0.0,
      α => 2.0]

p  = [lorenz1.σ => 10.0,
      lorenz1.ρ => 28.0,
      lorenz1.β => 8/3,
      lorenz2.σ => 10.0,
      lorenz2.ρ => 28.0,
      lorenz2.β => 8/3,
      γ => 2.0]

tspan = (0.0,100.0)
prob = ODEProblem(connected,u0,tspan,p)
sol = solve(prob,Rodas5())
@test maximum(sol[2,:] + sol[6,:] + 2sol[1,:]) < 1e-12
#using Plots; plot(sol,vars=(:α,Symbol(lorenz1.x),Symbol(lorenz2.y)))


text:
#=
add:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-06
=#

module Add
    function add(a, b)
        return a + b
    end

    function main()
        a = 5
        b = 9
        println("$a + $b = $(add(a, b))")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Add
    Add.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
const LITERATE_FENCER        = "julialit"
const LITERATE_JULIA_FENCE   = "```$LITERATE_FENCER"
const LITERATE_JULIA_FENCE_L = length(LITERATE_JULIA_FENCE)
const LITERATE_JULIA_FENCE_R = Regex(LITERATE_JULIA_FENCE)

"""
$SIGNATURES

Take a Literate.jl script and transform it into a Franklin-Markdown file.
"""
function literate_to_franklin(rpath::AS)::Tuple{String,Bool}
    startswith(rpath, "/") || throw(
        LiterateRelativePathError("Literate expects a path starting with '/'"))
    # rpath is of the form "/scripts/[path/]tutorial[.jl]"
    # split it so that when recombining it will lead to valid path inc windows
    srpath = split(rpath, '/')[2:end] # discard empty [1] since starts with "/"
    fpath  = joinpath(PATHS[:folder], srpath...)
    # append `.jl` if required
    endswith(fpath, ".jl") || (fpath *= ".jl")
    if !isfile(fpath)
        print_warning("""
            File not found when trying to  convert a literate file at '$fpath'.
            """)
        return "", true
    end
    outpath = joinpath(PATHS[:site], "assets", "literate",
                       srpath[2:end-1]...)
    isdir(outpath) || mkpath(outpath)
    # retrieve the file name
    fname = splitext(splitdir(fpath)[2])[1]
    spath = joinpath(outpath, fname * "_script.jl")
    prev  = ""
    if isfile(spath)
        prev = read(spath, String)
    end
    # don't show Literate's infos
    Logging.disable_logging(Logging.LogLevel(Logging.Info))
    # >> output the markdown
    Literate.markdown(fpath, outpath;
                      config=Dict("codefence" => (LITERATE_JULIA_FENCE => "```")),
                      documenter=false, postprocess=literate_post_process, credit=false)
    # >> output the script
    Literate.script(fpath, outpath; documenter=false,
                      postprocess=s->(MESSAGE_FILE_GEN_LIT * s),
                      name=fname * "_script", credit=false)
    # bring back logging
    Logging.disable_logging(Logging.LogLevel(Logging.Debug))
    # see if things have changed
    haschanged = (read(spath, String) != prev)
    # return path to md file
    return joinpath(outpath, fname * ".md"), haschanged
end

"""
$SIGNATURES

Take a markdown string generated by literate and post-process it to number each
code block and mark them as eval-ed ones.
"""
function literate_post_process(s::String)::String
    isempty(s) && return s
    em   = eachmatch(LITERATE_JULIA_FENCE_R, s)
    buf  = IOBuffer()
    write(buf, "<!--$MESSAGE_FILE_GEN-->\n")
    head = 1
    c    = 1
    for m in em
        write(buf, SubString(s, head, prevind(s, m.offset)))
        write(buf, "```julia:ex$c\n")
        head = nextind(s, m.offset + LITERATE_JULIA_FENCE_L)
        c   += 1
    end
    lis = lastindex(s)
    head < lis && write(buf, SubString(s, head, lis))
    return String(take!(buf))
end


text:
#=
Stores a single object of type `V` for each iteration.
`iterations` keep track of the indices (e.g. times) of the iterations.
=#
mutable struct History{I,V} <: UnivalueHistory{I}
    lastiter::I
    iterations::Vector{I}
    values::Vector{V}

    function History(::Type{V}, ::Type{I} = Int) where {I,V}
        new{I,V}(typemin(I), Array{I}(undef, 0), Array{V}(undef, 0))
    end
end

Base.length(history::History) = length(history.iterations)
Base.enumerate(history::History) = zip(history.iterations, history.values)
Base.first(history::History) = history.iterations[1], history.values[1]
Base.last(history::History) = history.iterations[end], history.values[end]
Base.get(history::History) = history.iterations, history.values
Base.eltype(history :: History) = eltype(history.values);
Base.values(h :: History) = h.values;
indices(h :: History) = h.iterations;

function Base.getindex(h :: History, idx)
    idxOut = findfirst(x -> isequal(x, idx),  indices(h));
    if isnothing(idxOut)
        return nothing
    else
        return values(h)[idxOut];
    end
end

"""
	$(SIGNATURES)

Make a history into a `Dict{String,ValueType}`.
The purpose is mainly to represent the history in a format that can be encoded as JSON.
"""
function history_to_dict(h :: History{I, V}) where {I, V}
    d = Dict{String, V}();
    n = length(h);
    idxV, valueV = get(h);
    for j = 1 : n
        d[string(idxV[j])] = valueV[j]
    end
    return d
end

"""
	$(SIGNATURES)

Make a `UnivalueHistory` from a vector of values.
"""
function make_history(idxV, valueV :: Vector)
    h = History(eltype(valueV));
    for (j, v) in enumerate(valueV)
        push!(h, idxV[j], v);
    end
    return h
end

function Base.push!(
        history::History{I,V},
        iteration::I,
        value::V) where {I,V}
    lastiter = history.lastiter
    iteration > lastiter || throw(ArgumentError("Iterations must increase over time"))
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

# Add a history entry. Simply increment iteration counter by 1.
function Base.push!(
        history::History{I,V},
        value::V) where {I,V}
    lastiter = history.lastiter == typemin(I) ? zero(I) : history.lastiter
    iteration = lastiter + one(history.lastiter)
    history.lastiter = iteration
    push!(history.iterations, iteration)
    push!(history.values, value)
    value
end

Base.print(io::IO, history::History{I,V}) where {I,V} = print(io, "$(length(history)) elements {$I,$V}")

function Base.show(io::IO, history::History{I,V}) where {I,V}
    println(io, "History")
    println(io, "  * types: $I, $V")
    print(io,   "  * length: $(length(history))")
end

"""
    increment!(trace, iter, val)

Increments the value for a given iteration if it exists, otherwise adds the iteration with an ordinary push.
"""
function increment!(trace::History{I,V}, iter::Number, val)  where {I,V}
    if !isempty(trace.iterations)
        if trace.lastiter == iter # Check most common case to make it faster
            i = length(trace.iterations)
        else
            i = findfirst(isequal(iter), trace.iterations)
        end
        if i != nothing
            return (trace.values[i] += val)
        end
    end
    push!(trace, iter, val)
end

# ----------------


text:
# Common univariate distributions.
distset = [
    Normal(1.23, 3.45),
    Exponential(2.12),
    Gamma(4.3),
    Beta(2.123),
    LogNormal(4.5),
    DiscreteUniform(1, 10),
    Binomial(101, 0.45),
    Categorical([0.2, 0.35, 0.15, 0.3])
]

for dist in distset
    μ = mean(dist)
    σ = std(dist)
    # No convenience call.
    E = expectation(dist)
    @test E(x -> x) ≈ μ
    @test E(x -> x^2) - μ^2 ≈ σ^2
    @test E(x -> ((x - μ)/σ)^3) + 1. ≈ skewness(dist) + 1. # To avoid comparisons to 0.0 exactly.
    # Convenience call.
    @test expectation(x -> x, dist) ≈ μ
    @test expectation(x -> x^2, dist) - μ^2 ≈ σ^2
    @test expectation(x -> ((x - μ)/σ)^3, dist) + 1. ≈ skewness(dist) + 1.
    println(dist)
end

# Linear operator behavior.
distset = [
    DiscreteUniform(1., 10.),
    Normal(1.45),
    Beta(2.1)
]

for dist in distset
    E = expectation(dist)
    h(x) = 2*x
    z = nodes(E)
    @test E * h.(z) ≈ E(x -> 2*x) # Right-multiplying.
    @test weights(2E) ≈ 2*weights(E) # Left-multiplying.
    @test 3E*z ≈ (3E) * z ≈ 3*(E * z) # Linearity.
end

# Error handling.
distset = [ # Noncompact dists
    Beta(Inf),
    LogNormal(Inf, Inf),
    Normal(0, Inf),
    Uniform(-Inf, Inf),
    Poisson(3)
]

for dist in distset
    @test_throws ArgumentError expectation(dist)
end

# Other errors.
E = expectation(DiscreteUniform(1, 10))
@test_throws Exception E(x -> dot(x, ones(7))) # Non-applicable functions.
@test_throws MethodError (x -> 2*x).(nodes(E)) * E # Non-commutativity.


text:
include("../basic_abstract_game.jl")

const NAME = "bigfish"

const COMPLETION_BONUS = 10
const POSITIVE_REWARD = 1

const FISH = 2

const FISH_MIN_R = Float32(0.25)
const FISH_MAX_R = Float32(2)

const FISH_QUOTA = 30

abstract type AbstractBigFish <: AbstractBasicAbstractGame end

mutable struct BigFish <: AbstractBigFish
    basic_game::BasicAbstractGame

    fish_eaten::Int
    r_inc::Float32

    function BigFish()

        basic_game = BasicAbstractGame(NAME)
        basic_game.game.timeout = 6000
        basic_game.main_width = 20
        basic_game.main_height = 20

        fish_eaten = 0
        r_inc = 0.0

        new(
            basic_game,
            fish_eaten,
            r_inc,
        )
    end
end

function get_big_fish(abstract_big_fish::AbstractBigFish)

    return abstract_big_fish
end

function get_basic_game(abstract_big_fish::AbstractBigFish)

    return get_big_fish(abstract_big_fish).basic_game
end

function load_background_images(abstract_big_fish::AbstractBigFish)

    big_fish = get_big_fish(abstract_big_fish)

    big_fish.basic_game.main_bg_images_ptr = water_backgrounds
    return nothing
end

function asset_for_type(abstract_big_fish::AbstractBigFish, type::Int, names::Vector{String})

    if type == Int(PLAYER)
        push!(names, "misc_assets/fishTile_072.png")
    elseif type == Int(FISH)
        push!(names, "misc_assets/fishTile_074.png")
        push!(names, "misc_assets/fishTile_078.png")
        push!(names, "misc_assets/fishTile_080.png")
    end
    return nothing
end

function handle_agent_collision(abstract_big_fish::AbstractBigFish, obj::Entity)

    big_fish = get_big_fish(abstract_big_fish)

    if obj.type == FISH
        if obj.rx > big_fish.basic_game.agent.rx
            big_fish.basic_game.game.step_data.done = true
        else
            big_fish.basic_game.game.step_data.reward += POSITIVE_REWARD
            obj.will_erase = true
            big_fish.basic_game.agent.rx += big_fish.r_inc
            big_fish.basic_game.agent.ry += big_fish.r_inc
            big_fish.fish_eaten += 1
        end
    end
    return nothing
end

function game_reset(abstract_big_fish::AbstractBigFish)

    big_fish = get_big_fish(abstract_big_fish)

    basic_game_reset(abstract_big_fish)

    big_fish.basic_game.game.options.center_agent = false
    big_fish.fish_eaten = 0

    start_r::Float32 = 0.5

    if big_fish.basic_game.game.options.distribution_mode == EasyMode
        start_r = 1
    end

    big_fish.r_inc = (FISH_MAX_R - start_r) / FISH_QUOTA

    big_fish.basic_game.agent.rx = start_r
    big_fish.basic_game.agent.ry = start_r
    big_fish.basic_game.agent.y = 1 + big_fish.basic_game.agent.ry
    return nothing
end

function game_step(abstract_big_fish::AbstractBigFish)

    big_fish = get_big_fish(abstract_big_fish)

    basic_game_step(abstract_big_fish)

    if randn(big_fish.basic_game.game.rand_gen, 10) == 1
        ent_r = (FISH_MAX_R - FISH_MIN_R) * rand01(big_fish.basic_game.game.rand_gen) ^ Float32(1.4) + FISH_MIN_R
        ent_y = rand01(big_fish.basic_game.game.rand_gen) * (big_fish.basic_game.main_height - 2 * ent_r)
        moves_right = rand01(big_fish.basic_game.game.rand_gen) < 0.5
        ent_vx = (Float32(0.15) + rand01(big_fish.basic_game.game.rand_gen) * Float32(0.25)) * (moves_right ? 1 : -1)
        ent_x = moves_right ? -1 * ent_r : big_fish.basic_game.main_width + ent_r
        type = FISH
        ent = add_entity(big_fish.basic_game, ent_x, ent_y, ent_vx, Float32(0), ent_r, type)
        choose_random_theme(abstract_big_fish, ent)
        match_aspect_ratio(abstract_big_fish, ent)
        ent.is_reflected = !moves_right
    end

    if big_fish.fish_eaten >= FISH_QUOTA
        big_fish.basic_game.game.step_data.done = true
        big_fish.basic_game.game.step_data.reward += COMPLETION_BONUS
        big_fish.basic_game.game.step_data.level_complete = true
    end

    if big_fish.basic_game.action_vx > 0
        big_fish.basic_game.agent.is_reflected = false
    end
    if big_fish.basic_game.action_vx < 0
        big_fish.basic_game.agent.is_reflected = true
    end

    return nothing
end