text:
"""
Convert a two-bit genotype (of ALT alleles) to a real number
of type `t` according to specified SNP model.
"""
function convert_gt(
    t::Type{T},
    a::NTuple{2, Bool},
    model::Symbol = :additive
    ) where T <: Real
    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

"""
    copy_gt!(A, reader; ...)

Fill the rows of matrix `A` by the GT data from VCF records in `reader`.
Records with missing genotypes are converted to `missing`.

# Input
- `A`: a matrix or vector such that `eltype(A) <: Union{Missing, Real}`. Each
    row is a person's genotype.
- `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`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.

# 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,
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real
    msg != "" && (pmeter = Progress(size(A, 2), 5, msg))
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    record = VCF.Record()
    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Reached end of reader; columns $j-$(size(A, 2)) are " *
                "set to missing values")
            fill!(view(A, :, j:size(A, 2)), missing)
            break
        else
            read!(reader, record)
        end
        gtkey = 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
        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))
        # second pass: impute, convert, center, scale
        _, _, _, _, _, alt_freq, _, _, _, _, _ = gtstats(record, nothing)
        ct = 2alt_freq
        wt = alt_freq == 0 ? 1.0 : 1.0 / √(2alt_freq * (1 - alt_freq))
        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() ≤ alt_freq, rand() ≤ alt_freq
                    A[i, j] = convert_gt(T, (a1, a2), 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), 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

        # update progress
        msg != "" && next!(pmeter)
    end
    return A
end

"""
    copy_gt_trans!(A, reader; ...)

Fill the rows of matrix `A` by the GT data from VCF records in `reader` where
the ALT allele in each record is interpreted as a `1`. Records with missing
genotypes are converted to `missing`.

# Input
- `A`: a matrix or vector such that `eltype(A) <: Union{Missing, Real}`. Each
    column is one person's genotype.
- `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`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.

# Output
- `A`: `ismissing(A[i, j]) == true` indicates missing genotype. If `impute=true`,
    `ismissing(A[i, j]) == false` for all entries.
"""
function copy_gt_trans!(
    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,
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real
    msg != "" && (pmeter = Progress(size(A, 1), 5, msg))
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    record = VCF.Record()

    for j in 1:size(A, 1)
        if eof(reader)
            @warn("Reached end of reader; rows $j-$(size(A, 1)) are set to " *
                "missing values")
            fill!(view(A, j:size(A, 1), :), missing)
            break
        else
            read!(reader, record)
        end
        gtkey = findgenokey(record, "GT")
        # if no GT field, fill by missing values
        if gtkey === nothing
            fill!(view(A, j, :), missing)
        end
        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))
        # second pass: impute, convert, center, scale
        _, _, _, _, _, alt_freq, _, _, _, _, _ = gtstats(record, nothing)
        ct = 2alt_freq
        wt = alt_freq == 0 ? 1.0 : 1.0 / √(2alt_freq * (1 - alt_freq))
        for i in 1:size(A, 2)
            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() ≤ alt_freq, rand() ≤ alt_freq
                    A[j, i] = convert_gt(T, (a1, a2), model)
                else
                    A[j, i] = 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[j, i] = convert_gt(T, (a1, a2), model)
            end
            # center and scale if asked
            center && !ismissing(A[j, i]) && (A[j, i] -= ct)
            scale && !ismissing(A[j, i]) && (A[j, i] *= wt)
        end
        # update progress
        msg != "" && next!(pmeter)
    end
    return nothing
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. Records with missing genotypes are converted to `missing`.

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

# 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`
- `trans`: whether to import data transposed so that each column is 1 genotype,
    default `false`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `save_snp_info`: Boolean. If true, will also output sample ID, chrom, pos,
    snp id, ref, and alt vectors

# Output
- `out`: The genotype matrix stored in column or rows depending on `trans`.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function convert_gt(
    t::Type{T},
    vcffile::AbstractString;
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false,
    trans::Bool = false,
    msg::String = "",
    save_snp_info::Bool = false
    ) where T <: Real
    reader = VCF.Reader(openvcf(vcffile, "r"))
    records = nrecords(vcffile)
    samples = nsamples(vcffile)
    sampleID   = (save_snp_info ? Vector{String}(undef, samples) : nothing)
    record_chr = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_pos = (save_snp_info ? zeros(Int, records) : nothing)
    record_ids = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)
    record_ref = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_alt = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)

    if trans
        out = Matrix{Union{t, Missing}}(undef, records, samples)
        copy_gt_trans!(out, reader; model = model, impute = impute,
            center = center, scale = scale, msg = msg, sampleID=sampleID,
            record_chr=record_chr, record_pos=record_pos, record_ids=record_ids,
            record_ref=record_ref, record_alt=record_alt)
    else
        out = Matrix{Union{t, Missing}}(undef, samples, records)
        copy_gt!(out, reader; model = model, impute = impute,
            center = center, scale = scale, msg = msg, sampleID=sampleID,
            record_chr=record_chr, record_pos=record_pos, record_ids=record_ids,
            record_ref=record_ref, record_alt=record_alt)
    end
    close(reader)
    if save_snp_info
        return out, sampleID, record_chr, record_pos, record_ids, record_ref,
            record_alt
    else
        return out
    end
end

"""
    convert_ht(t, vcffile)

Converts the GT data from a VCF file to a haplotype matrix of type `t`. One
VCF record will become 2 column/row in the resulting matrix. Records cannot
have missing genotypes and the target file must be phased.

# Input
- `t`: a type `t <: Real`. If `t` is `Bool`, output will be a BitMatrix.
- `vcffile`: a phased VCF file
- `trans`: whether to import data transposed so that each column is a
    haplotype, default `false`.
- `msg`: A message that will be printed to indicate import progress. Defaults to
    not printing.
- `save_snp_info`: Boolean. If true, will also output sample ID, chrom, pos,
    snp id, ref, and alt vectors

# Output
- `out`: The haplotype matrix stored in column or rows depending on `trans`.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function convert_ht(
    t::Type{T},
    vcffile::AbstractString;
    trans::Bool = false,
    msg::String = "",
    save_snp_info::Bool = false
    ) where T <: Real
    reader = VCF.Reader(openvcf(vcffile, "r"))
    records = nrecords(vcffile)
    samples = nsamples(vcffile)
    sampleID   = (save_snp_info ? Vector{String}(undef, samples) : nothing)
    record_chr = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_pos = (save_snp_info ? zeros(Int, records) : nothing)
    record_ids = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)
    record_ref = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_alt = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)

    M = (t == Bool ? BitArray{2} : Matrix{t})
    if trans
        out = M(undef, records, 2samples)
        copy_ht_trans!(out, reader, msg = msg, sampleID=sampleID,
            record_chr=record_chr, record_pos=record_pos, record_ids=record_ids,
            record_ref=record_ref, record_alt=record_alt)
    else
        out = M(undef, 2samples, records)
        copy_ht!(out, reader, msg = msg, sampleID=sampleID,
            record_chr=record_chr, record_pos=record_pos, record_ids=record_ids,
            record_ref=record_ref, record_alt=record_alt)
    end
    close(reader)
    if save_snp_info
        return out, sampleID, record_chr, record_pos, record_ids, record_ref,
            record_alt
    else
        return out
    end
end

"""
    copy_ht!(A, reader)

Fill 2 columns of `A` by the GT data from VCF records in `reader`, each record
filling 2 columns. Missing data is not allowed and genotypes must be phased.

# Input
- `A`: matrix where rows are haplotypes. Person `i`'s haplotype are filled in
    rows `2i - 1` and `2i`.
- `reader`: a VCF reader

# Optional inputs
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function copy_ht!(
    A::Union{AbstractMatrix{T}, AbstractVector{T}},
    reader::VCF.Reader;
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real

    n, p = size(A)
    nn   = n >>> 1
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    msg != "" && (pmeter = Progress(p, 5, msg)) # update every 5 seconds
    record = VCF.Record()

    for j in 1:p
        if eof(reader)
            @warn("Reached end of record! Columns $j-$p are filled " *
                "with $(zero(T))s and are NOT haplotypes!")
            fill!(view(A, :, j:size(A, 2)), zero(T))
            break
        else
            read!(reader, record)
        end
        gtkey = 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

        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))

        # 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 $(2i - 1). " *
                    "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

        # update progress
        msg != "" && next!(pmeter)
    end
    return A
end

"""
    copy_ht_trans!(A, reader)

Fill 2 rows of `A` by the GT data from VCF records in `reader`, each record
filling 2 columns. Missing data is not allowed and genotypes must be phased.

# Input
- `A`: matrix where columns are haplotypes. Person `i`'s haplotype are filled in
    rows `2i - 1` and `2i`.
- `reader`: a VCF reader

# Optional inputs
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function copy_ht_trans!(
    A::Union{AbstractMatrix{T}, AbstractVector{T}},
    reader::VCF.Reader;
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real

    p, n = size(A)
    nn   = n >>> 1
    msg != "" && (pmeter = Progress(p, 5, msg)) # update every 5 seconds
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    record = VCF.Record()

    for j in 1:p
        if eof(reader)
            @warn("Reached end of record! Rows $j-$p are filled with "
                * "$(zero(T))s and are NOT haplotypes!")
            fill!(view(A, j:size(A, 2)), zero(T))
            break
        else
            read!(reader, record)
        end
        gtkey = findgenokey(record, "GT")

        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))

        # 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 $(2i - 1). " *
                    "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[j, 2i - 1] = convert(T, a1)
                A[j, 2i] = convert(T, a2)
            end
        end

        # update progress
        msg != "" && next!(pmeter)
    end
    return nothing
end

"""
    convert_ds(t, vcffile; [key = "DS"], ...)

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").

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

# Optional argument
- `key`: The FIELD name of the VCF file that encodes dosages, defaults to `"DS"`
- `model`: genetic model `:additive` (default), `:dominant`, or `:recessive`
- `impute`: impute missing genotype with 2 times the ALT allele frequency 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`
- `trans`: whether to import data transposed so that each column is 1 genotype,
    default `false`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.

# Output
- `out`: Dosage matrix
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function convert_ds(
    t::Type{T},
    vcffile::AbstractString;
    key::String = "DS",
    model::Symbol = :additive,
    impute::Bool = false,
    center::Bool = false,
    scale::Bool = false,
    trans::Bool = false,
    save_snp_info::Bool=false,
    msg::String = ""
    ) where T <: Real
    reader = VCF.Reader(openvcf(vcffile, "r"))
    records = nrecords(vcffile)
    samples = nsamples(vcffile)
    sampleID   = (save_snp_info ? Vector{String}(undef, samples) : nothing)
    record_chr = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_pos = (save_snp_info ? zeros(Int, records) : nothing)
    record_ids = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)
    record_ref = (save_snp_info ? Vector{String}(undef, records) : nothing)
    record_alt = (save_snp_info ? Vector{Vector{String}}(undef, records) : nothing)

    if trans
        out = Matrix{Union{Missing, t}}(undef, records, samples)
        copy_ds_trans!(out, reader, key = key, impute = impute, center = center,
            scale = scale, msg = msg, sampleID=sampleID, record_chr=record_chr,
            record_pos=record_pos, record_ids=record_ids, record_ref=record_ref,
            record_alt=record_alt)
    else
        out = Matrix{Union{Missing, t}}(undef, samples, records)
        copy_ds!(out, reader, key = key, impute = impute, center = center,
            scale = scale, msg = msg, sampleID=sampleID, record_chr=record_chr,
            record_pos=record_pos, record_ids=record_ids, record_ref=record_ref,
            record_alt=record_alt)
    end
    close(reader)
    if save_snp_info
        return out, sampleID, record_chr, record_pos, record_ids, record_ref,
            record_alt
    else
        return out
    end
end

"""
    copy_ds!(A, reader; [key="DS"], ...)

Fill the columns of matrix `A` by the dosage data from VCF records in `reader`.
Records with missing fields (i.e. `key` field is `.`) are converted to
`missing`.

# Input
- `A`: matrix where rows are dosages.
- `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`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
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,
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real
    msg != "" && (pmeter = Progress(size(A, 2), 5, msg)) # update every 5 seconds
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    record = VCF.Record()
    for j in 1:size(A, 2)
        if eof(reader)
            @warn("Reached end of reader; columns $j-$(size(A, 2)) are set " *
                "to missing values")
            fill!(view(A, :, j:size(A, 2)), missing)
            break
        else
            read!(reader, record)
        end
        dskey = 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

        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))

        # loop over every marker in record
        _, _, _, _, _, alt_freq, _, _, _, _, _ = gtstats(record, nothing)
        ct = 2alt_freq
        wt = alt_freq == 0 ? 1.0 : 1.0 / √(2alt_freq * (1 - alt_freq))
        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

        # update progress
        msg != "" && next!(pmeter)
    end
end

"""
    copy_ds_trans!(A, reader; [key="DS"], ...)

Fill the rows of matrix `A` by the dosage data from VCF records in `reader`.
Records with missing fields (i.e. `key` field is `.`) are converted to `missing`.

# Input
- `A`: matrix where columns are dosages.
- `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`
- `msg`: A message that will be printed to indicate progress. Defaults to
    not printing.
- `sampleID`: The sample ID for each person
- `record_chr`: The chromosome number for each record
- `record_pos`: The position for each record
- `record_ids`: The record ID for each record. Only the first one is
    recorded when there are multiple IDs.
- `record_ref`: The ref allele for each record
- `record_alt`: The alternate allele for each record. Only the first one is
    recorded when there are multiple alternate alleles.
"""
function copy_ds_trans!(
    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,
    msg::String = "",
    sampleID::Union{AbstractVector, Nothing} = nothing,
    record_chr::Union{AbstractVector, Nothing} = nothing,
    record_pos::Union{AbstractVector, Nothing} = nothing,
    record_ids::Union{AbstractVector, Nothing} = nothing,
    record_ref::Union{AbstractVector, Nothing} = nothing,
    record_alt::Union{AbstractVector, Nothing} = nothing
    ) where T <: Real
    msg != "" && (pmeter = Progress(size(A, 1), 5, msg)) # update every 5 seconds
    sampleID === nothing || (sampleID .= header(reader).sampleID)
    record = VCF.Record()
    for j in 1:size(A, 1)
        if eof(reader)
            @warn("Reached end of reader; rows $j-$(size(A, 1)) " *
                "are set to missing values")
            fill!(view(A, j:size(A, 2), :), missing)
            break
        else
            read!(reader, record)
        end
        dskey = findgenokey(record, key)

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

        # save record's information
        record_chr === nothing || (record_chr[j] = VCF.chrom(record))
        record_pos === nothing || (record_pos[j] = VCF.pos(record))
        record_ids === nothing || (record_ids[j] = try VCF.id(record) catch; ["."] end)
        record_ref === nothing || (record_ref[j] = VCF.ref(record))
        record_alt === nothing || (record_alt[j] = VCF.alt(record))

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

        #update progress
        msg != "" && next!(pmeter)
    end
end


text:
"""
$(TYPEDEF)

Interpolant-basis weighted sum-of-squares polynomial (of dimension `U`) positive
semidefinite matrix (of side dimension `R`) cone, parametrized by vector of
matrices `Ps` derived from interpolant basis and polynomial domain constraints.

    $(FUNCTIONNAME){T}(R::Int, U::Int, Ps::Vector{Matrix{T}}, use_dual::Bool = false)
"""
mutable struct WSOSInterpPosSemidefTri{T <: Real} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    R::Int
    U::Int
    Ps::Vector{Matrix{T}}
    nu::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}
    use_hess_prod_slow::Bool
    use_hess_prod_slow_updated::Bool

    rt2::T
    rt2i::T
    tempU::Vector{T}
    tempLRLR::Vector{Symmetric{T, Matrix{T}}}
    tempLRLR2::Vector{Matrix{T}}
    tempLRUR::Vector{Matrix{T}}
    ΛFL::Vector
    ΛFLP::Vector{Matrix{T}}
    tempLU::Vector{Matrix{T}}
    PΛiP::Matrix{T}
    PΛiP_blocks_U
    Ps_times::Vector{Float64}
    Ps_order::Vector{Int}

    function WSOSInterpPosSemidefTri{T}(
        R::Int,
        U::Int,
        Ps::Vector{Matrix{T}};
        use_dual::Bool = false,
        ) where {T <: Real}
        for Pk in Ps
            @assert size(Pk, 1) == U
        end
        cone = new{T}()
        cone.use_dual_barrier = !use_dual # using dual barrier
        cone.dim = U * svec_length(R)
        cone.R = R
        cone.U = U
        cone.Ps = Ps
        cone.nu = R * sum(size(Pk, 2) for Pk in Ps)
        return cone
    end
end

reset_data(cone::WSOSInterpPosSemidefTri) = (cone.feas_updated =
    cone.grad_updated = cone.hess_updated = cone.inv_hess_updated =
    cone.hess_fact_updated = cone.use_hess_prod_slow =
    cone.use_hess_prod_slow_updated = false)

function setup_extra_data!(cone::WSOSInterpPosSemidefTri{T}) where {T <: Real}
    U = cone.U
    R = cone.R
    Ps = cone.Ps
    K = length(Ps)
    cone.rt2 = sqrt(T(2))
    cone.rt2i = inv(cone.rt2)
    cone.tempU = zeros(T, U)
    Ls = [size(Pk, 2) for Pk in cone.Ps]
    cone.tempLRLR = [Symmetric(zeros(T, L * R, L * R), :L) for L in Ls]
    cone.tempLRLR2 = [zeros(T, L * R, L * R) for L in Ls]
    cone.tempLRUR = [zeros(T, L * R, U * R) for L in Ls]
    cone.tempLU = [zeros(T, L, U) for L in Ls]
    cone.ΛFL = Vector{Any}(undef, K)
    cone.ΛFLP = [zeros(T, R * L, R * U) for L in Ls]
    cone.PΛiP = zeros(T, R * U, R * U)
    cone.PΛiP_blocks_U = [view(cone.PΛiP, block_idxs(U, r), block_idxs(U, s))
        for r in 1:R, s in 1:R]
    cone.Ps_times = zeros(K)
    cone.Ps_order = collect(1:K)
    return cone
end

function set_initial_point!(arr::AbstractVector, cone::WSOSInterpPosSemidefTri)
    arr .= 0
    block = 1
    @inbounds for i in 1:cone.R
        @views arr[block_idxs(cone.U, block)] .= 1
        block += i + 1
    end
    return arr
end

function update_feas(cone::WSOSInterpPosSemidefTri)
    @assert !cone.feas_updated

    # order the Ps by how long it takes to check feasibility, to improve efficiency
    sortperm!(cone.Ps_order, cone.Ps_times, initialized = true) # stochastic

    cone.is_feas = true
    for k in cone.Ps_order
        cone.Ps_times[k] = @elapsed @inbounds begin
            Pk = cone.Ps[k]
            LU = cone.tempLU[k]
            L = size(Pk, 2)
            Λ = cone.tempLRLR[k]

            @views for p in 1:cone.R, q in 1:p
                @. cone.tempU = cone.point[block_idxs(cone.U, svec_idx(p, q))]
                if p != q
                    cone.tempU .*= cone.rt2i
                end
                mul!(LU, Pk', Diagonal(cone.tempU)) # TODO check efficiency
                mul!(Λ.data[block_idxs(L, p), block_idxs(L, q)], LU, Pk)
            end

            ΛFLk = cone.ΛFL[k] = cholesky!(Λ, check = false)
            if !isposdef(ΛFLk)
                cone.is_feas = false
                break
            end
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::WSOSInterpPosSemidefTri{T}) where {T}
    # condition is necessary but not sufficient for dual feasibility
    block = 1
    @inbounds for i in 1:cone.R
        @views diag_i = cone.dual_point[block_idxs(cone.U, block)]
        all(>(eps(T)), diag_i) || return false
        block += i + 1
    end
    return true
end

function update_grad(cone::WSOSInterpPosSemidefTri)
    @assert is_feas(cone)
    U = cone.U
    R = cone.R
    grad = cone.grad
    cone.grad .= 0

    @inbounds for k in eachindex(cone.Ps)
        L = size(cone.Ps[k], 2)
        ΛFL = cone.ΛFL[k].L
        ΛFLP = cone.ΛFLP[k]

        # given cholesky L factor ΛFL, get ΛFLP = ΛFL \ kron(I, P')
        for p in 1:R
            block_U_p_idxs = block_idxs(U, p)
            block_L_p_idxs = block_idxs(L, p)
            @views ΛFLP_pp = ΛFLP[block_L_p_idxs, block_U_p_idxs]
            # ΛFLP_pp = ΛFL_pp \ P'
            @views ldiv!(ΛFLP_pp, LowerTriangular(
                ΛFL[block_L_p_idxs, block_L_p_idxs]), cone.Ps[k]')
            # to get off-diagonals in ΛFLP, subtract known blocks aggregated in ΛFLP_qp
            for q in (p + 1):R
                block_L_q_idxs = block_idxs(L, q)
                @views ΛFLP_qp = ΛFLP[block_L_q_idxs, block_U_p_idxs]
                ΛFLP_qp .= 0
                for p2 in p:(q - 1)
                    block_L_p2_idxs = block_idxs(L, p2)
                    @views mul!(ΛFLP_qp, ΛFL[block_L_q_idxs, block_L_p2_idxs],
                        ΛFLP[block_L_p2_idxs, block_U_p_idxs], -1, 1)
                end
                @views ldiv!(LowerTriangular(ΛFL[block_L_q_idxs,
                    block_L_q_idxs]), ΛFLP_qp)
            end
        end

        # update grad
        block_diag_prod!(grad, ΛFLP, ΛFLP, cone)
    end
    grad .*= -1

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::WSOSInterpPosSemidefTri{T}) where {T <: Real}
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    H = cone.hess.data
    R = cone.R
    U = cone.U
    PΛiP_blocks = cone.PΛiP_blocks_U
    H .= 0

    @inbounds for k in eachindex(cone.Ps)
        L = size(cone.Ps[k], 2)
        # PΛiP = ΛFLP' * ΛFLP
        ΛFLP = cone.ΛFLP[k]
        for p in 1:R, q in p:R
            # since ΛFLP is block lower triangular rows only from max(p,q)
            # start making a nonzero contribution to the product
            row_range = ((q - 1) * L + 1):(L * R)
            @views mul!(PΛiP_blocks[p, q], ΛFLP[row_range, block_idxs(U, p)]',
                ΛFLP[row_range, block_idxs(U, q)])
        end
        LinearAlgebra.copytri!(cone.PΛiP, 'U')

        for p in 1:R, q in 1:p
            block = svec_idx(p, q)
            idxs = block_idxs(U, block)

            for p2 in 1:R, q2 in 1:p2
                block2 = svec_idx(p2, q2)
                if block2 < block
                    continue
                end
                idxs2 = block_idxs(U, block2)

                @views Hview = H[idxs, idxs2]
                scal = (xor(p == q, p2 == q2) ? cone.rt2 : one(T))
                @. Hview += PΛiP_blocks[p, p2] * PΛiP_blocks[q, q2] * scal
                if (p != q) && (p2 != q2)
                    @inbounds @. Hview += PΛiP_blocks[p, q2] * PΛiP_blocks[q, p2]
                end
            end
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::WSOSInterpPosSemidefTri,
    )
    cone.use_hess_prod_slow_updated || update_use_hess_prod_slow(cone)
    @assert cone.hess_updated
    cone.use_hess_prod_slow || return hess_prod!(prod, arr, cone)
    return partial_prod!(prod, arr, false, cone)
end

function dder3(cone::WSOSInterpPosSemidefTri, dir::AbstractVector)
    @assert cone.grad_updated
    return partial_prod!(cone.dder3, dir, true, cone)
end

# diagonal from each (i, j) block in mat1' * mat2
function block_diag_prod!(
    vect::AbstractVector{T},
    mat1::Matrix{T},
    mat2::Matrix{T},
    cone::WSOSInterpPosSemidefTri{T},
    ) where T
    U = cone.U
    @inbounds for u in 1:U
        idx = u
        j_idx = u
        for j in 1:cone.R
            i_idx = u
            for i in 1:(j - 1)
                @views vect[idx] += dot(mat1[:, i_idx], mat2[:, j_idx]) * cone.rt2
                idx += U
                i_idx += U
            end
            @views vect[idx] += dot(mat1[:, j_idx], mat2[:, j_idx])
            j_idx += U
            idx += U
        end
    end
    return
end

function partial_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    use_symm_prod::Bool,
    cone::WSOSInterpPosSemidefTri,
    )
    @assert cone.grad_updated
    prod .= 0
    U = cone.U
    R = cone.R
    tempU = cone.tempU

    @inbounds for (k, Pk) in enumerate(cone.Ps)
        Lk = size(Pk, 2)
        LUk = cone.tempLU[k]
        LRURk = cone.tempLRUR[k]
        LRLRk = cone.tempLRLR2[k]
        ΛFLk = cone.ΛFL[k]
        ΛFLPk = cone.ΛFLP[k]
        left_prod = (use_symm_prod ? LRURk : ΛFLPk)

        @views for j in 1:size(arr, 2)
            delta = arr[:, j]

            for q in 1:R, p in 1:q
                @. tempU = delta[block_idxs(U, svec_idx(q, p))]
                if p != q
                    # svec scaling
                    tempU .*= cone.rt2i
                end
                mul!(LUk, Pk', Diagonal(tempU))
                mul!(LRLRk[block_idxs(Lk, p), block_idxs(Lk, q)], LUk, Pk)
            end

            LinearAlgebra.copytri!(LRLRk, 'U')
            ldiv!(ΛFLk.L, LRLRk)
            rdiv!(LRLRk, ΛFLk.L')
            mul!(LRURk, Symmetric(LRLRk), ΛFLPk)

            block_diag_prod!(prod[:, j], left_prod, LRURk, cone)
        end
    end

    return prod
end

ext: jl


text:
"""This module implements utility functions."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

function dispatch_index(N::Int, Nprocs::Int)
    r = N % Nprocs
    q = div(N, Nprocs)
    a = cumsum((q+1)*ones(r))
    b = a - q
    if length(a) > 0
        c = cumsum(q*ones(Nprocs-r))+a[end]
    else
        c = cumsum(q*ones(Nprocs-r))
    end
    d = c- (q-1)
    indmax::Array{Int,1} = vcat(a, c)
    indmin::Array{Int,1} = vcat(b, d)
    return indmin, indmax
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
using PyPlot
using LinearAlgebra
using Statistics
using Random
using NLopt
using UKrig

##
# data grid and interpolation grid
n₀ = 500 # grid evals
n₁ = 200  # obs points
x₁ = sort(rand(n₁))
x₀ = range(-0.1, 1.1, length=n₀)
x₀₀ = range(0, 1, length=n₀)

##
# set the simulation truth parameters
σεₒ = 0.05
σzₒ = 1.0
ρₒ  = 0.1
νₒ  = 1.25
mₒ  = 2
βₒ  = 2 * (rand(mₒ) .- 1/2)

##
# simulate data
Σ₁ = (σzₒ^2) .* Mnu.(UKrig.ℓ2_dist.(x₁,permutedims(x₁)) ./ ρₒ, νₒ)
L₁ = cholesky(Σ₁).L
Z₁ = L₁ * randn(n₁)
Y₁ = Z₁ .+ βₒ[1] .+ x₁ .* βₒ[2]
data = Y₁ + σεₒ .* randn(n₁)

## ===============================
#

musr = 2
ν    = 1.25

loglike_inv_nu1, loglike_chol_nu1, wgrad_loglike_inv_nu1 = loglike_Gnu(
	data, x₁;
	musr=musr,
	ν=ν, # ≈ differentiability
)

#=
loglike_inv_nu1(0.1, 2.0)
loglike_chol_nu1(0.1, 2.0)
wgrad_loglike_inv_nu1(0.1, 2.0)
ϵ = 1e-7
(loglike_inv_nu1(0.1+ϵ, 2.0) - loglike_inv_nu1(0.1-ϵ, 2.0))/(2ϵ)
(loglike_inv_nu1(0.1, 2.0+ϵ) - loglike_inv_nu1(0.1, 2.0-ϵ))/(2ϵ)
=#

#LN_algm = [:LN_BOBYQA, :LN_COBYLA, :LN_PRAXIS, :LN_NELDERMEAD, :LN_SBPLX]
#LD_algm = [:LD_MMA, :LD_SLSQP, :LD_LBFGS, :LD_TNEWTON]
opt1 = Opt(:LN_BOBYQA, 2) # second arg is the number of variables to optimize
opt2 = Opt(:LN_BOBYQA, 2) # second arg is the number of variables to optimize
opt3 = Opt(:LD_LBFGS, 2) # second arg is the number of variables to optimize
opt1.max_objective = (x,grad) -> loglike_inv_nu1(x[1],x[2])
opt2.max_objective = (x,grad) -> loglike_chol_nu1(x[1],x[2])
opt3.max_objective = function (x,grad)
	ll, grad1, grad2 = wgrad_loglike_inv_nu1(x[1],x[2])
	if length(grad)>0
		grad[1]=grad1
		grad[2]=grad2
	end
	return ll
end
opt1.lower_bounds = [0.0, 0.0]
opt2.lower_bounds = [0.0, 0.0]
opt3.lower_bounds = [0.0, 0.0]
opt1.maxtime = 10 # in seconds
opt2.maxtime = 10
opt3.maxtime = 10
# opt1.maxeval = 50
# opt2.maxeval = 50
# opt3.maxeval = 50
# opt1.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
# opt2.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
# opt3.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
xinit1 = [0.1, 2.0]
xinit2 = [0.1, 2.0]
xinit3 = [0.1, 2.0]
optf1, optx1, ret1 = optimize!(opt1, xinit1)
optf2, optx2, ret2 = optimize!(opt2, xinit2)
optf3, optx3, ret3 = optimize!(opt3, xinit3)
@show sqrt(length(data))/√2

σe = optx1[1]
σg = optx1[2]

rng_σe = σe |> x->range(x - .05 * x, x + .05 * x, length=50)
rng_σg = σg |> x->range(x - .15 * x, x + .15 * x, length=50)
llmat = loglike_chol_nu1.(rng_σe,rng_σg')
pcolor(rng_σe, rng_σg, llmat)

krig1 = krig_Gnu(
	data, x₁;
	musr=musr,
	ν=ν,
	σe=σe,
	σg=σg,
)

figure(figsize=(9,6))
plot(x₁, data, ".", label="data")
plot(x₀, krig1.(x₀))
legend()

## ===============================
# do the interpolation, vary signal/noise
figure(figsize=(9,6))
plot(x₁, data, ".", label="data")
for σg ∈ [1,5,10]
	krig1 = krig_Gnu(
		data, x₁;
		musr=2,
		ν=1.0,    # ≈ differentiability
		σe=σe, # ≈ noise/signal ratio
		σg=σg,    # ≈ signal/noise ratio
	)
	plot(x₀, krig1.(x₀), label = "signal/noise = $σg")
end
legend()

# do the interpolation, vary smoothness
figure(figsize=(9,6))
plot(x₁, data, ".", label="data")
for ν ∈ [0.25, .5, 1.0, 1.25, 2.0]
	krig1 = krig_Gnu(
		data, x₁;
		musr=2,
		ν=ν,   # ≈ differentiability
		σe=σe, # ≈ noise/signal ratio
		σg=σg, # ≈ signal/noise ratio
	)
	plot(x₀, krig1.(x₀), label = "nu = $ν")
end
legend()

# do the interpolation, interpolate
figure(figsize=(9,6))
plot(x₁, data, ".", label="data")
for ν ∈ [0.25, .5, 1.0]
	krig1 = krig_Gnu(
		data, x₁;
		musr=2,
		ν=ν,    # ≈ differentiability
		σe=1e-7, # interpolate
		σg=σg, # ≈ signal/noise ratio
	)
	plot(x₀₀, krig1.(x₀₀), label = "nu = $ν")
end
legend()


text:
# Physical constants in SI units
const qₑ = -1.60217662e-19  # electron charge, [C]
const mₑ = 9.10938356e-31   # electron mass, [kg]
const qᵢ = 1.60217662e-19   # proton mass, [C]
const mᵢ = 1.673557546e-27  # proton mass, [kg]
const c  = 3e8              # speed of light, [m/s]
const μ₀ = 4π*1e-7          # Vacuum permeability, [H/m]
const Rₑ = 6.38e6           # Earth radius, [m]
const BMoment_Earth = [0.0, 0.0, 7.94e22] # [V*s/(A*m)]


text:
using DiffEqJump
dims = (5,5)
num_nodes = prod(dims) # number of sites
grid = CartesianGrid(dims) # or use LightGraphs.grid(dims)

num_species = 3
starting_state = zeros(Int, num_species, num_nodes)
starting_state[1,1] = 25
starting_state[2,end] = 25
starting_state

tspan = (0.0, 3.0)
rates = [6.0, 0.05] # k_1 = rates[1], k_2 = rates[2]

prob = DiscreteProblem(starting_state, tspan, rates)

netstoch = [[1 => -1, 2 => -1, 3 => 1],[1 => 1, 2 => 1, 3 => -1]]
reactstoch = [[1 => 1, 2 => 1],[3 => 1]]
majumps = MassActionJump(rates, reactstoch, netstoch)

hopping_constants = ones(num_species, num_nodes)
hopping_constants[3, :] .= 0.0
hopping_constants

alg = NSM()
jump_prob = JumpProblem(prob, alg, majumps, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(true, false))

solution = solve(jump_prob, SSAStepper())

using Plots
is_static(spec) = (spec == 3) # true if spec does not hop
"get frame k"
function get_frame(k, sol, linear_size, labels, title)
    num_species = length(labels)
    h = 1/linear_size
    t = sol.t[k]
    state = sol.u[k]
    xlim=(0,1+3h/2); ylim=(0,1+3h/2);
    plt = plot(xlim=xlim, ylim=ylim, title = "$title, $(round(t, sigdigits=3)) seconds")

    species_seriess_x = [[] for i in 1:num_species]
    species_seriess_y = [[] for i in 1:num_species]
    CI = CartesianIndices((linear_size, linear_size))
    for ci in CartesianIndices(state)
        species, site = Tuple(ci)
        x,y = Tuple(CI[site])
        num_molecules = state[ci]
        sizehint!(species_seriess_x[species], num_molecules)
        sizehint!(species_seriess_y[species], num_molecules)
        if !is_static(species)
            randsx = rand(num_molecules)
            randsy = rand(num_molecules)
        else
            randsx = zeros(num_molecules)
            randsy = zeros(num_molecules)
        end
        for k in 1:num_molecules
            push!(species_seriess_x[species], x*h - h/4 + 0.5h*randsx[k])
            push!(species_seriess_y[species], y*h - h/4 + 0.5h*randsy[k])
        end
    end
    for species in 1:num_species
        scatter!(plt, species_seriess_x[species], species_seriess_y[species], label = labels[species], marker = 6)
    end
    xticks!(plt, range(xlim...,length = linear_size+1))
    yticks!(plt, range(ylim...,length = linear_size+1))
    xgrid!(plt, 1, 0.7)
    ygrid!(plt, 1, 0.7)
    return plt
end

"make an animation of solution sol in 2 dimensions"
function animate_2d(sol, linear_size; species_labels, title, verbose = true)
    num_frames = length(sol.t)
    anim = @animate for k=1:num_frames
        verbose && println("Making frame $k")
        get_frame(k, sol, linear_size, species_labels, title)
    end
    anim
end
# animate
anim=animate_2d(solution, 5, species_labels = ["A", "B", "C"], title = "A + B <--> C", verbose = false)
fps = 5
gif(anim, fps = fps)

dims = (2,3,4) # can pass in a 1-Tuple, a 2-Tuple or a 3-Tuple
num_nodes = prod(dims)
grid = CartesianGrid(dims)

using LightGraphs
graph = cycle_digraph(5) # directed cyclic graph on 5 nodes

hopping_constants = Matrix{Vector{Float64}}(undef, num_species, num_nodes)
for ci in CartesianIndices(hopping_constants)
    (species, site) = Tuple(ci)
    hopping_constants[species, site] = zeros(outdegree(grid, site))
    for (n, nb) in enumerate(neighbors(grid, site))
        if nb < site
            hopping_constants[species, site][n] = 1.0
        end
    end
end

species_hop_constants = ones(num_species)
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = ones(outdegree(grid, site))
end
hopping_constants=Pair(species_hop_constants, site_hop_constants)

species_hop_constants = ones(num_species, num_nodes)
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
for site in 1:num_nodes
    site_hop_constants[site] = ones(outdegree(grid, site))
end
hopping_constants=Pair(species_hop_constants, site_hop_constants)

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
using Word2Vec
using LinearAlgebra

function loadYears(fileName::String)::Vector{Int}
    file = open(fileName)
    years = Vector{Int}()
    for line in readlines(file)
        if line == ""
            continue
        end
        fields = split(line, ". ")
        try
        year = parse(Int, fields[1])
        push!(years, year)
        catch
            println(line)
        end
    end
    close(file)
    years
end

function generateTemporalTrainFile(inFileName::String, outFileName::String, years::Vector{Int}, toYear::Int)
    fileIn = open(inFileName)
    fileOut = open(outFileName, "w")
    iLine = 1
    for line in readlines(fileIn)
        year = years[iLine]
        if year <= toYear
            println(fileOut, line)
        end
        iLine += 1
    end
    close(fileIn)
    close(fileOut)
end

function findHighFreqPhrases(fileName::String, freqThres::Int)::Vector{String}
    file = open(fileName)
    phraseFreqs = Dict{String, Int}()
    for line in readlines(file)
        phrases = split(line)
        for phrase in phrases
            phraseFreqs[phrase] = get(phraseFreqs, phrase, 0) + 1
        end
    end
    close(file)

    highFreqPhrases = Vector{String}()
    for (phrase, freq) in phraseFreqs
        if freq >= freqThres
            push!(highFreqPhrases, phrase)
        end
    end

    highFreqPhrases
end

function loadNormalizedVectors(fileName::String, highFreqPhrases::Vector{String})::Array{Float64, 2}
    wordVectors = wordvectors(fileName)
    numDims, vocabSize = size(wordVectors)
    numVecs = length(highFreqPhrases)
    vectors = Array{Float64}(undef, (numDims, numVecs))
    for (id, phrase) in enumerate(highFreqPhrases)
        if in_vocabulary(wordVectors, phrase)
            vector = get_vector(wordVectors, phrase)
            w = 1.0 / norm(vector)
            vectors[:,i] = w .* vector
        else
            vectors[:,i] = 0.0
        end
    end
    vectors
end

years = loadYears("DualAspectTopicRepresentation/data/NIPS.txt")
for toYear in 2010:2020
    generateTemporalTrainFile("DualAspectTopicRepresentation/workspace/NIPS.w2v.training.txt", "DualAspectTopicRepresentation/workspace/NIPS.w2v.$toYear.txt", years, toYear)
    word2vec("DualAspectTopicRepresentation/workspace/NIPS.w2v.$toYear.txt", "DualAspectTopicRepresentation/workspace/NIPS.w2v.model-$toYear.txt", size=200, iter=100, min_count=3, verbose=true)
end

highFreqPhrases = findHighFreqPhrases("DualAspectTopicRepresentation/workspace/NIPS.w2v.2015.txt", 5)

prevRankMat = nothing
for toYear in [2010, 2020]
    vectors = loadNormalizedVectors("DualAspectTopicRepresentation/workspace/NIPS.w2v.$toYear.txt", highFreqPhrases)
    simMat = 0.5 .+ 0.5 .* (vectors' * vectors)
    simVec = simMat[:]
    simRank = sortperm(simVec, rev=true)
    rankMat = Array{Int,2}(undef, size(simMat))
    for (i, k) in enumerate(simRank)
        rankMat[k] = i
    end
    if prevRankMat !== nothing
        deltaRankMat = rankMat - prevRankMat
        deltaRankVec = deltaRankMat[:]
        jumpRank = sortperm(deltaRankVec, rev=true)
        saveTopJumps
    end
    prevRankMat = rankMat
end


text:
@enum XC_Type begin
  xc_type_lda
  xc_type_gga
  xc_type_mgga
  xc_type_hybrid
end

struct XC_Functional
  kind::XC_Type
  xc_func::Function
end

function xc_func(kind::XC_Type = xc_type_lda, spin_polarization::Bool = false,
                 exchange::Function = LDA_X, correlation::Function = LDA_C_VWN, threshold = eps())
  function xc_lda(ρ, vxc = similar(ρ), exc = similar(ρ))
    if !(length(ρ) == length(vxc) == length(exc))
      throw(DimensionMismatch("Dimensions of the arrays ρ, vxc, and exc doesn't match"))
    end
    @inbounds for i in eachindex(ρ)
      vxc[i], exc[i] = ρ[i] ≤ threshold ? (0, 0) : exchange(ρ[i]) .+ correlation(ρ[i])
    end
    return vxc, exc
  end

  if kind == xc_type_lda && spin_polarization == false
    return XC_Functional(kind, xc_lda)
  else
    return nothing
  end
end


text:
package java.lang;

public class IndexOutOfBoundsException extends RuntimeException {
    public IndexOutOfBoundsException() { }
    public IndexOutOfBoundsException(String s) { super(s); }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Test

m = page_model()
include("../src/components/LGcycle.jl")

add_comp!(m, LGcycle)

set_param!(m, :LGcycle, :e_globalLGemissions, readpagedata(m,"test/validationdata/e_globalLGemissions.csv"))
set_param!(m, :LGcycle, :y_year, [2009.,2010.,2020.,2030.,2040.,2050.,2075.,2100.,2150.,2200.]) #real value
set_param!(m, :LGcycle, :y_year_0, 2008.) #real value
set_param!(m, :LGcycle, :rtl_g_landtemperature, readpagedata(m,"test/validationdata/rtl_g_landtemperature.csv"))

p=load_parameters(m)
set_leftover_params!(m,p) #important for setting left over component values

# run Model
run(m)

conc=m[:LGcycle,  :c_LGconcentration]
conc_compare=readpagedata(m,"test/validationdata/c_LGconcentration.csv")

@test conc ≈ conc_compare rtol=1e-4


text:
include("../utils.jl")
include("group_utils.jl")

using Manifolds: invariant_metric_dispatch, default_metric_dispatch

@testset "Circle group" begin
    G = CircleGroup()
    @test repr(G) == "CircleGroup()"

    @test base_manifold(G) === Circle{ℂ}()

    @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
    @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(true)
    @test (@inferred Manifolds.biinvariant_metric_dispatch(G)) === Val(true)
    @test (@inferred default_metric_dispatch(MetricManifold(G, EuclideanMetric()))) ===
          Val(true)
    @test has_invariant_metric(G, LeftAction())
    @test has_invariant_metric(G, RightAction())
    @test has_biinvariant_metric(G)
    @test is_default_metric(MetricManifold(G, EuclideanMetric()))

    @testset "identity overloads" begin
        ig = Identity(G)
        @test inv(G, ig) === ig
        q = [1.0 * im]
        X = [Complex(0.5)]
        @test translate_diff(G, ig, q, X) === X

        @test identity_element(G) === 1.0
        @test identity_element(G, 1.0f0) === 1.0f0
        @test identity_element(G, [1.0f0]) == [1.0f0]
    end

    @testset "scalar points" begin
        pts = [1.0 + 0.0im, 0.0 + 1.0im, (1.0 + 1.0im) / √2]
        Xpts = [0.0 + 0.5im, 0.0 - 1.5im]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] * pts[1]
        @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ pts[2] * Xpts[1]
        test_group(
            G,
            pts,
            Xpts,
            Xpts;
            test_diff=true,
            test_mutating=false,
            test_invariance=true,
            test_lie_bracket=true,
            test_adjoint_action=true,
        )
    end

    @testset "vector points" begin
        pts = [[1.0 + 0.0im], [0.0 + 1.0im], [(1.0 + 1.0im) / √2]]
        Xpts = [[0.0 + 0.5im], [0.0 - 1.5im]]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] .* pts[1]
        @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ pts[2] .* Xpts[1]
        test_group(
            G,
            pts,
            Xpts,
            Xpts;
            test_diff=true,
            test_mutating=true,
            test_invariance=true,
            test_lie_bracket=true,
            test_adjoint_action=true,
        )
    end

    @testset "Group forwards to decorated" begin
        pts = [[1.0 + 0.0im], [0.0 + 1.0im], [(1.0 + 1.0im) / √2]]
        test_manifold(
            G,
            pts,
            basis_types_to_from=(Manifolds.VeeOrthogonalBasis(), DefaultOrthonormalBasis()),
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=false,
            test_project_tangent=true,
            test_musical_isomorphisms=false,
            test_default_vector_transport=true,
            is_mutating=true,
            exp_log_atol_multiplier=2.0,
            is_tangent_atol_multiplier=2.0,
        )
    end
end

@testset "Real circle group" begin
    G = RealCircleGroup()
    @test repr(G) == "RealCircleGroup()"

    @test base_manifold(G) === Circle{ℝ}()

    @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
    @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(true)
    @test (@inferred Manifolds.biinvariant_metric_dispatch(G)) === Val(true)
    @test (@inferred default_metric_dispatch(MetricManifold(G, EuclideanMetric()))) ===
          Val(true)
    @test has_invariant_metric(G, LeftAction())
    @test has_invariant_metric(G, RightAction())
    @test has_biinvariant_metric(G)
    @test is_default_metric(MetricManifold(G, EuclideanMetric()))

    @testset "identity overloads" begin
        ig = Identity(G)
        @test inv(G, ig) === ig
        q = [0.0]
        X = [0.5]
        @test translate_diff(G, ig, q, X) === X

        @test identity_element(G) === 0.0
        @test identity_element(G, 1.0f0) === 0.0f0
        @test identity_element(G, [0.0f0]) == [0.0f0]
    end

    @testset "scalar points" begin
        pts = [1.0, 0.5, -3.0]
        Xpts = [-2.0, 0.5, 2.0]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] + pts[1]
        @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ Xpts[1]
        test_group(
            G,
            pts,
            Xpts,
            Xpts;
            test_diff=true,
            test_mutating=false,
            test_invariance=true,
            test_lie_bracket=true,
            test_adjoint_action=true,
        )
    end

    @testset "vector points" begin
        pts = [[1.0], [0.5], [-3.0]]
        Xpts = [[-2.0], [0.5], [2.0]]
        @test compose(G, pts[2], pts[1]) ≈ pts[2] .+ pts[1]
        @test translate_diff(G, pts[2], pts[1], Xpts[1]) ≈ Xpts[1]
        test_group(
            G,
            pts,
            Xpts,
            Xpts;
            test_diff=true,
            test_mutating=true,
            test_invariance=true,
            test_lie_bracket=true,
            test_adjoint_action=true,
        )
    end

    @testset "Group forwards to decorated" begin
        pts = [[1.0], [0.5], [-3.0]]
        test_manifold(
            G,
            pts,
            basis_types_to_from=(Manifolds.VeeOrthogonalBasis(), DefaultOrthonormalBasis()),
            test_forward_diff=false,
            test_reverse_diff=false,
            test_vector_spaces=false,
            test_project_tangent=true,
            test_musical_isomorphisms=false,
            test_default_vector_transport=true,
            is_mutating=true,
            exp_log_atol_multiplier=2.0,
            is_tangent_atol_multiplier=2.0,
        )
    end
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 47bd5a91-c963-4fed-bbfd-2c3c28eaf366
using PlutoLinks: @use_task

# ╔═╡ b8b2c431-a18d-4b65-8207-8cc87e990ddf
using PlutoHooks: @use_state

# ╔═╡ 13946695-d509-4820-921d-9f1a23ea3ab0
using PlutoUI

# ╔═╡ 1ce83536-9746-11ec-342a-9db71b8eba53
md"""
# `@use_task`

```julia
@use_task(f::Function, dependencies::Vector{Any})::Task
```

`@use_task` is a macro to run a function asynchronously without blocking the Pluto cell execution. The created task can be used to provide updates to the notebook asynchronously. `@use_task` has several advantages over simply creating a `Task` in your cell because it manages the task lifecycle automatically.

## Explicit runs & dependencies

It terminates the previous task when a new task is created by sending an `InterruptException` to the previous task. The previous task will be automatically replaced or terminated in the following cases:

 - If the cell containing the `@use_task()` call is **explicitely** run. An explicit run means that the "Run cell" button was pressed for this cell or that the code for this cell was edited. If the cell is run because it is dependent on a cell which was explicitely run, this is not an explicit run and will not replace the task.

 - If one of the value in the `dependencies` array parameter has changed, the task will be replaced. This can be useful if the task has some initialization to do and must be reset when a parameter has changed.

Let's try using `@use_task` for simple function that sleeps for 3 seconds:
"""

# ╔═╡ 4f1b81f4-5262-4063-a8ce-1ce7da56a69b
@use_task(() -> sleep(3), [])

# ╔═╡ 360c71c5-5f4e-4118-93af-e8fd3c3f163f
md"""
When trying this yourself, you can see the task starting by being:

```
Task (runnable) @0x00007f7e61e136b0
```

and then after 3 seconds, the cell will update because the task is done:

```
Task (done) @0x00007f7e61e136b0
```

!!! note
	Notice that the cell is surrounded in blue! This means that the cell is using features from [PlutoHooks.jl](https://github.com/JuliaPluto/PlutoHooks.jl) and some of its state will be reset on explicit runs.

Now, let's try to bind the sleep duration to a julia variable instead of hardcoding 3 seconds:
"""

# ╔═╡ 994ab3fa-35e0-41a2-9f24-60b5024c8008
@bind sleeping_time PlutoUI.Slider(1:10, show_value=true)

# ╔═╡ f5b08153-58b2-45d6-8882-f593caba70f4
# We use the Julia do end syntax here to pass a function as the first argument
@use_task([]) do
	sleep(sleeping_time)
end

# ╔═╡ e7a2f49e-3551-4185-bc16-21ec5b2ac252
md"""
The task is not recreated when the `sleeping_time` value changes! This is because we are not doing an explicit run here. The cell is run because the value of `sleeping_time` changed. We need to explicitely add `sleeping_time` as a dependency to reset the task when it has a new value:
"""

# ╔═╡ af6b1697-d127-4865-8347-e80c73e27194
@use_task([sleeping_time]) do
	#      ^^^^^^^^^^^^^
	#      sleeping_time is added as an explicit dependency
	sleep(sleeping_time)
end

# ╔═╡ 0e8057e7-3f69-4ed2-9fe1-d5af8e7e50ef
md"""
With the added dependencies, the task is not re-created every time the value of `sleeping_time` changes!
"""

# ╔═╡ 285fbc70-aa43-44f7-af1d-689545945261
md"""
## Long running tasks & notebook updates

Very often, one wants to do more than sleeping (in a task at least) and run code until an event happens. Let's build a Julia version of `PlutoUI.Clock()` which is a widget that ticks counts from 1 to X.
"""

# ╔═╡ 562c408c-717e-4686-8c03-b315e72050fa
function wait_until_event()
	sleep(1.)
end

# ╔═╡ eb1ce3ae-47d9-436c-be84-8baa2d03879b
md"""
To trigger an update to the notebook from the current cell, the recommanded method is to use `@PlutoHooks.use_state()`. `@use_state` functions much like its React.js counterpart [`useState()`](https://reactjs.org/docs/hooks-state.html) (in Pluto, notebook cells play the role of React.js components) because it returns a tuple containing the state value and a callback to re-run the cell with an updated state. It takes only one argument which is the initial state value.
"""

# ╔═╡ 459f0704-5d59-4fde-bfa4-e3e2ff5d85c4
let
	# create a new state with a callback to update it
	count, set_count = @use_state(1)

	@use_task([]) do
		# capture the value of count in a new variable
		new_count = count

		# start the event loop
		while new_count < 10
			wait_until_event()
			new_count += 1
			set_count(new_count) # <- update the value of count (this will trigger a re-run)
		end
	end
	count
end

# ╔═╡ d9619966-e61a-40fa-bd54-f9efb512f14d
md"""
## Reusing hooks

For reusing code using hooks, the easiest thing is to create a macro returning the code:
"""

# ╔═╡ 63a79fa0-366c-4edc-a79d-6a4a78d74ef6
macro use_clock()
	quote
		let
			count, set_count = @use_state(1)
			@use_task([]) do
				# capture the value of count in a new variable
				new_count = count

				# start the event loop
				while new_count < 10
					wait_until_event()
					new_count += 1
					set_count(new_count) # <- update the value of count (this will trigger a re-run)
				end
			end
			count
		end
	end
end

# ╔═╡ b893ec96-90d4-4131-8380-f4f15bc5cf40
@use_clock()

# ╔═╡ 553818c9-2258-4993-a44d-17df498adcd1
md"""
## Common problems

##### 1. The hook does not work

Make sure that the macro call is not inside a function definition. Hooks only work if they are at the top level. To reuse hooks, package them as macros instead of functions.

##### 2. Some values become `nothing` in the task function after a while

Because of how Pluto works internally to redefine consts and structs, it may override the value of variables captured by the task by `nothing`. To prevent this from happening, the variables captured should be deepcopied into local variables before entering the task:

```julia
begin
	local new_variable = deepcopy(variable)
	@use_task([variable]) do
		# use only new_variable here...
		new_variable
	end
end
```

Have one yourself that is not on this list? [Open a discussion](https://github.com/JuliaPluto/PlutoLinks.jl/discussions) or don't hesitate to send a message on the pluto-development Zulip channel.
"""

# ╔═╡ 4bbd2db5-2a21-45e1-8f19-0a0f9a3b56fa
md"""
#### Packages
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoHooks = "0ff47ea0-7a50-410d-8455-4348d5de0774"
PlutoLinks = "0ff47ea0-7a50-410d-8455-4348d5de0420"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoHooks = "~0.0.4"
PlutoLinks = "~0.1.4"
PlutoUI = "~0.7.35"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "759a12cefe1cd1bb49e477bc3702287521797483"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.7"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "0a815f0060ab182f6c484b281107bfcd5bbb58dc"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.9.7"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "6b0440822974cab904c8b14d79743565140567f6"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.2.1"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "13468f237353112a01b2d6b32f3d0f80219944aa"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.2"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoHooks]]
deps = ["InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "bda062fe28bab89e96281ba3047afa515b06e8e1"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.4"

[[deps.PlutoLinks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "PlutoHooks", "Revise", "UUIDs"]
git-tree-sha1 = "147a4adcf0bf7715e2a120311b09ca35acdc4d64"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0420"
version = "0.1.4"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "85bf3e4bd279e405f91489ce518dedb1e32119cb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.35"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "606ddc4d3d098447a09c9337864c73d017476424"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.3.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─1ce83536-9746-11ec-342a-9db71b8eba53
# ╠═4f1b81f4-5262-4063-a8ce-1ce7da56a69b
# ╟─360c71c5-5f4e-4118-93af-e8fd3c3f163f
# ╠═994ab3fa-35e0-41a2-9f24-60b5024c8008
# ╠═f5b08153-58b2-45d6-8882-f593caba70f4
# ╟─e7a2f49e-3551-4185-bc16-21ec5b2ac252
# ╠═af6b1697-d127-4865-8347-e80c73e27194
# ╟─0e8057e7-3f69-4ed2-9fe1-d5af8e7e50ef
# ╟─285fbc70-aa43-44f7-af1d-689545945261
# ╠═562c408c-717e-4686-8c03-b315e72050fa
# ╟─eb1ce3ae-47d9-436c-be84-8baa2d03879b
# ╠═459f0704-5d59-4fde-bfa4-e3e2ff5d85c4
# ╟─d9619966-e61a-40fa-bd54-f9efb512f14d
# ╠═63a79fa0-366c-4edc-a79d-6a4a78d74ef6
# ╠═b893ec96-90d4-4131-8380-f4f15bc5cf40
# ╟─553818c9-2258-4993-a44d-17df498adcd1
# ╟─4bbd2db5-2a21-45e1-8f19-0a0f9a3b56fa
# ╠═47bd5a91-c963-4fed-bbfd-2c3c28eaf366
# ╠═b8b2c431-a18d-4b65-8207-8cc87e990ddf
# ╠═13946695-d509-4820-921d-9f1a23ea3ab0
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
Part of a measurement unit definition that corresponds to the SBML definition
of `Unit`. For example, the unit "per square megahour", Mh^(-2), is written as:

    SBML.UnitPart("second",  # base SI unit, this says we are measuring time
             -2,        # exponent, says "per square"
             6,         # log-10 scale of the unit, says "mega"
             1/3600)    # second-to-hour multiplier

Compound units (such as "volt-amperes" and "dozens of yards per ounce") are
built from multiple `UnitPart`s; see the definition of field `units` in
[`SBML.Model`](@ref).
"""
Base.@kwdef struct UnitPart
    kind::String
    exponent::Int
    scale::Int
    multiplier::Float64
end

"""
Abstract type for all kinds of gene product associations
"""
abstract type GeneProductAssociation end

"""
Gene product reference in the association expression
"""
struct GPARef <: GeneProductAssociation
    gene_product::String
end

"""
Boolean binary "and" in the association expression
"""
struct GPAAnd <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
Boolean binary "or" in the association expression
"""
struct GPAOr <: GeneProductAssociation
    terms::Vector{GeneProductAssociation}
end

"""
A simplified representation of MathML-specified math AST
"""
abstract type Math end

"""
A literal value (usually a numeric constant) in mathematical expression
"""
struct MathVal{T} <: Math where {T}
    val::T
end

"""
An identifier (usually a variable name) in mathematical expression
"""
struct MathIdent <: Math
    id::String
end

"""
A constant identified by name (usually something like `pi`, `e` or `true`) in
mathematical expression
"""
struct MathConst <: Math
    id::String
end

"""
A special value representing the current time of the simulation, with a special
name.
"""
struct MathTime <: Math
    id::String
end

"""
Function application ("call by name", no tricks allowed) in mathematical expression
"""
struct MathApply <: Math
    fn::String
    args::Vector{Math}
end

"""
Function definition (aka "lambda") in mathematical expression
"""
struct MathLambda <: Math
    args::Vector{String}
    body::Math
end

"""
$(TYPEDEF)

SBML Compartment with sizing information.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct Compartment
    name::Maybe{String} = nothing
    constant::Maybe{Bool} = nothing
    spatial_dimensions::Maybe{Int} = nothing
    size::Maybe{Float64} = nothing
    units::Maybe{String} = nothing
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end

"""
$(TYPEDEF)

Reaction with stoichiometry that assigns reactants and products their relative
consumption/production rates, lower/upper bounds (in tuples `lb` and `ub`, with
unit names), and objective coefficient (`oc`). Also may contains `notes` and
`annotation`.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct Reaction
    reactants::Dict{String,Float64} = Dict()
    products::Dict{String,Float64} = Dict()
    kinetic_parameters::Dict{String,Tuple{Float64,String}} = Dict()
    lower_bound::Maybe{String} = nothing
    upper_bound::Maybe{String} = nothing
    gene_product_association::Maybe{GeneProductAssociation} = nothing
    kinetic_math::Maybe{Math} = nothing
    reversible::Bool
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end

"""
$(TYPEDEF)

Abstract type representing SBML rules.
"""
abstract type Rule end

"""
$(TYPEDEF)

SBML algebraic rule.

# Fields
$(TYPEDFIELDS)
"""
struct AlgebraicRule <: Rule
    math::Math
end

"""
$(TYPEDEF)

SBML assignment rule.

# Fields
$(TYPEDFIELDS)
"""
struct AssignmentRule <: Rule
    id::String
    math::Math
end

"""
$(TYPEDEF)

SBML rate rule.

# Fields
$(TYPEDFIELDS)
"""
struct RateRule <: Rule
    id::String
    math::Math
end

"""
$(TYPEDEF)

SBML constraint.

# Fields
$(TYPEDFIELDS)
"""
struct Constraint
    math::Math
    message::String
end

"""
$(TYPEDEF)

Species metadata -- contains a human-readable `name`, a `compartment`
identifier, `formula`, `charge`, and additional `notes` and `annotation`.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct Species
    name::Maybe{String} = nothing
    compartment::String
    boundary_condition::Maybe{Bool} = nothing
    formula::Maybe{String} = nothing
    charge::Maybe{Int} = nothing
    initial_amount::Maybe{Tuple{Float64,Maybe{String}}} = nothing
    initial_concentration::Maybe{Tuple{Float64,Maybe{String}}} = nothing
    only_substance_units::Maybe{Bool} = nothing
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end

"""
$(TYPEDEF)

Gene product metadata.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct GeneProduct
    name::Maybe{String} = nothing
    label::Maybe{String} = nothing
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end

"""
Custom function definition.
"""
Base.@kwdef struct FunctionDefinition
    name::Maybe{String} = nothing
    body::Maybe{Math} = nothing
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end

"""
$(TYPEDEF)

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct EventAssignment
    variable::String
    math::Maybe{Math} = nothing
end

"""
$(TYPEDEF)

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct Event
    name::Maybe{String} = nothing
    trigger::Maybe{Math} = nothing
    event_assignments::Maybe{Vector{EventAssignment}} = nothing
end

"""
$(TYPEDEF)

Structure that collects the model-related data. Contains `parameters`, `units`,
`compartments`, `species` and `reactions` and `gene_products`, and additional
`notes` and `annotation` (also present internally in some of the data fields).
The contained dictionaries are indexed by identifiers of the corresponding
objects.

# Fields
$(TYPEDFIELDS)
"""
Base.@kwdef struct Model
    parameters::Dict{String,Tuple{Float64,String}} = Dict()
    units::Dict{String,Vector{UnitPart}} = Dict()
    compartments::Dict{String,Compartment} = Dict()
    species::Dict{String,Species} = Dict()
    initial_assignments::Dict{String,Math} = Dict()
    rules::Vector{Rule} = Rule[]
    constraints::Vector{Constraint} = Constraint[]
    reactions::Dict{String,Reaction} = Dict()
    objective::Dict{String,Float64} = Dict()
    gene_products::Dict{String,GeneProduct} = Dict()
    function_definitions::Dict{String,FunctionDefinition} = Dict()
    events::Dict{String,Event} = Dict()
    name::Maybe{String} = nothing
    id::Maybe{String} = nothing
    conversion_factor::Maybe{String} = nothing
    area_units::Maybe{String} = nothing
    extent_units::Maybe{String} = nothing
    length_units::Maybe{String} = nothing
    substance_units::Maybe{String} = nothing
    time_units::Maybe{String} = nothing
    volume_units::Maybe{String} = nothing
    notes::Maybe{String} = nothing
    annotation::Maybe{String} = nothing
end


text:
#=
Indexer struct for array type entities

An indexer has a symbol, range and current value
The range is an array of integers or any other Array subclass that can be used like "for i in I.range"
=#

"""
    Indexer

An entity representing an index to be applied to a variable or coefficient. It has a
symbol that can be used in expressions, a range of integer values, and a current
value that can be accessed with `I.value` for an index labeled I.

This should be built with the `index` function.
"""
mutable struct Indexer
    symbol::Symbol   # symbol used in expressions
    range::Array     # Range of integers.
    value::Int       # Current value
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
##### Moisture tendencies

#####
##### First order fluxes
#####

function flux(::Advect{TotalMoisture}, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_tot
end

function flux(::Advect{LiquidMoisture}, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_liq
end

function flux(::Advect{IceMoisture}, atmos, args)
    @unpack state = args
    u = state.ρu / state.ρ
    return u * state.moisture.ρq_ice
end

#####
##### Second order fluxes
#####

function flux(::MoistureDiffusion{TotalMoisture}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    return d_q_tot * state.ρ
end

function flux(::MoistureDiffusion{LiquidMoisture}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_liq = (-D_t) .* diffusive.moisture.∇q_liq
    return d_q_liq * state.ρ
end

function flux(::MoistureDiffusion{IceMoisture}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_ice = (-D_t) .* diffusive.moisture.∇q_ice
    return d_q_ice * state.ρ
end

#####
##### Sources
#####

function source(s::Subsidence{TotalMoisture}, m, args)
    @unpack state, aux, diffusive = args
    z = altitude(m, aux)
    w_sub = subsidence_velocity(s, z)
    k̂ = vertical_unit_vector(m, aux)
    return -state.ρ * w_sub * dot(k̂, diffusive.moisture.∇q_tot)
end

export CreateClouds
"""
    CreateClouds{PV <: Union{LiquidMoisture,IceMoisture}} <: TendencyDef{Source, PV}

A source/sink to `q_liq` and `q_ice` implemented as a relaxation towards
equilibrium in the Microphysics module.
The default relaxation timescales are defined in CLIMAParameters.jl.
"""
struct CreateClouds{PV <: Union{LiquidMoisture, IceMoisture}} <:
       TendencyDef{Source, PV} end

CreateClouds() = (CreateClouds{LiquidMoisture}(), CreateClouds{IceMoisture}())

function source(s::CreateClouds{LiquidMoisture}, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    # get current temperature and phase partition
    FT = eltype(state)
    q = PhasePartition(ts)
    T = air_temperature(ts)

    # phase partition corresponding to the current T and q.tot
    # (this is not the same as phase partition from saturation adjustment)
    ts_eq = PhaseEquil_ρTq(m.param_set, state.ρ, T, q.tot)
    q_eq = PhasePartition(ts_eq)

    # cloud condensate as relaxation source terms
    S_q_liq = conv_q_vap_to_q_liq_ice(m.param_set.microphys.liq, q_eq, q)

    return state.ρ * S_q_liq
end

function source(s::CreateClouds{IceMoisture}, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    # get current temperature and phase partition
    FT = eltype(state)
    q = PhasePartition(ts)
    T = air_temperature(ts)

    # phase partition corresponding to the current T and q.tot
    # (this is not the same as phase partition from saturation adjustment)
    ts_eq = PhaseEquil_ρTq(m.param_set, state.ρ, T, q.tot)
    q_eq = PhasePartition(ts_eq)

    # cloud condensate as relaxation source terms
    S_q_ice = conv_q_vap_to_q_liq_ice(m.param_set.microphys.ice, q_eq, q)

    return state.ρ * S_q_ice
end

function source(s::RemovePrecipitation{TotalMoisture}, m, args)
    @unpack state = args
    @unpack ts = args.precomputed
    if has_condensate(ts)
        nt = remove_precipitation_sources(s, m, args)
        return nt.S_ρ_qt
    else
        FT = eltype(state)
        return FT(0)
    end
end

function source(s::WarmRain_1M{TotalMoisture}, m, args)
    nt = warm_rain_sources(m, args)
    return nt.S_ρ_qt
end

function source(s::WarmRain_1M{LiquidMoisture}, m, args)
    nt = warm_rain_sources(m, args)
    return nt.S_ρ_ql
end

function source(s::RainSnow_1M{TotalMoisture}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_qt
end

function source(s::RainSnow_1M{LiquidMoisture}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_ql
end

function source(s::RainSnow_1M{IceMoisture}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_qi
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
module Basics

include("AbstractDataStructure.jl")
include("LinkedLists.jl")
include("SkipLists.jl")

end


text:
using Flux
using CUDA

Random.seed!(123)

@testset "train steps" begin
    crit = Chain(
        # Input: (64, 64, 1, :)
        Conv((4, 4), 1 => 64, selu, stride = 2, pad = SamePad()),   # -> (32, 32, 64, :)
        Conv((4, 4), 64 => 128, selu, stride = 2, pad = SamePad()), # -> (16, 16, 128, :)
        Conv((5, 5), 128 => 256, selu, stride = 2, pad = SamePad()),   # -> (8, 8, 256, :)
        Flux.flatten, # -> (8*8*256, :)
        Dropout(0.2),
        Dense(8 * 8 * 256, 1)
    ) |> gpu
    id = 64
    bs = 32
    opt_crit = ADAM(0.0002, (0.5, 0.9))
    x_t = cu(randn(Float32, id, id, 1, bs))
    x_f = cu(randn(Float32, id, id, 1, bs))
    @test step_critic!(opt_crit, crit, x_t, x_f, bs) isa Float32
end


text:
abstract type OrthantTrait end
abstract type Positive <: OrthantTrait end
abstract type Negative <: OrthantTrait end

struct Orthant{T <: OrthantTrait} <: AbstractProx
    function (::Type{Orthant})(::Type{T}) where T <: OrthantTrait
        return new{T}()
    end
end

function prox!(orthant::Orthant{Positive},step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    xcurr .= max.(xprev-step*gcurr, zeros(xprev))
end

function prox!(orthant::Orthant{Negative},step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    xcurr .= min.(xprev-step*gcurr, zeros(xprev))
end

PositiveOrthant() = Orthant(Positive)
NegativeOrthant() = Orthant(Negative)


text:
using GlobalSearchRegression, DataFrames, Distributions, CSV
# Insert here your working directory path
data = DataFrame(Array{Union{Missing,Float64}}(randn(1000,21)))
headers = [ :y ; [ Symbol("x$i") for i = 1:size(data,2) - 1 ] ]
names!(data, headers )
b=@elapsed gsreg("y x*", data; ttest=true, criteria=[:aic], parallel=16, csv="20x1000x16.csv")
t=DataFrame(A=Float64[])
push!(t,b)
CSV.write("juliatimmings.csv", t; append=true)
exit()


text:
using Documenter, Test, Weave, EconJobMarket
const DOCSPATH = realpath(joinpath(dirname(pathof(EconJobMarket)), "..", "docs"))
joinpath(DOCSPATH, "src") |>
    (dir -> isdir(dir) || mkdir(dir))
foreach(rm, joinpath(DOCSPATH, "src", file) for file ∈ readdir(joinpath(DOCSPATH, "src")))
foreach(filepath -> weave(joinpath(DOCSPATH, "jmd", filepath),
                          out_path = joinpath(DOCSPATH, "src"),
                          doctype = "github"),
        readdir(joinpath(DOCSPATH, "jmd")))
DocMeta.setdocmeta!(EconJobMarket, :DocTestSetup, :(using EconJobMarket), recursive = true)
makedocs(root = DOCSPATH,
         format = Documenter.HTML(),
         modules = [EconJobMarket],
         sitename = "EconJobMarket.jl",
         pages = ["Introduction" => "index.md",
                  "API" => "api.md"])
doctest(EconJobMarket)


text:
using PyPlot
using JLD

function main(code::AbstractString)
    log_path_record = "$(homedir())/open_loop/output-record-$(code).jld"
    d_rec = load(log_path_record)
    L_b = 0.125

    imu_meas    = d_rec["imu_meas"]
    gps_meas    = d_rec["gps_meas"]
    cmd_log     = d_rec["cmd_log"]
    cmd_pwm_log = d_rec["cmd_pwm_log"]
    vel_est     = d_rec["vel_est"]
    pos_info    = d_rec["pos_info"]

    t0 = max(cmd_pwm_log.t[1],vel_est.t[1],imu_meas.t[1])
    t_end = min(cmd_pwm_log.t[end],vel_est.t[end],imu_meas.t[end])

    t = t0+0.1:.02:t_end-0.1
    v = zeros(length(t))
    cmd = zeros(length(t))
    cmd_raw = zeros(length(t))

    for i=1:length(t)
        if cmd_log.z[t[i].>cmd_log.t,1][end] > 0 && vel_est.z[t[i].>vel_est.t,1][end] > 0.2
            v[i] = vel_est.z[t[i].>vel_est.t,1][end]
            cmd[i] = cmd_log.z[t[i].>cmd_log.t,1][end]
            cmd_raw[i] = cmd_pwm_log.z[t[i].>cmd_pwm_log.t,1][end]
        end
    end

    idx = (v.>1) & (cmd_raw .> 95)

    sz = size(v[idx],1)
    A = [cmd_raw[idx] ones(sz,1)]
    b = v[idx]
    coeff = A\b

    x = [90,110]
    y = [90 1;110 1]*coeff

    #plot(cmd,v)
    figure(1)
    ax1=subplot(211)
    plot(vel_est.t,vel_est.z)
    grid("on")
    legend(["v","v_fl","v_fr","v_bl","v_br"])
    subplot(212,sharex=ax1)
    plot(cmd_pwm_log.t,cmd_pwm_log.z[:,1])
    grid("on")

    figure(2)
    plot(cmd[cmd.>0],cmd[cmd.>0]./v[cmd.>0],"*")
    xlabel("u_a")
    ylabel("u_a/v")
    grid("on")

    figure(3)
    plot(cmd_raw[cmd_raw.>80],v[cmd_raw.>80],"*",x,y)
    grid("on")
    xlabel("PWM signal")
    ylabel("v [m/s]")

    println("c1 = $(1/coeff[1])")
    println("c2 = $(coeff[2]/coeff[1])")
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
function fillSonum8Tables()
    fillSonumTables(8)
end

function fillSonum16Tables()
    fillSonumTables(16)
end

function fillSonumTables(nbit::Int)

    # start with the lists - the fast bit
    fillListInv(nbit)
    fillListSqrt(nbit)
    fillTable(nbit,(+))
    fillTable(nbit,(-))
    fillTable(nbit,(*))
    fillTable(nbit,(/))
end

function fillTable(nbit::Int,operator)

    if nbit == 8
        sonum = sonum8
        Float2Sonum = Sonum8
    elseif nbit == 16
        sonum = sonum16
        Float2Sonum = Sonum16
    else
        throw(error("Only 8/16bit supported."))
    end

    # fill the underlying matrices not the symmetric views
    if operator == (+)
        if nbit == 8
            table = TableAdd8
        elseif nbit == 16
            table = TableAdd16
        end
    elseif operator == (-)
        if nbit == 8
            table = TableSub8
        elseif nbit == 16
            table = TableSub16
        end
    elseif operator == (*)
        if nbit == 8
            table = TableMul8
        elseif nbit == 16
            table = TableMul16
        end
    elseif operator == (/)
        if nbit == 8
            table = TableDiv8
        elseif nbit == 16
            table = TableDiv16
        end
    else
        throw(error("No ($operator) Table defined."))
    end

    n = length(sonum)
    if nbit == 16 print("Sonum16: Table lookup ($operator) ...") end
    t0 = time()
    if operator == (/) # not symmetric, nor anti-symmetric
        for i in 1:n
            for j in 1:n
                table[i,j] = Float2Sonum(operator(sonum[i],sonum[j]))
            end
        end

        # division by NaR should always yield NaR (and not 0 when interpreted as Inf)
        table[:,2^(nbit-1)+1] .= notareal(Float2Sonum)

    else
        for i in 1:n
            for j in 1:n
                if j >= i      # only upper triangle elements (symmetric or antisymmetric)
                    table[i,j] = Float2Sonum(operator(sonum[i],sonum[j]))
                end
            end
        end
    end
    t1 = time() - t0
    s = @sprintf "in %.1fs." t1
    if nbit == 16 print("\r\u1b[K"); println("Sonum16: Table lookup ($operator) $s") end
end

function fillListSqrt(nbit::Int)
    if nbit == 8
        sonum = sonum8
        ListSqrt = ListSqrt8
        Float2Sonum = Sonum8
    elseif nbit == 16
        sonum = sonum16
        ListSqrt = ListSqrt16
        Float2Sonum = Sonum16
    else
        throw(error("Only 8/16bit supported."))
    end

    for i in 1:length(sonum)
        ListSqrt[i] = Float2Sonum(sqrt(sonum[i]))
    end
end

function fillListInv(nbit::Int)
    if nbit == 8
        sonum = sonum8
        ListInv = ListInv8
        Float2Sonum = Sonum8
    elseif nbit == 16
        sonum = sonum16
        ListInv = ListInv16
        Float2Sonum = Sonum16
    else
        throw(error("Only 8/16bit supported."))
    end

    for i in 1:length(sonum)
        ListInv[i] = Float2Sonum(1/sonum[i])
    end
end


text:
using ODMXMLTools
using Test

path = dirname(@__FILE__)
cd(path)

@testset "ODMXMLTools.jl" begin

    io = IOBuffer();
    odm = ODMXMLTools.importxml(joinpath(path, "test.xml"))

    @test_nowarn ODMXMLTools.metadatalist(odm)
    st1 =  ODMXMLTools.findstudy(odm, "ST1")
    @test_nowarn ODMXMLTools.findelement(st1, :MetaDataVersion, "v2")
    @test_nowarn ODMXMLTools.findstudymetadata(odm, "ST1", "v2")

    @test_nowarn ODMXMLTools.buildmetadata(odm, "DEFS", "v1")
    mdb = ODMXMLTools.buildmetadata(odm, "ST1", "v2")

    @test_nowarn ODMXMLTools.studylist(odm)
    @test_nowarn ODMXMLTools.eventlist(mdb)
    @test_nowarn ODMXMLTools.formlist(mdb)
    @test_nowarn ODMXMLTools.itemgrouplist(mdb)
    @test_nowarn ODMXMLTools.itemlist(mdb)
    @test_nowarn ODMXMLTools.itemlist(mdb; optional = true)
    @test_nowarn ODMXMLTools.itemlist(mdb.el)
    @test_nowarn ODMXMLTools.itemgroupcontent(mdb, "IG_1")
    @test_nowarn ODMXMLTools.findelement(mdb, :ItemGroupDef, "IG_1")
    @test_nowarn ODMXMLTools.formcontent(mdb, "FORM_1")
    @test_nowarn ODMXMLTools.itemsformcontent(mdb, "FORM_1"; optional = true)

    @test_nowarn ODMXMLTools.validateodm(odm)

    cdat = ODMXMLTools.findclinicaldata(odm, "ST1", "v2")

    @test_nowarn ODMXMLTools.clinicaldatatable(cdat)

    @test_nowarn ODMXMLTools.subjectdatatable(odm; attrs = [:SubjectKey, :StudySubjectID])

    @test_nowarn ODMXMLTools.studyinfo(odm; io = io)
end