text:
"""
RosslerSRI

Holds the Butcher tableaus for a Roessler SRI method.
"""
struct RosslerSRI{T,T2} <: Tableau
  c₀::Vector{T2}
  c₁::Vector{T2}
  A₀::Matrix{T}
  A₁::Matrix{T}
  B₀::Matrix{T}
  B₁::Matrix{T}
  α::Vector{T}
  β₁::Vector{T}
  β₂::Vector{T}
  β₃::Vector{T}
  β₄::Vector{T}
  order::Rational{Int}
end

"""
RosslerSRA

Holds the Butcher tableaus for a Rosser SRA method.
"""
struct RosslerSRA{T,T2} <: Tableau
  c₀::Vector{T2}
  c₁::Vector{T2}
  A₀::Matrix{T}
  B₀::Matrix{T}
  α::Vector{T}
  β₁::Vector{T}
  β₂::Vector{T}
  order::Rational{Int}
end

"""
constructSRIW1()

Constructs the tableau type for the SRIW1 method.
"""
function constructSRIW1(T=Float64,T2=Float64)
  c₀ = [0;3//4;0;0]
  c₁ = [0;1//4;1;1//4]
  A₀ = [0 0 0 0
      3//4 0 0 0
      0 0 0 0
      0 0 0 0]
  A₁ = [0 0 0 0
      1//4 0 0 0
      1 0 0 0
      0 0 1//4 0]
  B₀ = [0 0 0 0
      3//2 0 0 0
      0 0 0 0
      0 0 0 0]
  B₁ = [0 0 0 0
      1//2 0 0 0
      -1 0 0 0
      -5 3 1//2 0]

  α = [1//3;2//3;0;0]

  β₁ = [-1;4//3;2//3;0]
  β₂ = -[1;-4//3;1//3;0]
  β₃ = [2;-4//3;-2//3;0]
  β₄ = [-2;5//3;-2//3;1]
  RosslerSRI(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,A₁),
             map(T,B₀),map(T,B₁),
             map(T,α),map(T,β₁),map(T,β₂),
             map(T,β₃),map(T,β₄),3//2)
end

"""
constructSRIW2()

Constructs the tableau type for the SRIW1 method.
"""
function constructSRIW2(T=Float64,T2=Float64)
  c₀ = [0;1;1//2;0]
  c₁ = [0;1//4;1;1//4]
  A₀ = [0 0 0 0
      1 0 0 0
      1//4 1//4 0 0
      0 0 0 0]
  A₁ = [0 0 0 0
      1//4 0 0 0
      1 0 0 0
      0 0 1//4 0]
  B₀ = [0 0 0 0
      0 0 0 0
      1 1//2 0 0
      0 0 0 0]
  B₁ = [0 0 0 0
      -1//2 0 0 0
      1 0 0 0
      2 -1 1//2 0]

  α = [1//6;1//6;2//3;0]

  β₁ = [-1;4//3;2//3;0]
  β₂ = [1;-4//3;1//3;0]
  β₃ = [2;-4//3;-2//3;0]
  β₄ = [-2;5//3;-2//3;1]
  RosslerSRI(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,A₁),
             map(T,B₀),map(T,B₁),
             map(T,α),map(T,β₁),map(T,β₂),
             map(T,β₃),map(T,β₄),3//2)
end

"""
constructSRIOpt1()

Opti6-12-11-10-01-47
"""
function constructSRIOpt1(T=Float64,T2=Float64)
  A₀ = [0.0 0.0 0.0 0.0; -0.04199224421316468 0.0 0.0 0.0; 2.842612915017106 -2.0527723684000727 0.0 0.0; 4.338237071435815 -2.8895936137439793 2.3017575594644466 0.0]
  A₁ = [0.0 0.0 0.0 0.0; 0.26204282091330466 0.0 0.0 0.0; 0.20903646383505375 -0.1502377115150361 0.0 0.0; 0.05836595312746999 0.6149440396332373 0.08535117634046772 0.0]
  B₀ = [0.0 0.0 0.0 0.0; -0.21641093549612528 0.0 0.0 0.0; 1.5336352863679572 0.26066223492647056 0.0 0.0; -1.0536037558179159 1.7015284721089472 -0.20725685784180017 0.0]
  B₁ = [0.0 0.0 0.0 0.0; -0.5119011827621657 0.0 0.0 0.0; 2.67767339866713 -4.9395031322250995 0.0 0.0; 0.15580956238299215 3.2361551006624674 -1.4223118283355949 0.0]
  α  = [1.140099274172029,-0.6401334255743456,0.4736296532772559,0.026404498125060714]
  β₁ = [-1.8453464565104432,2.688764531100726,-0.2523866501071323,0.40896857551684956]
  β₂ = [0.4969658141589478,-0.5771202869753592,-0.12919702470322217,0.2093514975196336]
  β₃ = [2.8453464565104425,-2.688764531100725,0.2523866501071322,-0.40896857551684945]
  β₄ = [0.11522663875443433,-0.57877086147738,0.2857851028163886,0.17775911990655704]
  e = ones(size(α))
  c₀ = A₀*e
  c₁ = A₁*e
  RosslerSRI(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,A₁),
             map(T,B₀),map(T,B₁),
             map(T,α),map(T,β₁),map(T,β₂),
             map(T,β₃),map(T,β₄),3//2)
end

"""
constructSRIOpt2()

Opti6-12-11-10-01-47
"""
function constructSRIOpt2(T=Float64,T2=Float64)
  c0 = [0.0,0.13804532298278663,0.9999999999999992,0.9999999999999994]
  c1 = [0.0,0.45605532163856893,0.999999999999996,0.9999999999999962]
  A0 = [0.0 0.0 0.0 0.0
        0.13804532298278663 0.0 0.0 0.0
        0.5818361298250374 0.4181638701749618 0.0 0.0
        0.4670018408674211 0.8046204792187386 -0.27162232008616016 0.0]
  A1 = [0.0 0.0 0.0 0.0
        0.45605532163856893 0.0 0.0 0.0
        0.7555807846451692 0.24441921535482677 0.0 0.0
        0.6981181143266059 0.3453277086024727 -0.04344582292908241 0.0]
  B0 = [0.0 0.0 0.0 0.0
        0.08852381537667678 0.0 0.0 0.0
        1.0317752458971061 0.4563552922077882 0.0 0.0
        1.73078280444124 -0.46089678470929774 -0.9637509618944188 0.0]
  B1 = [0.0 0.0 0.0 0.0
        0.6753186815412179 0.0 0.0 0.0
        -0.07452812525785148 -0.49783736486149366 0.0 0.0
        -0.5591906709928903 0.022696571806569924 -0.8984927888368557 0.0]
  α  = [-0.15036858140642623,0.7545275856696072,0.686995463807979,-0.2911544680711602]
  β1 = [-0.45315689727309133,0.8330937231303951,0.3792843195533544,0.24077885458934192]
  β2 = [-0.4994383733810986,0.9181786186154077,-0.25613778661003145,-0.16260245862427797]
  β3 = [1.4531568972730915,-0.8330937231303933,-0.3792843195533583,-0.24077885458934023]
  β4 = [-0.4976090683622265,0.9148155835648892,-1.4102107084476505,0.9930041932449877]
  RosslerSRI(map(T2,c0),map(T2,c1),
             map(T,A0),map(T,A1),
             map(T,B0),map(T,B1),
             map(T,α),map(T,β1),map(T,β2),
             map(T,β3),map(T,β4),3//2)
end

"""
constructSRA1()

Constructs the taleau type for the SRA1 method.
"""
function constructSRA1(T=Float64,T2=Float64)
  α  = [1//3;2//3]
  β₁ = [1;0]
  β₂ = [-1;1]
  A₀ = [0 0
       3//4 0]
  B₀ = [0 0
       3//2 0]
  c₀ = [0;3//4]
  c₁ = [1;0]
  RosslerSRA(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,B₀),
             map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
constructSRA2()

Constructs the taleau type for the SRA2 method.
"""
function constructSRA2(T=Float64,T2=Float64)
  α  = [1//3;2//3]
  β₁ = [0;1]
  β₂ = [3//2;-3//2]
  A₀ = [0 0
       3//4 0]
  B₀ = [0 0
       3//2 0]
  c₀ = [0;3//4]
  c₁ = [1//3;1]
  RosslerSRA(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,B₀),
             map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
constructSRA3()

Constructs the taleau type for the SRA3 method.
"""
function constructSRA3(T=Float64,T2=Float64)
  α  = [1//6;1//6;2//3]
  β₁ = [1;0;0]
  β₂ = [-1;1;0]
  A₀ = [0 0 0
        1 0 0
        1//4 1//4 0]
  B₀ = [0 0 0
        0 0 0
        1 1//2 0]
  c₀ = [0;1;1//2]
  c₁ = [1;0;0]
  RosslerSRA(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,B₀),
             map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
constructSOSRA()

Constructs the taleau type for the SOSRA method.
"""
function constructSOSRA(T=Float64,T2=Float64)
  a1 = 0.2889874966892885;
  a2 = 0.6859880440839937;
  a3 = 0.025024459226717772;
  c01 = 0;
  c02 = 0.6923962376159507;
  c03 = 1;
  c11 = 0;
  c12 = 0.041248171110700504;
  c13 = 1;
  b11 = -16.792534242221663;
  b12 = 17.514995785380226;
  b13 = 0.27753845684143835;
  b21 = 0.4237535769069274;
  b22 = 0.6010381474428539;
  b23 = -1.0247917243497813;
  A021 = 0.6923962376159507;
  A031 = -3.1609142252828395;
  A032 = 4.1609142252828395;
  B021 = 1.3371632704399763;
  B031 = 1.442371048468624;
  B032 = 1.8632741501139225;
  α  = [a1;a2;a3]
  β₁ = [b11;b12;b13]
  β₂ = [b21;b22;b23]
  A₀ = [0 0 0
        A021 0 0
        A031 A032 0]
  B₀ = [0 0 0
        B021 0 0
        B031 B032 0]
  c₀ = [c01;c02;c03]
  c₁ = [c11;c12;c13]
  RosslerSRA(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,B₀),
             map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
constructSOSRA2()

Constructs the taleau type for the SOSRA method.
"""
function constructSOSRA2(T=Float64,T2=Float64)
  a1 = 0.4999999999999998;
  a2 = -0.9683897375354181;
  a3 = 1.4683897375354185;
  c01 = 0;
  c02 = 1;
  c03 = 1;
  c11 = 0;
  c12 = 1.0;
  c13 = 1;
  b11 = 0.0;
  b12 = 0.92438032145683;
  b13 = 0.07561967854316998;
  b21 = 1.0;
  b22 = -0.8169981105823436;
  b23 = -0.18300188941765633;
  A021 = 1;
  A031 = 0.9511849235504364;
  A032 = 0.04881507644956362;
  B021 = 0.7686101171003622;
  B031 = 0.43886792994934987;
  B032 = 0.7490415909204886;
  α  = [a1;a2;a3]
  β₁ = [b11;b12;b13]
  β₂ = [b21;b22;b23]
  A₀ = [0 0 0
        A021 0 0
        A031 A032 0]
  B₀ = [0 0 0
        B021 0 0
        B031 B032 0]
  c₀ = [c01;c02;c03]
  c₁ = [c11;c12;c13]
  RosslerSRA(map(T2,c₀),map(T2,c₁),
             map(T,A₀),map(T,B₀),
             map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
checkSRIOrder(RosslerSRI)

Determines whether the order conditions are met via the tableaus of the SRI method.
"""
function checkSRIOrder(RosslerSRI;tol=1e-6)
  @unpack c₀,c₁,A₀,A₁,B₀,B₁,α,β₁,β₂,β₃,β₄ = RosslerSRI
  e = ones(size(α))
  conditions = Vector{Bool}(undef, 25)
  conditions[1] = abs(dot(α,e)-1)<tol
  conditions[2] = abs(dot(β₁,e)-1)<tol
  conditions[3] = abs(dot(β₂,e)-0)<tol
  conditions[4] = abs(dot(β₃,e)-0)<tol
  conditions[5] = abs(dot(β₄,e)-0)<tol
  conditions[6] = abs(sum(β₁'*B₁)-0)<tol
  conditions[7] = abs(sum(β₂'*B₁)-1)<tol
  conditions[8] = abs(sum(β₃'*B₁)-0)<tol
  conditions[9] = abs(sum(β₄'*B₁)-0)<tol
  conditions[10]= abs(sum(α'*A₀)- .5) < tol
  conditions[11]= abs(sum(α'*B₀)- 1) < tol
  conditions[12]= abs(sum(α'*(B₀*e).^2)- 3/2) < tol
  conditions[13]= abs(sum(β₁'*A₁)-1) < tol
  conditions[14]= abs(sum(β₂'*A₁)-0) < tol
  conditions[15]= abs(sum(β₃'*A₁)+1) < tol
  conditions[16]= abs(sum(β₄'*A₁)-0) < tol
  conditions[17]= abs(sum(β₁'*(B₁*e).^2)-1) < tol
  conditions[18]= abs(sum(β₂'*(B₁*e).^2)-0) < tol
  conditions[19]= abs(sum(β₃'*(B₁*e).^2)+1) < tol
  conditions[20]= abs(sum(β₄'*(B₁*e).^2)-2) < tol
  conditions[22]= abs(sum(β₂'*B₁*(B₁*e))-0) < tol
  conditions[21]= abs(sum(β₁'*B₁*(B₁*e))-0) < tol
  conditions[23]= abs(sum(β₃'*B₁*(B₁*e))-0) < tol
  conditions[24]= abs(sum(β₄'*B₁*(B₁*e))-1) < tol
  conditions[25]= abs.(.5*β₁'*(A₁*(B₀*e)) + (1/3)*β₃'*(A₁*(B₀*e)))[1] .< tol
  return(conditions)
end

"""
checkSRAOrder(RosslerSRI)

Determines whether the order conditions are met via the tableaus of the SRA method.
"""
function checkSRAOrder(SRA;tol=1e-6)
  @unpack c₀,c₁,A₀,B₀,α,β₁,β₂ = SRA
  e = ones(size(α))
  conditions = Vector{Bool}(undef, 8)
  conditions[1] = abs(dot(α,e)-1)<tol
  conditions[2] = abs(dot(β₁,e)-1)<tol
  conditions[3] = abs(dot(β₂,e)-0)<tol
  conditions[4] = abs(dot(α,B₀*e)-1).<tol
  conditions[5] = abs(dot(α,A₀*e)-1/2).<tol
  conditions[6] = abs(dot(α,(B₀*e).^2) - 3/2).<tol
  conditions[7] = abs(dot(β₁,c₁)-1)<tol
  conditions[8] = abs(dot(β₂,c₁)+1)<tol
  return(conditions)
end

"""
constructSOSRA2()

Constructs the taleau type for the SOSRA method.
"""
function constructSKenCarp(T=Float64,T2=Float64)
    γ  = convert(T,0.435866521508459)
    a31 = convert(T,0.2576482460664272)
    a32 = -convert(T,0.09351476757488625)
    a41 = convert(T,0.18764102434672383)
    a42 = -convert(T,0.595297473576955)
    a43 = convert(T,0.9717899277217721)
    # bhat1 = convert(T,2756255671327//12835298489170)
    # bhat2 = -convert(T,10771552573575//22201958757719)
    # bhat3 = convert(T,9247589265047//10645013368117)
    # bhat4 = convert(T,2193209047091//5459859503100)
    btilde1 = convert(T,0.027099261876665316) # bhat1-a41
    btilde2 = convert(T,0.11013520969201586) # bhat2-a42
    btilde3 = convert(T,-0.10306492520138458) # bhat3-a43
    btilde4 = convert(T,-0.0341695463672966) # bhat4-γ
    c3 = convert(T2,0.6)
    c2 = 2γ
    θ = c3/c2
    α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
    α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
    θ = 1/c2
    α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
    α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

    nb021 = convert(T,-12.246764387585056)
    nb043 = convert(T,-14.432096958608753)
    α  = [a41;a42;a43;γ]
    β₁ = [0,0,0,1]
    β₂ = [1,0,0,-1]
    A₀ = [0 0 0 0
          γ γ 0 0
          a31 a32 γ 0
          a41 a42 a43 γ]
    B₀ = [0 0 0 0
          nb021 0 0 0
          0 0 0 0
          0 0 nb043 0]
    c₀ = [0,c2,c3,1]
    c₁ = [0,0,0,1]
    RosslerSRA(map(T2,c₀),map(T2,c₁),
               map(T,A₀),map(T,B₀),
               map(T,α),map(T,β₁),map(T,β₂),4//2)
end

"""
constructSOSRA2()

Constructs the taleau type for the SOSRA method.
"""
function constructExplicitSKenCarp(T=Float64,T2=Float64)
    γ  = convert(T,0.435866521508459)
    a31 = convert(T,0.2576482460664272)
    a32 = -convert(T,0.09351476757488625)
    a41 = convert(T,0.18764102434672383)
    a42 = -convert(T,0.595297473576955)
    a43 = convert(T,0.9717899277217721)
    # bhat1 = convert(T,2756255671327//12835298489170)
    # bhat2 = -convert(T,10771552573575//22201958757719)
    # bhat3 = convert(T,9247589265047//10645013368117)
    # bhat4 = convert(T,2193209047091//5459859503100)
    btilde1 = convert(T,0.027099261876665316) # bhat1-a41
    btilde2 = convert(T,0.11013520969201586) # bhat2-a42
    btilde3 = convert(T,-0.10306492520138458) # bhat3-a43
    btilde4 = convert(T,-0.0341695463672966) # bhat4-γ
    c3 = convert(T2,0.6)
    c2 = 2γ
    θ = c3/c2
    α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
    α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
    θ = 1/c2
    α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
    α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

    ea21 = convert(T,1767732205903//2027836641118)
    ea31 = convert(T,5535828885825//10492691773637)
    ea32 = convert(T,788022342437//10882634858940)
    ea41 = convert(T,6485989280629//16251701735622)
    ea42 = -convert(T,4246266847089//9704473918619)
    ea43 = convert(T,10755448449292//10357097424841)
    eb1 = convert(T,1471266399579//7840856788654)
    eb2 = convert(T,-4482444167858//7529755066697 )
    eb3 = convert(T,11266239266428//11593286722821)
    eb4 = convert(T,1767732205903//4055673282236)

    nb021 = convert(T,-12.246764387585056)
    nb043 = convert(T,-14.432096958608753)

    α  = [eb1;eb2;eb3;eb4]
    β₁ = [0,0,0,1]
    β₂ = [1,0,0,-1]
    A₀ = [0 0 0 0
          ea21 0 0 0
          ea31 ea32 0 0
          ea41 ea42 ea43 0]
    B₀ = [0 0 0 0
          nb021 0 0 0
          0 0 0 0
          0 0 nb043 0]
    c₀ = [0,c2,c3,1]
    c₁ = [0,0,0,1]
    RosslerSRA(map(T2,c₀),map(T2,c₁),
               map(T,A₀),map(T,B₀),
               map(T,α),map(T,β₁),map(T,β₂),4//2)
end

struct SKenCarpTableau{T,T2}
    γ::T
    a31::T
    a32::T
    a41::T
    a42::T
    a43::T
    btilde1::T
    btilde2::T
    btilde3::T
    btilde4::T
    c3::T2
    α31::T
    α32::T
    α41::T
    α42::T
    ea21::T
    ea31::T
    ea32::T
    ea41::T
    ea42::T
    ea43::T
    eb1::T
    eb2::T
    eb3::T
    eb4::T
    ebtilde1::T
    ebtilde2::T
    ebtilde3::T
    ebtilde4 ::T
    nb021::T
    nb043::T
end

#=
# KenCarp3
# Predict z4 from Hermite z2 and z1
# Not z3 because c3 < c2 !

θ = c3/c2
dt = c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
θ = c4/c2
((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
=#
function SKenCarpTableau(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  γ  = convert(T,0.435866521508459)
  a31 = convert(T,0.2576482460664272)
  a32 = -convert(T,0.09351476757488625)
  a41 = convert(T,0.18764102434672383)
  a42 = -convert(T,0.595297473576955)
  a43 = convert(T,0.9717899277217721)
  # bhat1 = convert(T,2756255671327//12835298489170)
  # bhat2 = -convert(T,10771552573575//22201958757719)
  # bhat3 = convert(T,9247589265047//10645013368117)
  # bhat4 = convert(T,2193209047091//5459859503100)
  btilde1 = convert(T,0.027099261876665316) # bhat1-a41
  btilde2 = convert(T,0.11013520969201586) # bhat2-a42
  btilde3 = convert(T,-0.10306492520138458) # bhat3-a43
  btilde4 = convert(T,-0.0341695463672966) # bhat4-γ
  c3 = convert(T2,0.6)
  c2 = 2γ
  θ = c3/c2
  α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
  θ = 1/c2
  α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

  ea21 = convert(T,0.871733043016918)
  ea31 = convert(T,0.5275890119763004)
  ea32 = convert(T,0.0724109880236996)
  ea41 = convert(T,0.3990960076760701)
  ea42 = -convert(T,0.4375576546135194)
  ea43 = convert(T,1.0384616469374492)
  eb1 = convert(T,0.18764102434672383)
  eb2 = convert(T,-0.595297473576955)
  eb3 = convert(T,0.9717899277217721)
  eb4 = convert(T,0.435866521508459)
  ebtilde1 = convert(T,0.027099261876665316)
  ebtilde2 = convert(T,0.11013520969201586)
  ebtilde3 = -convert(T,0.10306492520138458)
  ebtilde4 = -convert(T,0.0341695463672966)

  nb021 = convert(T,-12.246764387585056)
  nb043 = convert(T,-14.432096958608753)
  SKenCarpTableau(γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,
                  α32,α41,α42,ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,
                  ebtilde1,ebtilde2,ebtilde3,ebtilde4,nb021,nb043)
end

function SKenCarpTableau(T,T2)
  γ  = convert(T,1767732205903//4055673282236)
  a31 = convert(T,2746238789719//10658868560708)
  a32 = -convert(T,640167445237//6845629431997)
  a41 = convert(T,1471266399579//7840856788654)
  a42 = -convert(T,4482444167858//7529755066697)
  a43 = convert(T,11266239266428//11593286722821)
  # bhat1 = convert(T,2756255671327//12835298489170)
  # bhat2 = -convert(T,10771552573575//22201958757719)
  # bhat3 = convert(T,9247589265047//10645013368117)
  # bhat4 = convert(T,2193209047091//5459859503100)
  btilde1 = convert(T,BigInt(681815649026867975666107)//BigInt(25159934323302256049469295)) # bhat1-a41
  btilde2 = convert(T,BigInt(18411887981491912264464127)//BigInt(167175311446532472108584143)) # bhat2-a42
  btilde3 = convert(T,BigInt(-12719313754959329011138489)//BigInt(123410692144842870217698057)) # bhat3-a43
  btilde4 = convert(T,BigInt(-47289384293135913063989)//BigInt(1383962894467812063558225)) # bhat4-γ
  c3 = convert(T2,3//5)
  c2 = 2γ
  θ = c3/c2
  α31 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α32 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)
  θ = 1/c2
  α41 = ((1 + (-4θ + 3θ^2)) + (6θ*(1-θ)/c2)*γ)
  α42 = ((-2θ + 3θ^2) + (6θ*(1-θ)/c2)*γ)

  # Explicit Tableau
  ea21 = convert(T,1767732205903//2027836641118)
  ea31 = convert(T,5535828885825//10492691773637)
  ea32 = convert(T,788022342437//10882634858940)
  ea41 = convert(T,6485989280629//16251701735622)
  ea42 = -convert(T,4246266847089//9704473918619)
  ea43 = convert(T,10755448449292//10357097424841)
  eb1 = convert(T,1471266399579//7840856788654)
  eb2 = convert(T,-4482444167858//7529755066697 )
  eb3 = convert(T,11266239266428//11593286722821)
  eb4 = convert(T,1767732205903//4055673282236)
  ebtilde1 = convert(T,BigInt(681815649026867975666107)//BigInt(25159934323302256049469295))
  ebtilde2 = convert(T,BigInt(18411887981491912264464127)//BigInt(167175311446532472108584143))
  ebtilde3 = -convert(T,BigInt(12719313754959329011138489)//BigInt(123410692144842870217698057))
  ebtilde4 = -convert(T,BigInt(47289384293135913063989)//BigInt(1383962894467812063558225))

  # Noise Tableau

  nb021 = convert(T,parse(BigFloat,"-12.246764387585055918338744103409192607986567514699471403397969732723452087723101"))
  nb043 = convert(T,parse(BigFloat,"-14.432096958608752822047165680776748797565142459789556194474191884258734697161106"))
  SKenCarpTableau(γ,a31,a32,a41,a42,a43,btilde1,btilde2,btilde3,btilde4,c3,α31,
                  α32,α41,α42,ea21,ea31,ea32,ea41,ea42,ea43,eb1,eb2,eb3,eb4,
                  ebtilde1,ebtilde2,ebtilde3,ebtilde4,nb021,nb043)
end

# Flip them all!

# ebtilde1 = big(1471266399579)//7840856788654 - big(2756255671327)//12835298489170
# ebtilde2 = -big(4482444167858)//7529755066697 + big(10771552573575)//22201958757719
# ebtilde3 = big(11266239266428)//11593286722821 - big(9247589265047)//10645013368117
# ebtilde4 = big(1767732205903)//4055673282236 - big(2193209047091)//5459859503100

"""

RoesslerRI

Holds the Butcher tableaus for a Roessler RI method. (high weak order)

"""

struct RoesslerRI{T,T2} <: Tableau
  c₀::Vector{T2}
  c₁::Vector{T2}
  c₂::Vector{T2}
  A₀::Matrix{T}
  A₁::Matrix{T}
  A₂::Matrix{T}
  B₀::Matrix{T}
  B₁::Matrix{T}
  B₂::Matrix{T}
  α::Vector{T}
  β₁::Vector{T}
  β₂::Vector{T}
  β₃::Vector{T}
  β₄::Vector{T}
  order::Rational{Int}
  quantile::T
end

function constructDRI1(T=Float64,T2=Float64)
  c₀ = [0;1//2;1]
  c₁ = [0;342//491;342//491]
  c₂ = [0;0;0]
  A₀ = [0 0 0
      1//2 0 0
      -1 2 0]
  A₁ = [0 0 0
      342//491 0 0
      342//491 0 0]
  A₂ = [0 0 0
      0 0 0
      0 0 0]
  B₀ = [0 0 0
       (6-sqrt(6))/10 0 0
       (3+2*sqrt(6))/5 0 0]
  B₁ = [0 0 0
       3*sqrt(38//491) 0 0
       -3*sqrt(38/491) 0 0]
  B₂ = [0 0 0
       -214//513*sqrt(1105//991) -491//513*sqrt(221//4955) -491//513*sqrt(221//4955)
       214//513*sqrt(1105//991) 491//513*sqrt(221//4955) 491//513*sqrt(221//4955)]
  α = [1//6;2//3;1//6]

  β₁ = [193//684;491//1368;491//1368]
  β₂ = [0;1//6*sqrt(491//38);-1//6*sqrt(491//38)]
  β₃ = [-4955//7072;4955//14144;4955//14144]
  β₄ = [0;-1//8*sqrt(4955//221);1//8*sqrt(4955//221)]
  RoesslerRI(map(T2,c₀),map(T2,c₁),map(T2,c₂),
             map(T,A₀),map(T,A₁),map(T,A₂),
             map(T,B₀),map(T,B₁),map(T,B₂),
             map(T,α),map(T,β₁),map(T,β₂),
             map(T,β₃),map(T,β₄),
             1//1,-0.9674215661017014)
end

function constructRI1(T=Float64,T2=Float64)
  c₀ = [0;2//3;2//3]
  c₁ = [0;1;1]
  c₂ = [0;0;0]
  A₀ = [0 0 0
      2//3 0 0
      -1//3 1 0]
  A₁ = [0 0 0
      1 0 0
      1 0 0]
  A₂ = [0 0 0
      0 0 0
      0 0 0]
  B₀ = [0 0 0
       1 0 0
       0 0 0]
  B₁ = [0 0 0
       1 0 0
       -1 0 0]
  B₂ = [0 0 0
       1 0 0
       -1 0 0]
  α = [1//4;1//2;1//4]

  β₁ = [1//2;1//4;1//4]
  β₂ = [0;1//2;-1//2]
  β₃ = [-1//2;1//4;1//4]
  β₄ = [0;1//2;-1//2]
  RoesslerRI(map(T2,c₀),map(T2,c₁),map(T2,c₂),
             map(T,A₀),map(T,A₁),map(T,A₂),
             map(T,B₀),map(T,B₁),map(T,B₂),
             map(T,α),map(T,β₁),map(T,β₂),
             map(T,β₃),map(T,β₄),
             1//1,-0.9674215661017014)
end

function checkRIOrder(RI;tol=1e-6)
  @unpack c₀,c₁,c₂,A₀,A₁,A₂,B₀,B₁,B₂,α,β₁,β₂,β₃,β₄ = RI
  e = ones(size(α))
  conditions = Vector{Bool}(undef, 59) # 9 conditions for first order, 59 in total
  #first order weak sense
  conditions[1] = abs(dot(α,e)-1)<tol
  conditions[2] = abs(dot(β₄,e)-0)<tol
  conditions[3] = abs(dot(β₃,e)-0)<tol
  conditions[4] = abs(dot(β₁,e)^2-1)<tol
  conditions[5] = abs(dot(β₂,e)-0)<tol
  conditions[6] = abs(dot(β₁,(B₁*e))-0)<tol
  conditions[7] = abs(dot(β₄,(A₂*e))-0)<tol
  conditions[8] = abs(dot(β₃,(B₂*e))-0)<tol
  conditions[9] = abs(dot(β₄,(B₂*e).^2)-0)<tol

  #second order weak sense
  conditions[10] = abs(dot(α,(A₀*e))-0.5)<tol
  conditions[11] = abs(dot(α,(B₀*e).^2)-0.5)<tol
  conditions[12] = abs(dot(β₁,e)*dot(α,(B₀*e))-0.5)<tol
  conditions[13] = abs(dot(β₁,e)*dot(β₁,(A₁*e))-0.5)<tol
  conditions[14] = abs(dot(β₃,(A₂*e))-0)<tol
  conditions[15] = abs(dot(β₂,(B₁*e))-1)<tol
  conditions[16] = abs(dot(β₄,(B₂*e))-1)<tol
  conditions[17] = abs(dot(β₁,e)*dot(β₁,(B₁*e).^2)-0.5)<tol
  conditions[18] = abs(dot(β₁,e)*dot(β₃,(B₂*e).^2)-0.5)<tol
  conditions[19] = abs(dot(β₁,B₁*(B₁*e))-0)<tol
  conditions[20] = abs(dot(β₃,B₂*(B₁*e))-0)<tol
  conditions[21] = abs(dot(β₃,B₂*(B₁*(B₁*e)))-0)<tol
  conditions[22] = abs(dot(β₁,A₁*(B₀*e))-0)<tol
  conditions[23] = abs(dot(β₃,A₂*(B₀*e))-0)<tol
  conditions[24] = abs(dot(β₄,(A₂*e).^2)-0)<tol
  conditions[25] = abs(dot(β₄,A₂*(A₀*e))-0)<tol
  conditions[26] = abs(dot(α,B₀*(B₁*e))-0)<tol
  conditions[27] = abs(dot(β₂,A₁*e)-0)<tol
  conditions[28] = abs(dot(β₁,(A₁*e).*(B₁*e))-0)<tol
  conditions[29] = abs(dot(β₃,(A₂*e).*(B₂*e))-0)<tol
  conditions[30] = abs(dot(β₄,A₂*(B₀*e))-0)<tol
  conditions[31] = abs(dot(β₂,A₁*(B₀*e))-0)<tol
  conditions[32] = abs(dot(β₄,((B₂*e).^2).*(A₂*e))-0)<tol
  conditions[33] = abs(dot(β₄,A₂*(B₀*e).^2)-0)<tol
  conditions[34] = abs(dot(β₂,A₁*(B₀*e).^2)-0)<tol
  conditions[35] = abs(dot(β₁,B₁*(A₁*e))-0)<tol
  conditions[36] = abs(dot(β₃,B₂*(A₁*e))-0)<tol
  conditions[37] = abs(dot(β₂,(B₁*e).^2)-0)<tol
  conditions[38] = abs(dot(β₄,B₂*(B₁*e))-0)<tol
  conditions[39] = abs(dot(β₂,B₁*(B₁*e))-0)<tol
  conditions[40] = abs(dot(β₁,(B₁*e).^3)-0)<tol
  conditions[41] = abs(dot(β₃,(B₂*e).^3)-0)<tol
  conditions[42] = abs(dot(β₁,B₁*((B₁*e).^2))-0)<tol
  conditions[43] = abs(dot(β₃,B₂*((B₁*e).^2))-0)<tol
  conditions[44] = abs(dot(β₄,(B₂*e).^4)-0)<tol
  conditions[45] = abs(dot(β₄,(B₂*(B₁*e)).^2)-0)<tol
  conditions[46] = abs(dot(β₄,(B₂*e).*(B₂*(B₁*e)))-0)<tol
  conditions[47] = abs(dot(α,(B₀*e).*(B₀*(B₁*e)))-0)<tol
  conditions[48] = abs(dot(β₁,(A₁*(B₀*e)).*(B₁*e))-0)<tol
  conditions[49] = abs(dot(β₃,(A₂*(B₀*e)).*(B₂*e))-0)<tol
  conditions[50] = abs(dot(β₁,A₁*(B₀*(B₁*e)))-0)<tol
  conditions[51] = abs(dot(β₃,A₂*(B₀*(B₁*e)))-0)<tol
  conditions[52] = abs(dot(β₄,(B₂*(A₁*e)).*(B₂*e))-0)<tol
  conditions[53] = abs(dot(β₁,B₁*(A₁*(B₀*e)))-0)<tol
  conditions[54] = abs(dot(β₃,B₂*(A₁*(B₀*e)))-0)<tol
  conditions[55] = abs(dot(β₁,(B₁*e).*(B₁*(B₁*e)))-0)<tol
  conditions[56] = abs(dot(β₃,(B₂*e).*(B₂*(B₁*e)))-0)<tol
  conditions[57] = abs(dot(β₁,B₁*(B₁*(B₁*e)))-0)<tol
  conditions[58] = abs(dot(β₄,(B₂*e).*(B₂*((B₁*e).^2)))-0)<tol
  conditions[59] = abs(dot(β₄,(B₂*e).*(B₂*(B₁*(B₁*e))))-0)<tol

  return(conditions)
end


text:
using LinearAlgebra

function corr_ss(κ, λ, ρ, μ)
    (t) -> κ / λ * exp(-abs(t) * λ)
end

function corr_xs(κ, λ, ρ, μ)
    function (t)
        if t >= 0
            ρ * κ / λ / (λ + μ) * exp(-λ * t)
        else
            corr_sx(κ, λ, ρ, μ)(-t)
        end
    end
end

function corr_sx(κ, λ, ρ, μ)
    function (t)
        if t >= 0
            a = ρ * κ / λ / (λ - μ)
            b1 = (1 + (λ - μ) / (λ + μ)) * exp(-μ * t)
            b2 = -exp(-λ * t)
            a * (b1 + b2)
        else
            corr_xs(κ, λ, ρ, μ)(-t)
        end
    end
end

function corr_xx(κ, λ, ρ, μ)
    function (t)
        c1 = exp(-μ * abs(t)) -  exp(-λ * abs(t))
        c2 = exp(-μ * abs(t))
        d1 = ρ^2 / (λ^2 - μ^2) * κ / λ
        d2 = (1 + ρ / (λ + μ)) * κ / λ * ρ / μ
        d1 * c1 + d2 * c2
    end
end

struct GaussianSystem{Signal <: MvNormal,Joint <: MvNormal,Likelihood <: MvNormal}
    prior::Signal
    joint::Joint
    likelihood::Likelihood
end

function GaussianSystem(
    κ::Real,
    λ::Real,
    ρ::Real,
    μ::Real;
    duration::Real=1.0,
    delta_t::Real=1.0,
    t::Matrix{<:Real}=time_matrix(duration=duration, delta_t=delta_t)
)
    c_ss = corr_ss(κ, λ, ρ, μ).(t)
    c_sx = corr_sx(κ, λ, ρ, μ).(t)
    c_xs = corr_xs(κ, λ, ρ, μ).(t)
    c_xx = corr_xx(κ, λ, ρ, μ).(t)
    c_z = [c_ss c_xs; c_sx c_xx]

    prior = MvNormal(LinearAlgebra.Symmetric(c_ss))
    joint = MvNormal(LinearAlgebra.Symmetric(c_z))

    regression_coef = c_sx * inv(c_ss)
    p_x_given_s_cov = LinearAlgebra.Symmetric(c_xx - regression_coef * c_xs)
    likelihood = MvNormal(p_x_given_s_cov)

    GaussianSystem(prior, joint, likelihood)
end

function GaussianSystem(;
    mean_s::Real=50,
    tau_s::Real=1.0,
    corr_time_ratio::Real=10,
    mean_x::Real=mean_s,
    tau_x::Real=tau_s / corr_time_ratio,
    duration::Real=1.0,
    delta_t::Real=1.0,
    t::Matrix{<:Real}=time_matrix(duration=duration, delta_t=delta_t)
)
    λ = 1 / tau_s
    κ = mean_s * λ
    μ = 1 / tau_x
    ρ = μ * mean_x / mean_s
    GaussianSystem(κ, λ, ρ, μ, t=t)
end

function time_matrix(n::Integer, delta_t::Real)
    t_range = range(0.0, length=n, step=delta_t)
    transpose(t_range) .- t_range
end

function time_matrix(; duration::Real, delta_t::Real)
    n = round(Int, duration / delta_t, RoundUp)
    time_matrix(n, duration / n)
end

function marginal(system::GaussianSystem, t)
    MvNormal(LinearAlgebra.Symmetric(corr_xx(system).(t)))
end

function joint(system::GaussianSystem, t)
    MvNormal(LinearAlgebra.Symmetric(corr_z(system, t)))
end

function likelihood(system::GaussianSystem, t; signal::AbstractVector)
    c_ss = corr_ss(system).(t)
    c_sx = corr_sx(system).(t)
    c_xs = corr_xs(system).(t)
    c_xx = corr_xx(system).(t)

    regression_coef = c_sx * inv(c_ss)
    p_x_given_s_cov = LinearAlgebra.Symmetric(c_xx - regression_coef * c_xs)

    MvNormal(regression_coef * signal, p_x_given_s_cov)
end

function posterior(system::GaussianSystem, t; response::AbstractVector)
    c_ss = corr_ss(system).(t)
    c_sx = corr_sx(system).(t)
    c_xs = corr_xs(system).(t)
    c_xx = corr_xx(system).(t)

    regression_coef = c_xs * inv(c_xx)
    p_s_given_x_cov = LinearAlgebra.Symmetric(c_ss - regression_coef * c_sx)

    MvNormal(regression_coef * response, p_s_given_x_cov)
end

function log_likelihood(system::GaussianSystem, t; signal::AbstractArray, response::AbstractArray)
    c_ss = corr_ss(system).(t)
    c_sx = corr_sx(system).(t)
    c_xs = corr_xs(system).(t)
    c_xx = corr_xx(system).(t)

    regression_coef = c_sx * inv(c_ss)
    p_x_given_s_cov = LinearAlgebra.Symmetric(c_xx - regression_coef * c_xs)

    scaled_normal = MvNormal(p_x_given_s_cov)

    logpdf(scaled_normal, response .- regression_coef * signal)
end

function log_joint(system::GaussianSystem, t; signal::AbstractArray, response::AbstractArray)
    c_z = corr_z(system, t)
    distr = MvNormal(c_z)
    logpdf(distr, vcat(signal, response))
end

function log_prior(system::GaussianSystem, t; signal::AbstractArray)
    distr = prior(system)
    logpdf(distr, signal)
end

mutable struct GaussianChain{System <: GaussianSystem} <: MarkovChain
    scale::Float64
    θ::Float64
    system::System
end

GaussianChain(system::GaussianSystem; θ=1.0, scale::Real=1.0) = GaussianChain(scale, θ, system)
chain(system::GaussianSystem; θ=1.0, scale::Real=1.0) = GaussianChain(scale, θ, system)

function generate_configuration(system::GaussianSystem)
    rand(system.joint)
end

function sample(old_conf::Vector{Float64}, system::GaussianSystem)
    n_dim = length(old_conf) ÷ 2
    new_conf = copy(old_conf)
    new_conf[1:n_dim] .= rand(system.prior)
    new_conf
end

function propose!(new_state::Vector{Float64}, old_state::Vector{Float64}, kernel::GaussianChain)
    n_dim = length(old_state) ÷ 2
    new_state .= old_state
    normal_distr = Normal(0.0, kernel.scale)
    for i in 1:n_dim
        new_state[i] += rand(normal_distr)
    end
    new_state
end

energy(state::Vector{Float64}, chain::GaussianChain) = energy(state, chain.system, chain.θ)

function energy(state::Array{Float64,1}, system::GaussianSystem, θ=1.0)
    result = 0.0
    n_dim = length(state) ÷ 2
    if θ > 0.0
        result += θ * logpdf(system.joint, state)
    end
    if (1.0 - θ) > 0.0
        result += (1.0 - θ) * logpdf(system.prior, @view state[1:n_dim])
    end
    -result
end

function energy_difference(state::Array{Float64,1}, system::GaussianSystem)
    energy(state, system, 1.0) - energy(state, system, 0.0)
end

ext: jl


text:
# This file was generated, do not modify it.

using RDatasets, DataFrames

boston = dataset("MASS", "Boston");

typeof(boston)

names(boston)

first(boston, 4)

boston.Crim[1:5]

boston[3, 5]

boston[1:5, [:Crim, :Zn]]

boston[1:5, 1:2]

b1 = select(boston, [:Crim, :Zn, :Indus])
first(b1, 2)

b2 = select(boston, Not(:NOx))
first(b2, 2)

select!(b1, Not(:Crim))
first(b1, 2)

using StatsBase
describe(boston, :min, :max, :mean, :median, :std)

foo(x) = sum(abs.(x)) / length(x)
d = describe(boston, :mean, :median, :foo => foo)
first(d, 3)

using Statistics

mat = convert(Matrix, boston)
mat[1:3, 1:3]

boston.Crim_x_Zn = boston.Crim .* boston.Zn;

mao = dataset("gap", "mao")
describe(mao, :nmissing)

std(mao.Age)

std(skipmissing(mao.Age))

iris = dataset("datasets", "iris")
first(iris, 3)

unique(iris.Species)

gdf = groupby(iris, :Species);

subdf_setosa = gdf[1]
describe(subdf_setosa, :min, :mean, :max)

by(iris, :Species, :PetalLength => mean)

by(iris, :Species, MPL = :PetalLength => mean, SPL = :PetalLength => std)

DataFrames.aggregate(iris, :Species, std)


text:
module ReferenceModels

using JLD2
using ..HelperFuncs
using Random
using DocStringExtensions

import JSON
import TurbulenceConvection
tc = pkgdir(TurbulenceConvection)
include(joinpath(tc, "driver", "main.jl"))
include(joinpath(tc, "driver", "generate_namelist.jl"))
export main1d

import NCDatasets
NC = NCDatasets

export ReferenceModel, ReferenceModelBatch
export get_t_start, get_t_end, get_t_start_Σ, get_t_end_Σ, get_z_obs
export y_dir, Σ_dir, scm_dir, num_vars, uuid
export y_nc_file, Σ_nc_file, get_scm_namelist
export data_directory, namelist_directory
export construct_reference_models
export get_minibatch!, reshuffle_on_epoch_end, write_ref_model_batch
export time_shift_reference_model, write_val_ref_model_batch

"""
    ReferenceModel{FT <: Real}

A structure containing information about the 'true' reference model
and the observation map used to compare the parameterized and
reference models.

# Fields

$(TYPEDFIELDS)

# Constructors

    ReferenceModel(
        y_names::Vector{String},
        y_dir::String,
        scm_dir::String,
        case_name::String,
        t_start::Real,
        t_end::Real;
        Σ_dir::Union{String, Nothing} = nothing,
        Σ_t_start::Union{Real, Nothing} = nothing,
        Σ_t_end::Union{Real, Nothing} = nothing,
        n_obs::Union{Integer, Nothing} = nothing,
        namelist_args = nothing,
    )

`ReferenceModel` constructor allowing for any or all of `Σ_dir`, `Σ_t_start`, `Σ_t_end` to be
unspecified, in which case they take their values from `y_dir`, `t_start` and `t_end`, respectively.

    ReferenceModel(
        y_names::Vector{String},
        y_dir::String,
        scm_parent_dir::String,
        scm_suffix::String,
        case_name::String,
        t_start::Real,
        t_end::Real;
        Σ_dir::Union{String, Nothing} = nothing,
        Σ_t_start::Union{Real, Nothing} = nothing,
        Σ_t_end::Union{Real, Nothing} = nothing,
        n_obs::Union{Integer, Nothing} = nothing,
        namelist_args = nothing,
    )

`ReferenceModel` constructor using `scm_parent_dir`, `case_name`, `scm_suffix`, to define `scm_dir`.
"""
Base.@kwdef struct ReferenceModel{FT <: Real}
    "Vector of reference variable names"
    y_names::Vector{String}
    "Directory for reference data to compute `y` mean vector"
    y_dir::String
    "Directory for reference data to compute `Σ` covariance matrix"
    Σ_dir::String
    "Directory for static data related to forward scm model (parameter file & vertical levels)"
    scm_dir::String
    "Name of case"
    case_name::String
    # TODO: Make t_start and t_end vectors for multiple time intervals per reference model.
    "Start time for computing mean statistics over"
    y_t_start::FT
    "End time for computing mean statistics over"
    y_t_end::FT
    "Start time for computing covariance statistics over"
    Σ_t_start::FT
    "End time for computing covariance statistics over"
    Σ_t_end::FT
    "Vector of observed vertical locations"
    z_obs::Vector{FT}
    "TurbulenceConvection namelist"
    namelist::Dict
end  # ReferenceModel struct

function ReferenceModel(
    y_names::Vector{String},
    y_dir::String,
    scm_dir::String,
    case_name::String,
    t_start::Real,
    t_end::Real;
    Σ_dir::Union{String, Nothing} = nothing,
    Σ_t_start::Union{Real, Nothing} = nothing,
    Σ_t_end::Union{Real, Nothing} = nothing,
    n_obs::Union{Integer, Nothing} = nothing,
    namelist_args = nothing,
)
    # Always create new namelist
    namelist = get_scm_namelist(scm_dir, case_name, y_dir = y_dir, namelist_args = namelist_args, overwrite = true)
    z_obs = construct_z_obs(namelist)
    z_obs = !isnothing(n_obs) ? Array(range(z_obs[1], z_obs[end], n_obs)) : z_obs
    FT = eltype(z_obs)

    Σ_dir = isnothing(Σ_dir) ? y_dir : Σ_dir
    Σ_t_start = isnothing(Σ_t_start) ? t_start : Σ_t_start
    Σ_t_end = isnothing(Σ_t_end) ? t_end : Σ_t_end

    ReferenceModel{FT}(
        y_names,
        y_dir,
        Σ_dir,
        scm_dir,
        case_name,
        FT(t_start),
        FT(t_end),
        FT(Σ_t_start),
        FT(Σ_t_end),
        z_obs,
        namelist,
    )
end

function ReferenceModel(
    y_names::Vector{String},
    y_dir::String,
    scm_parent_dir::String,
    scm_suffix::String,
    case_name::String,
    t_start::Real,
    t_end::Real;
    Σ_dir::Union{String, Nothing} = nothing,
    Σ_t_start::Union{Real, Nothing} = nothing,
    Σ_t_end::Union{Real, Nothing} = nothing,
    n_obs::Union{Integer, Nothing} = nothing,
    namelist_args = nothing,
)
    scm_dir = data_directory(scm_parent_dir, case_name, scm_suffix)
    args = (y_names, y_dir, scm_dir, case_name, t_start, t_end)
    return ReferenceModel(
        args...,
        Σ_dir = Σ_dir,
        Σ_t_start = Σ_t_start,
        Σ_t_end = Σ_t_end,
        n_obs = n_obs,
        namelist_args = namelist_args,
    )
end

get_t_start(m::ReferenceModel) = m.y_t_start
get_t_end(m::ReferenceModel) = m.y_t_end
get_t_start_Σ(m::ReferenceModel) = m.Σ_t_start
get_t_end_Σ(m::ReferenceModel) = m.Σ_t_end

"Returns the observed vertical locations for a reference model"
get_z_obs(m::ReferenceModel) = m.z_obs

y_dir(m::ReferenceModel) = m.y_dir
Σ_dir(m::ReferenceModel) = m.Σ_dir
scm_dir(m::ReferenceModel) = m.scm_dir
get_scm_namelist(m::ReferenceModel) = deepcopy(m.namelist)

# TODO: cache filename and move `get_stats_path` call to constructor.
y_nc_file(m::ReferenceModel) = get_stats_path(y_dir(m))
Σ_nc_file(m::ReferenceModel) = get_stats_path(Σ_dir(m))

data_directory(root::S, name::S, suffix::S) where {S <: AbstractString} = joinpath(root, "Output.$name.$suffix")
uuid(m::ReferenceModel) = String(split(scm_dir(m), ".")[end])

namelist_directory(root::String, m::ReferenceModel) = namelist_directory(root, m.case_name)
namelist_directory(root::S, casename::S) where {S <: AbstractString} = joinpath(root, "namelist_$casename.in")

num_vars(m::ReferenceModel) = length(m.y_names)

"""
    get_scm_namelist(
        output_dir::String,
        case_name::String;
        y_dir::Union{String, Nothing} = nothing,
        overwrite::Bool = false,
        namelist_args = nothing,
    )::Dict

Returns a TurbulenceConvection.jl namelist, given the case and a list of namelist arguments.

Inputs:

 - `output_dir`     :: Directory where the namelist will be stored or should be located.
 - `case_name`      :: Name of the TurbulenceConvection.jl case considered.
 - `y_dir`          :: Directory with LES data to drive the SCM with, if `case_name` is `LES_driven_SCM`.
 - `overwrite`      :: Whether to overwrite the namelist file, if it exists.
 - `namelist_args`  :: Vector of non-default arguments to be used in the namelist, defined as a vector of tuples.
Outputs:
 - `namelist`   :: The TurbulenceConvection.jl namelist.
"""
function get_scm_namelist(
    output_dir::String,
    case_name::String;
    y_dir::Union{String, Nothing} = nothing,
    overwrite::Bool = true,
    namelist_args = nothing,
)::Dict
    namelist_path = namelist_directory(output_dir, case_name)
    namelist = if ~isfile(namelist_path) | overwrite
        NameList.default_namelist(case_name, root = output_dir)
    else
        JSON.parsefile(namelist_path)
    end
    namelist["stats_io"]["calibrate_io"] = true
    if !isnothing(namelist_args)
        for namelist_arg in namelist_args
            change_entry!(namelist, namelist_arg)
        end
    end

    namelist["meta"]["uuid"] = String(split(output_dir, ".")[end])
    namelist["output"]["output_root"] = dirname(output_dir)
    # if `LES_driven_SCM` case, provide input LES stats file
    if case_name == "LES_driven_SCM"
        @assert !isnothing(y_dir) "lesfile must be specified in the construction of LES_driven_SCM namelist."
        namelist["meta"]["lesfile"] = get_stats_path(y_dir)
    end

    return namelist
end

"""
    construct_z_obs(namelist::Dict)

Constructs the vector of observed locations given a TurbulenceConvection.jl namelist.
"""
function construct_z_obs(namelist::Dict)
    grid = construct_grid(namelist)
    return vec(grid.zc)
end

"""
    construct_reference_models(kwarg_ld::Dict{Symbol, Vector{T} where T})::Vector{ReferenceModel}

Returns a vector of `ReferenceModel`s given a dictionary of keyword argument lists.

Inputs:

 - `kwarg_ld`     :: Dictionary of keyword argument lists

Outputs:

 - `ref_models`   :: Vector where the i-th ReferenceModel is constructed from the i-th element of every keyword argument list of the dictionary.
"""
function construct_reference_models(kwarg_ld::Dict{Symbol, Vector{T} where T})::Vector{ReferenceModel}
    n_RM = length(kwarg_ld[:case_name])
    ref_models = Vector{ReferenceModel}()
    for RM_i in 1:n_RM
        kw = Dict(k => v[RM_i] for (k, v) in pairs(kwarg_ld))  # unpack dict
        args = if haskey(kw, :scm_dir)
            # construct ref_models using `scm_dir` directly
            (kw[j] for j in (:y_names, :y_dir, :scm_dir, :case_name, :t_start, :t_end))
        elseif all(haskey.(Ref(kw), [:scm_parent_dir, :scm_suffix]))
            # construct ref_models using `scm_parent_dir` and `scm_suffix`
            (kw[j] for j in (:y_names, :y_dir, :scm_parent_dir, :scm_suffix, :case_name, :t_start, :t_end))
        else
            throw(
                ArgumentError(
                    "You need to specify either `scm_dir` or all of " *
                    "(`scm_parent_dir`, `case_name`, `scm_suffix`) to construct a `ReferenceModel`",
                ),
            )
        end
        push!(
            ref_models,
            ReferenceModel(
                args...,
                Σ_dir = get(kw, :Σ_dir, nothing),
                Σ_t_start = get(kw, :Σ_t_start, nothing),
                Σ_t_end = get(kw, :Σ_t_end, nothing),
                n_obs = get(kw, :n_obs, nothing),
                namelist_args = get(kw, :namelist_args, nothing),
            ),
        )
    end
    return ref_models
end

"""
    time_shift_reference_model(m::ReferenceModel, Δt::FT) where {FT <: Real}

Returns a time-shifted ReferenceModel, considering an interval relative to the last
available time step of the original model.

Inputs:

 - `m`     :: A ReferenceModel.
 - `Δt`  :: [LES last time - SCM start time (LES timeframe)]

Outputs:

 - The time-shifted ReferenceModel.
"""
function time_shift_reference_model(m::ReferenceModel, Δt::FT) where {FT <: Real}
    filename = y_nc_file(m)
    t = nc_fetch(filename, "t")
    t_start = t[end] - Δt + get_t_start(m)
    t_end = t[end] - Δt + get_t_end(m)
    Σ_t_start = t[end] - Δt + get_t_start_Σ(m)
    Σ_t_end = t[end] - Δt + get_t_end_Σ(m)

    @assert t_start >= 0 "t_start must be positive after time shift, but $t_start was given."
    @assert Σ_t_start >= 0 "Σ_t_start must be positive after time shift, but $Σ_t_start was given."
    @info string(
        "Shifting time windows for ReferenceModel $(m.case_name)",
        "to ty=($t_start, $t_end), tΣ=($Σ_t_start, $Σ_t_end).",
    )
    return ReferenceModel(
        m.y_names,
        m.y_dir,
        m.Σ_dir,
        m.scm_dir,
        m.case_name,
        t_start,
        t_end,
        Σ_t_start,
        Σ_t_end,
        m.z_obs,
        m.namelist,
    )
end

"""
    struct ReferenceModelBatch

A structure containing a batch of ReferenceModels and an evaluation
order for ReferenceModels within the current epoch.

# Fields

$(TYPEDFIELDS)

# Constructors

    ReferenceModelBatch(ref_models::Vector{ReferenceModel}, shuffling::Bool = true)

`ReferenceModelBatch` constructor given a vector of `ReferenceModel`s.

    ReferenceModelBatch(kwarg_ld::Dict{Symbol, Vector{T} where T}, shuffling::Bool = true)

`ReferenceModelBatch` constructor given a dictionary of keyword argument lists.

Inputs:

 - `kwarg_ld`     :: Dictionary of keyword argument lists
 - `shuffling`    :: Whether to shuffle the order of ReferenceModels.
"""
Base.@kwdef struct ReferenceModelBatch
    "Vector containing all reference models"
    ref_models::Vector{ReferenceModel}
    "Vector of indices defining the `ReferenceModel` evaluation order when batching"
    eval_order::Vector{Int}
end

function ReferenceModelBatch(kwarg_ld::Dict{Symbol, Vector{T} where T}, shuffling::Bool = true)
    ref_models = construct_reference_models(kwarg_ld)
    eval_order = shuffling ? shuffle(1:length(ref_models)) : 1:length(ref_models)
    return ReferenceModelBatch(ref_models, eval_order)
end

function ReferenceModelBatch(ref_models::Vector{ReferenceModel}, shuffling::Bool = true)
    eval_order = shuffling ? shuffle(1:length(ref_models)) : 1:length(ref_models)
    return ReferenceModelBatch(ref_models, eval_order)
end

"""
    get_minibatch!(ref_models::ReferenceModelBatch, batch_size::Int)

Returns a minibatch of `ReferenceModel`s from a ReferenceModelBatch and updates
the eval order.

The size of the minibatch is either the requested size, or the remainder of the
elements in the eval_order for this epoch.

Inputs:

 - `ref_model_batch` :: A ReferenceModelBatch.
 - `batch_size`      :: The number of `ReferenceModel`s to retrieve.

Outputs:

 - A vector of `ReferenceModel`s.
 - The indices of the returned `ReferenceModel`s.
"""
function get_minibatch!(ref_model_batch::ReferenceModelBatch, batch_size::Int)
    batch = min(batch_size, length(ref_model_batch.eval_order))
    indices = [pop!(ref_model_batch.eval_order) for i in 1:batch]
    return ref_model_batch.ref_models[indices], indices
end

"Restarts a shuffled evaluation order if the current epoch has finished."
function reshuffle_on_epoch_end(ref_model_batch::ReferenceModelBatch, shuffling::Bool = true)
    if isempty(ref_model_batch.eval_order)
        @info "Current epoch finished. Reshuffling dataset."
        return ReferenceModelBatch(ref_model_batch.ref_models, shuffling)
    else
        return ref_model_batch
    end
end

function write_ref_model_batch(ref_model_batch::ReferenceModelBatch; outdir_path::String = pwd())
    jldsave(joinpath(outdir_path, "ref_model_batch.jld2"); ref_model_batch)
end
function write_val_ref_model_batch(ref_model_batch::ReferenceModelBatch; outdir_path::String = pwd())
    jldsave(joinpath(outdir_path, "val_ref_model_batch.jld2"); ref_model_batch)
end

end # module


text:
# compound expressions
julia> volume = begin
           len = 10
           breadth = 20
           height = 30
       len * breadth * height
       end
6000

# Alternate method using ;
julia> volume = (length = 10; breadth = 20; height = 30; length * breadth * height)
6000

# conditional evaluation
# FizzBuzz algorithm
julia> for i in 1:30
           if i % 3 == 0 && i % 5 == 0
               println("FizzBuzz")
           elseif i % 3 == 0
               println("Fizz")
           elseif i % 5 == 0
               println("Buzz")
           else
               println(i)
           end
       end
1
2
Fizz
4
Buzz
Fizz
7
8
Fizz
Buzz
11
Fizz
13
14
FizzBuzz
16
17
Fizz
19
Buzz
Fizz
22
23
Fizz
Buzz
26
Fizz
28
29
FizzBuzz

# checking prime number
# first find factors
julia> function factors(num)
           factors = []
           for i in 1:num
               if rem(num, i) == 0
                   push!(factors, i)
               end
           end
           return factors
       end
factors (generic function with 1 method)

# check if the number was prime
julia> function is_prime(num)
           factors_array = factors(num)
           if length(factors_array) == 2
               true
           else
               false
           end
       end
is_prime (generic function with 1 method)

julia> is_prime(10)
false

julia> is_prime(11)
true

# if blocks are leacky
# here parity is being declared inside a if block, but is still accessible to the println function
julia> function f(n)
           if n % 2 == 0
               parity = "even"
           else
               parity = "odd"
           end
           println("the number is $parity")
       end
f (generic function with 1 method)

julia> f(2)
the number is even

julia> f(3)
the number is odd

# if blocks return values!
julia> switch = true
true

julia> y = if switch println("switch is ON") else println("switch is OFF") end
switch is ON

# although
julia> 1 == true
true

# BUT....If block does NOT translate 1/0 to true/false
# julia> if 1 println("true") end
# ERROR: TypeError: non-boolean (Int64) used in boolean context

### ternary operators
julia> name = "julia"
"julia"

julia> isa(name, String) ? "its a string" : "nopes, not a string"
"its a string"

### short - circuit evaluation

# this method is now deprecated
julia> isdigit("10") && println("reached me")
reached me

# newer method
julia> all(isdigit, "10") && println("reached me")
reached me

julia> all(isdigit, "ten") && println("reached me")
false

### Repeated evaluation

# over integers
julia> collection = [1,2,3,4,5,6,7,8,9]
9-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6
 7
 8
 9

julia> while length(collection) > 1
           pop!(collection)
           println(collection)
       end
[1, 2, 3, 4, 5, 6, 7, 8]
[1, 2, 3, 4, 5, 6, 7]
[1, 2, 3, 4, 5, 6]
[1, 2, 3, 4, 5]
[1, 2, 3, 4]
[1, 2, 3]
[1, 2]
[1]

# over strings
julia> statement = "This is a great example!"
"This is a great example!"

julia> for word in split(statement)
           println(word)
       end
This
is
a
great
example!

# slicing
julia> statement = "Yet another awesome example"
"Yet another awesome example"

julia> for word in split(statement)[2:4]
           println(word)
       end
another
awesome
example

julia> for word in split(statement)[2:3:4]
           println(word)
       end
another

# defining range of numbers
# range can be defined using a syntax similar to start:end example,
julia> for i in 1:5
           print("*",i,"\n", "\t"^i)
       end
*1
	*2
		*3
			*4
				*5

# playing around with loops for fun
julia> for i in 1:5
           for j in 1:i
               print("*")
           end
           println("")
       end
*
**
***
****
*****

julia> for i in 5:-1:1
           for j in 1:i
               print("*")
           end
           println("")
       end
*****
****
***
**
*

# reverse a string
julia> my_string = "This is a great thing to do!"
"This is a great thing to do!"

julia> my_string[length(my_string):-1:1]
"!od ot gniht taerg a si sihT"

# break/continue
julia> for letter in word
           println(letter)
           if letter == 'i'
               break
           end
       end
j
u
l
i

julia> for letter in word
           println(letter)
           if letter == 'l'
               continue
           end
       end
j
u
l
i
a

### Exception handling

# exception is a data type in julia
julia> typeof(Exception)
DataType

# using try catch along with throw()
julia> function say_hi(name)
           try
               if typeof(name) != String
                   throw(TypeError(:say_hi, "printing name", String, name))
               else
                   println("hi $name")
               end
           catch e
               println(typeof(e))
           end
       end
say_hi (generic function with 1 method)

# throws exception as `n` is a Char not a String type
julia> say_hi('n')
TypeError

julia> say_hi("rahul")
hi rahul

# error/warn/info
julia> info("This is an information")
INFO: This is an information

julia> warn("This is a warning")
WARNING: This is a warning

julia> error("This is an error")
ERROR: This is an error
Stacktrace:
 [1] error(::String) at ./error.jl:21

# try/catch/finally
julia> try
           exponent("alpha")
       catch e
           println(e)
       finally
           print("Goodbye")
       end
MethodError(exponent, ("alpha",), 0x0000000000005546)
Goodbye

### Tasks in Julia
julia> function add_one()
           a = 1
           # produce() is now deprecated
           produce(a)
           while true
               a += 1
               produce(a)
           end
       end

julia> generator = Task(add_one)
Task (runnable) @0x000000010eecb850

# consume is deprecated
julia> consume(generator)
1

julia> consume(generator)
2

julia> consume(generator)
3

# valid for julia 0.5
julia> function sample_task()
			for i in 1:3
				produce(i^2)
			end
			produce("END")
			end

julia> result = Task(sample_task)
Task(runnable) @0x0000000113bfc6d0

julia> consume(result)
1

julia> consume(result)
4

julia> consume(result)
9

julia> consume(result)
"END"

julia> consume(result)
()

## helpful external links:
# https://docs.julialang.org/en/latest/manual/parallel-computing/#Channels-1
# https://discourse.julialang.org/t/how-to-replace-consume-and-produce-with-channels/5125/9

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
export PCAICA, PCA #, FA

ica() = MultivariateStats.ICA(Array{Float64,1}(),
                              Array{Float64,2}(undef,0,0))
pca() = MultivariateStats.PCA(Array{Float64,1}(),
                              Array{Float64,2}(undef,0,0),
                              Array{Float64,1}(),
                              Float64(NaN),
                              Float64(NaN))
"""
mutable struct PCA <: PopulationModel
Principal Components Analysis
"""
mutable struct PCA <: PopulationModel
    k::Int
    pca::MultivariateStats.PCA
    function PCA(k=0::Int, pca=pca())
        new(k,pca)
    end
end

function fit!(model::PCA, Z::AbstractArray; use_marchenko=false)
    model.pca = fit(MultivariateStats.PCA, Z)
    if use_marchenko
        λmin,λmax = marchenko_thresh(size(Z)...)
        model.k = model.k == 0 ? sum(model.pca.prinvars.>λmax) : model.k
    else
        model.k = length(model.pca.prinvars)
    end
    return model
end

function weights(model::PCA)
    V = projection(model.pca)[:,1:model.k]
    normweightvectors!(V)
    return V
end

"""
mutable struct PCAICA <: PopulationModel
See: Detecting cell assemblies in large neuronal populations. Lopes dos Santos et al 2013.
"""
mutable struct PCAICA <: PopulationModel
    k::Int
    pca::MultivariateStats.PCA
    ica::MultivariateStats.ICA
    seed::AbstractRNG
    function PCAICA(k=0::Int, pca=pca(), ica=ica(), seed=Random.seed!(13))
        new(k, pca, ica, seed)
    end
end

function fit!(model::PCAICA, Z::AbstractArray)
    model.pca = fit(MultivariateStats.PCA, Z)
    λmin,λmax = marchenko_thresh(size(Z)...)
    model.k = model.k == 0 ? sum(model.pca.prinvars.>λmax) : model.k
    Psign = projection(model.pca)[:,1:model.k] # projection of first k pc's
    Zproj = transpose(Psign) * Z # project spikes onto pc subspace
    model.ica = fit(MultivariateStats.ICA, Zproj, model.k,
                    winit=rand(model.seed, size(Zproj,1), model.k), # for reproducibility
                    maxiter=10000, tol=1.0e-6)
end

function weights(model::PCAICA)
    Psign = projection(model.pca)[:,1:model.k]
    V = Psign * model.ica.W
    normweightvectors!(V)
    return V
end


text:
# Read files to set up a simulation
# See http://manual.gromacs.org/documentation/2016/user-guide/file-formats.html

export
    Atomtype,
    Bondtype,
    Angletype,
    Torsiontype,
    readinputs

"""
    Atomtype(mass, charge, σ, ϵ)

Gromacs atom type.
"""
struct Atomtype{T}
    mass::T
    charge::T
    σ::T
    ϵ::T
end

"""
    Bondtype(b0, kb)

Gromacs bond type.
"""
struct Bondtype{T}
    b0::T
    kb::T
end

"""
    Angletype(th0, cth)

Gromacs angle type.
"""
struct Angletype{T}
    th0::T
    cth::T
end

"""
    Torsiontype(f1, f2, f3, f4)

Gromacs torsion type.
"""
struct Torsiontype{T}
    f1::T
    f2::T
    f3::T
    f4::T
end

"""
    readinputs(topology_file, coordinate_file)
    readinputs(T, topology_file, coordinate_file)

Read a Gromacs topology flat file, i.e. all includes collapsed into one file,
and a Gromacs coordinate file.
Returns the atoms, specific interaction lists, general interaction lists,
non-bonded matrix, coordinates and box size.
"""
function readinputs(T::Type,
                    top_file::AbstractString,
                    coord_file::AbstractString)
    # Read forcefield and topology file
    atomtypes = Dict{String, Atomtype{T}}()
    bondtypes = Dict{String, Bondtype{T}}()
    angletypes = Dict{String, Angletype{T}}()
    torsiontypes = Dict{String, Torsiontype{T}}()
    atomnames = Dict{String, String}()

    name = "?"
    atoms = Atom{T}[]
    bonds = HarmonicBond{T}[]
    pairs = Tuple{Int, Int}[]
    angles = HarmonicAngle{T}[]
    possible_torsions = Tuple{Int, Int, Int, Int}[]
    torsions = Torsion{T}[]

    current_field = ""
    for l in eachline(top_file)
        sl = strip(l)
        if length(sl) == 0 || startswith(sl, ';')
            continue
        end
        if startswith(sl, '[') && endswith(sl, ']')
            current_field = strip(sl[2:end-1])
            continue
        end
        c = split(rstrip(first(split(sl, ";", limit=2))), r"\s+")
        if current_field == "bondtypes"
            bondtype = Bondtype(parse(T, c[4]), parse(T, c[5]))
            bondtypes["$(c[1])/$(c[2])"] = bondtype
            bondtypes["$(c[2])/$(c[1])"] = bondtype
        elseif current_field == "angletypes"
            # Convert th0 to radians
            angletype = Angletype(deg2rad(parse(T, c[5])), parse(T, c[6]))
            angletypes["$(c[1])/$(c[2])/$(c[3])"] = angletype
            angletypes["$(c[3])/$(c[2])/$(c[1])"] = angletype
        elseif current_field == "dihedraltypes" && c[1] != "#define"
            # Convert back to OPLS types
            f4 = parse(T, c[10]) / -4
            f3 = parse(T, c[9]) / -2
            f2 = 4 * f4 - parse(T, c[8])
            f1 = 3 * f3 - 2 * parse(T, c[7])
            torsiontypes["$(c[1])/$(c[2])/$(c[3])/$(c[4])"] = Torsiontype(f1, f2, f3, f4)
        elseif current_field == "atomtypes" && length(c) >= 8
            atomname = uppercase(c[2])
            atomnames[c[1]] = atomname
            # Take the first version of each atom type only
            if !haskey(atomtypes, atomname)
                atomtypes[atomname] = Atomtype(parse(T, c[4]), parse(T, c[5]),
                        parse(T, c[7]), parse(T, c[8]))
            end
        elseif current_field == "atoms"
            attype = atomnames[c[2]]
            push!(atoms, Atom(attype=attype, name=c[5], resnum=parse(Int, c[3]),
                resname=c[4], charge=parse(T, c[7]),
                mass=parse(T, c[8]), σ=atomtypes[attype].σ,
                ϵ=atomtypes[attype].ϵ))
        elseif current_field == "bonds"
            i, j = parse.(Int, c[1:2])
            bondtype = bondtypes["$(atoms[i].attype)/$(atoms[j].attype)"]
            push!(bonds, HarmonicBond(i, j, bondtype.b0, bondtype.kb))
        elseif current_field == "pairs"
            push!(pairs, (parse(Int, c[1]), parse(Int, c[2])))
        elseif current_field == "angles"
            i, j, k = parse.(Int, c[1:3])
            angletype = angletypes["$(atoms[i].attype)/$(atoms[j].attype)/$(atoms[k].attype)"]
            push!(angles, HarmonicAngle(i, j, k, angletype.th0, angletype.cth))
        elseif current_field == "dihedrals"
            i, j, k, l = parse.(Int, c[1:4])
            push!(possible_torsions, (i, j, k, l))
        elseif current_field == "system"
            name = rstrip(first(split(sl, ";", limit=2)))
        end
    end

    # Add torsions based on wildcard torsion types
    for inds in possible_torsions
        at_types = [atoms[x].attype for x in inds]
        desired_key = join(at_types, "/")
        if haskey(torsiontypes, desired_key)
            d = torsiontypes[desired_key]
            push!(torsions, Torsion(inds..., d.f1, d.f2, d.f3, d.f4))
        else
            best_score = 0
            best_key = ""
            for k in keys(torsiontypes)
                c = split(k, "/")
                for a in (c, reverse(c))
                    valid = true
                    score = 0
                    for (i, v) in enumerate(a)
                        if v == at_types[i]
                            score += 1
                        elseif v != "X"
                            valid = false
                            break
                        end
                    end
                    if valid && (score > best_score)
                        best_score = score
                        best_key = k
                    end
                end
            end
            # If a wildcard match is found, add a new specific torsion type
            if best_key != ""
                d = torsiontypes[best_key]
                push!(torsions, Torsion(inds..., d.f1, d.f2, d.f3, d.f4))
            end
        end
    end

    # Read coordinate file and add solvent atoms
    lines = readlines(coord_file)
    coords = SArray{Tuple{3}, T, 1, 3}[]
    for (i, l) in enumerate(lines[3:end-1])
        push!(coords, SVector(
            parse(T, l[21:28]),
            parse(T, l[29:36]),
            parse(T, l[37:44])
        ))

        # Some atoms are not specified explicitly in the topology so are added here
        if i > length(atoms)
            atname = strip(l[11:15])
            attype = replace(atname, r"\d+" => "")
            temp_charge = atomtypes[attype].charge
            if attype == "CL" # Temp hack to fix charges
                temp_charge = T(-1.0)
            end
            push!(atoms, Atom(attype=attype, name=atname,
                resnum=parse(Int, l[1:5]), resname=strip(l[6:10]),
                charge=temp_charge, mass=atomtypes[attype].mass,
                σ=atomtypes[attype].σ, ϵ=atomtypes[attype].ϵ))

            # Add O-H bonds and H-O-H angle in water
            if atname == "OW"
                bondtype = bondtypes["OW/HW"]
                push!(bonds, HarmonicBond(i, i + 1, bondtype.b0, bondtype.kb))
                push!(bonds, HarmonicBond(i, i + 2, bondtype.b0, bondtype.kb))
                angletype = angletypes["HW/OW/HW"]
                push!(angles, HarmonicAngle(i + 1, i, i + 2, angletype.th0, angletype.cth))
            end
        end
    end

    # Calculate matrix of pairs eligible for non-bonded interactions
    n_atoms = length(coords)
    nb_matrix = trues(n_atoms, n_atoms)
    for i in 1:n_atoms
        nb_matrix[i, i] = false
    end
    for b in bonds
        nb_matrix[b.i, b.j] = false
        nb_matrix[b.j, b.i] = false
    end
    for a in angles
        # Assume bonding is already specified
        nb_matrix[a.i, a.k] = false
        nb_matrix[a.k, a.i] = false
    end

    # Calculate matrix of pairs eligible for halved non-bonded interactions
    # This applies to specified pairs in the topology file, usually 1-4 bonded
    #for (i, j) in pairs
    #    nb_matrix[i, j] = T(0.5)
    #    nb_matrix[j, i] = T(0.5)
    #end

    lj = LennardJones(true)
    coulomb = Coulomb(true)

    # Bounding box for PBCs - box goes 0 to this value in 3 dimensions
    box_size = parse(T, first(split(strip(lines[end]), r"\s+")))

    specific_inter_lists = (bonds, angles, torsions)
    general_inters = (lj, coulomb)

    return atoms, specific_inter_lists, general_inters,
            nb_matrix, coords, box_size
end

function readinputs(top_file::AbstractString, coord_file::AbstractString)
    return readinputs(DefaultFloat, top_file, coord_file)
end


text:
using Test
using MyPkgTemplates

@testset "TestsInSrc" begin
    @test isfile(TestsInSrc.TEMPLATE_PATH)
end


text:
struct Serial <: ExecutionPolicy
    fval::Float64
    x::Vector{Float64}

    function (::Type{Serial})()
        return new(0., Vector{Float64}())
    end
end

function initialize!(serial::Serial, x₀::AbstractVector)
    resize!(serial.x,length(x₀))
    serial.x .= x₀
end
getx(serial::Serial) = serial.x

function initialize!(proxgrad::ProxGradient{Serial})
    init_c = @cfunction(init_wrapper, Nothing,
                        (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    boost = POLO.boosting(proxgrad)
    boosting_c = @cfunction(boost_wrapper, Ptr{Cdouble},
                              (Cint, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    step = POLO.stepsize(proxgrad)
    step_c = @cfunction(step_wrapper, Cdouble,
                          (Cint, Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    smooth = POLO.smoothing(proxgrad)
    smoothing_c = @cfunction(smooth_wrapper, Ptr{Cdouble},
                               (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    proxim = POLO.prox(proxgrad)
    prox_c = @cfunction(prox_wrapper, Ptr{Cdouble},
                          (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cvoid}))
    proxgrad.ptr = ccall(POLO.proxgradient_s, Ptr{Cvoid},
                         (Ptr{Cvoid},
                          Ptr{Cvoid}, Any,
                          Ptr{Cvoid}, Any,
                          Ptr{Cvoid}, Any,
                          Ptr{Cvoid}, Any),
                         init_c,
                         boosting_c, boost,
                         step_c, step,
                         smoothing_c, smooth,
                         prox_c, proxim)
    return proxgrad
end

execution_handle(::Serial) = POLO.run_serial
delete_handle(::Serial) = POLO.delete_proxgradient_s
getf_handle(::Serial) = POLO.getf_s
getx_handle(::Serial) = POLO.getx_s


text:
# Barycentric Dynamic Time Julian date of epoch J2000.0
const T0 = 2451545.0
# Speed of light in meters/second is a defining physical constant.
const C = 299792458.0
# Light-time for one astronomical unit (AU) in seconds, from DE-405.
const AU_SEC = 499.0047838061
# Speed of light in AU/day.  Value is 86400 / AU_SEC.
const C_AUDAY = 86400 / AU_SEC
# Astronomical unit in meters.  Value is AU_SEC * C.
const AU = AU_SEC * C
# Astronomical Unit in kilometers.
const AU_KM = 1.4959787069098932e+8
# Heliocentric gravitational constant in meters^3 / second^2, from
# DE-405.
const GS = 1.32712440017987e+20
# Geocentric gravitational constant in meters^3 / second^2, from
# DE-405.
const GE = 3.98600433e+14
# Radius of Earth in meters from IERS Conventions (2003).
const ERAD = 6378136.6
# Earth ellipsoid flattening from IERS Conventions (2003).
const F = 1 / 298.25642
# Rotational angular velocity of Earth in radians/sec from IERS
# Conventions (2003).
const ANGVEL = 7.2921150e-5
# Reciprocal masses of solar system bodies, from DE-405
# (Sun mass / body mass).
# MASS[0] = Earth/Moon barycenter, MASS[1] = Mercury, ...,
# MASS[9] = Pluto, MASS[10] = Sun, MASS[11] = Moon.
const RMASS = [328900.561400, 6023600.0, 408523.71, 332946.050895, 3098708.0, 1047.3486,
               3497.898, 22902.98, 19412.24, 135200000.0, 1.0, 27068700.387534]
# Number of arcseconds in 360 degrees.
const ASEC360 = 1296000.0
# Angle conversion constants
const ASEC2RAD = 4.848136811095359935899141e-6

# "Computed" constants from novas.c
const xi0 = -0.0166170 * ASEC2RAD
const eta0 = -0.0068192 * ASEC2RAD
const da0 = -0.01460 * ASEC2RAD

const frame_tie_rot = @SMatrix [1-0.5(da0^2 + xi0^2) -da0 xi0
                                da0 1-0.5(da0^2 + eta0^2) eta0
                                -xi0 -eta0 1-0.5(eta0^2 + xi0^2)]


text:
using Plots

# Simulation Input: Initial Condition / Design parameters
dim     = 3

p_M_0   = zeros(dim)

V_M_0   = 300
γ_M_0 	= deg2rad(45)
χ_M_0 	= deg2rad(90)

p_T_0   = [10E3; 0; 0]
v_T_0   = zeros(dim)

γ_f_d   = deg2rad(-90)
χ_f_d 	= deg2rad(180)

if dim == 2
    γ_M_0   = 0
    γ_f_d   = 0
end

v_M_0   = V_M_0*[cos(γ_M_0)*sin(χ_M_0); cos(γ_M_0)*cos(χ_M_0); sin(γ_M_0)]
v̂_f_d   = [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)]

if dim == 2
    (p_M_0, v_M_0, p_T_0, v_T_0) = (p_M_0, v_M_0, p_T_0, v_T_0) .|> x -> x[1:2]
end

σ_M_lim = deg2rad(90)
A_M_max = 300

N       = 3
s_Bias  = ComponentArray(α = 1, δ = 0.01, n = 1, r_ref = 10E3, k = 3, m = 0, k̂_d = [1; 0; 0], i_Ω_μ = 0, i_σ_M_lim = 1)

s_BPNG  = BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_zero, s_Bias)
# Bias options: Bias_zero, Bias_IACG_StationaryTarget

function sim_plot(p_M_0::Vector, v_M_0::Vector, p_T_0::Vector, v_T_0::Vector, s_BPNG::BPNG, A_M_max::Number)
    # Execute simulation
    df = main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max)

    ts = df.time
    p_Ms = df.sol |> Map(datum -> datum.pursuer.p) |> collect
    p_Ts = df.sol |> Map(datum -> datum.evador.p)  |> collect

    p_Ms = hcat(p_Ms...)'
    p_Ts = hcat(p_Ts...)'

    # Plotting
    legend_string = ["Missile" "Target"]

    if dim == 2
        f = fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [], legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0)

    elseif dim == 3
        f_2D= fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [], legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0)

        # f_3D = plot([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], label = legend_string,
        #         xlabel = "x [m]", ylabel = "y [m]", zlabel = "z [m]", aspect_ratio = :equal,
        #         camera =(45,30)
        #         )
        # f_3D = equal_AR_3D([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], [], f_3D)

        f_3D = equal_AR_3D([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], legend_string)
        plot!(f_3D, xlabel = "x [m]", ylabel = "y [m]", zlabel = "z [m]", camera = (45,30))

        f = plot(f_2D, f_3D, layout = (2,1), size = (600, 1200))
    end

    f = fig_print([], [], "Trajectory", [], [], [], f; ar_val = :equal) 	# To save the result
    display(f)

    return df, f
end


text:
module LearningHorse

using LinearAlgebra
using Statistics

export Regression, Classification, Preprocessing, LossFunction, Tree, NeuralNetwork

include("utils.jl")
include("./Regression/Regression.jl")
include("./Classification/Classification.jl")
include("./Tree/Tree.jl")
include("./NeuralNetwork/NN.jl")
include("./Preprocessing/Preprocessing.jl")
include("./LossFunction/LossFunctions.jl")

end


text:
function get_pull_requests(
    api::GitHub.GitHubAPI,
    repo::GitHub.Repo,
    state::String;
    per_page::Integer=100,
    page_limit::Integer=100,
)
    repo_name = repo.name
    repo_owner = repo.owner.login
    paginated_prs = GitForge.paginate(
        GitHub.get_pull_requests,
        api,
        repo_owner,
        repo_name;
        state=state,
        per_page=per_page,
        page_limit=page_limit,
    )

    return @mock unique(paginated_prs)
end

function get_pull_requests(
    api::GitLab.GitLabAPI,
    project::GitLab.Project,
    state::String;
    per_page::Integer=100,
    page::Integer=1,
)
    project_id = project.id
    paginated_prs = GitForge.paginate(
        GitLab.get_pull_requests, api, project_id; state=state, per_page=per_page, page=page
    )

    return @mock unique(paginated_prs)
end

function not_pr_fork(repo::GitHub.Repo, pr::GitHub.PullRequest)
    head_repo = pr.head.repo

    # [TODO] [FIXME] figure out why this can sometimes be `nothing`
    # For now, we'll assume that this means the PR was made from a fork, but that fork
    # has been deleted.
    head_repo === nothing && return false # "true" means "not a fork", so "false" means "yes a fork"

    return !head_repo.fork
end
not_pr_fork(repo::GitLab.Project, pr::GitLab.MergeRequest) = repo.id == pr.project_id

my_prs(username::String, pr::GitHub.PullRequest) = lower(pr.user.login) == lower(username)
function my_prs(username::String, pr::GitLab.MergeRequest)
    return lower(pr.author.username) == lower(username)
end

function get_pr_titles(
    forge::Forge, repo::Union{GitHub.Repo,GitLab.Project}, username::AbstractString
)
    state = repo isa GitHub.Repo ? "open" : "opened"
    all_open_prs = @mock get_pull_requests(forge, repo, state)

    return [
        convert(String, strip(pr.title)) for
        (pr, _) in all_open_prs if not_pr_fork(repo, pr) && my_prs(username, pr)
    ]
end


text:
#train.jl
using Flux
using Flux: @epochs
import Flux.Tracker: Params, gradient, data, update!, grad, back!
using ArgParse

include("util2.jl")
include("models.jl")

include("args.jl")

# Training Settings
FLAGS = parse_commandline()

# Load Data
adj, features, labels, idx_train, idx_val, idx_test = jlload_data("cora")

# process data
train_data = zip(features, labels)

# Model + optimizer
model = GCNModel(size(features, 2), FLAGS["hidden"], size(labels, 2), FLAGS["dropout"], adj)
optimizer = ADAMW(FLAGS["lr"], (0.9, 0.999), FLAGS["weight_decay"])

# Loss
function train_loss(x, y)
    # @show(size(x))
    output = model(x)
    return Flux.crossentropy(output[idx_train, :], y[idx_train, :])
end

function val_loss(x, y)
    output = model(x)
    return Flux.crossentropy(output[idx_val, :], y[idx_val, :])
end

train_accuracy(x, y) = jlaccuracy(model(x)[idx_train, :], y[idx_train, :])
val_accuracy(x, y) = jlaccuracy(model(x)[idx_val, :], y[idx_val, :])

ps = Flux.params(model)

function train(train_loss, val_loss, ps, data, labels, optim)
    t0 = time()
    ps = Params(ps)
    @show(train_loss(data, labels))
    @show(val_loss(data, labels))

    # @progress for i in 1:size(data,1)
        try
            l = train_loss(data, labels)
            @interrupts back!(l)
            _update_params!(optim, ps)
        catch ex
            rethrow(ex)
        end
    # end

    @show(time() - t0)
end

function test()
    testmode!(0, true)

end

# train model
# testmode!(0, false) #into training mode
@epochs(FLAGS["epochs"], train(train_loss, val_loss, ps, features, labels, optimizer))

max_stars_repo_name: mcbaron/fluxGCN

max_issues_repo_name: mcbaron/fluxGCN

max_forks_repo_name: mcbaron/fluxGCN


text:
function voigt(x, σ, γ)
    faddeeva(z) = erfcx(-im * z)
    z = (x + im * γ) / σ / sqrt(2)
    return real(faddeeva(z)) / (σ * sqrt(2pi))
end

# J.J.Olivero and R.L. Longbothum in Empirical fits to the Voigt line width: A brief review, JQSRT 17, P233, 1977
# http://snst-hu.lzu.edu.cn/zhangyi/ndata/Voigt_profile.html
function voigt_fwhm(σ, γ)
    @assert σ > 0
    @assert γ > 0
    fg = σ * 2.355
    fl = γ * 2
    return fl * 0.5346 + sqrt(0.2166 * fl^2 + fg^2)
    #=
    Introducing dampening parameter:
    a = fl / fg
    fwhm = fl * (0.5346 + sqrt(0.2166 + (1/a)^2))
    =#
end

function voigt_σγ(fwhm, log_a)
    a = 10. ^log_a
    fl = fwhm / (0.5346 + sqrt(0.2166 + (1/a)^2))
    fg = fl / a
    return (fg / 2.355, fl / 2)
end
# log_a = log10((γ * 2) / (σ * 2.355));  QSFit.voigt_σγ(QSFit.voigt_fwhm(σ, γ), log_a), (σ, γ)

function voigt_γ(fwhm, σ)
    A = 0.5346
    B = 0.2166
    C = (σ * 2.355)^2
    # fwhm = A * fl + sqrt(B * fl^2 + C)
    fl = A * fwhm / (A^2 - B) - sqrt((A^2 * C - B * C + B * fwhm^2) / ((A^2 - B)^2))
    γ = fl / 2
    @assert (fwhm - voigt_fwhm(σ, γ)) / fwhm < 1e-5
    return γ
end
# QSFit.voigt_γ(QSFit.voigt_fwhm(σ, γ), σ), γ


text:
function CentralGravity(Body::T, PosRsph::Vector{Float64}) where {T <: abstractCelestialBody}
  #This function computes the gravitational acceleration according to a Central
  #field model. Gravitational acceleration written in V-frame
  return g = [0, 0, mu(Body)/PosRsph[1]^2]
end


text:
module SkronTest

using Test
using DiffPrecTest: svec, skron
using SparseArrays, Random, LinearAlgebra

function Q(p::Int)
  _c = 1. / sqrt(2.)
  rp = div((p+1)*p,2)
  Q = spzeros( rp,  p*p )
  ind = 0
  for ci=1:p
    for ri=ci:p
      ind += 1
      if ri == ci
        Q[ind, (ci-1)*p + ri] = 1.
      else
        Q[ind, (ci-1)*p + ri] = _c
        Q[ind, (ri-1)*p + ci] = _c
      end
    end
  end
  Q
end

function Qt(p::Int)
  _c = 1. / sqrt(2.)
  rp = div((p+1)*p,2)
  Q = spzeros( p*p, rp )
  ind = 0
  for ci=1:p
    for ri=ci:p
      ind += 1
      if ri == ci
        Q[(ci-1)*p + ri, ind] = 1.
      else
        Q[(ci-1)*p + ri, ind] = _c
        Q[(ri-1)*p + ci, ind] = _c
      end
    end
  end
  Q
end

@testset "skron" begin
  p = 5
  x = randn(p)
  y = randn(p)

  num_row_reduced = div(p*(p+1),2)
  o = zeros( num_row_reduced )
  for i=1:num_row_reduced
      o[i] = svec(x, y, i)
  end

  @test o ≈ Q(p) * kron(x, y)

  o2 = zeros( num_row_reduced, num_row_reduced )
  for ri=1:num_row_reduced
      for ci=1:num_row_reduced
          o2[ri, ci] = skron(x, y, ri, ci)
      end
  end

  @test o2 ≈ Q(p) * ( kron(x*x', y*y') + kron(y*y', x*x') ) * Qt(p) / 2.

  @test skron(x, y, 1:4, 3:6) ≈ (Q(p) * ( kron(x*x', y*y') + kron(y*y', x*x') ) * Qt(p) / 2.)[1:4, 3:6]

  # Q_Sx_skron_Sy_Q
  X = randn(100, 5)
  Y = randn(100, 5)
  Sx = Symmetric( X'*X / 100 )
  Sy = Symmetric( Y'*Y / 100 )

  @test skron(Sx, Sy) ≈ Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.

  I = 1:num_row_reduced
  @test skron(Sx, Sy, I) ≈ Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.

  I = 3:10
  J = 4:7
  @test skron(Sx, Sy, I, J) ≈ (Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.)[I, J]

  I = [2, 5, 7, 8]
  @test skron(Sx, Sy, I) ≈ (Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.)[I, I]

  I = 10:15
  @test skron(Sx, Sy, I) ≈ (Q(p) * ( kron(Sx, Sy) + kron(Sy, Sx) ) * Qt(p) / 2.)[I, I]

  # Q Sx skron yy' Qt
  yy = view(Y, 1, :) * view(Y, 1, :)'
  @test skron(Sx, view(Y, 1, :), I, I) ≈ (Q(p) * ( kron(Sx, yy) + kron(yy, Sx) ) * Qt(p) / 2.)[I, I]
end

@testset "svec" begin
  p = 5
  X = randn(100, p)
  Y = randn(100, p)
  Sx = Symmetric( X'*X / 100 )
  Sy = Symmetric( Y'*Y / 100 )

  # Q_Sx_minus_Sy
  ox = svec(Sx)
  oy = svec(Sy)
  @test ox - oy ≈ Q(p) * vec( Sx - Sy )

  I = 10:15
  @test svec(Sx, I) ≈ (Q(p) * vec( Sx ))[I]

  I = [2, 5, 7, 8]
  @test svec(Sy, I) ≈ (Q(p) * vec( Sy ))[I]

  @test svec(view(X, 1, :), view(Y, 1, :), I) ≈ (Q(p) * kron(view(X, 1, :), view(Y, 1, :)))[I]
end

end


text:
using Pkg
using Base.Iterators
using BSON
using CSV
# using CuArrays
using DataFrames
using DataStructures: OrderedDict
using Dates
using FileIO
using Flux
# using JLD2
# using JLSO
using JSON
using LinearAlgebra
import MLDataUtils: shuffleobs, getobs, splitobs, stratifiedobs, RandomBatches, batchview
import MLDataPattern
using Printf
# using PyPlot
using Random
using Statistics

using LogGoodies
using MLGoodies
using FluxGoodies

# CuArrays.allowscalar(false)

function import_dataset(fname, split::AbstractFloat = 0.8)
    datafname = "$(fname).ser.gz"
    if isfile(datafname)
        @info "importing raw binary dataset: $(datafname)"
        # rawdata = JLSO.load(datafname)["data"]
        rawdata = deserialize_gzip(datafname)
    else
        df = CSV.read("$(fname).csv")
        n = size(df, 1)
        npassengers = (size(df, 2) - 4) ÷ 6
        @info "samples: $(n), #passengers: $npassengers"
        df = df[:, Symbol.(vcat(["onboard_count", "lat", "lon"], [["pickup_$(i)_lat",
                "pickup_$(i)_lon", "pickup_$(i)_maxtime",
                "dropoff_$(i)_lat",
                "dropoff_$(i)_lon", "dropoff_$(i)_maxtime"] for i in 1:npassengers]..., ["feasible"]))]
        desc = "VGA data for $npassengers passengers"
        Random.seed!(1234)
        @info "splitting"
        itrn, itst = first.(parentindices.(splitobs(shuffleobs(df), split))) #TODO there is method for this in  Dataset!
        @info "building raw dataset"
        rawdata = Dataset(df, desc, [:train => itrn, :test => itst])
        @info "saving raw binary dataset"
        # JLSO.save(datafname, "data" => rawdata; format=:julia_serialize)
        serialize_gzip(datafname, rawdata)

    end
    @info "#true: $(sum(rawdata[][!, :feasible])), #false: $(nobs(rawdata[]) - sum(rawdata[][!, :feasible]))"
    dtrans = tonumerical(rawdata, :, :T => Dict(:standardize => false))
    @info "transforming"
    data = MLGoodies.transform(rawdata, dtrans, Matrix{Float32})
    @info "storing the transform to JSON"
    open(f->JSON.print(f, dtrans, 3), "$(fname).transform.json", "w")
    data, dtrans
end

function train(Xtrn, Ttrn, Xtst, Ttst;
    model_dir = ".",
    epochs = 500,
    early_stop = 20,
    hidden = [1024, 512, 256],
    logepochs=1
    )
    mkpath(model_dir)

    Random.seed!()
    model = MLP(size(Xtrn, 1), hidden..., 1; outactivation=identity) |> gpu

    ps = params(model)
    lossforoutput(Y, T) = batchlogitbinarycrossentropy(Y, T)
    loss(X, T) = lossforoutput(model(getobs(X)), getobs(T))

    best_loss = Inf
    best_epoch = 0

    Random.seed!(1)
    total_start_time = time()

    function evaluatebatch(X, T, batchsize=4096; gpuloss=true)
        Y = similar(T)
        for idxs in partition(1:size(X, 2), batchsize)
            Y[:,idxs] .= model(X[:,idxs])
        end
        loss_ = lossforoutput(Y, T)
        Ycls = prob2class(round.(σ.(Y)))
        Tcls = prob2class(T)

        # serialize_gzip(joinpath(model_dir, "eval.ser.gz"), (Ycls, Tcls))
        acc_ = 100 * accuracy(Ycls, Tcls)
        cm_ = cmatrix(Ycls, Tcls)
        loss_, acc_, cm_
    end

    opt = ADAM()
    # bsize = 4096
    bsize = 1024

    nbatches = nobs(Xtrn) ÷ bsize
    if nbatches == 0
        nbatches = 1
        bsize = nobs(Xtrn)
    end
    @info "#batches: $nbatches, bsize: $bsize, #epochs: $epochs"
    ps = params(model)
    @info "#params: $(sum(length.(ps)))"

    test_loss, test_acc, test_cm = evaluatebatch(Xtst, Ttst)
    @info @sprintf "initial: test: %0.7f (%0.2f%%)\n" test_loss test_acc

    # dataset = flatten(repeated(RandomBatches((Xtrn, Ttrn), bsize, nbatches), epochs))
    # dataset = RandomBatches((Xtrn, Ttrn), bsize, epochs*nbatches)
    dataset = RandomBatches((Xtrn, Ttrn), bsize, nbatches)
    # @info "dataset ready: type = $(typeof(dataset))"

    for epoch in 1:epochs
        epoch_start_time = time()

        Flux.train!(loss, ps, dataset, opt; cb = () -> ())

        ste = time()
        train_loss, train_acc, train_cm = evaluatebatch(Xtrn, Ttrn)
        test_loss, test_acc, test_cm = evaluatebatch(Xtst, Ttst)
        este = time() - ste

        if test_loss < best_loss
            best_loss, best_epoch = test_loss, epoch
            cpumodel = model |> cpu
            BSON.@save joinpath(model_dir, "model.bson") cpumodel
            improvestr = "*"
        else
            improvestr = "-$(epoch-best_epoch)"
        end
        etime = time() - epoch_start_time
        ttime = time() - total_start_time

        if epoch % logepochs == 0
            @info @sprintf "E%d/%d: train: %0.7f (%0.2f%%) test: %0.7f (%0.2f%%) in %0.3fs (eval %0.3fs, total %0.3fs) %s\n" epoch epochs train_loss train_acc test_loss test_acc etime este ttime improvestr
        end
        if epoch >= best_epoch + early_stop
            @info "early stopping in epoch $epoch\n"
            break
        end
    end

    BSON.@load joinpath(model_dir, "model.bson") cpumodel
    model = cpumodel |> gpu

    train_loss, train_acc = NaN32, NaN32
    train_loss, train_acc, train_cm = evaluatebatch(Xtrn, Ttrn)
    test_loss, test_acc, test_cm = evaluatebatch(Xtst, Ttst)

    @info @sprintf "final: train: %0.5f/%0.2f%%, test: %0.5f/%0.2f%%\n" train_loss train_acc test_loss test_acc
    @info "train confusion matrix:\n$(train_cm)"
    @info "test confusion matrix:\n$(test_cm)"

    model
end

function main()
    hidden = [1024, 512, 256] # default
    # hidden = [512, 512, 256]
    # hidden = [512, 512, 256, 128]
    # hidden = [1024, 512, 256, 128]
    # hidden = [64, 64, 64, 64]
    # hidden = [1024, 1024, 1024]
    # hidden = [2048, 2048, 2048]

    # data_dir = "full/1/group_data_size"
    # data_dir = "full/2/group_data_size"
    # data_dir = "full/3/group_data_size"
    # data_dir = "full/4/group_data_size"
    # data_dir = "full/5/group_data_size"
    # data_dir = "full/6/group_data_size"
    # data_dir = "full/7/group_data_size"
    data_dir = "full/8/group_data_size"
    # data_dir = "full/9/group_data_size"
    # data_dir = "full/10/group_data_size"

    model_dir = "exp/$(data_dir)/$(join(hidden, "_"))"
    ispath(model_dir) && error("path $model_dir exists!")
    mkpath(model_dir)

    append_default_flushed_logger(joinpath(model_dir, "log.txt"))

    @info "importing dataset: $(data_dir)"
    data, dtrans = import_dataset("data/$(data_dir)", 0.98)
    @info "dataset imported"
    Xtrn = data[:train, :X] |> getobs |> gpu
    Ttrn = data[:train, :T] |> getobs |> gpu
    Xtst = data[:test, :X] |> getobs |> gpu
    Ttst = data[:test, :T] |> getobs |> gpu

    @info "copied to GPU"
    @info "Xtrn: $(summary(Xtrn)), Ttrn: $(summary(Ttrn)), #true: $(Int(sum(Ttrn))), #false: $(size(Ttrn, 2) - Int(sum(Ttrn)))"
    @info "Xtst: $(summary(Xtst)), Ttst: $(summary(Ttst)), #true: $(Int(sum(Ttst))), #false: $(size(Ttst, 2) - Int(sum(Ttst)))"

    train(Xtrn, Ttrn, Xtst, Ttst; model_dir=model_dir, hidden=hidden)
end

main()


text:
using Documenter, Gillespie

makedocs(modules=[Gillespie],
    doctest = false)

deploydocs(
    deps = Deps.pip("mkdocs","python-markdown-math"),
    repo = "github.com/sdwfrost/Gillespie.jl.git",
    julia = "release"
)


text:
include("periodize.jl")
include("kdtree.jl")
include("vxlcut.jl")
include("plotting.jl")


text:
# This file contains functions for the Generalized Distance Function DEA model
"""
    GeneralizedDFDEAModel
An data structure representing a generalized distance function DEA model.
"""
struct GeneralizedDFDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    alpha::Float64
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deagdf(X, Y, alpha)
Compute generalized distance function data envelopment analysis model for
inputs `X`, outputs `Y`, and `alpha`.

# Optional Arguments
- `alpha=0.5`: alpha value.
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `slack=true`: compute input and output slacks.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9];

julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6];

julia> deagdf(X, Y, alpha = 0.5, rts = :VRS)
Generalized DF DEA Model
DMUs = 5; Inputs = 2; Outputs = 2
alpha = 0.5; Returns to Scale = VRS
─────────────────────────────────────────────────────────────
   efficiency     slackX1     slackX2     slackY1     slackY2
─────────────────────────────────────────────────────────────
1     0.68185  0.605935    4.26672e-8  3.91163e-8  4.67865
2     1.0      5.34772e-8  6.70059e-8  2.7034e-8   4.8232e-8
3     1.0      4.82929e-8  7.26916e-8  6.00225e-8  1.66806e-8
4     0.25     4.6491e-8   9.94558e-8  5.39305e-8  9.75587e-8
5     0.36     0.2         3.4         3.0         8.07052e-8
─────────────────────────────────────────────────────────────
```
"""
function deagdf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    alpha::Float64 = 0.5, rts::Symbol = :CRS, slack::Bool = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::GeneralizedDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:NLP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, eff, start = 1.0)
        @variable(deamodel, lambda[1:nref] >= 0)

        @NLobjective(deamodel, Min, eff)

        @NLconstraint(deamodel, [j in 1:m], sum(Xref[t,j] * lambda[t] for t in 1:nref) <= eff^(1-alpha) * x0[j])
        @NLconstraint(deamodel, [j in 1:s], sum(Yref[t,j] * lambda[t] for t in 1:nref) >= y0[j] / (eff^alpha) )

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(lambda) == 1)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        lambdaeff[i,:] = JuMP.value.(lambda)

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    # Get first-stage X and Y targets
    Xtarget = X .* effi .^(1-alpha)
    Ytarget = Y ./ ( effi .^alpha )

    # Compute slacks
    if slack == true
        # Use additive model with X and Y targets to get slacks
        slacksmodel = deaadd(Xtarget, Ytarget, :Ones, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
        slackX = slacks(slacksmodel, :X)
        slackY = slacks(slacksmodel, :Y)

        # Get second-stage X and Y targets
        Xtarget = Xtarget - slackX
        Ytarget = Ytarget + slackY
    else
        if typeof(Xtarget) <: AbstractVector    Xtarget = Xtarget[:,:]  end
        if typeof(Ytarget) <: AbstractVector    Ytarget = Ytarget[:,:]  end

        slackX = Array{Float64}(undef, 0, 0)
        slackY = Array{Float64}(undef, 0, 0)
    end

    return GeneralizedDFDEAModel(n, m, s, alpha, rts, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::GeneralizedDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    dmunames = names(x)

    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Generalized DF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "alpha = ", x.alpha)
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], dmunames))
        end
    end

end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function arrow_test(fname)
    pagewidth, pageheight = 2000, 2000
    Drawing(pagewidth, pageheight, fname)
    origin() # move 0/0 to center
    background("ivory")
    setopacity(0.8)
    pagetiles = Tiler(pagewidth, pageheight, 40, 8, margin=50)
    theta = 0
    for theta in range(0, step=π/8, stop = 2π - π/8)
        @layer begin
            rotate(theta)
            translate(120, 0)
            sethue("magenta")
            line(O, O + (230, 230), :stroke)
            sethue("grey10")
            arrow(O, O + (230, 230), [rand(-80:80), rand(-80:80)],
                :stroke,
                linewidth=3,
                startarrow=rand(Bool),
                finisharrow=rand(Bool),
                arrowheadlength=30,
                decoration = rescale(theta, 0, 2π),
                decorate=() -> (randomhue(), circle(O, rescale(theta, 0, 2π, 5, 20), :fill)))
        end
    end
    @test finish() == true
end

Random.seed!(42)
fname = "arrow-bezier-test.pdf"
arrow_test(fname)
println("...finished test: output in $(fname)")


text:
#+FindTag,Extracted_Item_Base
struct Extract_Tag_Result <: Extracted_Item_Base
    _tokenized::Tokenized
    _idx::Int
    _tags::Array{String,1}
    _link::String
end

#+FindTag
tag_idx(o::Extract_Tag_Result)::Int = o._idx

#+FindTag
skip(o::Extract_Tag_Result)::Int = tag_idx(o)+1

#+FindTag
tags(o::Extract_Tag_Result)::Array{String,1} = o._tags

#+FindTag
link(o::Extract_Tag_Result)::String = o._link

#+FindTag
line(o::Extract_Tag_Result)::Int = first(startpos(o._tokenized[o._idx]))

#+FindTag L:extract_tag_link
function extract_tag(tok::Tokenized,idx::Int)::Union{Nothing,Extract_Tag_Result}

    if !is_comment(tok,idx)
        return nothing
    end

    s = untokenize(tok[idx])

    first_tag=match(r"^#[ ]?\+(\w+)",s)

    # Math tag?
    if first_tag==nothing
        return nothing
    end

    tag_list = Array{String,1}([first_tag[1]])

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

    regex_other_tags = r"\s*,\s*(\w+)"
    match_other_tags=first_tag

    @label more_tags
    offset = match_other_tags.offsets[end]+length(match_other_tags[1])
    match_other_tags=match(regex_other_tags,s,offset)

    if match_other_tags!=nothing
        push!(tag_list,match_other_tags[1])
        @goto more_tags
    end

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

    regex_link = r"\s*L:(\w+)\s*$"

    match_link=match(regex_link,s,offset)

    if match_link!=nothing
        link=match_link[1]
    else
        link=""
    end

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

    return Extract_Tag_Result(tok,idx,tag_list,link)
end

#+FindTag L:find_tag_master
#
# If some tags are found use the
# #+BEGIN_SRC julia
# predicate(r::Extract_Tag_Result)::Bool
# #+END_SRC
# to accept (or not) the tags
#
# *Returns:*
# - =(idx,Extract_Tag_Result)=: where =idx= is the position of the discovered tag
# - =(length(tok)+1,nothing)=: if no tag found
#
function find_tag(tok::Tokenized,idx::Int,predicate::Function)::Union{Nothing,Extract_Tag_Result}
    n = length(tok)

    for k in idx:n

        extracted = extract_tag(tok,k)

        if (extracted!=nothing)&&(predicate(extracted))
            return extracted
        end

    end

    return nothing
end

#+FindTag
#
# Convenience function that uses [[find_tag_master][]] with =x->true=
# predicate. It accepts all tags
function find_tag(tok::Tokenized,idx::Int)::Union{Nothing,Extract_Tag_Result}
    return find_tag(tok,idx,x->true)
end

#+FindTag,Obsolete
#
# Convenience function that uses  find_tag_master find_tag  with a
# predicate that checks for *Tag* existence.
#
function find_tag(tok::Tokenized,idx::Int,tag::String)::Union{Nothing,Extract_Tag_Result}
    return find_tag(tok,idx,x->(tag ∈ tags(x)))
end


text:
type DistMultiVec{T} <: ElementalMatrix{T}
    obj::Ptr{Void}
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:Complex64, :c),
                    (:Complex128, :z))
    @eval begin
        function DistMultiVec(::Type{$elty}, cm::ElComm = CommWorld)
            obj = Ref{Ptr{Void}}(C_NULL)
            err = ccall(($(string("ElDistMultiVecCreate_", ext)), libEl), Cuint,
                (Ref{Ptr{Void}}, ElComm),
                obj, cm)
            err == 0 || throw(ElError(err))
            return DistMultiVec{$elty}(obj[])
        end

        function comm(A::DistMultiVec{$elty})
            cm = Ref{ElComm}()
            err = ccall(($(string("ElDistMultiVecComm_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElComm}),
                A.obj, cm)
            err == 0 || throw(ElError(err))
            return cm[]
        end

        function get(x::DistMultiVec{$elty}, i::Integer = size(x, 1), j::Integer = 1)
            v = Ref{$elty}()
            err = ccall(($(string("ElDistMultiVecGet_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, Ref{$elty}),
                x.obj, i - 1, j - 1, v)
            err == 0 || throw(ElError(err))
            return v[]
        end

        function getLocal(A::DistMultiVec{$elty}, i::Integer, j::Integer)
            rv = Ref{$elty}(0)
            err = ccall(($(string("ElDistMultiVecGetLocal_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, Ref{$elty}),
                A.obj, i - 1, j - 1, rv)
            err == 0 || throw(ElError(err))
            return rv[]
        end

        function globalRow(A::DistMultiVec{$elty}, i::Integer)
            rv = Ref{ElInt}(0)
            err = ccall(($(string("ElDistMultiVecGlobalRow_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, Ref{ElInt}),
                A.obj, i - 1, rv)
            err == 0 || throw(ElError(err))
            return rv[] + 1
        end

        function height(x::DistMultiVec{$elty})
            i = Ref{ElInt}()
            err = ccall(($(string("ElDistMultiVecHeight_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                x.obj, i)
            err == 0 || throw(ElError(err))
            return i[]
        end

        function localHeight(A::DistMultiVec{$elty})
            rv = Ref{ElInt}(0)
            err = ccall(($(string("ElDistMultiVecLocalHeight_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                A.obj, rv)
            err == 0 || throw(ElError(err))
            return rv[]
        end

        function processQueues(A::DistMultiVec{$elty})
            err = ccall(($(string("ElDistMultiVecProcessQueues_", ext)), libEl), Cuint,
                (Ptr{Void},), A.obj)
            err == 0 || throw(ElError(err))
          return nothing
        end

        function queueUpdate(A::DistMultiVec{$elty}, i::Integer, j::Integer, value::$elty)
            err = ccall(($(string("ElDistMultiVecQueueUpdate_", ext)), libEl), Cuint,
              (Ptr{Void}, ElInt, ElInt, $elty),
              A.obj, i - 1, j - 1, value)
            err == 0 || throw(ElError(err))
            return nothing
        end

        function reserve(A::DistMultiVec{$elty}, numEntries::Integer)
            err = ccall(($(string("ElDistMultiVecReserve_", ext)), libEl), Cuint,
              (Ptr{Void}, ElInt),
              A.obj, numEntries)
            err == 0 || throw(ElError(err))
            return nothing
        end

        function resize!(A::DistMultiVec{$elty}, m::Integer, n::Integer = 1) # to mimic vector behavior
            err = ccall(($(string("ElDistMultiVecResize_", ext)), libEl), Cuint,
              (Ptr{Void}, ElInt, ElInt),
              A.obj, ElInt(m), ElInt(n))
            err == 0 || throw(ElError(err))
            return A
        end

        function width(x::DistMultiVec{$elty})
            i = Ref{ElInt}()
            err = ccall(($(string("ElDistMultiVecWidth_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
                x.obj, i)
            err == 0 || throw(ElError(err))
            return i[]
        end
    end
end

#########################
### Julia convenience ###
#########################

getindex(x::DistMultiVec, i, j) = get(x, i, j)

function similar{T}(::DistMultiVec, ::Type{T}, sz::Dims, cm::ElComm = CommWorld)
    A = DistMultiVec(T, cm)
    resize!(A, sz...)
    return A
end

# FixMe! Should this one handle vectors of matrices?
function hcat{T}(x::Vector{DistMultiVec{T}})
    l    = length(x)
    if l == 0
        throw(ArgumentError("cannot flatten empty vector"))
    else
        x1   = x[1]
        m, n = size(x1, 1), size(x1, 2)
        if n != 1
            throw(ArgumentError("elements has to be vectors, i.e. the second dimension has to have size one"))
        end
        A    = DistMultiVec(T)
        zeros!(A, m, l*n)
        for j = 1:l
            xj = x[j]
            for k = 1:width(xj)
                for i = 1:localHeight(xj)
                    xji = getLocal(xj, i, 1)
                    queueUpdate(A, globalRow(xj, i), j, xji)
                end
            end
        end
        processQueues(A)
        return A
    end
end