text:
"""
Convert a two-bit genotype to a real number (minor allele count) of type
`t` according to specified SNP model. Missing genotype is converted
to null. `minor_allele==true` indicates `REF` is the minor allele;
`minor_allele==false` indicates `ALT` is the minor allele.
"""
function convert_gt(
    t::Type{T},
    a::NTuple{2, Bool},
    minor_allele::Bool,
    model::Symbol = :additive
    ) where T <: Real
    if minor_allele # REF is the minor allele
        if model == :additive
            return convert(T, a[1] + a[2])
        elseif model == :dominant
            return convert(T, a[1] | a[2])
        elseif model == :recessive
            return convert(T, a[1] & a[2])
        else
            throw(ArgumentError("un-recognized SNP model: $model"))
        end
    else # ALT is the minor allele
        if model == :additive
            return convert(T, !a[1] + !a[2])
        elseif model == :dominant
            return convert(T, !a[1] | !a[2])
        elseif model == :recessive
            return convert(T, !a[1] & !a[2])
        else
            throw(ArgumentError("un-recognized SNP model: $model"))
        end
    end
end

"""
    copy_gt!(A, reader; [model=:additive], [impute=false], [center=false], [scale=false])

Fill the columns of matrix `A` by the GT data from VCF records in `reader` where the minor
allele for each record is computed on-the-fly. Each column of `A` corresponds to one record.
Record without GT field is converted to `missing`.

# Input
- `A`: a matrix or vector such that `eltype(A) <: Union{Missing, Real}`
- `reader`: a VCF reader

# Optional argument
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype or not, default `false`
- `center`: center gentoype by 2maf or not, default `false`
- `scale`: scale genotype by 1/√2maf(1-maf) or not, default `false`

# Output
- `A`: `ismissing(A[i, j]) == true` indicates missing genotype. If `impute=true`,
    `ismissing(A[i, j]) == false` for all entries.
"""
function copy_gt!(
    A::Union{AbstractMatrix{Union{Missing, T}}, AbstractVector{Union{Missing, T}}},
    reader::VCF.Reader;
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false
    ) where T <: Real
    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Only $j records left in reader; columns $(j+1)-$(size(A, 2)) are set to missing values")
            fill!(view(A, :, (j + 1):size(A, 2)), missing)
            break
        else
            record = read(reader)
        end
        gtkey = VCF.findgenokey(record, "GT")
        # if no GT field, fill by missing values
        if gtkey == nothing
            @inbounds @simd for i in 1:size(A, 1)
                A[i, j] = missing
            end
        end
        # convert GT field to numbers according to specified genetic model
        _, _, _, _, _, _, _, _, minor_allele, maf, _ = gtstats(record, nothing)
        # second pass: impute, convert, center, scale
        ct = 2maf
        wt = maf == 0 ? 1.0 : 1.0 / √(2maf * (1 - maf))
        for i in 1:size(A, 1)
            geno = record.genotype[i]
            # Missing genotype: dropped field or when either haplotype contains "."
            if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
                if impute
                    if minor_allele # REF is the minor allele
                        a1, a2 = rand() ≤ maf, rand() ≤ maf
                    else # ALT is the minor allele
                        a1, a2 = rand() > maf, rand() > maf
                    end
                    A[i, j] = convert_gt(T, (a1, a2), minor_allele, model)
                else
                    A[i, j] = missing
                end
            else # not missing
                # "0" (REF) => 0x30, "1" (ALT) => 0x31
                a1 = record.data[geno[gtkey][1]] == 0x31
                a2 = record.data[geno[gtkey][3]] == 0x31
                A[i, j] = convert_gt(T, (a1, a2), minor_allele, model)
            end
            # center and scale if asked
            center && !ismissing(A[i, j]) && (A[i, j] -= ct)
            scale && !ismissing(A[i, j]) && (A[i, j] *= wt)
        end
    end
    A
end

"""
    copy_gt_as_is!(A, reader; [model=:additive], [impute=false], [center=false], [scale=false])

Fill the columns of matrix `A` by the GT data from VCF records in `reader` where the ALT allele
in each record is interpreted as a `1`. Each column of `A` corresponds to one record.
Record without GT field is converted to `missing`.

# Input
- `A`: a matrix or vector such that `eltype(A) <: Union{Missing, Real}`
- `reader`: a VCF reader

# Optional argument
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype or not, default `false`
- `center`: center gentoype by 2maf or not, default `false`
- `scale`: scale genotype by 1/√2maf(1-maf) or not, default `false`

# Output
- `A`: `ismissing(A[i, j]) == true` indicates missing genotype. If `impute=true`,
    `ismissing(A[i, j]) == false` for all entries.
"""
function copy_gt_as_is!(
    A::Union{AbstractMatrix{Union{Missing, T}}, AbstractVector{Union{Missing, T}}},
    reader::VCF.Reader;
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false
    ) where T <: Real
    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Only $j records left in reader; columns $(j+1)-$(size(A, 2)) are set to missing values")
            fill!(view(A, :, (j + 1):size(A, 2)), missing)
            break
        else
            record = read(reader)
        end
        gtkey = VCF.findgenokey(record, "GT")
        # if no GT field, fill by missing values
        if gtkey == nothing
            @inbounds @simd for i in 1:size(A, 1)
                A[i, j] = missing
            end
        end
        # second pass: impute, convert, center, scale
        _, _, _, _, _, maf, _, _, _, _, _ = gtstats(record, nothing) #reffreq in gtstats actually count number of "1"
        ct = 2maf
        wt = maf == 0 ? 1.0 : 1.0 / √(2maf * (1 - maf))
        for i in 1:size(A, 1)
            geno = record.genotype[i]
            # Missing genotype: dropped field or when either haplotype contains "."
            if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
                if impute
                    a1, a2 = rand() ≤ maf, rand() ≤ maf
                    A[i, j] = convert_gt(T, (a1, a2), true, model)
                else
                    A[i, j] = missing
                end
            else # not missing
                # "0" (REF) => 0x30, "1" (ALT) => 0x31
                a1 = record.data[geno[gtkey][1]] == 0x31
                a2 = record.data[geno[gtkey][3]] == 0x31
                A[i, j] = convert(T, a1 + a2)
            end
            # center and scale if asked
            center && !ismissing(A[i, j]) && (A[i, j] -= ct)
            scale && !ismissing(A[i, j]) && (A[i, j] *= wt)
        end
    end
    return A
end

function geno_ismissing(record::VCF.Record, range::UnitRange{Int})
    return record.data[first(range)] == UInt8('.') || record.data[last(range)] == UInt8('.')
end

"""
    convert_gt!(t, vcffile; [impute=false], [center=false], [scale=false])

Convert the GT data from a VCF file to a matrix of type `t`, allowing for
missing values. Each column of the matrix corresponds to one VCF record.
Record without GT field is converted to equivalent of missing genotypes.

# Input
- `t`: a type `t <: Real`
- `vcffile`: VCF file path

# Optional argument
- `as_minorallele`: convert VCF data (1 indicating ALT allele) to minor allele count or not. If `false`, 0 and 1 will be read as stored in VCF file, default `false`.
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype or not, default `false`
- `center`: center gentoype by 2maf or not, default `false`
- `scale`: scale genotype by 1/√2maf(1-maf) or not, default `false`

# Output
- `A`: matrix where `eltype(A) <: Union{missing, Real}`. `ismissing(A[i, j]) == true`
    indicates missing genotype.
"""
function convert_gt(
    t::Type{T},
    vcffile::AbstractString;
    as_minorallele::Bool = false,
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false,
    ) where T <: Real
    out = Matrix{Union{t, Missing}}(undef, nsamples(vcffile), nrecords(vcffile))
    reader = VCF.Reader(openvcf(vcffile, "r"))
    if as_minorallele
        copy_gt!(out, reader; model = model, impute = impute,
            center = center, scale = scale)
    else
        copy_gt_as_is!(out, reader; model = model, impute = impute,
            center = center, scale = scale)
    end
    close(reader)
    out
end

"""
    convert_ht(t, vcffile; [as_minorallele=false])

Converts the GT data from a VCF file to a haplotype matrix of type `t`. One
column of the VCF record will become 2 column in the resulting matrix. If haplotypes
are not phased, 1 will be on the left column (i.e. `1/0`).

# Input
- `t`: a type `t <: Real`
- `vcffile`: VCF file path

# Optional argument
- `as_minorallele`: convert VCF data (1 indicating ALT allele) to minor allele count or not. If `false`, 0 and 1 will be read as stored in VCF file, default `false`.
"""
function convert_ht(
    t::Type{T},
    vcffile::AbstractString;
    as_minorallele::Bool = false
    ) where T <: Real
    out = Matrix{t}(undef, 2nsamples(vcffile), nrecords(vcffile))
    reader = VCF.Reader(openvcf(vcffile, "r"))
    if as_minorallele
        copy_ht!(out, reader)
    else
        copy_ht_as_is!(out, reader)
    end
    close(reader)
    return out
end

"""
Convert a one-bit haplotype to a real number (minor allele count) of type
`t`. Here `minor_allele==true` indicates `REF` is the minor allele;
`minor_allele==false` indicates `ALT` is the minor allele.
"""
function convert_ht(
    t::Type{T},
    a::Bool,
    minor_allele::Bool
    ) where T <: Real
    if minor_allele # REF is the minor allele
        return convert(T, a)
    else # ALT is the minor allele
        return convert(T, !a)
    end
end

"""
    copy_ht!(A, reader)

Fill 2 columns of BitArray `A` by the GT data from VCF records in `reader`, each record filling 2 columns.
The minor allele for each record is computed on-the-fly. Record without GT field is converted to `missing`.
Missing GT field is NOT allowed.

# Input
- `A`: a BitArray{2} or BitArray{1}.
- `reader`: a VCF reader
"""
function copy_ht!(
    A::Union{AbstractMatrix{T}, AbstractVector{T}},
    reader::VCF.Reader
    ) where T <: Real

    n, p = size(A)
    nn   = Int(n / 2)

    # loop over every record, each filling 2 columns of A
    for j in 1:p
        if eof(reader)
            @warn("Reached end! Check if output is correct!")
            A = A[:, 1:2j]
            break
        else
            record = read(reader)
        end
        gtkey = VCF.findgenokey(record, "GT")

        # haplotype reference files must have GT field
        if gtkey == nothing
            error("Missing GT field for record $j. Reference panels cannot have missing data!")
        end

        # check whether allele is REF or ALT
        _, _, _, _, _, _, _, _, minor_allele, _, _ = gtstats(record, nothing)

        # loop over every marker in record
        for i in 1:nn
            geno = record.genotype[i]
            # Missing genotype: dropped field or when either haplotype contains "."
            if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
                error("Missing GT field for record $j entry $i. Reference panels cannot have missing data!")
            else # not missing
                # "0" (REF) => 0x30, "1" (ALT) => 0x31
                a1 = record.data[geno[gtkey][1]] == 0x31
                a2 = record.data[geno[gtkey][3]] == 0x31
                A[2i - 1, j] = convert_ht(T, a1, minor_allele)
                A[2i    , j] = convert_ht(T, a2, minor_allele)
            end
        end
    end
end

"""
    copy_ht_as_is!(A, reader)

Fill 2 columns of BitArray `A` by the GT data from VCF records in `reader`, each record filling 2 columns.
The minor allele rfor each ecord is the ALT allele (i.e. will read 0s and 1s of the vcffile as-is). Record
without GT field is converted to `missing`. Missing GT field is NOT allowed.

# Input
- `A`: a BitArray{2} or BitArray{1}.
- `reader`: a VCF reader
"""
function copy_ht_as_is!(
    A::Union{AbstractMatrix{T}, AbstractVector{T}},
    reader::VCF.Reader;
    ) where T <: Real

    n, p = size(A)
    nn   = Int(n / 2)

    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Reached end! Check if output is correct!")
            A = A[:, 1:2j]
            break
        else
            record = read(reader)
        end
        gtkey = VCF.findgenokey(record, "GT")

        # haplotype reference files must have GT field
        if gtkey == nothing
            error("Missing GT field for record $j. Reference panels cannot have missing data!")
        end

        # second pass: convert
        for i in 1:nn
            geno = record.genotype[i]
            # Missing genotype: dropped field or when either haplotype contains "."
            if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
                error("Missing GT field for record $j entry $i. Reference panels cannot have missing data!")
            else # not missing
                # "0" (REF) => 0x30, "1" (ALT) => 0x31
                a1 = record.data[geno[gtkey][1]] == 0x31
                a2 = record.data[geno[gtkey][3]] == 0x31
                A[2i - 1, j] = convert(T, a1)
                A[2i    , j] = convert(T, a2)
            end
        end
    end
    return A
end

"""
    convert_ds(t, vcffile; [key = "DS"], [impute=false], [center=false], [scale=false])

Converts dosage data from a VCF file to a numeric matrix of type `t`. Here `key` specifies
the FORMAT field of the VCF file that encodes the dosage (default = "DS"). If `impute=true`,
we fill missing entries with 2 times the ALT allele frequency.

# Arguments
- `t`: type of output matrix.
- `vcffile`: VCF file path

# Optional argument
- `key`: The FIELD name if the VCF file that encodes dosages, defaults to `"DS"`
- `as_minorallele`: convert VCF data (1 indicating ALT allele) to minor allele count or not. If `false`, 0 and 1 will be read as stored in VCF file, default `false`.
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype or not, default `false`
- `center`: center gentoype by 2maf or not, default `false`
- `scale`: scale genotype by 1/√2maf(1-maf) or not, default `false`
"""
function convert_ds(
    t::Type{T},
    vcffile::AbstractString;
    key::String = "DS",
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false
    ) where T <: Real
    out = Matrix{Union{Missing, t}}(undef, nsamples(vcffile), nrecords(vcffile))
    reader = VCF.Reader(openvcf(vcffile, "r"))
    copy_ds!(out, reader, key = key, impute = impute, center = center, scale = scale)
    close(reader)
    return out
end

"""
    copy_ds!(A, reader; [key="DS"], [model=:additive], [impute=false], [center=false], [scale=false])

Fill the columns of matrix `A` by the dosage data from VCF records in `reader`. Record without GT field
is converted to `missing`.

# Input
- `A`: a matrix or vector such that `eltype(A) <: Union{Missing, Real}`
- `reader`: a VCF reader

# Optional argument
- `key`: The field key for accessing dosage values, default `"DS"`
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype or not, default `false`
- `center`: center gentoype by 2maf or not, default `false`
- `scale`: scale genotype by 1/√2maf(1-maf) or not, default `false`

# Output
- `A`: `ismissing(A[i, j]) == true` indicates missing genotype. If `impute=true`,
    `ismissing(A[i, j]) == false` for all entries.
"""
function copy_ds!(
    A::Union{AbstractMatrix{Union{Missing, T}}, AbstractVector{Union{Missing, T}}},
    reader::VCF.Reader;
    key::String = "DS",
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false
    ) where T <: Real
    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Only $j records left in reader; columns $(j+1)-$(size(A, 2)) are set to missing values")
            fill!(view(A, :, (j + 1):size(A, 2)), missing)
            break
        else
            record = read(reader)
        end
        dskey = VCF.findgenokey(record, key)

        # if no dosage field, fill by missing values
        if dskey == nothing
            @inbounds @simd for i in 1:size(A, 1)
                A[i, j] = missing
            end
        end

        # loop over every marker in record
        _, _, _, _, _, maf, _, _, _, _, _ = gtstats(record, nothing) #reffreq in gtstats actually count number of "1"
        ct = 2maf
        wt = maf == 0 ? 1.0 : 1.0 / √(2maf * (1 - maf))
        for i in 1:size(A, 1)
            geno = record.genotype[i]
            # Missing genotype: dropped field or "." => 0x2e
            if dskey > lastindex(geno) || record.data[geno[dskey]] == [0x2e]
                A[i, j] = (impute ? ct : missing)
            else # not missing
                A[i, j] = parse(T, String(record.data[geno[dskey]]))
            end
            # center and scale if asked
            center && !ismissing(A[i, j]) && (A[i, j] -= ct)
            scale && !ismissing(A[i, j]) && (A[i, j] *= wt)
        end
    end
end


text:
facts("File formats") do

    context("avr") do
        include("avr.jl")
    end
    context("bsa") do
        include("bsa.jl")
    end
    context("dat") do
        include("dat.jl")
    end
    context("elp") do
        include("elp.jl")
    end
    context("evt") do
        include("evt.jl")
    end
    context("rba") do
        include("rba.jl")
    end
    context("sfp") do
        include("rba.jl")
    end
end

ext: jl


text:
abstract type RosenbrockMutableCache <: OrdinaryDiffEqMutableCache end
################################################################################

# Shampine's Low-order Rosenbrocks

mutable struct Rosenbrock23Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::du2Type
  f₁::rateType
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock23Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock23Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock23Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock23Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

mutable struct Rosenbrock32Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du1::rateType
  du2::du2Type
  f₁::rateType
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock32Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock32Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock32Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock32Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rosenbrock23ConstantCache(uEltypeNoUnits,identity,identity)
  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)

  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)

  Rosenbrock23Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  # f₀ = similar(u) fsalfirst
  f₁ = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rosenbrock32ConstantCache(uEltypeNoUnits,identity,identity)

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock32Cache(u,uprev,k₁,k₂,k₃,du1,du2,f₁,vectmp,vectmp2,vectmp3,fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

struct Rosenbrock23ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock23ConstantCache(T::Type,tf,uf)
  c₃₂ = T(6 + sqrt(2))
  d = T(1/(2+sqrt(2)))
  Rosenbrock23ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock23,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock23ConstantCache(uEltypeNoUnits,tf,uf)
end

struct Rosenbrock32ConstantCache{T,TF,UF} <: OrdinaryDiffEqConstantCache
  c₃₂::T
  d::T
  tf::TF
  uf::UF
end

function Rosenbrock32ConstantCache(T::Type,tf,uf)
  c₃₂ = T(6 + sqrt(2))
  d = T(1/(2+sqrt(2)))
  Rosenbrock32ConstantCache(c₃₂,d,tf,uf)
end

function alg_cache(alg::Rosenbrock32,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock32ConstantCache(uEltypeNoUnits,tf,uf)
end

################################################################################

### 3rd order specialized Rosenbrocks

struct Rosenbrock33ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

mutable struct Rosenbrock33Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::du2Type
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  vectmp4::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock33Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock33Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock33Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock33Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = ROS3PConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock33Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::ROS3P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock33ConstantCache(tf,uf,ROS3PConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

mutable struct Rosenbrock34Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::du2Type
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  vectmp4::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock34Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock34Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock34Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock34Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rodas3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock34Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

struct Rosenbrock34ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

function alg_cache(alg::Rodas3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock34ConstantCache(tf,uf,Rodas3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

################################################################################

struct Rosenbrock4ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

mutable struct Rosenbrock4Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  du::rateType
  du1::rateType
  du2::du2Type
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  vectmp4::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock4Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock4Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock4Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock4Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::RosShamp4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = RosShamp4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::RosShamp4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,RosShamp4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::Veldd4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Veldd4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Veldd4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,Veldd4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::Velds4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Velds4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Velds4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,Velds4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::GRK4T,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = GRK4TConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::GRK4T,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,GRK4TConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::GRK4A,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = GRK4AConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::GRK4A,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,GRK4AConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::Ros4LStab,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Ros4LStabConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock4Cache(u,uprev,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Ros4LStab,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock4ConstantCache(tf,uf,Ros4LStabConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

###############################################################################

### Rodas methods

struct Rodas4ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

mutable struct Rodas4Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::du2Type
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  vectmp4::vecuType
  vectmp5::vecuType
  vectmp6::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rodas4Cache) = (c.dT,c.tmp)
du_cache(c::Rodas4Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rodas4Cache) = (c.J,c.W)
vecu_cache(c::Rodas4Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  vectmp5 = vec(similar(u,axes(u)))
  vectmp6 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rodas4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    vectmp5,vectmp6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  vectmp5 = vec(similar(u,axes(u)))
  vectmp6 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rodas42ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    vectmp5,vectmp6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas42,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas42ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  vectmp5 = vec(similar(u,axes(u)))
  vectmp6 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rodas4PConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rodas4Cache(u,uprev,dense1,dense2,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    vectmp5,vectmp6,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas4P,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rodas4ConstantCache(tf,uf,Rodas4PConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end

################################################################################

### Rosenbrock5

struct Rosenbrock5ConstantCache{TF,UF,Tab} <: OrdinaryDiffEqConstantCache
  tf::TF
  uf::UF
  tab::Tab
end

mutable struct Rosenbrock5Cache{uType,uArrayType,rateType,du2Type,LinuType,vecuType,JType,WType,TabType,TFType,UFType,F,JCType,GCType} <: RosenbrockMutableCache
  u::uType
  uprev::uType
  dense1::rateType
  dense2::rateType
  du::rateType
  du1::rateType
  du2::du2Type
  vectmp::vecuType
  vectmp2::vecuType
  vectmp3::vecuType
  vectmp4::vecuType
  vectmp5::vecuType
  vectmp6::vecuType
  vectmp7::vecuType
  vectmp8::vecuType
  fsalfirst::rateType
  fsallast::rateType
  dT::uArrayType
  J::JType
  W::WType
  tmp::uArrayType
  tab::TabType
  tf::TFType
  uf::UFType
  linsolve_tmp::LinuType
  linsolve_tmp_vec::vecuType
  linsolve::F
  jac_config::JCType
  grad_config::GCType
end

u_cache(c::Rosenbrock5Cache) = (c.dT,c.tmp)
du_cache(c::Rosenbrock5Cache) = (c.k₁,c.k₂,c.k₃,c.du1,c.du2,c.f₁,c.fsalfirst,c.fsallast,c.linsolve_tmp)
jac_cache(c::Rosenbrock5Cache) = (c.J,c.W)
vecu_cache(c::Rosenbrock5Cache) = (c.vectmp,c.vectmp2,c.vectmp3)

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  dense1 = zero(rate_prototype)
  dense2 = zero(rate_prototype)
  du = zero(rate_prototype)
  du1 = zero(rate_prototype)
  du2 = zero(rate_prototype)
  vectmp = vec(similar(u,axes(u)))
  vectmp2 = vec(similar(u,axes(u)))
  vectmp3 = vec(similar(u,axes(u)))
  vectmp4 = vec(similar(u,axes(u)))
  vectmp5 = vec(similar(u,axes(u)))
  vectmp6 = vec(similar(u,axes(u)))
  vectmp7 = vec(similar(u,axes(u)))
  vectmp8 = vec(similar(u,axes(u)))
  fsalfirst = zero(rate_prototype)
  fsallast = zero(rate_prototype)
  dT = similar(u,axes(u))
  if DiffEqBase.has_jac(f) && !DiffEqBase.has_invW(f) && f.jac_prototype != nothing
    W = WOperator(f, dt)
    J = nothing # is J = W.J better?
  else
    J = fill(zero(uEltypeNoUnits),length(u),length(u)) # uEltype?
    W = similar(J)
  end
  tmp = similar(u,axes(u))
  tab = Rodas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

  tf = DiffEqDiffTools.TimeGradientWrapper(f,uprev,p)
  uf = DiffEqDiffTools.UJacobianWrapper(f,t,p)
  linsolve_tmp = similar(u,axes(u))
  linsolve_tmp_vec = vec(linsolve_tmp)
  linsolve = alg.linsolve(Val{:init},uf,u)
  grad_config = build_grad_config(alg,f,tf,du1,t)
  jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,du2)
  Rosenbrock5Cache(u,uprev,dense1,dense2,du,du1,du2,vectmp,vectmp2,vectmp3,vectmp4,
                    vectmp5,vectmp6,vectmp7,vectmp8,
                    fsalfirst,fsallast,dT,J,W,tmp,tab,tf,uf,linsolve_tmp,
                    linsolve_tmp_vec,linsolve,jac_config,grad_config)
end

function alg_cache(alg::Rodas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  tf = DiffEqDiffTools.TimeDerivativeWrapper(f,u,p)
  uf = DiffEqDiffTools.UDerivativeWrapper(f,t,p)
  Rosenbrock5ConstantCache(tf,uf,Rodas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits)))
end


text:
export HLAAllele, HLAType, parse_allele, is_valid_allele, limit_hla_accuracy, hla_accuracy

struct HLAAllele
    gene::String
    field_1::String
    field_2::Union{Missing, String}
    field_3::Union{Missing, String}
    field_4::Union{Missing, String}
    suffix::Union{Missing, String}
end

struct HLAType
    alleles::NTuple{6, Union{HLAAllele, Missing}}

    function HLAType(alleles::NTuple{6, Union{HLAAllele, Missing}})
        genes = tuple(ismissing(allele) ? missing : allele.gene for allele in alleles)

        for gene in ("A", "B", "C")
            c = count(genes .== gene)
            if c > 2
                error("HLA type must contain at  2 $gene alleles, got $(c)")
            end
        end

        new(alleles)
    end
end

function HLAType(::Missing)
    HLAType((missing, missing, missing, missing, missing, missing))
end

function hla_matrix(data::Vector{HLAType}; allele_depth::Int = 1)
    alleles = sort(unique_alleles(data, allele_depth = allele_depth))
    m = zeros(Int, length(data), length(alleles))

    for i in eachindex(data)
        hla_type = data[i]
        for allele in hla_type.alleles
            a = limit_hla_accuracy(allele, allele_depth = allele_depth)
            (ismissing(a) || hla_accuracy(a) != allele_depth) && continue

            m[i, findfirst(x -> x == a, alleles)] = 1
        end
    end

    return m
end

function hla_matrix_standardized(data::Vector{HLAType}; allele_depth::Int = 1)
    alleles = sort(unique_alleles(data, allele_depth = allele_depth))
    m = zeros(Float64, length(data), length(alleles))

    for (i, allele) in enumerate(alleles)
        frequency = count(x -> allele in x, data) / length(data)
        for (j, hla_type) in enumerate(data)
            if allele in hla_type
                m[j, i] = 1 - frequency
            else
                m[j, i] = -frequency
            end
        end
    end

    return m
end

function Base.in(::Missing, hla_type::HLAType)
    if any(ismissing.(hla_type.alleles))
        return true
    end

    return false
end

function Base.in(allele::HLAAllele, hla_type::HLAType)
    for a in hla_type.alleles
        ismissing(a) && continue

        if allele == limit_hla_accuracy(a, allele_depth = hla_accuracy(allele))
            return true
        end
    end

    return false
end

function HLAAllele(s::AbstractString)
    if is_valid_allele(s)
        regex = raw"HLA-([ABC])\*(\d{2}(?!\d)):?((?<=:)\d{2,3}(?!\d))?:?" *
                raw"((?<=:)\d\d(?!\d))?:?((?<=:)\d\d(?!\d))?([NLSCAQ])?" |> Regex

        captures = match(regex, s).captures
        return HLAAllele((i === nothing ? missing : String(i) for i in captures)...)
    else
        throw(error("$s does not follow HLA nomenclature " *
              "(http://hla.alleles.org/nomenclature/naming.html)"))
    end
end

function is_valid_allele(s::AbstractString)
    regex = raw"HLA-([ABC])\*(\d{2}(?!\d)):?((?<=:)\d{2,3}(?!\d))?:?" *
            raw"((?<=:)\d\d(?!\d))?:?((?<=:)\d\d(?!\d))?([NLSCAQ])?" |> Regex

    if match(regex, s) === nothing || length(match(regex, s).captures) < 3
        return false
    end

    return true
end

function Base.isless(x::HLAAllele, y::HLAAllele)
    for name in fieldnames(typeof(x))
        x_value = getfield(x, name)
        y_value = getfield(y, name)

        xor(ismissing(x_value), ismissing(y_value)) && return missing
        ismissing(x_value) && ismissing(y_value) && return false

        if x_value < y_value
            return true
        elseif x_value > y_value
            return false
        end
    end

    return false
end

function parse_allele(::Missing)
    return missing
end

function parse_allele(x::AbstractString)
    if occursin("non", x)
        return HLAAllele(missing)
    end

    s = replace(x, r"[()]" => s"")

    # place * after gene, remove w in Cw, add 'HLA-' prefix
    s = replace(s, r"([ABC])w?\*?(\d.*)" => s"HLA-\1*\2")

    # append 0 for single-digit allele groups
    s = replace(s, r"HLA-([ABC])w?\*?([1-9])(\D|$)(.*)" => s"HLA-\1*0\2")

    # place field seperator after allele group
    s = replace(s, r"HLA-([ABC]\*)(\d{2}):?(.+)" => s"HLA-\1\2:\3")

    # place field seperator for 3 digit field_1
    s = replace(s, r"HLA-([ABC]\*\d{2}:)(\d{3}):?(\d{2}):?(\d{2})(.*)" =>
                s"HLA-\1\2:\3:\4\5")

    # place field seperator for 2 digit field_1
    s = replace(s, r"HLA-([ABC]\*\d{2}:)(\d{2}):?(\d{2}):?(\d{2})(.*)" =>
                s"HLA-\1\2:\3:\4\5")

    # place field seperator for 6 digit HLA alleles
    s = replace(s , r"HLA-([ABC]\*\d{2}:)(\d{2}):?(\d{2})(.*)" => s"HLA-\1\2:\3\4")

    if is_valid_allele(s)
        return HLAAllele(s)
    else
        return missing
    end
end

function parse_allele(x...)
    return NTuple{length(x), HLAAllele}(parse_allele(allele) for allele in x)
end

function convert(::Type{String}, x::HLAAllele)
    s = "HLA-"

    s = s * x.gene * "*" * x.field_1

    s = !ismissing(x.field_2) ? s * ":" * x.field_2 : s
    s = !ismissing(x.field_3) ? s * ":" * x.field_3 : s
    s = !ismissing(x.field_4) ? s * ":" * x.field_4 : s
    s = !ismissing(x.suffix) ? s * x.suffix : s

    return s
end

function Base.show(io::IO, x::HLAAllele)
    print(io, convert(String, x))
end

function Base.rand(S::Type{HLAAllele}, gene::Symbol)
    gene ∉ [:A, :B, :C] && error("gene must be either :A, :B, or :C")
    fasta_path = joinpath(@__DIR__, "..", "data", "HLA_alleles",
                         "$(string(gene))_nuc.fasta")
    fasta_reader = FASTA.Reader(open(fasta_path, "r"))
    allele_strings = Vector{String}()

    for record in fasta_reader
        allele_string = split(FASTA.description(record))[1]
        push!(allele_strings, allele_string)
    end

    return parse_allele(rand(allele_strings))
end

function Base.rand(::Type{HLAType})
    hla_alleles = Tuple(rand(HLAAllele, s) for s in (:A, :A, :B, :B, :C, :C))

    return HLAType(hla_alleles)
end

function Base.rand(::Type{HLAType}, d::Int)
    return [rand(HLAType) for i in 1:d]
end

function unique_alleles(hla_types::Vector{HLAType}; allele_depth::Int = 1)
    alleles = Vector{HLAAllele}()

    for hla_type in hla_types
        for allele in hla_type.alleles
            a = limit_hla_accuracy(allele, allele_depth = allele_depth)
            (ismissing(a) || hla_accuracy(a) != allele_depth) && continue

            if a ∉ alleles
                push!(alleles, a)
            end
        end
    end

    return sort(alleles)
end

function limit_hla_accuracy(::Missing; allele_depth::Int = 1)
    return missing
end

function limit_hla_accuracy(s::HLAAllele; allele_depth::Int = 1)
    1 <= allele_depth <= 5 || error("allele_depth must be between 1 and 5.")

    HLA_components = Vector{Union{Missing, String}}(undef, 6)
    fill!(HLA_components, missing)

    HLA_components[1] = s.gene
    HLA_components[2] = s.field_1

    if allele_depth >= 2
        HLA_components[3] = s.field_2
    elseif allele_depth >= 3
        HLA_components[4] = s.field_3
    elseif allele_depth >= 4
        HLA_components[5] = s.field_4
    elseif allele_depth == 6
        HLA_components[6] = s.suffix
    end

    return HLAAllele(HLA_components[1], HLA_components[2], HLA_components[3],
        HLA_components[4], HLA_components[5], HLA_components[6])
end

function hla_accuracy(s::HLAAllele)
    for (i, field) in enumerate(fieldnames(HLAAllele)[3:end])
        if ismissing(getfield(s, field))
            return i
        end
    end

    return 5
end


text:
include("main.jl")

using JLD
using PyPlot
using CompEcon

evaluate(coefs, k::Vector, x::Vector, v::Vector, order=[0 0 0]) =
    funeval(coefs, m.basis, [k x v], order)

evaluate(coefs, k::Vector, x::Real, v::Real, order=[0 0 0]) =
    funeval(coefs, m.basis, [k fill(x, length(k)) fill(v, length(k))], order)

evaluate(coefs, k::Real, x::Vector, v::Real, order=[0 0 0]) =
    funeval(coefs, m.basis, [fill(k, length(x)) x fill(v, length(x))], order)

evaluate(coefs, k::Real, x::Real, v::Vector, order=[0 0 0]) =
    funeval(coefs, m.basis, [fill(k, length(v)) fill(x, length(v)) v], order)

add_params!(ax, ρ, α, β) =
    ax[:annotate](s="rho:$(ρ)\nalpha:$(α)\nbeta:$(β)", xy=(0.05, 0.85),
                  xycoords="axes fraction")

function plot(m::X.BFZ, var::Symbol; trans=(x,y,z)->(x.^(m.agent.ρ-1).*y),
              lab=L"log(J^{\rho-1}J_k)", fn_suffix=1)
    j_coefs, j, c, kp, grp_key = X.solve(m, false)

    grid, (kgrid, xgrid, vgrid) = nodes(m.basis)
    fig, ax = subplots()

    k, x, v = var == :k ? (kgrid, mean(xgrid), mean(vgrid)) :
              var == :x ? (mean(kgrid), xgrid, mean(vgrid)) :
              var == :v ? (mean(kgrid), mean(xgrid), vgrid) :
              error("var must be one of `:k`, `:x`, or `:v`")

    j = evaluate(j_coefs, k, x, v, [0 0 0])
    j_k = evaluate(j_coefs, k, x, v, [1 0 0])

    # get c
    c_coefs = get_coefs(m.basis, m.basis_struc, c)
    c = evaluate(c_coefs, k, x, v, [0 0 0])

    y = log(trans(j, j_k, c))

    x = var == :k ? log(kgrid) : var == :x ? xgrid : var == :v ? vgrid :
               error("var must be one of `:k`, `:x`, or `:v`")

    ax[:plot](x, y)
    ax[:set_xlabel](var ==:k ? "log(k)" : "$(var)")
    ax[:set_ylabel](lab)

    approx_slope = (y[end] - y[1]) / (x[end] - x[1])
    ax[:set_title]("slope: $approx_slope")

    add_params!(ax, m.agent.ρ, m.agent.α, m.agent.β)

    fig[:savefig]("figures/model_$(grp_key)/plot_$(var)_$(fn_suffix).png")

end

function plot_kapitalratio(m::X.BFZ)
    j_coefs, j, c, kp, grp_key = X.solve(m, false)

    # Take the solution to the model and get relevant info
    c_coefs, x, v, c_sim, k_sim, y_sim = X.simulate_solution(m, j_coefs, c, T=5000)

    fig, ax = subplots()
    ax[:hist](k_sim[1:end-1] ./ y_sim)
    ax[:set_title]("Capital to Income Ratio")
    add_params!(ax, m.agent.ρ, m.agent.α, m.agent.β)
    fig[:savefig]("figures/model_$(grp_key)/capital_output.png")

end

function main(m::X.BFZ)
    grp_key = X.solve(m, false)[end];
    mkpath("figures/model_$(grp_key)")
    for var in (:k, :x, :v)
        plot(m, var, fn_suffix="1")
        plot(m, var, trans=(x,y,z)->x, lab=L"log(J)", fn_suffix="lJ")
        plot(m, var, trans=(x,y,z)->y, lab=L"log(J_k)", fn_suffix="lJk")
        plot(m, var, trans=(x,y,z)->z, lab=L"log(c)", fn_suffix="lc")
    end
    plot_kapitalratio(m)
end

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
#
#	Process
#

function draw_realization!(
	result::Array{Complex{T}} where T<:AbstractFloat,
	g::Function,
	V::AbstractArray)

	@. result = randn(ComplexF64) * sqrt(g(V))
end

#
#	Field
#
function draw_realization!(Y::Cscalar,var::Function)
	for i=1:Npoints(Y)
		Y.values[i] = randn(ComplexF64) * sqrt(var(Y.cl.points[i,1:end]))
	end
end


text:
export iauTaitt
"""
Time scale transformation:  International Atomic Time, TAI, to
Terrestrial Time, TT.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tai1,tai2  double    TAI as a 2-part Julian Date

Returned:
   tt1,tt2    double    TT as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Note:

   tai1+tai2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tai1 is the Julian
   Day Number and tai2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   Explanatory Supplement to the Astronomical Almanac,
   P. Kenneth Seidelmann (ed), University Science Books (1992)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTaitt(double tai1, double tai2, double *tt1, double *tt2)

function iauTaitt(tai1::Real, tai2::Real)
   ref_tt1 = Ref{Float64}(0.0)
   ref_tt2 = Ref{Float64}(0.0)

   status = ccall((:iauTaitt, libsofa_c), Cint,
                   (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tai1), convert(Float64, tai2), ref_tt1, ref_tt2)

   return status, ref_tt1[], ref_tt2[]
end


text:
using FlameGraphs
using Base.StackTraces: StackFrame

function stackframe(func, file, line; C=false, inlined=false)
    if func === :eval
        mi = first(Base.method_instances(eval, Tuple{Expr}))
        mi.def.module = Core
    else
        mi = nothing
    end
    StackFrame(Symbol(func), Symbol(file), line, mi, C, inlined, 0)
end

backtraces = UInt64[
    0,                                     39, 30, 32, 33, 21, 13,
    0,     12, 15, 38, 26, 25, 24, 25, 34, 31, 30, 32, 33, 21, 13,
    0,     16, 19, 38, 26, 25, 24, 25, 34, 31, 30, 32, 33, 21, 13,
    0,     16, 19, 38, 26, 25, 24, 25, 34, 31, 30, 32, 33, 21, 13,
    0,  7, 20, 28, 29, 27, 25, 24, 25, 34, 31, 30, 32, 33, 21, 13,
    0,  7, 20, 28, 29, 27, 25, 24, 25, 34, 31, 30, 32, 33, 21, 13,
    0,                          4,  5, 35, 31, 30, 32, 33, 21, 13,
    0,                      6, 14, 37, 36, 31, 30, 32, 33, 21, 13,
    0,                                     31, 30, 32, 33, 21, 13,
    0,                                         17,  2,  3, 22, 13,
    0,                         99, 18, 11, 10, 17,  2,  3, 22, 13,
    0,                             18, 11, 10, 17,  2,  3, 22, 13,
    0,                             98, 11, 10, 17,  2,  3, 22, 13,
    0,                                 11, 10, 17,  2,  3, 22, 13,
    0,                                 98, 10, 17,  2,  3, 22, 13,
    0,                                     10, 17,  2,  3, 22, 13,
    0,                                     10, 17,  2,  3, 22, 13,
    0,                                     10, 17,  2,  3, 22, 13,
    0,                                     98, 17,  2,  3, 22, 13,
    0,                         99,  8,  9, 41, 40,  1,  3, 23, 13,
    0,                              8,  9, 41, 40,  1,  3, 23, 13,
    0,                              8,  9, 41, 40,  1,  3, 23, 13,
    0,                              8,  9, 41, 40,  1,  3, 23, 13,
    0,                              8,  9, 41, 40,  1,  3, 23, 13,
    0]

lidict = Dict{UInt64,StackFrame}(
    1  => stackframe("#mapslices#115", ".\\abstractarray.jl", 2018),
    2  => stackframe("#mapslices#115", ".\\abstractarray.jl", 2029),
    3  => stackframe("mapslices##kw", ".\\abstractarray.jl", 1972),
    4  => stackframe(:*, ".\\float.jl", 405),
    5  => stackframe(:*, ".\\promotion.jl", 312),
    6  => stackframe(:-, ".\\float.jl", 403),
    7  => stackframe(:+, ".\\int.jl", 53),
    8  => stackframe(:Array, ".\\boot.jl", 407),
    9  => stackframe(:Array, ".\\boot.jl", 415),
    10 => stackframe(:concatenate_setindex!, ".\\abstractarray.jl", 2058),
    11 => stackframe(:dotview, ".\\broadcast.jl", 1138),
    12 => stackframe(:dsfmt_fill_array_close1_open2!, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\DSFMT.jl", 86),
    13 => stackframe(:eval, ".\\boot.jl", 331),
    14 => stackframe(:exp, ".\\special\\exp.jl", 136),
    15 => stackframe(:gen_rand, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 186, inlined=true),
    16 => stackframe(:getproperty, ".\\Base.jl", 33, inlined=true),
    17 => stackframe(:inner_mapslices!, ".\\abstractarray.jl", 2039),
    18 => stackframe(:maybeview, ".\\views.jl", 124),
    19 => stackframe(:mt_empty, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 180, inlined=true),
    20 => stackframe(:mt_pop!, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 183),
    21 => stackframe(:profile_test, ".\\REPL[2]", 3),
    22 => stackframe(:profile_test, ".\\REPL[2]", 5),
    23 => stackframe(:profile_test, ".\\REPL[2]", 7),
    24 => stackframe(:rand, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\generation.jl", 119, inlined=true),
    25 => stackframe(:rand, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\Random.jl", 253, inlined=true),
    26 => stackframe(:rand, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 370, inlined=true),
    27 => stackframe(:rand, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 371, inlined=true),
    28 => stackframe(:rand_inbounds, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 362),
    29 => stackframe(:rand_inbounds, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 366, inlined=true),
    30 => stackframe(:randn!, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 173),
    31 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 167),
    32 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 184),
    33 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 190, inlined=true),
    34 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 40, inlined=true),
    35 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 43),
    36 => stackframe(:randn, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 45),
    37 => stackframe(:randn_unlikely, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\normal.jl", 58),
    38 => stackframe(:reserve_1, "D:\\buildbot\\worker\\package_win64\\build\\usr\\share\\julia\\stdlib\\v1.4\\Random\\src\\RNGs.jl", 190, inlined=true),
    39 => stackframe(:setindex!, ".\\array.jl", 826, inlined=true),
    40 => stackframe(:similar, ".\\abstractarray.jl", 626),
    41 => stackframe(:similar, ".\\array.jl", 361),
    98 => stackframe(:jl_apply_generic, "gf.c", 2318, C=true),
    99 => stackframe(:jl_gc_collect, "gc.c", 3105, C=true))

g = flamegraph(backtraces, lidict=lidict)


text:
abstract PropertyPrimitive

immutable Property{P <: PropertyPrimitive} <: ComposeNode
	primitives::Vector{P}
end

function isempty(p::Property)
    return isempty(p.primitives)
end

function isscalar(p::Property)
    return length(p.primitives) == 1
end

# Some properties can be applied multiple times, most cannot.
function isrepeatable(p::Property)
    return false
end

function absolute_units{T}(p::Property{T}, t::Transform, units::UnitBox,
                           box::AbsoluteBoundingBox)
    return Property{T}([absolute_units(primitive, t, units, box)
                        for primitive in p.primitives])
end

# Property primitive catchall: most properties don't need measure transforms
function absolute_units{T <: PropertyPrimitive}(primitive::T, t::Transform,
                                                units::UnitBox,
                                                box::AbsoluteBoundingBox)
    return primitive
end

# Stroke
# ------

immutable StrokePrimitive <: PropertyPrimitive
	color::Maybe(ColorValue)
end

typealias Stroke Property{StrokePrimitive}

function stroke(c::Nothing)
    return Stroke([StrokePrimitive(c)])
end

function stroke(c::Union(ColorValue, String))
	return Stroke([StrokePrimitive(color(c))])
end

function stroke(cs::AbstractArray)
	return Stroke([StrokePrimitive(c == nothing ? c : color(c)) for c in cs])
end

# Fill
# ----

immutable FillPrimitive <: PropertyPrimitive
	color::Union(ColorValue, AlphaColorValue, Nothing)
end

typealias Fill Property{FillPrimitive}

function fill(c::Nothing)
    return Fill([FillPrimitive(c)])
end

function fill(c::Union(ColorValue, AlphaColorValue, String))
	return Fill([FillPrimitive(color(c))])
end

function fill(cs::AbstractArray)
	return Fill([FillPrimitive(c == nothing ? c : color(c)) for c in cs])
end

# StrokeDash
# ----------

immutable StrokeDashPrimitive <: PropertyPrimitive
    value::Vector{Measure}
end

typealias StrokeDash Property{StrokeDashPrimitive}

function strokedash(values::AbstractArray)
    return StrokeDash([StrokeDashPrimitive(collect(Measure, values))])
end

function strokedash(values::AbstractArray{AbstractArray})
    return StrokeDash([StrokeDashPrimitive(collect(Measure, value)) for value in values])
end

function absolute_units(primitive::StrokeDashPrimitive, t::Transform,
                        units::UnitBox, box::AbsoluteBoundingBox)
    return StrokeDashPrimitive([Measure(abs=absolute_units(v, t, units, box))
                                for v in primitive.value])
end

# StrokeLineCap
# -------------

abstract LineCap
immutable LineCapButt <: LineCap end
immutable LineCapSquare <: LineCap end
immutable LineCapRound <: LineCap end

immutable StrokeLineCapPrimitive <: PropertyPrimitive
    value::LineCap

    function StrokeLineCapPrimitive(value::LineCap)
        return new(value)
    end

    function StrokeLineCapPrimitive(value::Type{LineCap})
        return new(value())
    end
end

typealias StrokeLineCap Property{StrokeLineCapPrimitive}

function strokelinecap(value::Union(LineCap, Type{LineCap}))
    return StrokeLineCap([StrokeLineCapPrimitive(value)])
end

function strokelinecap(values::AbstractArray)
    return StrokeLineCap([StrokeLineCapPrimitive(value) for value in values])
end

# StrokeLineJoin
# --------------

abstract LineJoin
immutable LineJoinMiter <: LineJoin end
immutable LineJoinRound <: LineJoin end
immutable LineJoinBevel <: LineJoin end

immutable StrokeLineJoinPrimitive <: PropertyPrimitive
    value::LineJoin

    function StrokeLineJoinPrimitive(value::LineJoin)
        return new(value)
    end

    function StrokeLineCapPrimitive(value::Type{LineJoin})
        return new(value())
    end
end

typealias StrokeLineJoin Property{StrokeLineJoinPrimitive}

function strokelinejoin(value::Union(LineJoin, Type{LineJoin}))
    return StrokeLineJoin([StrokeLineJoinPrimitive(value)])
end

function strokelinejoin(values::AbstractArray)
    return StrokeLineJoin([StrokeLineJoinPrimitive(value) for value in values])
end

# LineWidth
# ---------

immutable LineWidthPrimitive <: PropertyPrimitive
    value::Measure

    function LineWidthPrimitive(value)
        return new(size_measure(value))
    end
end

typealias LineWidth Property{LineWidthPrimitive}

function linewidth(value::Union(Measure, Number))
    return LineWidth([LineWidthPrimitive(value)])
end

function linewidth(values::AbstractArray)
    return LineWidth([LineWidthPrimitive(value) for value in values])
end

function absolute_units(primitive::LineWidthPrimitive, t::Transform,
                        units::UnitBox, box::AbsoluteBoundingBox)
    return LineWidthPrimitive(Measure(abs=absolute_units(primitive.value, t, units, box)))
end

# Visible
# -------

immutable VisiblePrimitive <: PropertyPrimitive
    value::Bool
end

typealias Visible Property{VisiblePrimitive}

function visible(value::Bool)
    return Visible([VisiblePrimitive(value)])
end

function visible(values::AbstractArray)
    return Visible([VisiblePrimitive(value) for value in values])
end

# FillOpacity
# -----------

immutable FillOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function FillOpacityPrimitive(value_::Number)
        value = float64(value_)
        if value < 0.0 || value > 1.0
            error("Opacity must be between 0 and 1.")
        end
        return new(value)
    end
end

typealias FillOpacity Property{FillOpacityPrimitive}

function fillopacity(value::Float64)
    return FillOpacity([FillOpacityPrimitive(value)])
end

function fillopacity(values::AbstractArray)
    return FillOpacity([FillOpacityPrimitive(value) for value in values])
end

# StrokeOpacity
# -------------

immutable StrokeOpacityPrimitive <: PropertyPrimitive
    value::Float64

    function StrokeOpacityPrimitive(value_::Number)
        value = float64(value_)
        if value < 0.0 || value > 1.0
            error("Opacity must be between 0 and 1.")
        end
        return new(value)
    end
end

typealias StrokeOpacity Property{StrokeOpacityPrimitive}

function strokeopacity(value::Float64)
    return StrokeOpacity([StrokeOpacityPrimitive(value)])
end

function strokeopacity(values::AbstractArray)
    return StrokeOpacity([StrokeOpacityPrimitive(value) for value in values])
end

# Clip
# ----

immutable ClipPrimitive <: PropertyPrimitive
    points::Vector{Point}
end

typealias Clip Property{ClipPrimitive}

function clip()
    return Clip([ClipPrimitive(Array(Point, 0))])
end

function clip(points::XYTupleOrPoint...)
    return Clip([ClipPrimitive([convert(Point, point) for point in points])])
end

function clip(point_arrays::AbstractArray...)
    clipprims = Array(ClipPrimitive, length(point_arrays))
    for (i, point_array) in enumerate(point_arrays)
        clipprims[i] = ClipPrimitive([convert(Point, point)
                                      for point in point_array])
    end
    return Clip(clipprims)
end

function absolute_units(primitive::ClipPrimitive, t::Transform, units::UnitBox,
                        box::AbsoluteBoundingBox)
    return ClipPrimitive([absolute_units(point, t, units, box)
                          for point in primitive.points])
end

# Font
# ----

immutable FontPrimitive <: PropertyPrimitive
    family::String
end

typealias Font Property{FontPrimitive}

function font(family::String)
    return Font([FontPrimitive(family)])
end

function font(families::AbstractArray)
    return Font([FontPrimitive(family) for family in families])
end

# FontSize
# --------

immutable FontSizePrimitive <: PropertyPrimitive
    value::Measure

    function FontSizePrimitive(value)
        return new(size_measure(value))
    end
end

typealias FontSize Property{FontSizePrimitive}

function fontsize(value::Union(Number, Measure))
    return FontSize([FontSizePrimitive(value)])
end

function fontsize(values::AbstractArray)
    return FontSize([FontSizePrimitive(value) for value in values])
end

function absolue_units(primitive::FontSizePrimitive, t::Transform,
                       units::UnitBox, box::AbsoluteBoundingBox)
    return FontSizePrimitive(Measure(abs=absolute_units(primitive.value, t, units, box)))
end

# SVGID
# -----

immutable SVGIDPrimitive <: PropertyPrimitive
    value::String
end

typealias SVGID Property{SVGIDPrimitive}

function svgid(value::String)
    return SVGID([SVGIDPrimitive(value)])
end

function svgid(values::AbstractArray)
    return SVGID([SVGIDPrimitive(value) for value in values])
end

# SVGClass
# --------

immutable SVGClassPrimitive <: PropertyPrimitive
    value::String
end

typealias SVGClass Property{SVGClassPrimitive}

function svgclass(value::String)
    return SVGClass([SVGClassPrimitive(value)])
end

function svgclass(values::AbstractArray)
    return SVGClass([SVGClassPrimitive(value) for value in values])
end

# SVGAttribute
# ------------

immutable SVGAttributePrimitive <: PropertyPrimitive
    attribute::String
    value::String
end

typealias SVGAttribute Property{SVGAttributePrimitive}

function svgattribute(attribute::String, value)
    return SVGAttribute([SVGAttributePrimitive(attribute, string(value))])
end

function svgattribute(attribute::String, values::AbstractArray)
    return SVGAttribute([SVGAttributePrimitive(attribute, string(value))
                         for value in values])
end

function svgattribute(attributes::AbstractArray, values::AbstractArray)
    return SVGAttribute([SVGAttributePrimitive(attribute, string(value))
                         for (attribute, value) in cyclezip(attributes, values)])
end

# JSInclude
# ---------

immutable JSIncludePrimitive <: PropertyPrimitive
    value::String
end

typealias JSInclude Property{JSIncludePrimitive}

function jsinclude(value::String)
    return JSInclude([JSIncludePrimitive(value)])
end

# Don't bother with a vectorized version of this. It wouldn't really make #
# sense.

# JSCall
# ------

immutable JSCallPrimitive <: PropertyPrimitive
    code::String
    args::Vector{Measure}
end

typealias JSCall Property{JSCallPrimitive}

function jscall(code::String, arg::Vector{Measure}=Measure[])
    return JSCall([JSCallPrimitive(code, arg)])
end

function jscall(codes::AbstractArray,
                args::AbstractArray{Vector{Measure}}=Vector{Measure}[Measure[]])
    return JSCall([JSCallPrimitive(code, arg)
                   for (code, arg) in cyclezip(codes, args)])
end

function absolute_units(primitive::JSCallPrimitive, t::Transform,
                        units::UnitBox, box::AbsoluteBoundingBox)
    # we are going to build a new string by scanning across "code" and
    # replacing %x with translated x values, %y with translated y values
    # and %s with translated size values.
    newcode = IOBuffer()

    i = 1
    validx = 1
    while true
        j = search(primitive.code, '%', i)

        if j == 0
            write(newcode, primitive.code[i:end])
            break
        end

        write(newcode, primitive.code[i:j-1])
        if j == length(primitive.code)
            write(newcode, '%')
            break
        elseif primitive.code[j+1] == '%'
            write(newcode, '%')
        elseif primitive.code[j+1] == 'x'
            val = absolute_x_position(primitive.args[validx], t, units, box)
            write(newcode, svg_fmt_float(val))
            validx += 1
        elseif primitive.code[j+1] == 'y'
            val = absolute_y_position(primitive.args[validx], t, units, box)
            write(newcode, svg_fmt_float(val))
            validx += 1
        elseif primitive.code[j+1] == 's'
            val = absolute_units(primitive.args[validx], t, units, box)
            write(newcode, svg_fmt_float(val.abs))
            validx += 1
        else
            write(newcode, '%', primitive.code[j+1])
        end

        i = j + 2
    end

    return JSCallPrimitive(takebuf_string(newcode), Measure[])
end

function isrepeatable(p::JSCall)
    return true
end


text:
using JLProject
using Test

@testset "JLProject.jl" begin
    # Write your tests here.
end


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

kernelabstractions = Pkg.PackageSpec(path = root_directory)

BACKEND = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", "all")
CI = parse(Bool, get(ENV, "CI", "false"))

@info "Develop..." CI

Pkg.develop(kernelabstractions)
if !CI || BACKEND == "ROCM"
    rockernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "ROCKernels"))
    Pkg.develop(rockernels)
end

if !CI || BACKEND == "CUDA"
    cudakernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "CUDAKernels"))
    Pkg.develop(cudakernels)
end

if VERSION < v"1.8"
    kernelgradients = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "KernelGradients"))
    Pkg.develop(kernelgradients)
end
Pkg.build()
Pkg.precompile()


text:
# This file was generated, do not modify it. # hide
forest.n = 150;


text:
using Raytracing

using GeometryBasics
using Random
using CoordinateTransformations
using ColorTypes
using ColorVectorSpace
using StaticArrays
using Images
using ImageIO
using FileIO
using ImageMagick

function main()
    T = Float64
    origin = Point3{T}(0,0,1.5)
    lookat = Point3{T}(0,0,0)
    vup = Vec3{T}(0,1,0)
    aspectratio = T(16/9)
    vfov = T(40)
    camera = Raytracing.make_camera(origin, lookat, vup, vfov, aspectratio)

    #spheresdf = Raytracing.SphereSignedDistanceField{T}(0.5)

    mandelsdf = Raytracing.MandelBulbSignedDistanceField(50)
    boundingmandelspheresdf = Raytracing.SphereSignedDistanceField{T}(1.25)
    overshootfactor = T(1.5)
    boundingmandelsdf = Raytracing.BoundingVolumeSignedDistanceField(overshootfactor, boundingmandelspheresdf, mandelsdf)

    point = Point3{T}(0, 0, 1)
    mandelsdf(point)
    Raytracing.normalforwarddiff(mandelsdf, point)

    normalmat = Raytracing.NormalShading()

    id = IdentityTransformation()
    mandelalbedo = Vec{3, T}(0.8, 0.25, 0.25)

    #lambertianmat = Raytracing.Lambertian(mandelalbedo)

    #mandelshape = Raytracing.Shape(id, mandelsdf, lambertianmat)
    mandelshape = Raytracing.Shape(id, boundingmandelsdf, normalmat)
    #mandelshape = Raytracing.Shape(id, mandelsdf, normalmat)

    skyboxspheresdf = Raytracing.SphereSignedDistanceField{T}(100)
    skyboxsphereshape = Raytracing.Shape(id, skyboxspheresdf, normalmat)

    sceneshapes = (skyboxsphereshape, mandelshape)
    scene = Raytracing.Scene(sceneshapes)

    image_height = 100
    samples_per_pixel = 3
    max_bounces_per_ray = 1
    max_steps_per_bounce = 400
    distance_tolerance = 1e-4
    rng = MersenneTwister(1234)
    rat = 2.0 / 0.70
    len = 1.2
    camera3 = Raytracing.make_camera(Point3{T}(0.0, rat * len, 1/rat * len), lookat, vup, vfov, aspectratio)
    @time image = Raytracing.raytrace_image(scene, camera3, image_height, samples_per_pixel, max_bounces_per_ray, max_steps_per_bounce, distance_tolerance, rng)

    savename = "lofi_mandelbulb_50_normalshading.png"
    lofibig = Raytracing.supersample_crisp(image, 10)
    Raytracing.saveimage("$(savename)", lofibig)

end


text:
module SegmentationModels

using Reexport
using Flux
@reexport using Metalhead

const Classifier = Union{ResNet, VGG}

include("encoders.jl")
include("unet.jl")

export UNet

end


text:
using Stipple
using Test

@testset "Stipple.jl" begin

    @test 1 == 1

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ADIOS2"
version = v"2.7.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ornladios/ADIOS2/archive/refs/tags/v2.7.1.tar.gz", "c8e237fd51f49d8a62a0660db12b72ea5067512aa7970f3fcf80b70e3f87ca3e"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd ADIOS2-2.7.1
# See <https://github.com/ornladios/ADIOS2/issues/2705>
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/gettid.patch
# PR <https://github.com/ornladios/ADIOS2/pull/2712>
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ndims.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/shlwapi.patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/sockaddr_in.patch
mkdir build
cd build
if [[ "$target" == *-apple-* ]]; then
    # Set up a wrapper script for the assembler. GCC's assembly output
    # isn't accepted by the LLVM assembler.
    as=$(which as)
    mv "$as" "$as.old"
    export AS="$as.old"
    ln -s "$WORKSPACE/srcdir/scripts/as.llvm" "$as"
fi
archopts=
if [[ "$target" == x86_64-w64-mingw32 ]]; then
    # - The MSMPI Fortran bindings are missing a function; see
    #   <https://github.com/microsoft/Microsoft-MPI/issues/7>
    echo 'void __guard_check_icall_fptr(unsigned long ptr) {}' >cfg_stub.c
    gcc -c cfg_stub.c
    ar -crs libcfg_stub.a cfg_stub.o
    cp libcfg_stub.a $prefix/lib
    # - cmake's auto-detection for MPI doesn't work on Windows.
    # - The SST and Table ADIOS2 components don't build on Windows
    #   (reported in <https://github.com/ornladios/ADIOS2/issues/2705>)
    export FFLAGS="-I$prefix/src -I$prefix/include -fno-range-check"
    archopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI -DMPI_C_LIBRARIES=msmpi64 -DMPI_CXX_LIBRARIES=msmpi64 -DMPI_Fortran_LIBRARIES='msmpifec64;msmpi64;cfg_stub' -DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
elif [[ "$target" == *-mingw* ]]; then
    archopts="-DMPI_HOME=$prefix -DMPI_GUESS_LIBRARY_NAME=MSMPI -DADIOS2_USE_SST=OFF -DADIOS2_USE_Table=OFF"
fi
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN%.*}_gcc.cmake -DCMAKE_BUILD_TYPE=Release -DADIOS2_BUILD_EXAMPLES=OFF -DBUILD_TESTING=OFF ${archopts} ..
make -j${nproc}
make -j${nproc} install
install_license ../Copyright.txt ../LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# 32-bit architectures are not supported; see
# <https://github.com/ornladios/ADIOS2/issues/2704>
platforms = filter(p -> nbits(p) ≠ 32, platforms)
# Apparently, macOS doesn't use different C++ string APIs
platforms = expand_cxxstring_abis(platforms; skip=Sys.isapple)
# TODO: Windows doesn't build with libcxx="cxx03"
platforms = expand_gfortran_versions(platforms)
# x86_64-apple-darwin-libgfortran3 encounters an ICE in GCC
platforms = filter(p -> !(Sys.isapple(p) && libgfortran_version(p) == v"3"), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libadios2_c", :libadios2_c),
    LibraryProduct("libadios2_c_mpi", :libadios2_c_mpi),
    LibraryProduct("libadios2_core", :libadios2_core),
    LibraryProduct("libadios2_core_mpi", :libadios2_core_mpi),
    LibraryProduct("libadios2_cxx11", :libadios2_cxx11),
    LibraryProduct("libadios2_cxx11_mpi", :libadios2_cxx11_mpi),
    LibraryProduct("libadios2_fortran", :libadios2_fortran),
    LibraryProduct("libadios2_fortran_mpi", :libadios2_fortran_mpi),
    LibraryProduct("libadios2_taustubs", :libadios2_taustubs),

    # Missing on Windows:
    # LibraryProduct("libadios2_atl", :libadios2_atl),
    # LibraryProduct("libadios2_dill", :libadios2_dill),
    # LibraryProduct("libadios2_evpath", :libadios2_evpath),
    # LibraryProduct("libadios2_ffs", :libadios2_ffs),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Blosc_jll")),
    # We don't want to use Bzip2 because this would lock us into Julia ≥1.6
    # Dependency(PackageSpec(name="Bzip2_jll")),
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae")),
    Dependency(PackageSpec(name="MPICH_jll")),
    Dependency(PackageSpec(name="MicrosoftMPI_jll")),
    Dependency(PackageSpec(name="ZeroMQ_jll")),
    Dependency(PackageSpec(name="libpng_jll")),
    Dependency(PackageSpec(name="zfp_jll")),
    # We cannot use HDF5 because we need a HDF5 configuration with MPI support
    # Dependency(PackageSpec(name="HDF5_jll")),
]

# Build the tarballs, and possibly a `build.jl` as well.
# GCC 4 is too old for Windows; it doesn't have <regex.h>
# GCC 5 is too old for FreeBSD; it doesn't have `std::to_string`
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")

max_stars_repo_name: K-JW/Yggdrasil

max_issues_repo_name: K-JW/Yggdrasil

max_forks_repo_name: K-JW/Yggdrasil


text:
using Distributions
using Statistics

@testset "Radial Uniform" begin
    dist = RadialUniform(0, 1)

    @testset "Interface" begin
        @test minimum(dist) == 0
        @test maximum(dist) == 1

        t = RadialUniform(0, 1)
        @test typeof(t) === RadialUniform{Float64}
        t = RadialUniform(0, 1.0)
        @test typeof(t) === RadialUniform{Float64}
        t = RadialUniform(BigInt(0), BigInt(1))
        @test typeof(t) === RadialUniform{BigFloat}
    end

    @testset "Accuracy" begin
        @test pdf(dist, -1) ≈ 0
        @test pdf(dist, 0.5) ≈ 1
        @test pdf(dist, 2) ≈ 0

        @test logpdf(dist, -1) ≈ -Inf
        @test logpdf(dist, 0.5) ≈ 0
        @test logpdf(dist, 2) ≈ -Inf

        @test cdf(dist, -1) ≈ 0
        @test cdf(dist, 0.5) ≈ 1 / 4
        @test cdf(dist, 2) ≈ 1

        @test logcdf(dist, -1) ≈ -Inf
        @test logcdf(dist, 0.5) ≈ -log(4)
        @test logcdf(dist, 2) ≈ 0

        @test quantile(dist, 1) ≈ 1
        @test quantile(dist, 0) ≈ 0

        for i in range(0, 1, length = 10)
            @test cdf(dist, quantile(dist, i)) ≈ i
        end
    end
end

@testset "Poisson Invariant" begin
    dist = PoissonInvariant(0, 1)

    @testset "Interface" begin
        @test minimum(dist) == 0
        @test maximum(dist) == 1

        t = PoissonInvariant(0, 1)
        @test typeof(t) === PoissonInvariant{Float64}
        t = PoissonInvariant(0, 1.0)
        @test typeof(t) === PoissonInvariant{Float64}
        t = PoissonInvariant(BigInt(0), BigInt(1))
        @test typeof(t) === PoissonInvariant{BigFloat}
    end

    @testset "Accuracy" begin
        @test pdf(dist, -1) ≈ 0
        @test pdf(dist, 0.5) ≈ 1 / 2 / √0.5
        @test pdf(dist, 2) ≈ 0

        @test logpdf(dist, -1) ≈ -Inf
        @test logpdf(dist, 0.5) ≈ -log(2 * √0.5)
        @test logpdf(dist, 2) ≈ -Inf

        @test cdf(dist, -1) ≈ 0
        @test cdf(dist, 0.5) ≈ √0.5
        @test cdf(dist, 2) ≈ 1

        @test logcdf(dist, -1) ≈ -Inf
        @test logcdf(dist, 0.5) ≈ -log(2) / 2
        @test logcdf(dist, 2) ≈ 0

        @test quantile(dist, 1) ≈ 1
        @test quantile(dist, 0) ≈ 0

        for i in range(0, 1, length = 10)
            @test cdf(dist, quantile(dist, i)) ≈ i
        end
    end
end


text:
using DiscreteElasticRods
DER = DiscreteElasticRods
using Zygote

# Basic linalg
n = 40
a = randn(3,n)
b = randn(3,n)
c = randn(3,n)

linalg3_test(a,b,c) = sum(DER.dot3(c, DER.cross3(a,b))) + sum(DER.norm3(b))

linalg3_test(a,b,c)
out, back = Zygote.pullback(linalg3_test, a, b, c)
back(rand())

# Parallel transport
t1 = rand(3,n); DER.norm3!(t1)
t2 = rand(3,n); DER.norm3!(t2)
v = rand(3,n);
out, back = Zygote.pullback(DER.ptransport, v, t1, t2)
back(rand(3,n))

# Rotations
θ = rand(1,n)
out, back = Zygote.pullback(DER.rotate_orthogonal_unit, t1, t2, θ)
back(rand(3,n))

# Kinematics
r = DER.random_rod(n)

DER.full_kinematics(r.x, r.d)
out, back = Zygote.pullback(DER.full_kinematics, r.x, r.d)
out, back = Zygote.pullback(DER.full_kinematics, r)
back(out)

# Test case: the entire energy function
ref_rod = DER.copy(r)
ref_strains = DER.full_kinematics(r)
rod_props = DER.elliptical_stiffness(10, 0.3, 0.01, 0.04)

function total_energy(rod_coords::DER.rod_delta)
    new_rod = DER.rod_update(ref_rod, rod_coords)
    strains = DER.full_kinematics(new_rod)
    DER.elastic_energy(ref_strains, strains, rod_props)
end

rod_coord = DER.rod_delta(rand(3,n), rand(1,n-1))
out, back = Zygote.pullback(total_energy, rod_coord)

grad_total_energy(rod_coords::DER.rod_delta) = Zygote.gradient(total_energy, rod_coords)

using BenchmarkTools
@btime grad_total_energy($rod_coord)


text:
using BasisFunctions, GridArrays, DomainSets
using BasisFunctions: FFreq, DiscreteMeasure, DiscreteProductMeasure,
    innerproduct_fourier_part, default_dict_innerproduct, AbstractMeasure,
	OuterProductArray
import BasisFunctions: quadweights, stencilarray, name, element, elements,
    iscomposite, weights, grid, subindices, isnormalized, supermeasure,
    support, unsafe_weight, strings, innerproduct_native, restrict, discretemeasure
export submeasure, DiscreteTensorSubMeasure, SubMeasure, DiscreteSubMeasure

"""
    struct SubMeasure{M,D,T} <: Measure{T}

Submeasure of a Measure.
"""
struct SubMeasure{M<:Measure,D<:Domain,T} <: Measure{T}
    measure     ::  M
    domain      ::  D

    function SubMeasure(measure::Measure, domain::Domain)
        @assert eltype(domain) == eltype(support(measure))
        new{typeof(measure), typeof(domain), eltype(domain)}(measure, domain)
    end
end

"""
    submeasure(measure::Measure, domain::Domain)

Restrict a measure to a subset of its support
"""
submeasure(measure::Measure, domain::Domain) = SubMeasure(measure, domain)
name(m::SubMeasure) = "Restriction of a measure"
supermeasure(measure::SubMeasure) = measure.measure
support(measure::SubMeasure) = measure.domain
unsafe_weight(m::SubMeasure, x) = unsafe_weight(supermeasure(m), x)

"""
    restrict(measure::Measure, domain::Domain)

Restrict a measure to a subset of its support
"""
restrict(measure::Measure, domain::Domain) = submeasure(measure, domain)
strings(m::SubMeasure) = (name(m), (string(support(m)),), strings(supermeasure(m)))
submeasure(measure::ProductMeasure, domain::ProductDomain) = ProductMeasure(map(submeasure, elements(measure), elements(domain))...)

quadweights(grid::AbstractSubGrid, measure::SubMeasure) =
    quadweights(supergrid(grid), supermeasure(measure))[subindices(grid)]

function innerproduct_native(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq, m::SubMeasure{<:FourierMeasure}; options...)
	d = support(m)
	if d isa AbstractInterval
		innerproduct_fourier_part(b1, i, b2, j, infimum(d), supremum(d))
	else
		default_dict_innerproduct(b1, i, b2, j, m; options...)
	end
end

"""
    struct DiscreteSubMeasure{T,M<:DiscreteMeasure{T},G<:AbstractGrid} <: DiscreteMeasure{T}

Submeasure of a discrete measure.
"""
struct DiscreteSubMeasure{T,M<:DiscreteMeasure{T},G<:AbstractGrid} <: DiscreteMeasure{T}
    supermeasure   :: M
    subgrid        :: G
    DiscreteSubMeasure(measure::DiscreteMeasure{T}, grid::AbstractGrid) where T = new{T,typeof(measure),typeof(grid)}(measure, grid)
    DiscreteSubMeasure(measure::DiscreteSubMeasure{T}, grid::AbstractGrid) where T = error()
end

submeasure(measure::DiscreteMeasure, domain::Domain) = _discretesubmeasure(subgrid(points(measure), domain), weights(measure))
subindices(measure::DiscreteSubMeasure) = subindices(measure.subgrid)
supermeasure(measure::DiscreteSubMeasure) = measure.supermeasure
discretemeasure(grid::Union{AbstractSubGrid,TensorSubGrid}) = DiscreteSubMeasure(discretemeasure(supergrid(grid)), grid)
_discretesubmeasure(grid::Union{AbstractSubGrid,TensorSubGrid},weights) = DiscreteSubMeasure(discretemeasure(supergrid(grid),weights), grid)
restrict(measure::DiscreteMeasure, domain::Domain) = DiscreteSubMeasure(measure, subgrid(points(measure), domain))
weights(measure::DiscreteSubMeasure) = subweights(measure, subindices(measure), weights(supermeasure(measure)))
subweights(_, subindices, w) = w[subindices]
points(measure::DiscreteSubMeasure) = measure.subgrid
isnormalized(::DiscreteSubMeasure) = false
name(m::DiscreteSubMeasure) = "Restriction of a "*name(supermeasure(m))

"""
    const DiscreteTensorSubMeasure{T,G,W} = DiscreteSubMeasure{T,M,G} where {T,M<:DiscreteProductMeasure,G<:TensorSubGrid}

A tensor product of discrete submeasures.
"""
const DiscreteTensorSubMeasure{T,G,W} = DiscreteSubMeasure{T,M,G} where {T,M<:DiscreteProductMeasure,G<:TensorSubGrid}
name(m::DiscreteTensorSubMeasure) = "Tensor of submeasures (supermeasure:"*name(supermeasure(m))
elements(m::DiscreteTensorSubMeasure) = map(_discretesubmeasure,elements(points(m)),elements(weights(supermeasure(m))))
element(m::DiscreteTensorSubMeasure, i) = _discretesubmeasure(element(points(m),i),element(weights(supermeasure(m)),i))
iscomposite(m::DiscreteTensorSubMeasure) = true
weights(m::DiscreteTensorSubMeasure) = OuterProductArray(map(weights, elements(m))...)

function stencilarray(m::DiscreteTensorSubMeasure)
    A = Any[]
    push!(A, element(m,1))
    for i = 2:length(elements(m))
        push!(A," ⊗ ")
        push!(A, element(m,i))
    end
    A
end