text:
using UnitSystems, Test

for S ∈ UnitSystems.Systems
    U = eval(S)
    S ≠ :IAU && @testset "UnitSystem: $S" begin
        @testset "Dimensionless constants" begin
            @test μₑᵤ == electronmass(U)/atomicmass(U)
            @test μₚᵤ == protonmass(U)/atomicmass(U)
            @test μₚₑ ≈ protonmass(U)/electronmass(U)
            @test 1/αinv ≈ (charge(U)/charge(PlanckGauss,U))^2
            @test αG ≈ (electronmass(U)/mass(PlanckGauss,U))^2
            @test 1/αinv ≈ charge(U)^2*rationalization(U)/4π/permittivity(U)/planckreduced(U)/lightspeed(U)
            @test 1/αinv ≈ permeability(U)*lightspeed(U)*(charge(U)*lorentz(U))^2*rationalization(U)/4π/planckreduced(U)
            @test 1/αinv ≈ coulomb(U)*charge(U)^2/planckreduced(U)/lightspeed(U)
            @test 1/αinv ≈ charge(U)^2*rationalization(U)/2permittivity(U)/lightspeed(U)/planck(U)
            @test 1/αinv ≈ lightspeed(U)*permeability(U)*rationalization(U)*lorentz(U)^2/2klitzing(U)
            @test 1/αinv ≈ charge(U)^2*impedance(U)/2planck(U)
        end
        @testset "Fundamental constants" begin
            @testset "lightspeed" begin
                @test lightspeed(U) ≈ 1/lorentz(U)/sqrt(permeability(U)*permittivity(U))
                @test lightspeed(U) ≈ αinv*sqrt(hartree(U)/electronmass(U))
                @test lightspeed(U) ≈ planckreduced(U)/αinv/electronmass(U)/electronradius(U)
                @test lightspeed(U) ≈ charge(U)^2*coulomb(U)/planckreduced(U)*αinv
                @test lightspeed(U) ≈ electronmass(U)^2*newton(U)/planckreduced(U)/αG
            end
            @testset "planck" begin
                @test planck(U) == 2π*planckreduced(U)
                @test planck(U) == 2charge(U)*lorentz(U)/josephson(U)
                @test planck(U) ≈ 8/αinv/rationalization(U)/lightspeed(U)/permeability(U)/josephson(U)^2
                @test planck(U) ≈ 4lorentz(U)^2/josephson(U)^2/klitzing(U)
            end
            @testset "planckmass" begin
                @test planckmass(U) ≈ sqrt(planckreduced(U)*lightspeed(U)/newton(U))
                @test planckmass(U) ≈ electronmass(U)/sqrt(αG)
                @test planckmass(U) ≈ 2rydberg(U)*planck(U)/lightspeed(U)*αinv^2/sqrt(αG)
            end
            @testset "newton" begin
                @test newton(U) == planckreduced(U)*lightspeed(U)/planckmass(U)^2
                @test newton(U) ≈ planckreduced(U)*lightspeed(U)*αG/electronmass(U)^2
                @test newton(U) ≈ lightspeed(U)^3/αinv^4*αG/8π/rydberg(U)^2/planck(U)
                @test newton(U) ≈ einstein(U)*lightspeed(U)^4/8π
            end
            @testset "einstein" begin
                @test einstein(U) ≈ 8π*newton(U)/lightspeed(U)^4
                @test einstein(U) ≈ 8π*planckreduced(U)/lightspeed(U)^3/planckmass(U)^2
                @test einstein(U) ≈ 8π*planckreduced(U)*αG/lightspeed(U)^3/electronmass(U)^2
                @test einstein(U) ≈ αG/αinv^4/rydberg(U)^2/planck(U)/lightspeed(U)
            end
        end
        @testset "Atomic constants" begin
            @testset "atomicmass" begin
                @test atomicmass(U) ≈ molarmass(U)/avogadro(U)
                @test atomicmass(U) == electronmass(U)/μₑᵤ
                @test atomicmass(U) == protonmass(U)/μₚᵤ
                @test atomicmass(U) ≈ 2rydberg(U)*planck(U)/μₑᵤ/lightspeed(U)*αinv^2
                @test atomicmass(U) ≈ planckmass(U)*sqrt(αG)/μₑᵤ
            end
            @testset "protonmass" begin
                @test protonmass(U) == μₚᵤ*atomicmass(U)
                @test protonmass(U) ≈ μₚᵤ*molarmass(U)/avogadro(U)
                @test protonmass(U) == μₚₑ*electronmass(U)
                @test protonmass(U) ≈ μₚₑ*2rydberg(U)*planck(U)/lightspeed(U)*αinv^2
                @test protonmass(U) ≈ planckmass(U)*μₚₑ*sqrt(αG)
            end
            @testset "electronmass" begin
                @test electronmass(U) == μₑᵤ*atomicmass(U)
                @test electronmass(U) == μₑᵤ*molarmass(U)/avogadro(U)
                @test electronmass(U) ≈ protonmass(U)/μₚₑ
                @test electronmass(U) ≈ 2rydberg(U)*planck(U)/lightspeed(U)*αinv^2
                @test electronmass(U) ≈ planckmass(U)*sqrt(αG)
            end
            @testset "hartree" begin
                @test hartree(U) ≈ electronmass(U)*(lightspeed(U)/αinv)^2
                @test hartree(U) ≈ planckreduced(U)*lightspeed(U)/αinv/bohr(U)
                @test hartree(U) ≈ planckreduced(U)^2/electronmass(U)/bohr(U)^2
                @test hartree(U) ≈ 2rydberg(U)*planck(U)*lightspeed(U)
                @test hartree(U) ≈ planckmass(U)*sqrt(αG)*(lightspeed(U)/αinv)^2
            end
            @testset "rydberg" begin
                @test rydberg(U) == hartree(U)/2planck(U)/lightspeed(U)
                @test rydberg(U) ≈ electronmass(U)*lightspeed(U)/αinv^2/2planck(U)
                @test rydberg(U) ≈ 1/αinv/4π/bohr(U)
                @test rydberg(U) ≈ electronmass(U)*electronradius(U)*lightspeed(U)/2planck(U)/bohr(U)
                @test rydberg(U) ≈ electronmass(U)*lightspeed(U)/αinv^2/4π/planckreduced(U)
                @test rydberg(U) ≈ planckmass(U)*lightspeed(U)*sqrt(αG)/αinv^2/2planck(U)
            end
            @testset "bohr" begin
                @test bohr(U) ≈ planckreduced(U)/electronmass(U)/lightspeed(U)*αinv
                @test bohr(U) ≈ planckreduced(U)^2/coulomb(U)/electronmass(U)/charge(U)^2
                @test bohr(U) ≈ μₚₑ*bohrreduced(U)/(μₚₑ+1)
                @test bohr(U) ≈ electronradius(U)*αinv^2
                @test bohr(U) ≈ 1/αinv/4π/rydberg(U)
                @test bohrreduced(U) == (1+1/μₚₑ)*bohr(U)
            end
            @testset "electronradius" begin
                @test electronradius(U) ≈ planckreduced(U)/αinv/electronmass(U)/lightspeed(U)
                @test electronradius(U) ≈ bohr(U)/αinv^2
                @test electronradius(U) ≈ charge(U)^2*coulomb(U)/electronmass(U)/lightspeed(U)^2
                @test electronradius(U) ≈ 2planck(U)*rydberg(U)*bohr(U)/electronmass(U)/lightspeed(U)
                @test electronradius(U) ≈ 1/αinv^3/4π/rydberg(U)
            end
        end
        @testset "Thermodynamic constants" begin
            @testset "molarmass" begin
                @test molarmass(U) ≈ atomicmass(U)*avogadro(U)
                @test molarmass(U) ≈ avogadro(U)*electronmass(U)/μₑᵤ
                @test molarmass(U) ≈ avogadro(U)*protonmass(U)/μₚᵤ
                @test molarmass(U) ≈ avogadro(U)*2rydberg(U)*planck(U)/μₑᵤ/lightspeed(U)*αinv^2
            end
            @testset "avogadro" begin
                @test avogadro(U) ≈ universal(U)/boltzmann(U)
                @test avogadro(U) ≈ molarmass(U)/atomicmass(U)
                @test avogadro(U) == molarmass(U)*μₑᵤ/electronmass(U)
                @test avogadro(U) ≈ molarmass(U)*μₑᵤ*lightspeed(U)/αinv^2/2rydberg(U)/planck(U)
            end
            @testset "boltzmann" begin
                @test boltzmann(U) == universal(U)/avogadro(U)
                @test boltzmann(U) ≈ atomicmass(U)*universal(U)/molarmass(U)
                @test boltzmann(U) ≈ electronmass(U)*universal(U)/μₑᵤ/molarmass(U)
                @test boltzmann(U) ≈ 2universal(U)*rydberg(U)*planck(U)/molarmass(U)/μₑᵤ/lightspeed(U)*αinv^2
            end
            @testset "universal" begin
                @test universal(U) == boltzmann(U)*avogadro(U)
                @test universal(U) ≈ boltzmann(U)*molarmass(U)/atomicmass(U)
                @test universal(U) ≈ boltzmann(U)*molarmass(U)*μₑᵤ/electronmass(U)
                @test universal(U) ≈ boltzmann(U)*molarmass(U)*μₑᵤ*lightspeed(U)/αinv^2/2planck(U)/rydberg(U)
            end
            @testset "stefan" begin
                @test stefan(U) ≈ 2π^5*boltzmann(U)^4/15planck(U)^3/lightspeed(U)^2
                @test stefan(U) ≈ π^2*boltzmann(U)^4/60planckreduced(U)^3/lightspeed(U)^2
                @test stefan(U) ≈ 32π^5*planck(U)/15lightspeed(U)^6*αinv^8*(universal(U)*rydberg(U)/μₑᵤ/molarmass(U))^4
            end
            @testset "radiationdensity" begin
                @test radiationdensity(U) ≈ 4stefan(U)/lightspeed(U)
                @test radiationdensity(U) ≈ 8π^5*boltzmann(U)^4/15planck(U)^3/lightspeed(U)^3
                @test radiationdensity(U) ≈ π^2*boltzmann(U)^4/15planckreduced(U)^3/lightspeed(U)^3
                @test radiationdensity(U) ≈ 2^7*π^5*planck(U)/15lightspeed(U)^7*αinv^8*(universal(U)*rydberg(U)/μₑᵤ/molarmass(U))^4
            end
        end
        @testset "Electromagnetic constants" begin
            @testset "rationalization" begin
                @test rationalization(U) ≈ 4π*biotsavart(U)/permeability(U)/lorentz(U)
                @test rationalization(U) ≈ 4π*coulomb(U)*permittivity(U)
                @test rationalization(U) ≈ impedance(U)*permittivity(U)*lightspeed(U)
            end
            @testset "permeability" begin
                @test permeability(U) ≈ 1/permittivity(U)/(lightspeed(U)*lorentz(U))^2
                @test permeability(U) ≈ 4π*coulomb(U)/rationalization(U)/(lightspeed(U)*lorentz(U))^2
                @test permeability(U) ≈ 2planck(U)/αinv/rationalization(U)/lightspeed(U)/(charge(U)*lorentz(U))^2
                @test permeability(U) ≈ 2klitzing(U)/αinv/rationalization(U)/lightspeed(U)/lorentz(U)^2
            end
            @testset "permittivity" begin
                @test permittivity(U) ≈ 1/permeability(U)/(lightspeed(U)*lorentz(U))^2
                @test permittivity(U) ≈ rationalization(U)/4π/coulomb(U)
                @test permittivity(U) ≈ rationalization(U)*charge(U)^2*αinv/2planck(U)/lightspeed(U)
                @test permittivity(U) ≈ rationalization(U)/2klitzing(U)/lightspeed(U)*αinv
            end
            @testset "coulomb" begin
                @test coulomb(U) == rationalization(U)/4π/permittivity(U)
                @test coulomb(U) ≈ permeability(U)*rationalization(U)*(lorentz(U)*lightspeed(U))^2/4π
                @test coulomb(U) ≈ planckreduced(U)*lightspeed(U)/αinv/charge(U)^2
                @test coulomb(U) ≈ klitzing(U)*lightspeed(U)/αinv/2π
                @test coulomb(U) ≈ biotsavart(U)/lorentz(U)/permeability(U)/permittivity(U)
                @test coulomb(U) ≈ ampere(U)*lightspeed(U)^2
            end
            @testset "ampere" begin
                @test ampere(U) == lorentz(U)*biotsavart(U)
                @test ampere(U) ≈ permeability(U)*lorentz(U)^2*rationalization(U)/4π
                @test ampere(U) ≈ coulomb(U)/lightspeed(U)^2
                @test ampere(U) ≈ planckreduced(U)/αinv/lightspeed(U)/charge(U)^2
                @test ampere(U) ≈ klitzing(U)/αinv/2π/lightspeed(U)
            end
            @testset "lorentz" begin
                @test lorentz(U) ≈ 1/lightspeed(U)/sqrt(permeability(U)*permittivity(U))
                @test lorentz(U) ≈ biotsavart(U)/permeability(U)/permittivity(U)/coulomb(U)
                @test lorentz(U) ≈ 4π*biotsavart(U)/rationalization(U)/permeability(U)
                @test lorentz(U) == ampere(U)/biotsavart(U)
            end
            @testset "biotsavart" begin
                @test biotsavart(U) ≈ permeability(U)*lorentz(U)*rationalization(U)/4π
                @test biotsavart(U) ≈ lorentz(U)*permeability(U)*permittivity(U)*coulomb(U)
                @test biotsavart(U) == ampere(U)/lorentz(U)
                @test biotsavart(U) ≈ coulomb(U)*sqrt(permeability(U)*permittivity(U))/lightspeed(U)
            end
            @testset "charge" begin
                @test charge(U) ≈ sqrt(2planck(U)/αinv/impedance(U))
                @test charge(U) ≈ 2lorentz(U)/josephson(U)/klitzing(U)
                @test charge(U) == sqrt(planck(U)/klitzing(U))
                @test charge(U) ≈ planck(U)*josephson(U)/2lorentz(U)
                @test charge(U) == faraday(U)/avogadro(U)
            end
            @testset "faraday" begin
                @test faraday(U) == charge(U)*avogadro(U)
                @test faraday(U) ≈ avogadro(U)*sqrt(2planck(U)/αinv/impedance(U))
                @test faraday(U) ≈ 2avogadro(U)*lorentz(U)/josephson(U)/klitzing(U)
                @test faraday(U) == avogadro(U)*sqrt(planck(U)/klitzing(U))
                @test faraday(U) ≈ planck(U)*josephson(U)*avogadro(U)/2lorentz(U)
            end
            @testset "impedance" begin
                @test impedance(U) ≈ permeability(U)*rationalization(U)*lightspeed(U)*lorentz(U)^2
                @test impedance(U) ≈ rationalization(U)/permittivity(U)/lightspeed(U)
                @test impedance(U) ≈ rationalization(U)*lorentz(U)*sqrt(permeability(U)/permittivity(U))
                @test impedance(U) ≈ 2planck(U)/αinv/charge(U)^2
                @test impedance(U) ≈ 2klitzing(U)/αinv
            end
            @testset "conductance" begin
                @test conductance(U) == 2charge(U)^2/planck(U)
                @test conductance(U) ≈ 4/αinv/impedance(U)
                @test conductance(U) ≈ 2/klitzing(U)
                @test conductance(U) ≈ planck(U)*josephson(U)^2/2lorentz(U)^2
                @test conductance(U) ≈ 2faraday(U)^2/planck(U)/avogadro(U)^2
            end
            @testset "klitzing" begin
                @test klitzing(U) == planck(U)/charge(U)^2
                @test klitzing(U) ≈ impedance(U)*αinv/2
                @test klitzing(U) ≈ 2/conductance(U)
                @test klitzing(U) ≈ 4lorentz(U)^2/planck(U)/josephson(U)^2
                @test klitzing(U) ≈ planck(U)*avogadro(U)^2/faraday(U)^2
            end
            @testset "josephson" begin
                @test josephson(U) == 2charge(U)*lorentz(U)/planck(U)
                @test josephson(U) ≈ lorentz(U)*sqrt(8/αinv/planck(U)/impedance(U))
                @test josephson(U) ≈ lorentz(U)*sqrt(4/planck(U)/klitzing(U))
                @test josephson(U) ≈ 1/magneticflux(U)
                @test josephson(U) ≈ 2faraday(U)*lorentz(U)/planck(U)/avogadro(U)
            end
            @testset "magneticflux" begin
                @test magneticflux(U) ≈ planck(U)/2charge(U)/lorentz(U)
                @test magneticflux(U) ≈ sqrt(planck(U)*impedance(U)*αinv/8)/lorentz(U)
                @test magneticflux(U) ≈ sqrt(planck(U)*klitzing(U)/4)/lorentz(U)
                @test magneticflux(U) == 1/josephson(U)
                @test magneticflux(U) ≈ planck(U)*avogadro(U)/2faraday(U)/lorentz(U)
            end
            @testset "magneton" begin
                @test magneton(U) == charge(U)*planckreduced(U)*lorentz(U)/2electronmass(U)
                @test magneton(U) ≈ planckreduced(U)*lorentz(U)^2/electronmass(U)/josephson(U)/klitzing(U)
                @test magneton(U) ≈ planck(U)^2*josephson(U)/8π/electronmass(U)
                @test magneton(U) ≈ lorentz(U)*planckreduced(U)*faraday(U)/2electronmass(U)/avogadro(U)
                @test magneton(U) ≈ charge(U)*lightspeed(U)*lorentz(U)/αinv^2/8π/rydberg(U)
            end
        end
    end
end

@testset "CGS conversions" begin
    @test molarmass(Natural) == molarmass(CGS) == 1000molarmass(Metric)
    @test molarmass(CGS2019) == 1000molarmass(SI2019)

    C = 100𝘤

    @test charge(Metric,ESU) ≈ C/10
    @test charge(Metric,EMU) ≈ 1/10
    @test charge(Metric,Gauss) ≈ C/10

    @test current(Metric,ESU) ≈ C/10
    @test current(Metric,EMU) ≈ 1/10
    @test current(Metric,Gauss) ≈ C/10

    @test electricpotential(Metric,ESU) ≈ 1e8/C
    @test electricpotential(Metric,EMU) ≈ 1e8
    @test electricpotential(Metric,Gauss) ≈ 1e8/C

    @test electricfield(Metric,ESU) ≈ 1e6/C
    @test electricfield(Metric,EMU) ≈ 1e6
    @test electricfield(Metric,Gauss) ≈ 1e6/C

    @test electricfluxdensity(Metric,ESU) ≈ 4π*C/1e5
    @test electricfluxdensity(Metric,EMU) ≈ 4π/1e5
    @test electricfluxdensity(Metric,Gauss) ≈ 4π*C/1e5

    @test electricdipolemoment(Metric,ESU) ≈ 10C
    @test electricdipolemoment(Metric,EMU) ≈ 10
    @test electricdipolemoment(Metric,Gauss) ≈ 10C

    @test magneticdipolemoment(Metric,ESU) ≈ 1e3*C
    @test magneticdipolemoment(Metric,EMU) ≈ 1e3
    @test magneticdipolemoment(Metric,Gauss) ≈ 1e3

    @test magneticfield(Metric,ESU) ≈ 4π*C/1e3
    @test magneticfield(Metric,EMU) ≈ 4π/1e3
    @test magneticfield(Metric,Gauss) ≈ 4π/1e3

    @test magneticfluxdensity(Metric,ESU) ≈ 1e4/C
    @test magneticfluxdensity(Metric,EMU) ≈ 1e4
    @test magneticfluxdensity(Metric,Gauss) ≈ 1e4

    @test magneticflux(Metric,ESU) ≈ 1e8/C
    @test magneticflux(Metric,EMU) ≈ 1e8
    @test magneticflux(Metric,Gauss) ≈ 1e8

    @test resistance(Metric,ESU) ≈ 1e9/C^2
    @test resistance(Metric,EMU) ≈ 1e9
    @test resistance(Metric,Gauss) ≈ 1e9/C^2

    @test resistivity(Metric,ESU) ≈ 1e11/C^2
    @test resistivity(Metric,EMU) ≈ 1e11
    @test resistivity(Metric,Gauss) ≈ 1e11/C^2

    @test capacitance(Metric,ESU) ≈ C^2/1e9
    @test capacitance(Metric,EMU) ≈ 1e-9
    @test capacitance(Metric,Gauss) ≈ C^2/1e9

    @test inductance(Metric,ESU) ≈ 1e9/C^2
    @test inductance(Metric,EMU) ≈ 1e9
    @test inductance(Metric,Gauss) ≈ 1e9/C^2

    # extra

    @test conductance(Metric,ESU) ≈ C^2/1e9
    @test conductance(Metric,EMU) ≈ 1e-9
    @test conductance(Metric,Gauss) ≈ C^2/1e9

    @test chargedensity(Metric,ESU) ≈ C/1e7
    @test chargedensity(Metric,EMU) ≈ 1e-7
    @test chargedensity(Metric,Gauss) ≈ C/1e7

    @test magneticpotential(Metric,ESU) ≈ 4π/10*C
    @test magneticpotential(Metric,EMU) ≈ 4π/10
    @test magneticpotential(Metric,Gauss) ≈ 4π/10

    @test susceptibility(Metric,ESU) ≈ 1/4π
    @test susceptibility(Metric,EMU) ≈ 1/4π
    @test susceptibility(Metric,Gauss) ≈ 1/4π

    # magnetisation ?

    #@test magneticpolarization(Metric,ESU) ≈ 1e2/4π/𝘤
    #@test magneticpolarization(Metric,EMU) ≈ 1e4/4π
    #@test magneticpolarization(Metric,Gauss) ≈ 1e4/4π

    @test polestrength(Metric,ESU) ≈ 10C
    @test polestrength(Metric,EMU) ≈ 10
    @test polestrength(Metric,Gauss) ≈ 10

    #@test reluctance(Metric,ESU) ≈ 4π/1e9 ?N/A?
    @test reluctance(Metric,EMU) ≈ 4π/1e9
    @test reluctance(Metric,Gauss) ≈ 4π/1e9

    # other

    @test currentdensity(Metric,ESU) ≈ C/1e5
    @test currentdensity(Metric,EMU) ≈ 1e-5
    @test currentdensity(Metric,Gauss) ≈ C/1e5

    @test permittivity(ESU,Metric) ≈ ε₀
    @test permeability(EMU,Metric) ≈ μ₀
    @test permeability(Gauss,Metric) ≈ μ₀

    @test specificsusceptibility(ESU,Metric) ≈ 4π/1e3
    @test specificsusceptibility(EMU,Metric) ≈ 4π/1e3
    @test specificsusceptibility(Gauss,Metric) ≈ 4π/1e3

    @test demagnetizingfactor(ESU,Metric) ≈ 1/4π
    @test demagnetizingfactor(EMU,Metric) ≈ 1/4π
    @test demagnetizingfactor(Gauss,Metric) ≈ 1/4π

    @test electricpolarizability(Metric,EMU) ≈ 1e-5 # ??
    @test electricpolarizability(Metric,ESU) ≈ 1e6/4π/ε₀
    @test electricpolarizability(Metric,Gauss) ≈ 1e6/4π/ε₀

    @test magneticpolarizability(Metric,ESU) ≈ 1e6/4π
    @test magneticpolarizability(Metric,EMU) ≈ 1e6/4π
    @test magneticpolarizability(Metric,Gauss) ≈ 1e6/4π

    # uncertain:
    @test electricflux(Metric,ESU) ≈ 1e10/C
    @test electricflux(Metric,EMU) ≈ 1e10
    @test electricflux(Metric,Gauss) ≈ 1e10/C
    @test magneticmoment(Metric,ESU) ≈ 1e10/C # prefer: 10/𝘤
    @test magneticmoment(Metric,EMU) ≈ 1e10 # prefer: 1e3
    @test magneticmoment(Metric,Gauss) ≈ 1e10 # prefer: 1e3
    @test specificmagnetization(ESU,Metric) ≈ 1e7/C # prefer: 1
    @test specificmagnetization(EMU,Metric) ≈ 1e7 # prefer: 1
    @test specificmagnetization(Gauss,Metric) ≈ 1e7 # prefer: 1
end


text:
#=
See also:
https://github.com/FluxML/Flux.jl/blob/master/src/onehot.jl
https://github.com/cossio/OneHot.jl
=#

"""
    onehot_encode(A, code)

Given an array `A` of `N` dimensions, returns a one-hot encoded `BitArray` of
`N + 1` dimensions where single entries of the first dimension are one.
"""
function onehot_encode(A::AbstractArray, code = sort(unique(A)))
    return reshape(A, 1, size(A)...) .== code
end

"""
    onehot_decode(X)

Given a onehot encoded array `X` of `N + 1` dimensions, returns the
equivalent categorical array of `N` dimensions.
"""
function onehot_decode(X::AbstractArray)
    idx = dropdims(argmax(X; dims=1); dims=1)
    return first.(Tuple.(idx))
end

"""
    categorical_sample_from_logits(logits)

Given a logits array `logits` of size `(q, *)` (where `q` is the number of classes),
returns an array `X` of size `(*)`, such that `X[i]` is a categorical random sample
from the distribution with logits `logits[:,i]`.
"""
function categorical_sample_from_logits(logits::AbstractArray)
    p = softmax(logits; dims=1)
    return categorical_sample(p)
end

"""
    categorical_sample_from_logits_gumbel(logits)

Like categorical_sample_from_logits, but using the Gumbel trick.
"""
function categorical_sample_from_logits_gumbel(logits::AbstractArray)
    z = logits .+ randgumbel.(eltype(logits))
    idx = dropdims(argmax(z; dims=1); dims=1)
    c = first.(Tuple.(idx))
    return c
end

"""
    categorical_sample(P)

Given a probability array `P` of size `(q, *)`, returns an array
`C` of size `(*)`, such that `C[i] ∈ 1:q` is a random sample from the
categorical distribution `P[:,i]`.
You must ensure that `P` defines a proper probability distribution.
"""
function categorical_sample(P::AbstractArray)
    idx = CartesianIndices(tail(size(P)))
    C = Array{Int}(undef, tail(size(P)))
    @inbounds for i in idx
        ps = @view P[:,i]
        C[i] = categorical_rand(ps)
    end
    return C
end

"""
    categorical_rand(ps)

Randomly draw `i` with probability `ps[i]`.
You must ensure that `ps` defines a proper probability distribution.
"""
function categorical_rand(ps::AbstractVector)
    i = 0
    u = rand(eltype(ps))
    for p in ps
        u -= p
        i += 1
        u ≤ 0 && break
    end
    return i
end

"""
    randgumbel(T = Float64)

Generates a random Gumbel variate.
"""
randgumbel(::Type{T} = Float64) where {T} = -log(randexp(T))

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ b78a501c-fda3-11ea-06bb-5175ada30398
using Pkg, DrWatson

# ╔═╡ b78a9662-fda3-11ea-25be-871b1258e163
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ b0b04904-fda3-11ea-2a13-efbb8da05490
md"## Clip-05-28-34s.jl"

# ╔═╡ b78b1934-fda3-11ea-119c-5398009c1521
md"### snippet 5.28 - 5.31"

# ╔═╡ b7997a6a-fda3-11ea-325a-7d0fc706cdff
begin
	df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';')
	df.lmass = log.(df.mass)
	df = filter(row -> !(row[:neocortex_perc] == "NA"), df)
	df.neocortex_perc = parse.(Float64, df.neocortex_perc)
	scale!(df, [:kcal_per_g, :neocortex_perc, :lmass])
end;

# ╔═╡ 9fe858dc-fda3-11ea-38c8-df4b1d6e8419
stan_5_5_draft = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] K; // Outcome
 vector[N] NC; // Predictor
}
parameters {
 real a; // Intercept
 real bN; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}
model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 1);           //Priors
  bN ~ normal(0, 1);
  sigma ~ exponential(1);
  mu = a + bN * NC;
  K ~ normal(mu , sigma);     // Likelihood
}
";

# ╔═╡ d65e8eea-fda3-11ea-0f78-6d8a183f5818
md"##### Define the SampleModel, etc."

# ╔═╡ d65ec40a-fda3-11ea-2e15-2f55dd308dd9
begin
	m5_5_drafts = SampleModel("m5.5.draft", stan_5_5_draft);
	m5_5_data = Dict("N" => size(df, 1), "NC" => df.neocortex_perc_s,
		"K" => df.kcal_per_g_s);
	rc5_5_drafts = stan_sample(m5_5_drafts, data=m5_5_data)
end;

# ╔═╡ 407845b4-fd4a-11ea-1ed5-65581b68ccda
if success(rc5_5_drafts)
  post5_5_drafts_df = read_samples(m5_5_drafts, :dataframe)
end;

# ╔═╡ d8d01f36-fda3-11ea-2b26-59b86a4e23ca
md"## Result rethinking."

# ╔═╡ d8d1b594-fda3-11ea-2bae-e1be39ef751e
rethinking = "
        mean   sd  5.5% 94.5%
  a     0.09 0.24 -0.28  0.47
  bN    0.16 0.24 -0.23  0.54
  sigma 1.00 0.16  0.74  1.26
";

# ╔═╡ d8d8a4e4-fda3-11ea-3bda-0b8c90e63d65
Particles(post5_5_drafts_df)

# ╔═╡ d8db84ca-fda3-11ea-35c2-1f671cea8a32
if success(rc5_5_drafts)
  p = plot(title="m5.5.drafts: a ~ Normal(0, 1), bN ~ Normal(0, 1)")
  x = -2:0.01:2
  for j in 1:100
    y = post5_5_drafts_df[j, :a] .+ post5_5_drafts_df[j, :bN]*x
    plot!(p, x, y, color=:lightgrey, leg=false)
  end
	plot(p)
end

# ╔═╡ d8e7b710-fda3-11ea-13e3-e7141d1245d8
md"## End of clip-05-28-34s.jl"

# ╔═╡ Cell order:
# ╟─b0b04904-fda3-11ea-2a13-efbb8da05490
# ╠═b78a501c-fda3-11ea-06bb-5175ada30398
# ╠═b78a9662-fda3-11ea-25be-871b1258e163
# ╟─b78b1934-fda3-11ea-119c-5398009c1521
# ╠═b7997a6a-fda3-11ea-325a-7d0fc706cdff
# ╠═9fe858dc-fda3-11ea-38c8-df4b1d6e8419
# ╟─d65e8eea-fda3-11ea-0f78-6d8a183f5818
# ╠═d65ec40a-fda3-11ea-2e15-2f55dd308dd9
# ╠═407845b4-fd4a-11ea-1ed5-65581b68ccda
# ╟─d8d01f36-fda3-11ea-2b26-59b86a4e23ca
# ╠═d8d1b594-fda3-11ea-2bae-e1be39ef751e
# ╠═d8d8a4e4-fda3-11ea-3bda-0b8c90e63d65
# ╠═d8db84ca-fda3-11ea-35c2-1f671cea8a32
# ╟─d8e7b710-fda3-11ea-13e3-e7141d1245d8


text:
include("TwoPartCost.jl")
include("ThreePartCost.jl")
include("StorageManagementCost.jl")
include("MarketBidCost.jl")
include("MultiStartCost.jl")
include("Area.jl")
include("LoadZone.jl")
include("Bus.jl")
include("Arc.jl")
include("Line.jl")
include("MonitoredLine.jl")
include("PhaseShiftingTransformer.jl")
include("TapTransformer.jl")
include("Transformer2W.jl")
include("HVDCLine.jl")
include("VSCDCLine.jl")
include("InterruptibleLoad.jl")
include("FixedAdmittance.jl")
include("PowerLoad.jl")
include("ExponentialLoad.jl")
include("SingleCageInductionMachine.jl")
include("SimplifiedSingleCageInductionMachine.jl")
include("HydroEnergyReservoir.jl")
include("HydroDispatch.jl")
include("HydroPumpedStorage.jl")
include("RenewableDispatch.jl")
include("RenewableFix.jl")
include("ThermalStandard.jl")
include("ThermalMultiStart.jl")
include("GenericBattery.jl")
include("BatteryEMS.jl")
include("StaticReserve.jl")
include("StaticReserveNonSpinning.jl")
include("StaticReserveGroup.jl")
include("ReserveDemandCurve.jl")
include("VariableReserve.jl")
include("VariableReserveNonSpinning.jl")
include("AGC.jl")
include("Transfer.jl")
include("AVRFixed.jl")
include("AVRSimple.jl")
include("SEXS.jl")
include("ESDC1A.jl")
include("ESDC2A.jl")
include("IEEET1.jl")
include("AVRTypeI.jl")
include("AVRTypeII.jl")
include("SCRX.jl")
include("ESAC1A.jl")
include("EXAC1A.jl")
include("EXAC1.jl")
include("EXAC2.jl")
include("ESAC6A.jl")
include("ESST1A.jl")
include("EXPIC1.jl")
include("ESST4B.jl")
include("EXST1.jl")
include("EX4VSA.jl")
include("BaseMachine.jl")
include("RoundRotorMachine.jl")
include("SalientPoleMachine.jl")
include("AndersonFouadMachine.jl")
include("FullMachine.jl")
include("MarconatoMachine.jl")
include("OneDOneQMachine.jl")
include("SimpleAFMachine.jl")
include("SimpleFullMachine.jl")
include("SimpleMarconatoMachine.jl")
include("PSSFixed.jl")
include("PSSSimple.jl")
include("IEEEST.jl")
include("STAB1.jl")
include("SingleMass.jl")
include("FiveMassShaft.jl")
include("TGFixed.jl")
include("GasTG.jl")
include("GeneralGovModel.jl")
include("SteamTurbineGov1.jl")
include("HydroTurbineGov.jl")
include("IEEETurbineGov1.jl")
include("TGTypeI.jl")
include("TGTypeII.jl")
include("AverageConverter.jl")
include("RenewableEnergyConverterTypeA.jl")
include("FixedDCSource.jl")
include("ZeroOrderBESS.jl")
include("LCLFilter.jl")
include("LCFilter.jl")
include("RLFilter.jl")
include("KauraPLL.jl")
include("ReducedOrderPLL.jl")
include("FixedFrequency.jl")
include("VirtualInertia.jl")
include("ActivePowerDroop.jl")
include("ActivePowerPI.jl")
include("ActiveRenewableControllerAB.jl")
include("ReactiveRenewableControllerAB.jl")
include("ReactivePowerDroop.jl")
include("ReactivePowerPI.jl")
include("VoltageModeControl.jl")
include("CurrentModeControl.jl")
include("RECurrentControlB.jl")
include("AggregateDistributedGenerationA.jl")
include("Source.jl")
include("PeriodicVariableSource.jl")
include("GenericDER.jl")

export get_A
export get_A1
export get_A2
export get_A3
export get_A4
export get_A5
export get_A6
export get_AT
export get_A_set
export get_Accel
export get_Ae
export get_At
export get_B
export get_B_shunt
export get_Be
export get_Brkpt
export get_C
export get_D
export get_DB_h
export get_DB_l
export get_D_12
export get_D_23
export get_D_34
export get_D_45
export get_D_T
export get_D_dn
export get_D_ex
export get_D_hp
export get_D_ip
export get_D_lp
export get_D_turb
export get_D_up
export get_Dm
export get_E_lim
export get_E_sat
export get_Efd_lim
export get_FES_lim
export get_FRT_pnts
export get_Freq_Flag
export get_Ftrip_Flag
export get_G
export get_Gen_Flag
export get_H
export get_H_ex
export get_H_hp
export get_H_ip
export get_H_lim
export get_H_lp
export get_I_lr
export get_I_max
export get_Iflim
export get_Io_lim
export get_Iq_lim
export get_Iqinj_lim
export get_Iqr_lim
export get_Iqr_lims
export get_K
export get_K0
export get_K1
export get_K2
export get_K3
export get_K4
export get_K5
export get_K6
export get_K7
export get_K8
export get_KT
export get_K_c
export get_K_d
export get_K_ex
export get_K_hp
export get_K_hv
export get_K_i
export get_K_ig
export get_K_im
export get_K_ip
export get_K_ir
export get_K_lp
export get_K_lr
export get_K_p
export get_K_pg
export get_K_pm
export get_K_pr
export get_K_qi
export get_K_qp
export get_K_qv
export get_K_turb
export get_K_vi
export get_K_vp
export get_K_ω
export get_Ka
export get_Kb
export get_Kc
export get_Kd
export get_Kd_gov
export get_Ke
export get_Kf
export get_Kg
export get_Kh
export get_Ki
export get_Ki_gov
export get_Ki_load
export get_Ki_mw
export get_Ki_p
export get_Ki_q
export get_Kig
export get_Kip
export get_Kiq
export get_Kl
export get_Kp
export get_Kp_gov
export get_Kp_load
export get_Kp_p
export get_Kp_q
export get_Kpg
export get_Kpp
export get_Kpq
export get_Ks
export get_Kt
export get_Kv
export get_L_1d
export get_L_1q
export get_L_ad
export get_L_aq
export get_L_d
export get_L_f1d
export get_L_ff
export get_L_q
export get_Ld_ref
export get_Ls_lim
export get_Lv_pnts
export get_Lvpl1
export get_Lvpl_sw
export get_Oel_lim
export get_PF_Flag
export get_PQ_Flag
export get_PSS_flags
export get_P_lim
export get_P_lim_inner
export get_P_ref
export get_PerOp_Flag
export get_Pf_Flag
export get_Pfa_ref
export get_Q_Flag
export get_Q_lim
export get_Q_lim_inner
export get_Q_ref
export get_Qref_Flag
export get_R
export get_R_1d
export get_R_1q
export get_R_c
export get_R_close
export get_R_f
export get_R_lim
export get_R_open
export get_R_r
export get_R_s
export get_R_source
export get_R_th
export get_Recon_Flag
export get_Ref_Flag
export get_Rp
export get_Rrpwr
export get_Rselect
export get_SOC_ini
export get_SOC_lim
export get_Se
export get_Spar
export get_T
export get_T1
export get_T1T3
export get_T2
export get_T2T4
export get_T3
export get_T4
export get_T5
export get_T6
export get_T7
export get_TFRT_pnts
export get_TVRT_pnts
export get_T_AA
export get_T_act
export get_T_eng
export get_T_fltr
export get_T_ft
export get_T_fv
export get_T_g
export get_T_iq
export get_T_p
export get_T_pord
export get_T_rate
export get_T_rv
export get_Ta
export get_Ta_2
export get_Ta_3
export get_Ta_4
export get_Ta_Tb
export get_Tb
export get_Tb1
export get_Tc
export get_Tc1
export get_Td0_p
export get_Td0_pp
export get_Td_gov
export get_Te
export get_Tf
export get_Tf_1
export get_Tf_2
export get_Tf_load
export get_Tg
export get_Th
export get_Tj
export get_Tk
export get_Tp
export get_Tpelec
export get_Tpord
export get_Tq0_p
export get_Tq0_pp
export get_Tr
export get_Trf
export get_Trv
export get_Ts
export get_Tsa
export get_Tsb
export get_Tv
export get_Tw
export get_U0
export get_UEL_flags
export get_U_c
export get_VB_max
export get_VC_Flag
export get_VELM
export get_VES_lim
export get_VFE_lim
export get_VH_max
export get_VRT_pnts
export get_VV_pnts
export get_V_Flag
export get_V_frz
export get_V_lim
export get_V_lr
export get_V_pss
export get_V_ref
export get_V_ref0
export get_Va_lim
export get_Vcl
export get_Vcu
export get_Vdc_ref
export get_Vdip_lim
export get_Vf
export get_Vi_lim
export get_Vm_lim
export get_Vo_lim
export get_Vpr
export get_Vr_lim
export get_Vrfrac
export get_Vtrip_Flag
export get_Wf_nl
export get_X_ad
export get_X_aq
export get_X_c
export get_X_lr
export get_X_ls
export get_X_m
export get_X_p
export get_X_rr
export get_X_source
export get_X_ss
export get_X_th
export get_Xd
export get_Xd_p
export get_Xd_pp
export get_Xl
export get_Xq
export get_Xq_p
export get_Xq_pp
export get_Y
export get_Zerox
export get_active_power
export get_active_power_coefficient
export get_active_power_flow
export get_active_power_limits
export get_active_power_limits_from
export get_active_power_limits_pump
export get_active_power_limits_to
export get_ancillary_services
export get_angle
export get_angle_limits
export get_arc
export get_area
export get_available
export get_b
export get_base_power
export get_base_voltage
export get_battery_resistance
export get_battery_voltage
export get_bias
export get_branch_id_control
export get_bus
export get_bus_control
export get_bustype
export get_cf
export get_contributing_services
export get_conversion_factor
export get_d
export get_dP_lim
export get_db
export get_dbd_pnts
export get_dc_dc_inductor
export get_dc_link_capacitance
export get_delta_t
export get_dynamic_injector
export get_e_lim
export get_efficiency
export get_energy_shortage_cost
export get_energy_surplus_cost
export get_eq_p
export get_ext
export get_f
export get_fdbd_pnts
export get_fe_lim
export get_fh
export get_fixed
export get_fl
export get_flow_limits
export get_frequency
export get_from
export get_from_branch_control
export get_fs
export get_fuel
export get_fuel_flag
export get_gate_position_limits
export get_inflow
export get_initial_ace
export get_initial_energy
export get_initial_storage
export get_input_active_power_limits
export get_input_code
export get_internal_angle
export get_internal_angle_bias
export get_internal_angle_coefficients
export get_internal_angle_frequencies
export get_internal_voltage
export get_internal_voltage_bias
export get_internal_voltage_coefficients
export get_internal_voltage_frequencies
export get_inv_d_fluxlink
export get_inv_q_fluxlink
export get_inverter_firing_angle
export get_inverter_tap_limits
export get_inverter_xrc
export get_kWh_Cap
export get_kad
export get_kd
export get_kffi
export get_kffv
export get_ki_pll
export get_kic
export get_kii
export get_kiv
export get_kp_pll
export get_kpc
export get_kpi
export get_kpv
export get_kq
export get_kω
export get_lf
export get_lg
export get_load_response
export get_load_zone
export get_loss
export get_lv
export get_magnitude
export get_max_active_power
export get_max_reactive_power
export get_model
export get_must_run
export get_n_states
export get_name
export get_no_load
export get_number
export get_operation_cost
export get_outflow
export get_output_active_power_limits
export get_peak_active_power
export get_peak_reactive_power
export get_power_factor
export get_power_trajectory
export get_primary_shunt
export get_prime_mover
export get_pump_efficiency
export get_q_nl
export get_r
export get_ramp_limits
export get_ramp_limits_pump
export get_rate
export get_rated_current
export get_rated_voltage
export get_rating
export get_rating_pump
export get_rc_rfd
export get_reactive_power
export get_reactive_power_coefficient
export get_reactive_power_flow
export get_reactive_power_limits
export get_reactive_power_limits_from
export get_reactive_power_limits_pump
export get_reactive_power_limits_to
export get_rectifier_firing_angle
export get_rectifier_tap_limits
export get_rectifier_xrc
export get_remote_bus_control
export get_requirement
export get_rf
export get_rg
export get_rrpwr
export get_rv
export get_saturation_coeffs
export get_services
export get_shut_down
export get_speed_error_signal
export get_start_time_limits
export get_start_types
export get_start_up
export get_state_of_charge_limits
export get_states
export get_states_types
export get_status
export get_storage_capacity
export get_storage_target
export get_switch
export get_tF_delay
export get_tV_delay
export get_tap
export get_tfh
export get_tfl
export get_time_at_status
export get_time_frame
export get_time_limits
export get_time_limits_pump
export get_time_series_container
export get_to
export get_to_branch_control
export get_valve_position_limits
export get_variable
export get_vh_pnts
export get_vl_pnts
export get_voltage
export get_voltage_limits
export get_x
export get_α
export get_γ_d1
export get_γ_d2
export get_γ_q1
export get_γ_q2
export get_γ_qd
export get_γd
export get_γq
export get_θp
export get_θp_rad
export get_τ_limits
export get_τ_ref
export get_ω_lp
export get_ωad
export get_ωf
export get_ωz
export set_A!
export set_A1!
export set_A2!
export set_A3!
export set_A4!
export set_A5!
export set_A6!
export set_AT!
export set_A_set!
export set_Accel!
export set_Ae!
export set_At!
export set_B!
export set_B_shunt!
export set_Be!
export set_Brkpt!
export set_C!
export set_D!
export set_DB_h!
export set_DB_l!
export set_D_12!
export set_D_23!
export set_D_34!
export set_D_45!
export set_D_T!
export set_D_dn!
export set_D_ex!
export set_D_hp!
export set_D_ip!
export set_D_lp!
export set_D_turb!
export set_D_up!
export set_Dm!
export set_E_lim!
export set_E_sat!
export set_Efd_lim!
export set_FES_lim!
export set_FRT_pnts!
export set_Freq_Flag!
export set_Ftrip_Flag!
export set_G!
export set_Gen_Flag!
export set_H!
export set_H_ex!
export set_H_hp!
export set_H_ip!
export set_H_lim!
export set_H_lp!
export set_I_lr!
export set_I_max!
export set_Iflim!
export set_Io_lim!
export set_Iq_lim!
export set_Iqinj_lim!
export set_Iqr_lim!
export set_Iqr_lims!
export set_K!
export set_K0!
export set_K1!
export set_K2!
export set_K3!
export set_K4!
export set_K5!
export set_K6!
export set_K7!
export set_K8!
export set_KT!
export set_K_c!
export set_K_d!
export set_K_ex!
export set_K_hp!
export set_K_hv!
export set_K_i!
export set_K_ig!
export set_K_im!
export set_K_ip!
export set_K_ir!
export set_K_lp!
export set_K_lr!
export set_K_p!
export set_K_pg!
export set_K_pm!
export set_K_pr!
export set_K_qi!
export set_K_qp!
export set_K_qv!
export set_K_turb!
export set_K_vi!
export set_K_vp!
export set_K_ω!
export set_Ka!
export set_Kb!
export set_Kc!
export set_Kd!
export set_Kd_gov!
export set_Ke!
export set_Kf!
export set_Kg!
export set_Kh!
export set_Ki!
export set_Ki_gov!
export set_Ki_load!
export set_Ki_mw!
export set_Ki_p!
export set_Ki_q!
export set_Kig!
export set_Kip!
export set_Kiq!
export set_Kl!
export set_Kp!
export set_Kp_gov!
export set_Kp_load!
export set_Kp_p!
export set_Kp_q!
export set_Kpg!
export set_Kpp!
export set_Kpq!
export set_Ks!
export set_Kt!
export set_Kv!
export set_L_1d!
export set_L_1q!
export set_L_ad!
export set_L_aq!
export set_L_d!
export set_L_f1d!
export set_L_ff!
export set_L_q!
export set_Ld_ref!
export set_Ls_lim!
export set_Lv_pnts!
export set_Lvpl1!
export set_Lvpl_sw!
export set_Oel_lim!
export set_PF_Flag!
export set_PQ_Flag!
export set_PSS_flags!
export set_P_lim!
export set_P_lim_inner!
export set_P_ref!
export set_PerOp_Flag!
export set_Pf_Flag!
export set_Pfa_ref!
export set_Q_Flag!
export set_Q_lim!
export set_Q_lim_inner!
export set_Q_ref!
export set_Qref_Flag!
export set_R!
export set_R_1d!
export set_R_1q!
export set_R_c!
export set_R_close!
export set_R_f!
export set_R_lim!
export set_R_open!
export set_R_r!
export set_R_s!
export set_R_source!
export set_R_th!
export set_Recon_Flag!
export set_Ref_Flag!
export set_Rp!
export set_Rrpwr!
export set_Rselect!
export set_SOC_ini!
export set_SOC_lim!
export set_Se!
export set_Spar!
export set_T!
export set_T1!
export set_T1T3!
export set_T2!
export set_T2T4!
export set_T3!
export set_T4!
export set_T5!
export set_T6!
export set_T7!
export set_TFRT_pnts!
export set_TVRT_pnts!
export set_T_AA!
export set_T_act!
export set_T_eng!
export set_T_fltr!
export set_T_ft!
export set_T_fv!
export set_T_g!
export set_T_iq!
export set_T_p!
export set_T_pord!
export set_T_rate!
export set_T_rv!
export set_Ta!
export set_Ta_2!
export set_Ta_3!
export set_Ta_4!
export set_Ta_Tb!
export set_Tb!
export set_Tb1!
export set_Tc!
export set_Tc1!
export set_Td0_p!
export set_Td0_pp!
export set_Td_gov!
export set_Te!
export set_Tf!
export set_Tf_1!
export set_Tf_2!
export set_Tf_load!
export set_Tg!
export set_Th!
export set_Tj!
export set_Tk!
export set_Tp!
export set_Tpelec!
export set_Tpord!
export set_Tq0_p!
export set_Tq0_pp!
export set_Tr!
export set_Trf!
export set_Trv!
export set_Ts!
export set_Tsa!
export set_Tsb!
export set_Tv!
export set_Tw!
export set_U0!
export set_UEL_flags!
export set_U_c!
export set_VB_max!
export set_VC_Flag!
export set_VELM!
export set_VES_lim!
export set_VFE_lim!
export set_VH_max!
export set_VRT_pnts!
export set_VV_pnts!
export set_V_Flag!
export set_V_frz!
export set_V_lim!
export set_V_lr!
export set_V_pss!
export set_V_ref!
export set_V_ref0!
export set_Va_lim!
export set_Vcl!
export set_Vcu!
export set_Vdc_ref!
export set_Vdip_lim!
export set_Vf!
export set_Vi_lim!
export set_Vm_lim!
export set_Vo_lim!
export set_Vpr!
export set_Vr_lim!
export set_Vrfrac!
export set_Vtrip_Flag!
export set_Wf_nl!
export set_X_ad!
export set_X_aq!
export set_X_c!
export set_X_lr!
export set_X_ls!
export set_X_m!
export set_X_p!
export set_X_rr!
export set_X_source!
export set_X_ss!
export set_X_th!
export set_Xd!
export set_Xd_p!
export set_Xd_pp!
export set_Xl!
export set_Xq!
export set_Xq_p!
export set_Xq_pp!
export set_Y!
export set_Zerox!
export set_active_power!
export set_active_power_coefficient!
export set_active_power_flow!
export set_active_power_limits!
export set_active_power_limits_from!
export set_active_power_limits_pump!
export set_active_power_limits_to!
export set_ancillary_services!
export set_angle!
export set_angle_limits!
export set_arc!
export set_area!
export set_available!
export set_b!
export set_base_power!
export set_base_voltage!
export set_battery_resistance!
export set_battery_voltage!
export set_bias!
export set_branch_id_control!
export set_bus!
export set_bus_control!
export set_bustype!
export set_cf!
export set_contributing_services!
export set_conversion_factor!
export set_d!
export set_dP_lim!
export set_db!
export set_dbd_pnts!
export set_dc_dc_inductor!
export set_dc_link_capacitance!
export set_delta_t!
export set_dynamic_injector!
export set_e_lim!
export set_efficiency!
export set_energy_shortage_cost!
export set_energy_surplus_cost!
export set_eq_p!
export set_ext!
export set_f!
export set_fdbd_pnts!
export set_fe_lim!
export set_fh!
export set_fixed!
export set_fl!
export set_flow_limits!
export set_frequency!
export set_from!
export set_from_branch_control!
export set_fs!
export set_fuel!
export set_fuel_flag!
export set_gate_position_limits!
export set_inflow!
export set_initial_ace!
export set_initial_energy!
export set_initial_storage!
export set_input_active_power_limits!
export set_input_code!
export set_internal_angle!
export set_internal_angle_bias!
export set_internal_angle_coefficients!
export set_internal_angle_frequencies!
export set_internal_voltage!
export set_internal_voltage_bias!
export set_internal_voltage_coefficients!
export set_internal_voltage_frequencies!
export set_inv_d_fluxlink!
export set_inv_q_fluxlink!
export set_inverter_firing_angle!
export set_inverter_tap_limits!
export set_inverter_xrc!
export set_kWh_Cap!
export set_kad!
export set_kd!
export set_kffi!
export set_kffv!
export set_ki_pll!
export set_kic!
export set_kii!
export set_kiv!
export set_kp_pll!
export set_kpc!
export set_kpi!
export set_kpv!
export set_kq!
export set_kω!
export set_lf!
export set_lg!
export set_load_response!
export set_load_zone!
export set_loss!
export set_lv!
export set_magnitude!
export set_max_active_power!
export set_max_reactive_power!
export set_model!
export set_must_run!
export set_n_states!
export set_name!
export set_no_load!
export set_number!
export set_operation_cost!
export set_outflow!
export set_output_active_power_limits!
export set_peak_active_power!
export set_peak_reactive_power!
export set_power_factor!
export set_power_trajectory!
export set_primary_shunt!
export set_prime_mover!
export set_pump_efficiency!
export set_q_nl!
export set_r!
export set_ramp_limits!
export set_ramp_limits_pump!
export set_rate!
export set_rated_current!
export set_rated_voltage!
export set_rating!
export set_rating_pump!
export set_rc_rfd!
export set_reactive_power!
export set_reactive_power_coefficient!
export set_reactive_power_flow!
export set_reactive_power_limits!
export set_reactive_power_limits_from!
export set_reactive_power_limits_pump!
export set_reactive_power_limits_to!
export set_rectifier_firing_angle!
export set_rectifier_tap_limits!
export set_rectifier_xrc!
export set_remote_bus_control!
export set_requirement!
export set_rf!
export set_rg!
export set_rrpwr!
export set_rv!
export set_saturation_coeffs!
export set_services!
export set_shut_down!
export set_speed_error_signal!
export set_start_time_limits!
export set_start_types!
export set_start_up!
export set_state_of_charge_limits!
export set_states!
export set_states_types!
export set_status!
export set_storage_capacity!
export set_storage_target!
export set_switch!
export set_tF_delay!
export set_tV_delay!
export set_tap!
export set_tfh!
export set_tfl!
export set_time_at_status!
export set_time_frame!
export set_time_limits!
export set_time_limits_pump!
export set_time_series_container!
export set_to!
export set_to_branch_control!
export set_valve_position_limits!
export set_variable!
export set_vh_pnts!
export set_vl_pnts!
export set_voltage!
export set_voltage_limits!
export set_x!
export set_α!
export set_γ_d1!
export set_γ_d2!
export set_γ_q1!
export set_γ_q2!
export set_γ_qd!
export set_γd!
export set_γq!
export set_θp!
export set_θp_rad!
export set_τ_limits!
export set_τ_ref!
export set_ω_lp!
export set_ωad!
export set_ωf!
export set_ωz!


text:
using EngEconomics, Roots, Plots

# Given
aInit = -200000
aAnnual = 20000
aAnnualG = -5000
aLife = 5

MARR = 0.12

# Find IRR
IRR_eq(x) = aInit - 40000 * seriesPresentAmountFactor(x, 6) +  40000 / x + (20000 - 5000 * sinkingFundPaymentArithmetic(x, 5)) * seriesPresentAmountFactor(x, 5)
xVec = collect(0:0.01:1)
yVec = zeros(size(xVec))
plot(xVec, IRR_eq.(xVec))
plot!(xVec, yVec)

find_zero(IRR_eq, MARR)


text:
# # Measuring overheads of some threading API

# You can find the code for benchmarking and plotting at
# <https://github.com/tkf/ThreadsAPIBenchmarks.jl>.

using DisplayAs
using ThreadsAPIBenchmarksAnalyzers: TaskOverheads, NotifyOverheads

# ## Overhead of `@spawn` and `@sync`
#
# Measuring the overhead of `@spawn` and `@sync` is not entirely trivial.
# `@btime wait(Threads.@spawn nothing)` does not reflect the actual overhead
# since `@spawn` may use the current thread to run the task. That is to say,
# the following check likely to pass even if `nthreads() > 1`:

let n = 0
    for _ in 1:1000
        n += fetch(Threads.@spawn Threads.threadid()) == Threads.threadid()
    end
    @assert n > 0
end

# ThreadsAPIBenchmarks.jl contains a benchmark for measuring the overheads of
# `@sync` and `@spawn`. That is to say, it measures `@spawn` overhead `t1 - t0`
# and `wait` (sync) overhead `t2 - t1` of
#
# ```julia
# t0 = time_ns()
# wait(@spawn t1 = time_ns())
# t2 = time_ns()
# ```
#
# while making sure `@spawn`ed task is executed in a different OS threads.
#
# The distribution of the above measurement is shown below for different number
# `N = nthreads()` of threads:

begin
    task_overheads = TaskOverheads.analyze()
    task_overheads.hists
end |> DisplayAs.PNG
#-

# Here is a more direct plot of the raw data to show that the multimodality of
# the distribution comes from the temporal proximity of the measurements (some
# state transitions in Julia runtime, OS, or hardware?):

begin
    ## task_overheads = TaskOverheads.analyze()
    task_overheads.seq
end |> DisplayAs.PNG

# ## Overhead of `notify` (`schedule`)
#
# Following distribution describes the overheads of `wait` and `notify` on
# `Threads.Condition`. It is essentially the overhead of re-`schedule`.
#
# The first two rows (`max` and `min`) shows the distribution of the time took
# for `notify` to end the `wait` of the peer threads. The first row `max` shows
# the slowest among all threads and second row `min` is for the fastest. The
# last row `notify` is `@time notify(...)`.

begin
    notify_overheads = NotifyOverheads.analyze()
    notify_overheads.hists
end
#-

# Compared to `@sync` and `@spawn`, there are not much temporal irregularities
# especially for `N > 2`:

begin
    ## notify_overheads = NotifyOverheads.analyze()
    notify_overheads.seq
end |> DisplayAs.PNG


text:
"""
    Code for LineList types used by CCfs
Author: Eric Ford
Created: August 2020
"""

"""A struct implementing a line list should be a subtype of AbstractLineList. """
abstract type AbstractLineList end

""" A basic line list for passing to compute CCFs.
Contains (views into) arrays specifying the minimum and maximum wavelength range and weight for each line. """
struct BasicLineList{T<:Real, AA<:AbstractArray{T,1} } <: AbstractLineList
    λ::AA
    weight::AA
end

""" `BasicLineList( λ, weight )` """
function BasicLineList{T}(λ::AA, w::AA) where { T<:Real, AA<:AbstractArray{T,1} }
    @assert length(λ) == length(w)
    @assert length(λ) >= 1
    @assert 0.0 .<= w .<= 1.0
    BasicLineList{eltype(w),typeof(w)}(λ,w)
end

function BasicLineList( df::DataFrame )
    @assert hasproperty(df,:lambda)
    @assert hasproperty(df,:weight)
    BasicLineList{eltype(df[:,:lambda]),typeof(df[:,:weight])}(df[:,:lambda],df[:,:weight])
end

function EmptyBasicLineList()
    return BasicLineList{Float64,Array{Float64,1}}(zeros(0),zeros(0))
end

""" A basic line list for passing to compute CCFs.
Contains (views into) arrays specifying the minimum and maximum wavelength range and weight for each line. """
struct BasicLineList2D{T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Integer,  AA2<:AbstractArray{T2,1} } <: AbstractLineList
    λ::AA1
    weight::AA1
    order::AA2
end

""" `BasicLineList( λ, weight )` """
function BasicLineList2D{T1,T2}(λ::AA1, w::AA1, o::AA2) where {T1<:Real, AA1<:AbstractArray{T1,1}, T2<:Integer,  AA2<:AbstractArray{T2,1} }
    @assert length(λ) == length(w) == length(o)
    @assert length(λ) >= 1
    @assert all(0.0 .<= w .<= 1.0)
    @assert all(1 .<= o .<= 200)           # arbitrary limit for maximum number of orders
    BasicLineList2D{eltype(w),typeof(w),eltype{o},typeof(o)}(λ,w,o)
end

function BasicLineList2D( df::DataFrame )
    @assert hasproperty(df,:lambda)
    @assert hasproperty(df,:weight)
    @assert hasproperty(df,:order)
    @assert all(0.0 .<= df[!,:weight] .<= 1.0)
    @assert all(1 .<= df[!,:order] .<= 200)           # arbitrary limit for maximum number of orders
    BasicLineList2D{eltype(df[:,:lambda]),typeof(df[:,:weight]),eltype(df[:,:order]),typeof(df[:,:order])}(df[!,:lambda],df[!,:weight],df[!,:order])
end

function EmptyBasicLineList2D()
    return BasicLineList2D{Float64,Array{Float64,1},Int32,Array{Int32,1}}(zeros(0),zeros(0), zeros(Int32,0) )
end

#=  Not fully implemented/tested yet
""" A line list for passing to compute CCFs with variable line widths.
Contains (views into) arrays specifying the minimum and maximum wavelength range and weight for each line. """
struct VarWidthLineList{T<:Real, AA<:AbstractArray{T,1} } <: AbstractLineList
    λ_lo::AA
    λ_hi::AA
    weight::AA
end

""" VarWidthLineList( λ_lo, λ_hi, weight ) """
function VarWidthLineList{T}(lo::AA, hi::AA, w::AA) where { T<:Real, AA<:AbstractArray{T,1} }
    @assert length(lo) == length(hi) == length(w)
    @assert length(lo) >= 1
    @assert 0.0 .<= λ_hi.-λ_lo .<= 2.0
    @assert 0.0 .<= w .<= 1.0
    VarWidthLineList{eltype(w),typeof(w)}(lo,hi,w)
end
=#

import Base.length
""" Return length of line linst. """
length(ll::AbstractLineList) = length(ll.weight)

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
module CachePath

using Base: PkgId, root_module, root_module_exists, PkgOrigin,
    pkgorigins, package_locks, toplevel_load,
    locate_package, JLOptions, PrecompilableError, @logmsg,
    cache_file_entry, isfile_casesensitive, stale_cachefile,
    _include_from_serialized, _concrete_dependencies,
    loaded_modules, module_build_id, CoreLogging, create_expr_cache, _crc32c,
    preferences_hash, slug, MAX_NUM_PRECOMPILE_FILES, rename,
    _tryrequire_from_serialized, _require_from_serialized,
    package_callbacks
#, @constprop

export @cpimport

const CACHE_PATH_DEBUG = false

function cachepathdebug(args...)
    if CACHE_PATH_DEBUG
        println(args...)
    end
end

const CACHE_PATH_COMPAT =
if VERSION < v"1.7.999"
    true
else
    false
end

# This could cause a problem, as there will be two locks when running v1.8
# _require_lock and require_lock
# Only _require_lock is used in this file
const _require_lock = ReentrantLock()

if CACHE_PATH_COMPAT
    _Condition() = Condition()
else
    _Condition() = Threads.Condition(_require_lock)
end

# run_package_callbacks not present in v1.6.5
function _run_package_callbacks(modkey::PkgId)
    unlock(_require_lock)
    try
        for callback in package_callbacks
            Base.invokelatest(callback, modkey)
        end
    catch
        # Try to continue loading if a callback errors
        errs = current_exceptions()
        @error "Error during package callback" exception=errs
    finally
        lock(_require_lock)
    end
    nothing
end

# Macro @lock not present in v1.6.5
macro _lock(l, expr)
    quote
        temp = $(esc(l))
        lock(temp)
        try
            $(esc(expr))
        finally
            unlock(temp)
        end
    end
end

"""
    @cpimport module depot_path::AbstractString

Import `module` using `depot_path` to store and retrieve the precompile
cache. Precompile caches existing elsewhere are ignored.

The semantics of `@cpimport` probably differ from those of `import` in
other ways.

# Examples:
Import the module `Example` using "./newdepot" for storing and retrieving
the precompile cache.

    julia> @cpimport Example "./newdepot"
"""
macro cpimport(mod, depot_path)
    qmod = QuoteNode(mod)
    :(const $(esc(mod)) = CachePath.require(Main, $(esc(qmod)), $depot_path); nothing)
end

"""
    CachePath.require(into::Module, module::Symbol, depot_path::AbstractString)

This function is modified from `Base.require`. The precompiled cache (`.ji` file) for
`module` will only be searched for and stored in `depot_path`.

Loads a source file, in the context of the `Main` module, on every active node, searching
standard locations for files. `require` is considered a top-level operation, so it sets the
current `include` path but does not use it to search for files (see help for [`include`](@ref)).
This function is typically used to load library code, and is implicitly called by `using` to
load packages.

When searching for files, `require` first looks for package code in the global array
[`LOAD_PATH`](@ref). `require` is case-sensitive on all platforms, including those with
case-insensitive filesystems like macOS and Windows.

For more details regarding code loading, see the manual sections on [modules](@ref modules) and
[parallel computing](@ref code-availability).
"""
function require(into::Module, mod::Symbol, depot_path::AbstractString)
    @_lock _require_lock begin
    # LOADING_CACHE[] = LoadingCache() # Does not exist in v1.6.5
    try
        uuidkey = Base.identify_package(into, String(mod))
        # Core.println("require($(PkgId(into)), $mod) -> $uuidkey")
        if uuidkey === nothing
            where = PkgId(into)
            if where.uuid === nothing
                hint, dots = begin
                    if isdefined(into, mod) && getfield(into, mod) isa Module
                        true, "."
                    elseif isdefined(parentmodule(into), mod) && getfield(parentmodule(into), mod) isa Module
                        true, ".."
                    else
                        false, ""
                    end
                end
                hint_message = hint ? ", maybe you meant `import/using $(dots)$(mod)`" : ""
                start_sentence = hint ? "Otherwise, run" : "Run"
                throw(ArgumentError("""
                    Package $mod not found in current path$hint_message.
                    - $start_sentence `import Pkg; Pkg.add($(repr(String(mod))))` to install the $mod package."""))
            else
                throw(ArgumentError("""
                Package $(where.name) does not have $mod in its dependencies:
                - You may have a partially installed environment. Try `Pkg.instantiate()`
                  to ensure all packages in the environment are installed.
                - Or, if you have $(where.name) checked out for development and have
                  added $mod as a dependency but haven't updated your primary
                  environment's manifest file, try `Pkg.resolve()`.
                - Otherwise you may need to report an issue with $(where.name)"""))
            end
        end
        if Base._track_dependencies[]
            push!(Base._require_dependencies, (into, binpack(uuidkey), 0.0))
        end
        return _require_prelocked(uuidkey, depot_path)
    finally
#        LOADING_CACHE[] = nothing
    end
    end
end

require(uuidkey::PkgId, depot_path::AbstractString) =
    @_lock _require_lock _require_prelocked(uuidkey, depot_path)

"""
    Base.require(package::AbstractString, depot_path::AbstractString)
    Base.require(package::PkgId, depot_path::AbstractString)

Load `package` and store the cached precompile file in the depot specified
by `depot_path`. If it does not exist, `depot_path` is created. If the
cached precompiled file is found in `depot_path`, then it is loaded.

# Examples
```julia
const Example = Base.require("Example", "./a_depot")
```
"""
require(package::AbstractString, depot_path::AbstractString) =
    require(Base.identify_package(package), depot_path)

function _require_prelocked(uuidkey::PkgId, depot_path::AbstractString)
    just_loaded_pkg = false
    if !root_module_exists(uuidkey)
        cachefile = _require(uuidkey, depot_path)
        if cachefile !== nothing
            get!(PkgOrigin, pkgorigins, uuidkey).cachepath = cachefile
        end
        # After successfully loading, notify downstream consumers
        _run_package_callbacks(uuidkey)
        just_loaded_pkg = true
    end
    if just_loaded_pkg && !root_module_exists(uuidkey)
        # error("package `$(uuidkey.name)` did not define the expected \
        #       module `$(uuidkey.name)`, check for typos in package module name")
        error("package `$(uuidkey.name)` did not define the expected module `$(uuidkey.name)`, check for typos in package module name")
    end
    return root_module(uuidkey)
end

# Returns `nothing` or the name of the newly-created cachefile
function _require(pkg::PkgId, depot_path::AbstractString)
    # handle recursive calls to require
    loading = get(package_locks, pkg, false)
    if loading !== false
        # load already in progress for this module
        wait(loading)
        return
    end
    package_locks[pkg] = _Condition()

    last = toplevel_load[]
    try
        toplevel_load[] = false
        # perform the search operation to select the module file require intends to load
        path = locate_package(pkg)
        get!(PkgOrigin, pkgorigins, pkg).path = path
        if path === nothing
            throw(ArgumentError("""
                Package $pkg is required but does not seem to be installed:
                 - Run `Pkg.instantiate()` to install all recorded dependencies.
                """))
        end

        # attempt to load the module file via the precompile cache locations
        if JLOptions().use_compiled_modules != 0
            m = _require_search_from_serialized(pkg, depot_path, path)
            if !isa(m, Bool)
                return
            end
        end

        # if the module being required was supposed to have a particular version
        # but it was not handled by the precompile loader, complain
        for (concrete_pkg, concrete_build_id) in _concrete_dependencies
            if pkg == concrete_pkg
                @warn """Module $(pkg.name) with build ID $concrete_build_id is missing from the cache.
                     This may mean $pkg does not support precompilation but is imported by a module that does."""
                if JLOptions().incremental != 0
                    # during incremental precompilation, this should be fail-fast
                    throw(PrecompilableError())
                end
            end
        end

        if JLOptions().use_compiled_modules != 0
            if (0 == ccall(:jl_generating_output, Cint, ())) || (JLOptions().incremental != 0)
                # spawn off a new incremental pre-compile task for recursive `require` calls
                # or if the require search declared it was pre-compiled before (and therefore is expected to still be pre-compilable)
                cachefile = compilecache(pkg, path, depot_path)
                if isa(cachefile, Exception)
                    if precompilableerror(cachefile)
                        verbosity = isinteractive() ? CoreLogging.Info : CoreLogging.Debug
                        @logmsg verbosity "Skipping precompilation since __precompile__(false). Importing $pkg."
                    else
                        @warn "The call to compilecache failed to create a usable precompiled cache file for $pkg" exception=m
                    end
                    # fall-through to loading the file locally
                else
                    m = _require_from_serialized(cachefile)
                    if isa(m, Exception)
                        @warn "The call to compilecache failed to create a usable precompiled cache file for $pkg" exception=m
                    else
                        return cachefile
                    end
                end
            end
        end

        # just load the file normally via include
        # for unknown dependencies
        uuid = pkg.uuid
        uuid = (uuid === nothing ? (UInt64(0), UInt64(0)) : convert(NTuple{2, UInt64}, uuid))
        old_uuid = ccall(:jl_module_uuid, NTuple{2, UInt64}, (Any,), __toplevel__)
        if uuid !== old_uuid
            ccall(:jl_set_module_uuid, Cvoid, (Any, NTuple{2, UInt64}), __toplevel__, uuid)
        end
        if ! CACHE_PATH_COMPAT
            unlock(_require_lock)
        end
        try
            include(__toplevel__, path)
            return
        finally
            if ! CACHE_PATH_COMPAT
                lock(_require_lock)
            end
            if uuid !== old_uuid
                ccall(:jl_set_module_uuid, Cvoid, (Any, NTuple{2, UInt64}), __toplevel__, old_uuid)
            end
        end
    finally
        toplevel_load[] = last
        loading = pop!(package_locks, pkg)
        notify(loading, all=true)
    end
    nothing
end

const _TIMING_IMPORTS =
    if CACHE_PATH_COMPAT
        Base.Threads.Atomic{Int}(0)
    else
        Base.TIMING_IMPORTS
    end

macro myconstprop(setting, ex)
    if isa(setting, QuoteNode)
        setting = setting.value
    end
    setting === :aggressive && return esc(isa(ex, Expr) ? Base.pushmeta!(ex, :aggressive_constprop) : ex)
    setting === :none && return esc(isa(ex, Expr) ? Base.pushmeta!(ex, :no_constprop) : ex)
    throw(ArgumentError("@myconstprop $setting not supported"))
end

# returns `true` if require found a precompile cache for this sourcepath, but couldn't load it
# returns `false` if the module isn't known to be precompilable
# returns the set of modules restored if the cache load succeeded
@myconstprop :none function _require_search_from_serialized(pkg::PkgId, depot_path::String, sourcepath::String, depth::Int = 0)
    t_before = time_ns()
    paths = find_all_in_cache_path(pkg, depot_path)
    for path_to_try in paths::Vector{String}
        cachepathdebug("path_to_try ", path_to_try)
        cachepathdebug("sourcepath ", sourcepath)
        staledeps = stale_cachefile(sourcepath, path_to_try)
        if staledeps === true
            continue
        end
        try
            touch(path_to_try) # update timestamp of precompilation file
        catch # file might be read-only and then we fail to update timestamp, which is fine
        end
        # finish loading module graph into staledeps
        for i in 1:length(staledeps)
            dep = staledeps[i]
            dep isa Module && continue
            modpath, modkey, build_id = dep::Tuple{String, PkgId, UInt64}
            dep = _tryrequire_from_serialized(modkey, build_id, modpath, depth + 1)
            if dep === nothing
                @debug "Required dependency $modkey failed to load from cache file for $modpath."
                staledeps = true
                break
            end
            staledeps[i] = dep::Module
        end
        if staledeps === true
            continue
        end
        cachepathdebug("_include_from_serialized($path_to_try, staledeps)")
        restored = _include_from_serialized(path_to_try, staledeps)
        if isa(restored, Exception)
            @debug "Deserialization checks failed while attempting to load cache from $path_to_try" exception=restored
        else
            if _TIMING_IMPORTS[] > 0
                elapsed = round((time_ns() - t_before) / 1e6, digits = 1)
                tree_prefix = depth == 0 ? "" : "  "^(depth-1)*"┌ "
                print(lpad(elapsed, 9), " ms  ")
                printstyled(tree_prefix, color = :light_black)
                println(pkg.name)
            end
            cachepathdebug("Restored ", restored)
            return restored
        end
    end
    return !isempty(paths)
end

function find_all_in_cache_path(pkg::PkgId, depot_path::String)
    paths = String[]
    entrypath, entryfile = cache_file_entry(pkg)
    # Use the second line to require that the cache file is in depot_path
    # The first line will search in other depots if cache not found in depot_path
#   for path in (joinpath(depot_path, entrypath), joinpath.(DEPOT_PATH, entrypath)...)
    for path in (joinpath(depot_path, entrypath), )
        isdir(path) || continue
        for file in readdir(path, sort = false) # no sort given we sort later
            if !((pkg.uuid === nothing && file == entryfile * ".ji") ||
                 (pkg.uuid !== nothing && startswith(file, entryfile * "_")))
                 continue
            end
            filepath = joinpath(path, file)
            isfile_casesensitive(filepath) && push!(paths, filepath)
        end
    end
    if length(paths) > 1
        # allocating the sort vector is less expensive than using sort!(.. by=mtime), which would
        # call the relatively slow mtime multiple times per path
        p = sortperm(mtime.(paths), rev = true)
        return paths[p]
    else
        return paths
    end
end

function compilecache_dir(pkg::PkgId, depot_path::String)
    entrypath, entryfile = cache_file_entry(pkg)
    cachepathdebug("compilecache_dir($pkg ::PkgId, $depot_path ::String)")
    cachepathdebug("joinpath($depot_path, $entrypath)")
    return joinpath(depot_path, entrypath)
end

function compilecache(pkg::PkgId, path::String, depot_path::String, internal_stderr::IO = stderr, internal_stdout::IO = stdout,
                      ignore_loaded_modules::Bool = true)

    @nospecialize internal_stderr internal_stdout
    # decide where to put the resulting cache file
    cachepathdebug("compilecache: depot_path $depot_path")
    cachepath = compilecache_dir(pkg, depot_path)

    # build up the list of modules that we want the precompile process to preserve
    concrete_deps = copy(_concrete_dependencies)
    if ignore_loaded_modules
        for (key, mod) in loaded_modules
            if !(mod === Main || mod === Core || mod === Base)
                push!(concrete_deps, key => module_build_id(mod))
            end
        end
    end
    # run the expression and cache the result
    verbosity = isinteractive() ? CoreLogging.Info : CoreLogging.Debug
    @logmsg verbosity "Precompiling $pkg"

    # create a temporary file in `cachepath` directory, write the cache in it,
    # write the checksum, _and then_ atomically move the file to `cachefile`.
    mkpath(cachepath)
    tmppath, tmpio = mktemp(cachepath)
    local p
    try
        close(tmpio)
        p = create_expr_cache(pkg, path, tmppath, concrete_deps, internal_stderr, internal_stdout)
        if success(p)
            # append checksum to the end of the .ji file:
            open(tmppath, "a+") do f
                write(f, _crc32c(seekstart(f)))
            end
            # inherit permission from the source file (and make them writable)
            chmod(tmppath, filemode(path) & 0o777 | 0o200)

            # Read preferences hash back from .ji file (we can't precompute because
            # we don't actually know what the list of compile-time preferences are without compiling)
            prefs_hash = preferences_hash(tmppath)
            cachefile = compilecache_path(pkg, depot_path, prefs_hash)

            # prune the directory with cache files
            if pkg.uuid !== nothing
                entrypath, entryfile = cache_file_entry(pkg)
                cachefiles = filter!(x -> startswith(x, entryfile * "_"), readdir(cachepath))
                if length(cachefiles) >= MAX_NUM_PRECOMPILE_FILES[]
                    idx = findmin(mtime.(joinpath.(cachepath, cachefiles)))[2]
                    rm(joinpath(cachepath, cachefiles[idx]))
                end
            end

            # this is atomic according to POSIX:
            cachepathdebug("rename($tmppath, $cachefile; force=true)")
            rename(tmppath, cachefile; force=true)
            return cachefile
        end
    finally
        rm(tmppath, force=true)
    end
    if p.exitcode == 125
        return PrecompilableError()
    else
        error("Failed to precompile $pkg to $tmppath.")
    end
end

function compilecache_path(pkg::PkgId, depot_path::String, prefs_hash::UInt64)::String
    entrypath, entryfile = cache_file_entry(pkg)
    cachepath = joinpath(depot_path, entrypath)
    isdir(cachepath) || mkpath(cachepath)
    if pkg.uuid === nothing
         return_path = abspath(cachepath, entryfile) * ".ji"
    else
        crc = _crc32c(something(Base.active_project(), ""))
        crc = _crc32c(unsafe_string(JLOptions().image_file), crc)
        crc = _crc32c(unsafe_string(JLOptions().julia_bin), crc)
        crc = _crc32c(prefs_hash, crc)
        project_precompile_slug = slug(crc, 5)
        return_path = abspath(cachepath, string(entryfile, "_", project_precompile_slug, ".ji"))
    end
    cachepathdebug("compile_cache_path: $return_path")
    return return_path
end

# Careful to not make this a method of compilecache, else the second arg will
# be interpreted as the package path (source?)
function compilecache_depot(pkg::PkgId, depot_path::String, internal_stderr::IO = stderr, internal_stdout::IO = stdout)
    @nospecialize internal_stderr internal_stdout
    path = locate_package(pkg)
    path === nothing && throw(ArgumentError("$pkg not found during precompilation"))
    return compilecache(pkg, path, depot_path, internal_stderr, internal_stdout)
end

end # module CachePath


text:
# Calculate the Kullback-Leibler divergences between pairs of well-known parameteric
# distributions, and derivatives with respect to the parameters of the first distributions.

######################
# Beta Distributions #
######################

"""
Returns the KL divergence between a pair of beta distributions

Args:
- α₁: the shape parameter for the first beta distribution
- β₁: the scale parameter for the first beta distribution
- α₂: the shape parameter for the second beta distribution
- β₂: the scale parameter for the second beta distribution
"""
function beta_kl(α₁, β₁, α₂, β₂)
    α₁_plus_β₁ = α₁ + β₁
    α₂_plus_β₂ = α₂ + β₂
    α₁_minus_α₂ = α₁ - α₂
    β₁_minus_β₂ = β₁ - β₂
    log_term = lgamma(α₁_plus_β₁) - lgamma(α₁) - lgamma(β₁)
    log_term -= lgamma(α₂_plus_β₂) - lgamma(α₂) - lgamma(β₂)
    apart_term = α₁_minus_α₂ * digamma(α₁) + β₁_minus_β₂ * digamma(β₁)
    together_term = -(α₁_minus_α₂ + β₁_minus_β₂) * digamma(α₁_plus_β₁)
    return log_term + apart_term + together_term
end

#############################
# Categorical Distributions #
#############################

"""
Returns the KL divergence between a pair of categorical distributions
"""
categorical_kl(p₁, p₂) = sum(a * (log(a) - log(b)) for (a, b) in zip(p₁, p₂))

##########################
# Gaussian Distributions #
##########################

"""
Returns the KL divergence between a pair of univariate Gaussian distributions

Args:
- μ₁: the mean for the first Gaussian distribution
- var₁: the variance for the first Gaussian distribution
- μ₂: the mean for the second Gaussian distribution
- var₂: the variance for the second Gaussian distribution
"""
gaussian_kl(μ₁, var₁, μ₂, var₂) = .5 * (log(var₂) - log(var₁) + (var₁ + (μ₁ - μ₂)^2) / var₂ - 1)

#####################################
# Multivariate Normal Distributions #
#####################################

"""
Returns the KL divergence between a pair of multivariate normal distributions,
the first having a diagonal covariance matrix

Args:
- μ₁: the mean for the first normal distribution
- var₁: the variance parameters for the first normal distribution
- μ₂: the mean for the second normal distribution
- Σ₂: the covariance matrix for the second normal distribution
"""
function diagmvn_mvn_kl(μ₁, var₁, μ₂, Σ₂, inv_Σ₂ = inv(Σ₂), logdet_Σ₂ = logdet(Σ₂))
    μ₂_minus_μ₁ = μ₂ - μ₁
    kl = sum(diag(inv_Σ₂) .* var₁) - length(μ₂)
    kl += dot(μ₂_minus_μ₁, inv_Σ₂ * μ₂_minus_μ₁)
    kl += logdet_Σ₂ - sum(log.(var₁))
    return 0.5 * kl
end

###########################################
# Subtracting KL divergences from sources #
###########################################

subtract_kl_a(vs) = -(categorical_kl(vs[ids.a[:, 1]], prior.a))

function subtract_kl_r(vs)
    kl = zero(eltype(vs))
    for i in 1:Ia
        μ₁, var₁ = vs[ids.r1[i]], vs[ids.r2[i]]
        μ₂, var₂ = prior.r_mean[i], prior.r_var[i]
        kl -= vs[ids.a[i, 1]] * gaussian_kl(μ₁, var₁, μ₂, var₂)
    end
    return kl
end

function subtract_kl_k(vs)
    kl = zero(eltype(vs))
    for i in 1:Ia
        kl -= vs[ids.a[i, 1]] * categorical_kl(vs[ids.k[:, i]], prior.k[:, i])
    end
    return kl
end

function subtract_kl_c(vs)
    kl = zero(eltype(vs))
    for i in 1:Ia
        μ₁, var₁ = vs[ids.c1[:, i]], vs[ids.c2[:, i]]
        a = vs[ids.a[i, 1]]
        for d in 1:D
            μ₂, Σ₂ = prior.c_mean[:, d, i], prior.c_cov[:, :, d, i]
            kl -= a * vs[ids.k[d, i]] * diagmvn_mvn_kl(μ₁, var₁, μ₂, Σ₂)
        end
    end
    return kl
end

"""
Subtract the KL divergences for a single source.
"""
function subtract_kl(vs)
    kl = zero(eltype(vs))
    kl += subtract_kl_a(vs)
    kl += subtract_kl_k(vs)
    kl += subtract_kl_r(vs)
    kl += subtract_kl_c(vs)
    return kl
end

###################
# Differentiation #
###################

using ForwardDiff: Dual, JacobianConfig, pickchunksize
using ReverseDiff: compile_gradient

immutable KLHelper{G,H}
    gradient!::G
    hessian!::H
end

const PARAM_LENGTH = length(CanonicalParams)
const DEFAULT_DUAL_TYPE = Dual{pickchunksize(PARAM_LENGTH),Float64}

function KLHelper{N,T}(::Type{Dual{N,T}} = DEFAULT_DUAL_TYPE)
    dual_buffer = zeros(Dual{N,T}, PARAM_LENGTH)
    jacobian_config = JacobianConfig{N}(rand(PARAM_LENGTH))
    gradient! = compile_gradient(subtract_kl, rand(PARAM_LENGTH))
    nested_gradient! = compile_gradient(subtract_kl, rand(Dual{N,T}, PARAM_LENGTH))
    nested_gradient = x -> nested_gradient!(dual_buffer, x)
    hessian! = (out, x) -> ForwardDiff.jacobian!(out, nested_gradient, x, jacobian_config)
    return KLHelper(gradient!, hessian!)
end

const KL_HELPER_POOL = ntuple(n -> KLHelper(), nthreads())

###############
# Entry Point #
###############

function subtract_kl_source!(kl_source::SensitiveFloat, result::DiffBase.DiffResult,
                             vs, helper::KLHelper)
    if kl_source.has_gradient
        helper.gradient!(result, vs)
        kl_source.v[] = DiffBase.value(result)
    else
        kl_source.v[] = subtract_kl(vs)
    end
    if kl_source.has_hessian
        helper.hessian!(kl_source.h, vs)
    end
    return kl_source
end

function subtract_kl_all_sources!{T}(ea::ElboArgs,
                                     accum::SensitiveFloat,
                                     kl_source::SensitiveFloat{T},
                                     helper::KLHelper)
    result = DiffBase.DiffResult(zero(T), kl_source.d)
    for sa in 1:length(ea.active_sources)
        subtract_kl_source!(kl_source, result, ea.vp[ea.active_sources[sa]], helper)
        add_sources_sf!(accum, kl_source, sa)
    end
    return accum
end


text:
include("experiment_parallel.jl")

reps = parse(Int, ARGS[1])
limit = parse(Int, ARGS[2])

function run(n,lim)
    for k = 0:lim
        println()
        println("Running: jonah_parallel($n,$k)")
        println("Average turns per game = ", jonah_parallel(n, k))
    end
end

@time run(reps,limit)


text:
using Sundials
using Grid ## for interpolating grids

##
## Solving the 1D cable equation with IDA.
##
## The DAE system solved is a spatial discretization of the PDE
##
##  C_m dV/dt = d/R d^2V/x^2 - G_m V - R_m J
##
## The PDE is treated with central differences on a uniform M x 1
## grid. The values of u at the interior points satisfy ODEs, and
## equations u = 0 at the boundaries are appended, to form a DAE
## system of size M.
##

## Length of cable in um
L = 4000

## Number of space steps to simulate
dx = 10
xsteps = int(L/dx)

## Number of timesteps to simulate
tf = 10.0
dt = 0.1
timesteps = int(tf/dt)
t = 0.0:dt:(dt*(timesteps-1))

d = 4.0 * 1e-4 ## Cable diameter in cm
R_m = 2.5e11 ## Membrane resistance in Ohms/cm^2
G_m = 1.0/R_m ## Membrane conductance
R_i = 30.0  ## Longitudinal resistivity in Ohms*cm
C_m = 1e-6  ## Membrane capacitance in F/cm^2

## Define the injected current J
J = zeros(timesteps,xsteps)
I_in = 0.1
## When to start injecting current
I_delay = 1.0
## Duration of injected current
I_dur = 5.0
J[max(1,int((I_delay)/dt)):int((I_delay+I_dur)/dt),int(1*xsteps/2)] = I_in

G_J = map(i -> CoordInterpGrid(t, vec(J[:,i]), 0.0, InterpQuadratic),1:xsteps)

##
## cableres: cable equation system residual function. This uses
## central differencing on the interior points, and includes algebraic
## equations for the boundary values.
##

function cableres(t, u, up, r)

    r[:] = u ## Initialize r to u, to take care of boundary equations.

    ## Loop over segments; set res = up - (central difference).
    for i = 2:(xsteps-2)
        loc = i
        r[loc] = C_m * up[loc] - (d/(4*R_i)) * (u[loc-1] + u[loc+1] -
                                                2.0 * u[loc]) +
                 G_m * u[loc] - R_m * (G_J[loc])[t]

    end

    return (0)
end

function initial ()

    u = zeros(xsteps)

    u[2:xsteps-2] = -60.0 ## initial value -60 mV

    id = ones(xsteps)

    up = zeros (xsteps)
    r  = zeros (xsteps)

    cableres(0.0, u, up, r)

    ## Copy -res into up to get correct interior initial up values.
    up[:] = -1.0 * r

    return (u,up,id)

end

nvector = Sundials.nvector
function idabandsol(f::Function, y0::Vector{Float64}, yp0::Vector{Float64},
                    id::Vector{Float64}, t::Vector{Float64};
                    reltol::Float64=1e-4, abstol::Float64=1e-6)

    neq = length(y0)
    mem = Sundials.IDACreate()
    flag = Sundials.IDAInit(mem, cfunction(Sundials.idasolfun, Int32,
                                           (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Function)),
                            t[1], nvector(y0), nvector(yp0))
    assert (flag == 0)
    flag = Sundials.IDASetId(mem,nvector(id))
    assert (flag == 0)
    flag = Sundials.IDASetUserData(mem, f)
    assert (flag == 0)
    flag = Sundials.IDASStolerances(mem, reltol, abstol)
    assert (flag == 0)
    mu = xsteps-1
    ml = xsteps-1
    flag = Sundials.IDABand(mem, neq, mu, ml)
    ##flag = Sundials.IDADense(mem, neq)
    assert (flag == 0)
    rtest = zeros(neq)
    flag = Sundials.IDACalcIC(mem, Sundials.IDA_YA_YDP_INIT, t[2])
    assert (flag == 0)
    yres = zeros(length(t), length(y0))
    ypres = zeros(length(t), length(y0))
    yres[1,:] = y0
    ypres[1,:] = yp0
    y = copy(y0)
    yp = copy(yp0)
    tout = [0.0]
    for k in 2:length(t)
        retval = Sundials.IDASolve(mem, t[k], tout, y, yp, Sundials.IDA_NORMAL)
        yres[k,:] = y[:]
        ypres[k,:] = yp[:]
    end
    return yres, ypres
end

u0, up0, id = initial ()

yout, ypout = @time idabandsol(cableres, u0, up0, id, map(x -> x, t),
                               reltol = 1e-3, abstol = 1e-4)


text:
# module Jive

using FileWatching

const watched_folders = Dict{String, Tuple{FolderMonitor,Task}}()

"""
    watch(callback::Function, dir::String; targets=ARGS, sources::Union{Vector{Any},Vector{String}}=[])

watch the folders.
"""
function watch(callback::Function, dir::String; targets=ARGS, sources::Union{Vector{Any},Vector{String}}=[])
    (all_files, start_idx) = get_all_files(dir, [], targets)
    for src in sources
        for (root, dirs, files) in walkdir(isfile(src) ? dirname(src) : src)
            for filename in files
                !endswith(filename, ".jl") && continue
                subpath = path_separator_to_slash(relpath(normpath(root, filename), dir))
                push!(all_files, subpath)
            end
        end
    end

    function run_callback(cb, path)
        local oldstderr = stderr
        local errread, errwrite, errstream
        try
            (errread, errwrite) = redirect_stderr()
            cb(path)
        catch err
            @info :watch_catch err
        finally
            errstream = @async read(errread, String)
            redirect_stderr(oldstderr)
            close(errwrite)
        end
        errmsg = fetch(errstream)
        if !isempty(errmsg)
            for line in split(errmsg, '\n')
                !startswith(line, "WARNING: replacing module test") && println(stderr, line)
            end
        end
    end

    folders = Set(dirname.(normpath.(dir, all_files)))
    for folder in folders
        fm = FolderMonitor(folder)
        loop = Task() do
            last_time = 0
            while isopen(fm.notify)
                (fname, events) = wait(fm)::Pair
                if splitext(fname)[2] == ".jl" && time() - last_time > 0.010
                    filepath = normpath(folder, fname)
                    run_callback(callback, relpath(filepath, dir))
                    last_time = time()
                end
            end
        end
        watched_folders[folder] = (fm, loop)
    end
    for (folder, (fm, loop)) in watched_folders
        schedule(loop)
    end

    printstyled("watching folders ...\n", color=:green)
    for folder in (sort ∘ collect ∘ keys)(watched_folders)
        println("  - ", relpath(folder, dir))
    end
end

"""
    Jive.stop(::typeof(watch))

stop watching folders.
"""
function stop(::typeof(watch))
    for (folder, (fm, loop)) in watched_folders
        close(fm)
    end
    empty!(watched_folders)
    println("stopped watching folders.")
end


text:
using LinearAlgebra, ProgressMeter, JLD2
import KitBase

cd(@__DIR__)
D = KitBase.read_dict("naca0012.txt")
set = KitBase.set_setup(D)

ps = KitBase.set_geometry(D)
for i in eachindex(ps.cellType)
    if ps.cellType[i] == 1 &&
       -0.5 < ps.cellCenter[i, 1] < 1.5 &&
       -0.1 < ps.cellCenter[i, 2] < 0.1
        ps.cellType[i] = 2

        for j = 1:3
            if ps.faceType[ps.cellFaces[i, j]] == 1
                ps.faceType[ps.cellFaces[i, j]] = 2
            end
        end
    elseif ps.cellType[i] == 1 && ps.cellCenter[i, 1] > 0.5
        ps.cellType[i] = 3

        for j = 1:3
            if ps.faceType[ps.cellFaces[i, j]] == 1
                ps.faceType[ps.cellFaces[i, j]] = 3
            end
        end
    end
end

vs = KitBase.set_velocity(D)
gas = KitBase.set_property(D)
c0 = KitBase.sound_speed(1.0, gas.γ) * gas.Ma
α = π / 180 * 7.0

begin
    primL = [1.0, c0 * cos(α), c0 * sin(α), 1.0]
    #primL = [1.0, KitBase.sound_speed(1.0, gas.γ) * gas.Ma, 0.0, 1.0]
    wL = KitBase.prim_conserve(primL, gas.γ)
    hL = KitBase.maxwellian(vs.u, vs.v, primL)
    bL = @. hL * gas.K / 2 / primL[end]
    primR = [1.0, 0.0, 0.0, 1.0]
    wR = KitBase.prim_conserve(primR, gas.γ)
    hR = KitBase.maxwellian(vs.u, vs.v, primR)
    bR = @. hR * gas.K / 2 / primR[end]
    ib = KitBase.IB2F(wL, primL, hL, bL, primL, wL, primL, hL, bL, primR)
end

ks = KitBase.SolverSet(set, ps, vs, gas, ib, @__DIR__)

ctr, face = KitBase.init_fvm(ks, ks.pSpace)
#@load "ctr.jld2" ctr

dt = KitBase.timestep(ks, ctr, 0.0)
nt = ks.set.maxTime ÷ dt |> Int
res = zeros(4)
@showprogress for iter = 1:nt
    @inbounds Threads.@threads for i in eachindex(face)
        vn = ks.vSpace.u .* face[i].n[1] .+ ks.vSpace.v .* face[i].n[2]
        vt = ks.vSpace.v .* face[i].n[1] .- ks.vSpace.u .* face[i].n[2]

        if !(-1 in ps.faceCells[i, :])
            KitBase.flux_kfvs!(
                face[i].fw,
                face[i].fh,
                face[i].fb,
                ctr[ps.faceCells[i, 1]].h,
                ctr[ps.faceCells[i, 1]].b,
                ctr[ps.faceCells[i, 2]].h,
                ctr[ps.faceCells[i, 2]].b,
                vn,
                vt,
                ks.vSpace.weights,
                dt,
                face[i].len,
            )
            face[i].fw .= KitBase.global_frame(face[i].fw, face[i].n[1], face[i].n[2])
        else
            idx = ifelse(ps.faceCells[i, 1] != -1, 1, 2)

            if ps.cellType[ps.faceCells[i, idx]] == 2
                bc = KitBase.local_frame(ks.ib.primR, face[i].n[1], face[i].n[2])

                KitBase.flux_boundary_maxwell!(
                    face[i].fw,
                    face[i].fh,
                    face[i].fb,
                    bc,
                    ctr[ps.faceCells[i, idx]].h,
                    ctr[ps.faceCells[i, idx]].b,
                    vn,
                    vt,
                    ks.vSpace.weights,
                    ks.gas.K,
                    dt,
                    face[i].len,
                )

                face[i].fw .= KitBase.global_frame(face[i].fw, face[i].n[1], face[i].n[2])
            end
        end
    end

    sumres = zeros(4)
    sumavg = zeros(4)
    @inbounds Threads.@threads for i in eachindex(ctr)
        if ps.cellType[i] in (0, 2)
            dirc = [sign(dot(ctr[i].n[j], face[ps.cellFaces[i, j]].n)) for j = 1:3]

            KitBase.step!(
                ctr[i].w,
                ctr[i].prim,
                ctr[i].h,
                ctr[i].b,
                face[ps.cellFaces[i, 1]].fw,
                face[ps.cellFaces[i, 1]].fh,
                face[ps.cellFaces[i, 1]].fb,
                face[ps.cellFaces[i, 2]].fw,
                face[ps.cellFaces[i, 2]].fh,
                face[ps.cellFaces[i, 2]].fb,
                face[ps.cellFaces[i, 3]].fw,
                face[ps.cellFaces[i, 3]].fh,
                face[ps.cellFaces[i, 3]].fb,
                ks.vSpace.u,
                ks.vSpace.v,
                ks.vSpace.weights,
                ks.gas.K,
                ks.gas.γ,
                ks.gas.μᵣ,
                ks.gas.ω,
                ks.gas.Pr,
                ks.pSpace.cellArea[i],
                dirc,
                dt,
                sumres,
                sumavg,
                :bgk,
            )
        end
    end

    for i in eachindex(ps.cellType)
        if ps.cellType[i] == 3
            ids = ps.cellNeighbors[i, :]
            deleteat!(ids, findall(x -> x == -1, ids))
            id1, id2 = ids
            ctr[i].w .= 0.5 .* (ctr[id1].w .+ ctr[id2].w)
            ctr[i].h .= 0.5 .* (ctr[id1].h .+ ctr[id2].h)
            ctr[i].b .= 0.5 .* (ctr[id1].b .+ ctr[id2].b)
            ctr[i].prim .= KitBase.conserve_prim(ctr[i].w, ks.gas.γ)
        end
    end

    res = sqrt.(length(ctr) .* sumres) ./ (sumavg .+ 1e-6)
    if iter % 1000 == 0
        @show res
        @save "ctr.jld2" ctr
        KitBase.write_vtk(ks, ctr)
    end
end

KitBase.write_vtk(ks, ctr)
@save "ctr.jld2" ctr


text:
#console.jl: Console tools.

#==Output
===============================================================================#

#Print with underline
function println_ul(text::String, ulchar = "─")
	underline = repeat(ulchar, length(text))
	println(text)
	println(underline)
end

#==Input
===============================================================================#

#Wrapper functions:
_parse(::Type{String}, value::AbstractString) = string(value)
_parse(::Type{T}, value::AbstractString) where T<:Real = parse(T, value)

#Prompt user for a value
function input(::Type{T}, prompt::String) where T
	print(prompt)
	return _parse(T, strip(readline(stdin)))
end

#Prompt user for a value, using default
function input(::Type{T}, prompt::String, default) where T
	println(prompt)
	print("[$default]: ")
	result = _parse(T, strip(readline(stdin)))
	if "" == result
		result = default
	end
	return result
end

#Prompt user for string:
input(prompt::String) = input(String, prompt)

pause() = input("Press enter to continue...")

#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick