text:
import DataStructures:OrderedDict

"""
Variant type is typically used to represent a record of VCF
See VCF spec:
VCF 4.1 http://samtools.github.io/hts-specs/VCFv4.1.pdf
VCF 4.2 http://samtools.github.io/hts-specs/VCFv4.2.pdf
"""
type Variant
    chrom::ASCIIString
    pos::Int64
    id::ASCIIString
    ref::ASCIIString
    alt::ASCIIString
    qual::ASCIIString
    filter::ASCIIString
    info::ASCIIString
    format::ASCIIString
    samples::Array{ASCIIString, 1}
end

"""
Initialize a variant
requires: #CHROM, POS, ID, REF, ALT
optionals: QUAL, FILTER, INFO, FORMAT
"""
function Variant(chrom::ASCIIString,pos::Int64,id::ASCIIString= ".",ref::ASCIIString=".",alt::ASCIIString=".";
        qual::ASCIIString="",filter::ASCIIString="",info::ASCIIString="",format::ASCIIString="", samples::Array{ASCIIString, 1}=Array{ASCIIString, 1}())
    return Variant(chrom, pos, id, ref, alt, qual, filter, info, format, samples)
end

type VcfHeader
    metas::OrderedDict{ASCIIString, Array{Any, 1}}
    columns::Array{ASCIIString, 1}
    VcfHeader() = new(OrderedDict{ASCIIString, Array{Any, 1}}(), Array{ASCIIString, 1}())
end

type Vcf
    header::VcfHeader
    data::Array{Variant, 1}
end

function vcf_has_format_column(hd::VcfHeader)
    return haskey(hd.metas, "FORMAT")
end

function vcf_has_format_column(vcf::Vcf)
    return vcf_has_format_column(vcf.header)
end

function vcf_samples(vcf::Vcf)
    samples = Array{ASCIIString, 1}()
    columns = vcf.header.columns
    const tags = [ "CHROM", "CHROM", "POS", "ID","REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT"]
    for c in columns
        if !(c in tags)
            push!(samples, c)
        end
    end
    return samples
end

function version(vcf::Vcf)
    if "fileformat" in keys(vcf.header.metas)
        vstr = replace(vcf.header.metas["fileformat"][1], "VCFv", "")
        return VersionNumber(vstr)
    else
        return VersionNumber(0)
    end
end

function var_lt(v1::Variant, v2::Variant)
    if v1.chrom < v2.chrom || (v1.chrom == v2.chrom && v1.pos < v2.pos)# || (v1.chrom == v2.chrom && v1.pos == v2.pos && v1.id < v2.id)
        return true
    end
    return false
end

function var_gt(v1::Variant, v2::Variant)
    if v1.chrom > v2.chrom || (v1.chrom == v2.chrom && v1.pos > v2.pos)# || (v1.chrom == v2.chrom && v1.pos == v2.pos && v1.id > v2.id)
        return true
    end
    return false
end

<(v1::Variant, v2::Variant) = var_lt(v1, v2)
>(v1::Variant, v2::Variant) = var_gt(v1, v2)

function sort!(obj::Vcf)
    data = sort(obj.data, lt = <)
    obj.data = data
    return obj
end

function sort(obj::Vcf)
    data = sort(obj.data, lt=var_comp)
    header = deepcopy(obj.header)
    return Vcf(header, data)
end

function vcf_issorted(obj::Vcf)
    data = obj.data
    for i in 1:length(data)-1
        v1 = data[i]
        v2 = data[i+1]
        if v1 > v2
            return false
        end
    end
    return true
end

function vcf_merge(v1::Vcf, v2::Vcf)
    if !issorted(v1)
        #info("The first vcf is not sorted, sort it now")
        sort!(v1)
    end
    if !issorted(v2)
        #info("The second vcf is not sorted, sort it now")
        sort!(v2)
    end
    result = Array{Variant, 1}()
    i = 1
    j = 1
    last_push = nothing
    while i<=length(v1.data) || j<=length(v2.data)
        take_v1 = true
        if i>length(v1.data)
            take_v1 = false
        elseif j>length(v2.data)
            take_v1 = true
        else
            if v1.data[i] < v2.data[j]
                take_v1 = true
            elseif v1.data[i] > v2.data[j]
                take_v1 = false
            else
                # same position, check if they are identical
                if v1.data[i].ref != v2.data[j].ref || v1.data[i].alt != v2.data[j].alt
                    warn("mismatch discarded at ", v1.data[i].chrom, ":", v1.data[i].pos, " ", v1.data[i].ref, "-->", v1.data[i].alt, " vs ", v2.data[j].ref, "-->", v2.data[j].alt)
                    i += 1
                    j += 1
                    continue
                else
                    # if identical, only take the one from v1, and skip this one from v2
                    take_v1 = true
                    j += 1
                end
            end
        end
        if take_v1
            if last_push == nothing || last_push < v1.data[i]
                push!(result, v1.data[i])
                last_push = v1.data[i]
            end
            i += 1
        else
            if last_push == nothing || last_push < v2.data[j]
                push!(result, v2.data[j])
                last_push = v2.data[j]
            end
            j += 1
        end
    end
    # use header only from v1
    header = deepcopy(v1.header)
    return Vcf(header, result)
end

function vcf_intersect(v1::Vcf, v2::Vcf)
    if !issorted(v1)
        #info("The first vcf is not sorted, sort it now")
        sort!(v1)
    end
    if !issorted(v2)
        #info("The second vcf is not sorted, sort it now")
        sort!(v2)
    end
    result = Array{Variant, 1}()
    i = 1
    j = 1
    while i<=length(v1.data) && j<=length(v2.data)
        if v1.data[i] < v2.data[j]
            i += 1
        elseif v1.data[i] > v2.data[j]
            j += 1
        else
            if v1.data[i].ref == v2.data[j].ref && v1.data[i].alt == v2.data[j].alt
                push!(result, v1.data[i])
            end
            i += 1
            j += 1
        end
    end
    # use header only from v1
    header = deepcopy(v1.header)
    return Vcf(header, result)
end

function vcf_minus(v1::Vcf, v2::Vcf)
    common = v1 * v2
    result = Array{Variant, 1}()
    i = 1
    j = 1
    while i<=length(v1.data)
        if j > length(common.data)
            push!(result, v1.data[i])
            i += 1
            continue
        end
        if v1.data[i] < common.data[j]
            push!(result, v1.data[i])
            i += 1
        elseif v1.data[i] > common.data[j]
            j += 1
        else
            if v1.data[i].ref != common.data[j].ref || v1.data[i].alt != common.data[j].alt
                push!(result, v1.data[i])
            end
            i += 1
        end
    end
    # use header only from v1
    header = deepcopy(v1.header)
    return Vcf(header, result)
end

function parse_genotype(sample_str::ASCIIString, format::ASCIIString)
    metas = split(format, ":")
    gt_pos = -1
    ad_pos = -1
    rd_pos = -1
    for i in 1:length(metas)
        if metas[i] == "GT"
            gt_pos = i
        elseif metas[i] == "AD"
            ad_pos = i
        elseif metas[i] == "RD"
            rd_pos = i
        end
    end
    if gt_pos == -1 || ad_pos == -1
        error("GT or AD is not present in FORMAT")
    end
    values = split(sample_str, ":")
    gt = values[gt_pos]
    if gt == "./."
        return gt, -1, -1
    end
    ad = split(values[ad_pos], ",")
    ref_num = 0
    alt_num = 0
    if length(ad)<2
        # GT:GQ:DP:RD:AD:FREQ:DP4
        if rd_pos == -1
            error("ref depth is not given neithger by AD nor RD")
        end
        ref_num = parse(Int64, strip(values[rd_pos]))
        alt_num = parse(Int64, strip(ad[1]))
    else
        # "GT:AD:DP:GQ:PL"
        ref_num = parse(Int64, strip(ad[1]))
        alt_num = parse(Int64, strip(ad[2]))
    end
    return gt, ref_num, alt_num
end

function vcf_filter(vcf::Vcf; total_depth=0, allele_depth=0, allele_freq=0.0, sample="")
    if !vcf_has_format_column(vcf)
        error("This VCF doesn't have a FORMAT column")
    end
    all_samples = vcf_samples(vcf)
    if length(all_samples)==0
        error("This VCF doesn't have any sample column")
    end
    if sample==""
        sample=all_samples[1]
    elseif !(sample in all_samples)
        error("This VCF doesn't have sample $sample")
    end
    # sample index
    s = 1
    for i=1:length(all_samples)
        if all_samples[i] == sample
            s = i
            break
        end
    end
    good = Variant[]
    bad = Variant[]
    # GT:AD:DP:GQ:PL    0/1:37,14:51:99:111,0,614
    ad_pos = 2
    for v in vcf.data
        gt, ref_num, alt_num = parse_genotype(v.samples[s], v.format)
        total_num = ref_num + alt_num
        alt_freq = float(alt_num)/float(total_num)
        if total_num >= total_depth && alt_num >= allele_depth && alt_freq >= allele_freq
            push!(good, deepcopy(v))
        else
            push!(bad, deepcopy(v))
        end
    end
    good_vcf = Vcf(deepcopy(vcf.header), good)
    bad_vcf = Vcf(deepcopy(vcf.header), bad)
    return good_vcf, bad_vcf
end

# a region is a string like chr1:13579-2312332
function vcf_filter_region(vcf::Vcf, regions::Array{ASCIIString, 1})
    parsed_regions = []
    for r in regions
        items = split(r, ":")
        if length(items) != 2
            warn("error format of region, should be like chr1:13579-2312332")
            continue
        end
        chr = items[1]
        pos = split(items[2], "-")
        start_pos = parse(Int, pos[1])
        end_pos = start_pos
        if length(pos) > 1
            end_pos = parse(Int, pos[2])
        end
        push!(parsed_regions, (chr, start_pos, end_pos))
    end
    good = Variant[]
    bad = Variant[]
    for v in vcf.data
        found = false
        for (chr, start_pos, end_pos) in parsed_regions
            if v.chrom == chr && v.pos >= start_pos && v.pos <= end_pos
                found = true
                break
            end
        end
        new_var = deepcopy(v)
        if found
            push!(good, new_var)
        else
            push!(bad, new_var)
        end
    end
    good_vcf = Vcf(deepcopy(vcf.header), good)
    bad_vcf = Vcf(deepcopy(vcf.header), bad)
    return good_vcf, bad_vcf
end

# a region is like chr1:13579-2312332;chr3:132229-4533233
function vcf_filter_region(vcf::Vcf, regions::ASCIIString)
    splitted_regions = ASCIIString[]
    for r in split(regions, ";")
        push!(splitted_regions, ASCIIString(r))
    end
    return vcf_filter_region(vcf, splitted_regions)
end

function vcf_diff_genotype(v1::Vcf, v2::Vcf, v1_sample_id = 1, v2_sample_id = 1)
    if !issorted(v1)
        #info("The first vcf is not sorted, sort it now")
        sort!(v1)
    end
    if !issorted(v2)
        #info("The second vcf is not sorted, sort it now")
        sort!(v2)
    end
    if v1_sample_id > length(vcf_samples(v1))
        error("Wrong v1_sample_id")
    elseif v2_sample_id > length(vcf_samples(v2))
        error("Wrong v2_sample_id")
    end
    v1_diff = Variant[]
    v2_diff = Variant[]
    common = v1 * v2
    v1_common = v1 * common
    v2_common = v2 * common
    if length(v1_common) != length(v2_common)
        error("length(v1_common) != length(v2_common)")
    end
    for i in 1:length(v1_common)
        gt1, ref_num1, alt_num1 = parse_genotype(v1_common.data[i].samples[v1_sample_id], v1_common.data[i].format)
        gt2, ref_num2, alt_num2 = parse_genotype(v2_common.data[i].samples[v2_sample_id], v2_common.data[i].format)
        if gt1 != gt2
            push!(v1_diff, deepcopy(v1_common.data[i]))
            push!(v2_diff, deepcopy(v2_common.data[i]))
        end
    end
    return Vcf(deepcopy(v1.header), v1_diff), Vcf(deepcopy(v2.header), v2_diff)
end

function vcf_split_sample(vcf::Vcf)
    samples = vcf_samples(vcf)
    const tags = [ "CHROM", "CHROM", "POS", "ID","REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT"]
    columns_without_sample = ASCIIString[]
    for c in vcf.header.columns
        if (c in tags)
            push!(columns_without_sample, c)
        end
    end
    splitted = Vcf[]
    for s in samples
        header = deepcopy(vcf.header)
        columns = deepcopy(columns_without_sample)
        push!(columns, s)
        header.columns = columns
        data = Variant[]
        push!(splitted, Vcf(header, data))
    end
    for v in vcf.data
        for s in 1:length(samples)
            gt, ref_num, alt_num = parse_genotype(v.samples[s], v.format)
            if gt != "./." && gt != "0/0"
                var = deepcopy(v)
                var.samples = ASCIIString[]
                push!(var.samples, v.samples[s])
                push!(splitted[s].data, var)
            end
        end
    end
    return splitted
end

+(v1::Vcf, v2::Vcf) = vcf_merge(v1, v2)
-(v1::Vcf, v2::Vcf) = vcf_minus(v1, v2)
*(v1::Vcf, v2::Vcf) = vcf_intersect(v1, v2)

Base.issorted(obj::Vcf) = vcf_issorted(obj)
Base.merge(v1::Vcf, v2::Vcf) = vcf_merge(v1, v2)
Base.length(v::Vcf) = length(v.data)
Base.getindex(obj::Vcf, i::Int64) = getindex(obj.data, i)


text:
export MvNormalWeightedMeanPrecision

import Distributions: logdetcov, distrname, sqmahal, sqmahal!, AbstractMvNormal
import LinearAlgebra: diag, Diagonal, dot
import Base: ndims, precision, length, size, prod

struct MvNormalWeightedMeanPrecision{ T <: Real, M <: AbstractVector{T}, P <: AbstractMatrix{T} } <: AbstractMvNormal
    xi :: M
    Λ  :: P
end

function MvNormalWeightedMeanPrecision(xi::AbstractVector{ <: Real }, Λ::AbstractMatrix{ <: Real })
    T = promote_type(eltype(xi), eltype(Λ))
    return MvNormalWeightedMeanPrecision(convert(AbstractArray{T}, xi), convert(AbstractArray{T}, Λ))
end

function MvNormalWeightedMeanPrecision(xi::AbstractVector{ <: Integer}, Λ::AbstractMatrix{ <: Integer })
    return MvNormalWeightedMeanPrecision(float.(xi), float.(Λ))
end

function MvNormalWeightedMeanPrecision(xi::AbstractVector, λ::AbstractVector)
    return MvNormalWeightedMeanPrecision(xi, matrix_from_diagonal(promote_type(eltype(xi), eltype(λ)), λ))
end

function MvNormalWeightedMeanPrecision(xi::AbstractVector{T}) where T
    return MvNormalWeightedMeanPrecision(xi, convert(AbstractArray{T}, ones(length(xi))))
end

Distributions.distrname(::MvNormalWeightedMeanPrecision) = "MvNormalWeightedMeanPrecision"

weightedmean(dist::MvNormalWeightedMeanPrecision) = dist.xi

function mean_cov(dist::MvNormalWeightedMeanPrecision)
    Σ = cov(dist)
    μ = Σ * weightedmean(dist)
    return (μ, Σ)
end

function Distributions.mean(dist::MvNormalWeightedMeanPrecision)
    z = fastcholesky(precision(dist))
    return z \ weightedmean(dist)
end
Distributions.mode(dist::MvNormalWeightedMeanPrecision)      = mean(dist)
Distributions.var(dist::MvNormalWeightedMeanPrecision)       = diag(cov(dist))
Distributions.cov(dist::MvNormalWeightedMeanPrecision)       = cholinv(dist.Λ)
Distributions.invcov(dist::MvNormalWeightedMeanPrecision)    = dist.Λ
Distributions.std(dist::MvNormalWeightedMeanPrecision)       = cholsqrt(cov(dist))
Distributions.logdetcov(dist::MvNormalWeightedMeanPrecision) = -chollogdet(invcov(dist))

Distributions.sqmahal(dist::MvNormalWeightedMeanPrecision, x::AbstractVector) = sqmahal!(similar(x), dist, x)

function Distributions.sqmahal!(r, dist::MvNormalWeightedMeanPrecision, x::AbstractVector)
    μ = mean(dist)
    for i in 1:length(r)
        @inbounds r[i] = μ[i] - x[i]
    end
    return dot(r, invcov(dist), r)
end

Base.eltype(::MvNormalWeightedMeanPrecision{T})     where T = T
Base.precision(dist::MvNormalWeightedMeanPrecision)         = invcov(dist)
Base.length(dist::MvNormalWeightedMeanPrecision)            = length(weightedmean(dist))
Base.ndims(dist::MvNormalWeightedMeanPrecision)             = length(dist)
Base.size(dist::MvNormalWeightedMeanPrecision)              = (length(dist), )

Base.convert(::Type{ <: MvNormalWeightedMeanPrecision }, xi::AbstractVector, Λ::AbstractMatrix) = MvNormalWeightedMeanPrecision(xi, Λ)

function Base.convert(::Type{ <: MvNormalWeightedMeanPrecision{T} }, xi::AbstractVector, Λ::AbstractMatrix) where { T <: Real }
    MvNormalWeightedMeanPrecision(convert(AbstractArray{T}, xi), convert(AbstractArray{T}, Λ))
end

vague(::Type{ <: MvNormalWeightedMeanPrecision }, dims::Int) = MvNormalWeightedMeanPrecision(zeros(dims), fill(tiny, dims))

prod_analytical_rule(::Type{ <: MvNormalWeightedMeanPrecision }, ::Type{ <: MvNormalWeightedMeanPrecision }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdPreserveType, left::MvNormalWeightedMeanPrecision, right::MvNormalWeightedMeanPrecision)
    xi = weightedmean(left) + weightedmean(right)
    Λ  = invcov(left) + invcov(right)
    return MvNormalWeightedMeanPrecision(xi, Λ)
end

function Base.prod(::ProdAnalytical, left::MvNormalWeightedMeanPrecision, right::MvNormalWeightedMeanPrecision)
    xi = weightedmean(left) + weightedmean(right)
    Λ  = invcov(left) + invcov(right)
    return MvNormalWeightedMeanPrecision(xi, Λ)
end

ext: jl


text:
# In this example, we will demonstrate using Rimu with MPI.  A runnable script is located at
# [here](../../scripts/BHM-example-mpi.jl). Run it as `mpirun julia BHM-example.jl`.

# We start by importing `Rimu` and `Rimu.RMPI`, which contains MPI-replated
# functionality.
using Rimu
using Rimu.RMPI

# We will compute the ground-state of a Bose-Hubbard model in momentum space with 10 particles
# in 10 sites.

# First, we define the Hamiltonian. We want to start from an address with zero momentum.

address = BoseFS((0, 0, 0, 0, 10, 0, 0, 0, 0, 0))

# We will set the interaction strength `u` to `6`. The hopping strength `t` defaults to `1.0`.

hamiltonian = HubbardMom1D(address; u=6.0)

# Next, we construct the starting vector. Wrap a vector in `MPIData` to make it MPI
# distributed. We set the vector's style to [`IsDynamicSemistochastic`](@ref), which
# improves statistics and reduces the sign problem.

dvec = MPIData(DVec(address => 1.0; style=IsDynamicSemistochastic()))

# We set a reporting strategy. We will use [`ReportToFile`](@ref), which writes the reports
# directly to a file. This is useful for reducing memory use in long-running jobs, as we
# don't need to keep the results in memory. Setting `save_if=is_mpi_root()` will ensure only
# the root MPI rank will write to the file. The `chunk_size` parameter determines how often
# the data is saved to the file.

r_strat = ReportToFile(filename="result.arrow", save_if=is_mpi_root(), chunk_size=1000)

# Now, we can set other parameters as usual. We will perform the computation with 10_000
# walkers. We will also compute the projected energy.

s_strat = DoubleLogUpdate(targetwalkers=10_000)
post_step = ProjectedEnergy(hamiltonian, dvec)

# Finally, we can run the computation. The `@mpi_root` macro performs an action on the root
# rank only, which is useful for printing.

lomc!(hamiltonian, dvec; r_strat, s_strat, post_step, dτ=1e-4, laststep=10_000)
@mpi_root println("Finished!")


text:
# BAM File Format
# ===============

module BAM

import BGZFStreams
import BioAlignments
import XAM: XAM, SAM
import GenomicFeatures: GenomicFeatures, Interval
import BioSequences
import BioCore: BioCore, isfilled, header

include("bai.jl")
include("auxdata.jl")
include("reader.jl")
include("record.jl")
include("writer.jl")
include("overlap.jl")

end


text:
function ground(args...)
    background("white")
    sethue("black")
end

@testset "Javis Viewer" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    render(vid; pathname = "")

    action_list = [back, star_obj]

    viewer_win, frame_dims, r_slide, tbox, canvas, actions, total_frames, video =
        Javis._javis_viewer(vid, 100, action_list, false)
    visible(viewer_win, false)

    @test get_gtk_property(viewer_win, :title, String) == "Javis Viewer"

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    second_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 2

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test Reactive.value(r_slide) == 1

    @test first_frame != second_frame

    Javis._decrement(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    last_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == total_frames

    Javis._increment(video, [r_slide, tbox], actions, frame_dims, canvas, total_frames)
    sleep(0.1)
    curr_frame = Reactive.value(r_slide)
    first_frame = Javis.get_javis_frame(video, actions, curr_frame)
    @test curr_frame == 1

    @test last_frame != first_frame
end

@testset "Jupyter Viewer" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    objects = vid.objects
    frames = Javis.preprocess_frames!(objects)

    img = Javis._jupyter_viewer(vid, length(frames), objects, 30)
    @test img.output.val == Javis.get_javis_frame(vid, objects, 1)

    txt = Interact.textbox(1:length(frames), typ = "Frame", value = 2)
    frm = Interact.slider(1:length(frames), label = "Frame", value = txt[] + 1)
    @test Javis.get_javis_frame(vid, objects, 2) ==
          Javis.get_javis_frame(vid, objects, txt[])
    @test Javis.get_javis_frame(vid, objects, 3) ==
          Javis.get_javis_frame(vid, objects, frm[])

    for i in 4:length(frames)
        output = Javis.get_javis_frame(vid, objects, i)
        wdg = Widget(["frm" => frm, "txt" => txt], output = output)
        img = @layout! wdg vbox(hbox(:frm, :txt), output)
        @test img.output.val == output
    end
end

@testset "Pluto Viewer" begin
    v = Javis.PlutoViewer("foo.png")
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    objects = vid.objects
    frames = Javis.preprocess_frames!(objects)

    @test v.filename === "foo.png"
    img = Javis._pluto_viewer(vid, length(frames), objects)
    for i in 1:length(img)
        @test img[i] == Javis.get_javis_frame(vid, objects, i)
    end
end

@testset "Livestreaming" begin
    astar(args...; do_action = :stroke) = star(O, 50, 5, 0.5, 0, do_action)
    acirc(args...; do_action = :stroke) = circle(Point(100, 100), 50, do_action)

    vid = Video(500, 500)
    back = Background(1:100, ground)
    star_obj = Object(1:100, astar)
    act!(star_obj, Action(morph_to(acirc; do_action = :fill)))

    conf_local = setup_stream(:local, address = "0.0.0.0", port = 8081)
    @test conf_local isa Javis.StreamConfig
    @test conf_local.livestreamto == :local
    @test conf_local.protocol == "udp"
    @test conf_local.address == "0.0.0.0"
    @test conf_local.port == 8081

    conf_twitch_err = setup_stream(:twitch)
    conf_twitch = setup_stream(:twitch, twitch_key = "foo")
    @test conf_twitch_err isa Javis.StreamConfig
    @test conf_twitch_err.livestreamto == :twitch
    @test isempty(conf_twitch_err.twitch_key)
    @test conf_twitch.twitch_key == "foo"

    render(vid, streamconfig = conf_local)

    # errors with macos; a good test to have
    # test_local = run(pipeline(`lsof -i -P -n`, `grep ffmpeg`))
    # @test test_local isa Base.ProcessChain
    # @test test_local.processes isa Vector{Base.Process}

    cancel_stream()
    @test_throws ProcessFailedException run(
        pipeline(
            `ps aux`,
            pipeline(`grep ffmpeg`, pipeline(`grep stream_loop`, `awk '{print $2}'`)),
        ),
    )

    @test_throws ErrorException render(vid, streamconfig = conf_twitch_err)
end


text:
# __precompile__(true)

module Expectations

# Load external dependencies.
using Compat, FastGaussQuadrature
using Compat.LinearAlgebra
using SpecialFunctions
using Distributions

# Load internal files.
include("types.jl"),
include("iterable.jl")

# Export
export expectation, Expectation, IterableExpectation, Gaussian, Trapezoidal, FiniteDiscrete, QuadratureAlgorithm, ExplicitQuadratureAlgorithm, nodes, weights, QuantileRange

end # module


text:
struct PointValue{P} <: Dof
  point::P
end

"""
    struct LagrangianDofBasis{P,V} <: AbstractArray{<:Dof}
      nodes::Vector{P}
      dof_to_node::Vector{Int}
      dof_to_comp::Vector{Int}
      node_and_comp_to_dof::Vector{V}
    end

Type that implements a Lagrangian dof basis.

Fields:

- `nodes::Vector{P}` vector of points (`P<:Point`) storing the nodal coordinates
- `node_and_comp_to_dof::Vector{V}` vector such that `node_and_comp_to_dof[node][comp]` returns the dof associated with node `node` and the component `comp` in the type `V`.
- `dof_to_node::Vector{Int}` vector of integers such that `dof_to_node[dof]` returns the node id associated with dof id `dof`.
- `dof_to_comp::Vector{Int}` vector of integers such that `dof_to_comp[dof]` returns the component id associated with dof id `dof`.

"""
struct LagrangianDofBasis{P,V} <: AbstractVector{PointValue{P}}
  nodes::Vector{P}
  dof_to_node::Vector{Int}
  dof_to_comp::Vector{Int}
  node_and_comp_to_dof::Vector{V}
end

@inline Base.size(a::LagrangianDofBasis) = (length(a.nodes),)
@inline Base.axes(a::LagrangianDofBasis) = (axes(a.nodes,1),)
# @santiagobadia : Not sure we want to create the monomial machinery
@inline Base.getindex(a::LagrangianDofBasis,i::Integer) = PointValue(a.nodes[i])
@inline Base.IndexStyle(::LagrangianDofBasis) = IndexLinear()

# This one takes a basis and replaces the nodes
function LagrangianDofBasis(dofs::LagrangianDofBasis{P},nodes::Vector{P}) where P
  @check length(nodes) == length(dofs.nodes)
  LagrangianDofBasis(
    nodes,
    dofs.dof_to_node,
    dofs.dof_to_comp,
    dofs.node_and_comp_to_dof)
end

"""
    LagrangianDofBasis(::Type{T},nodes::Vector{<:Point}) where T

Creates a `LagrangianDofBasis` for fields of value type `T` associated
with the vector of nodal coordinates `nodes`.
"""
function LagrangianDofBasis(::Type{T},nodes::Vector{<:Point}) where T
  r = _generate_dof_layout_node_major(T,length(nodes))
  LagrangianDofBasis(nodes,r...)
end

get_nodes(b::LagrangianDofBasis) = b.nodes

function _generate_dof_layout_node_major(::Type{<:Real},nnodes::Integer)
  ndofs = nnodes
  dof_to_comp = ones(Int,ndofs)
  dof_to_node = collect(1:nnodes)
  node_and_comp_to_dof = collect(1:ndofs)
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

# Node major implementation
function _generate_dof_layout_node_major(::Type{T},nnodes::Integer) where T<:MultiValue
  ncomps = num_components(T)
  V = change_eltype(T,Int)
  ndofs = ncomps*nnodes
  dof_to_comp = zeros(Int,ndofs)
  dof_to_node = zeros(Int,ndofs)
  node_and_comp_to_dof = zeros(V,nnodes)
  m = zero(Mutable(V))
  for node in 1:nnodes
    for comp in 1:ncomps
      o = nnodes*(comp-1)
      dof = node+o
      dof_to_comp[dof] = comp
      dof_to_node[dof] = node
      m[comp] = dof
    end
    node_and_comp_to_dof[node] = m
  end
  (dof_to_node, dof_to_comp, node_and_comp_to_dof)
end

function return_cache(b::LagrangianDofBasis,field)
  cf = return_cache(field,b.nodes)
  vals = evaluate!(cf,field,b.nodes)
  ndofs = length(b.dof_to_node)
  r = _lagr_dof_cache(vals,ndofs)
  c = CachedArray(r)
  (c, cf)
end

function _lagr_dof_cache(node_comp_to_val::AbstractVector,ndofs)
  T = eltype(node_comp_to_val)
  r = zeros(eltype(T),ndofs)
end

function _lagr_dof_cache(node_pdof_comp_to_val::AbstractMatrix,ndofs)
  _, npdofs = size(node_pdof_comp_to_val)
  T = eltype(node_pdof_comp_to_val)
  r = zeros(eltype(T),ndofs,npdofs)
end

@inline function evaluate!(cache,b::LagrangianDofBasis,field)
  c, cf = cache
  vals = evaluate!(cf,field,b.nodes)
  ndofs = length(b.dof_to_node)
  T = eltype(vals)
  ncomps = num_components(T)
  @check ncomps == num_components(eltype(b.node_and_comp_to_dof)) """\n
  Unable to evaluate LagrangianDofBasis. The number of components of the
  given Field does not match with the LagrangianDofBasis.

  If you are trying to interpolate a function on a FESpace make sure that
  both objects have the same value type.

  For instance, trying to interpolate a vector-valued funciton on a scalar-valued FE space
  would raise this error.
  """
  _evaluate_lagr_dof!(c,vals,b.node_and_comp_to_dof,ndofs,ncomps)
end

function _evaluate_lagr_dof!(c::AbstractVector,node_comp_to_val,node_and_comp_to_dof,ndofs,ncomps)
  setsize!(c,(ndofs,))
  r = c.array
  for node in LinearIndices(node_and_comp_to_dof)
    comp_to_dof = node_and_comp_to_dof[node]
    comp_to_val = node_comp_to_val[node]
    for comp in 1:ncomps
      dof = comp_to_dof[comp]
      val = comp_to_val[comp]
      r[dof] = val
    end
  end
  r
end

function _evaluate_lagr_dof!(c::AbstractMatrix,node_pdof_comp_to_val,node_and_comp_to_dof,ndofs,ncomps)
  _, npdofs = size(node_pdof_comp_to_val)
  setsize!(c,(ndofs,npdofs))
  r = c.array
  for node in LinearIndices(node_and_comp_to_dof)
    comp_to_dof = node_and_comp_to_dof[node]
    for pdof in 1:npdofs
      comp_to_val = node_pdof_comp_to_val[node,pdof]
      for comp in 1:ncomps
        dof = comp_to_dof[comp]
        val = comp_to_val[comp]
        r[dof,pdof] = val
      end
    end
  end
  r
end


text:
using Documenter

deploydocs(repo = "github.com/uva-bi-sdad/COVID_Tracking_Clean.git",
           push_preview = true)


text:
using Random

const rseed = Ref(Random.GLOBAL_RNG)  # to get new random directories (see julia #24445)
function randtmp()
    Random.seed!(rseed[])
    dirname = joinpath(tempdir(), randstring(10))
    rseed[] = Random.GLOBAL_RNG
    dirname
end

@static if Sys.isapple()
    yry() = (sleep(1.1); revise(); sleep(1.1))
else
    yry() = (sleep(0.1); revise(); sleep(0.1))
end


text:
struct Contiguous{N} end
Base.@pure Contiguous(N::Int) = Contiguous{N}()
_get(::Contiguous{N}) where {N} = N
"""
contiguous_axis(::Type{T}) -> Contiguous{N}

Returns the axis of an array of type `T` containing contiguous data.
If no axis is contiguous, it returns `Contiguous{-1}`.
If unknown, it returns `nothing`.
"""
contiguous_axis(x) = contiguous_axis(typeof(x))
contiguous_axis(::Type) = nothing
contiguous_axis(::Type{<:Array}) = Contiguous{1}()
contiguous_axis(::Type{<:Tuple}) = Contiguous{1}()
function contiguous_axis(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}) where {T,A<:AbstractVector{T}}
    c = contiguous_axis(A)
    isnothing(c) && return nothing
    c === Contiguous{1}() ? Contiguous{2}() : Contiguous{-1}()
end
function contiguous_axis(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}) where {T,A<:AbstractMatrix{T}}
    c = contiguous_axis(A)
    isnothing(c) && return nothing
    contig = _get(c)
    new_contig = contig == -1 ? -1 : 3 - contig
    Contiguous{new_contig}()
end
function contiguous_axis(::Type{<:PermutedDimsArray{T,N,I1,I2,A}}) where {T,N,I1,I2,A<:AbstractArray{T,N}}
    c = contiguous_axis(A)
    isnothing(c) && return nothing
    new_contig = I2[_get(c)]
    Contiguous{new_contig}()
end
function contiguous_axis(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    _contiguous_axis(S, contiguous_axis(A))
end
_contiguous_axis(::Any, ::Nothing) = nothing
@generated function _contiguous_axis(::Type{S}, ::Contiguous{C}) where {C,N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    n = 0
    new_contig = contig = C
    for np in 1:NP
        if I.parameters[np] <: AbstractUnitRange
            n += 1
            if np == contig
                new_contig = n
            end
        else
            if np == contig
                new_contig = -1
            end
        end
    end
    # If n != N, then an axis was indexed by something other than an integer or `AbstractUnitRange`, so we return `nothing`.
    n == N || return nothing
    Expr(:call, Expr(:curly, :Contiguous, new_contig))
end

"""
contiguous_axis_indicator(::Type{T}) -> Tuple{Vararg{<:Val}}

Returns a tuple boolean `Val`s indicating whether that axis is contiguous.
"""
contiguous_axis_indicator(::Type{A}) where {D, A <: AbstractArray{<:Any,D}} = contiguous_axis_indicator(contiguous_axis(A), Val(D))
contiguous_axis_indicator(::A) where {A <: AbstractArray} = contiguous_axis_indicator(A)
Base.@pure contiguous_axis_indicator(::Contiguous{N}, ::Val{D}) where {N,D} = ntuple(d -> Val{d == N}(), Val{D}())

"""
If the contiguous dimension is not the dimension with `Stride_rank{1}`:
"""
struct ContiguousBatch{N} end
Base.@pure ContiguousBatch(N::Int) = ContiguousBatch{N}()
_get(::ContiguousBatch{N}) where {N} = N

"""
contiguous_batch_size(::Type{T}) -> ContiguousBatch{N}

Returns the Base.size of contiguous batches if `!isone(stride_rank(T, contiguous_axis(T)))`.
If `isone(stride_rank(T, contiguous_axis(T)))`, then it will return `ContiguousBatch{0}()`.
If `contiguous_axis(T) == -1`, it will return `ContiguousBatch{-1}()`.
If unknown, it will return `nothing`.
"""
contiguous_batch_size(x) = contiguous_batch_size(typeof(x))
contiguous_batch_size(::Type) = nothing
contiguous_batch_size(::Type{Array{T,N}}) where {T,N} = ContiguousBatch{0}()
contiguous_batch_size(::Type{<:Tuple}) = ContiguousBatch{0}()
contiguous_batch_size(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}) where {T,A<:AbstractVecOrMat{T}} = contiguous_batch_size(A)
contiguous_batch_size(::Type{<:PermutedDimsArray{T,N,I1,I2,A}}) where {T,N,I1,I2,A<:AbstractArray{T,N}} = contiguous_batch_size(A)
function contiguous_batch_size(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    _contiguous_batch_size(S, contiguous_batch_size(A), contiguous_axis(A))
end
_contiguous_batch_size(::Any, ::Any, ::Any) = nothing
@generated function _contiguous_batch_size(::Type{S}, ::ContiguousBatch{B}, ::Contiguous{C}) where {B,C,N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    if I.parameters[C] <: AbstractUnitRange
        Expr(:call, Expr(:curly, :ContiguousBatch, B))
    else
        Expr(:call, Expr(:curly, :ContiguousBatch, -1))
    end
end

struct StrideRank{R} end
Base.@pure StrideRank(R::NTuple{<:Any,Int}) = StrideRank{R}()
_get(::StrideRank{R}) where {R} = R
Base.collect(::StrideRank{R}) where {R} = collect(R)
@inline Base.getindex(::StrideRank{R}, i::Integer) where {R} = R[i]
@inline Base.getindex(::StrideRank{R}, ::Val{I}) where {R,I} = StrideRank{permute(R, I)}()

"""
rank_to_sortperm(::StrideRank) -> NTuple{N,Int}

Returns the `sortperm` of the stride ranks.
"""
function rank_to_sortperm(R::NTuple{N,Int}) where {N}
    sp = ntuple(zero, Val{N}())
    r = ntuple(n -> sum(R[n] .≥ R), Val{N}())
    @inbounds for n in 1:N
        sp = Base.setindex(sp, n, r[n])
    end
    sp
end
@generated Base.sortperm(::StrideRank{R}) where {R} = rank_to_sortperm(R)

stride_rank(x) = stride_rank(typeof(x))
stride_rank(::Type) = nothing
stride_rank(::Type{Array{T,N}}) where {T,N} = StrideRank{ntuple(identity, Val{N}())}()
stride_rank(::Type{<:Tuple}) = StrideRank{(1,)}()

stride_rank(::Type{B}) where {T, A <: AbstractVector{T}, B <: Union{Transpose{T,A},Adjoint{T,A}}} = StrideRank{(2, 1)}()
stride_rank(::Type{B}) where {T, A <: AbstractMatrix{T}, B <: Union{Transpose{T,A},Adjoint{T,A}}} = _stride_rank(B, stride_rank(A))
_stride_rank(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}, ::Nothing) where {T,A<:AbstractMatrix{T}} = nothing
_stride_rank(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}, rank) where {T,A<:AbstractMatrix{T}} = rank[Val{(2,1)}()]

stride_rank(::Type{B}) where {T,N,I1,I2,A<:AbstractArray{T,N},B<:PermutedDimsArray{T,N,I1,I2,A}} = _stride_rank(B, stride_rank(A))
_stride_rank(::Type{B}, ::Nothing) where {T,N,I1,I2,A<:AbstractArray{T,N},B<:PermutedDimsArray{T,N,I1,I2,A}} = nothing
_stride_rank(::Type{B}, rank) where {T,N,I1,I2,A<:AbstractArray{T,N},B<:PermutedDimsArray{T,N,I1,I2,A}} = rank[Val{I1}()]
function stride_rank(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    _stride_rank(S, stride_rank(A))
end
_stride_rank(::Any, ::Any) = nothing
@generated function _stride_rank(::Type{S}, ::StrideRank{R}) where {R,N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    rankv = collect(R)
    rank_new = Int[]
    n = 0
    for np in 1:NP
        r = rankv[np]
        if I.parameters[np] <: AbstractUnitRange
            n += 1
            push!(rank_new, r)
        end
    end
    # If n != N, then an axis was indexed by something other than an integer or `AbstractUnitRange`, so we return `nothing`.
    n == N || return nothing
    ranktup = Expr(:tuple); append!(ranktup.args, rank_new) # dynamic splats bad
    Expr(:call, Expr(:curly, :StrideRank, ranktup))
end
stride_rank(x, i) = stride_rank(x)[i]

"""
is_column_major(A) -> Val{true/false}()
"""
is_column_major(A) = is_column_major(stride_rank(A))
@generated function is_column_major(::StrideRank{R}) where {R}
    N = length(R)
    for n ∈ 2:N
        if R[n] ≤ R[n-1]
            return :(Val{false}())
        end
    end
    :(Val{true}())
end
is_column_major(::Nothing) = Val{false}()

struct DenseDims{D} end
Base.@pure DenseDims(D::NTuple{<:Any,Bool}) = DenseDims{D}()
@inline Base.getindex(::DenseDims{D}, i::Integer) where {D} = D[i]
@inline Base.getindex(::DenseDims{D}, ::Val{I}) where {D,I} = DenseDims{permute(D, I)}()
"""
dense_dims(::Type{T}) -> NTuple{N,Bool}

Returns a tuple of indicators for whether each axis is dense.
An axis `i` of array `A` is dense if `stride(A, i) * Base.size(A, i) == stride(A, j)` where `stride_rank(A)[i] + 1 == stride_rank(A)[j]`.
"""
dense_dims(x) = dense_dims(typeof(x))
dense_dims(::Type) = nothing
dense_dims(::Type{Array{T,N}}) where {T,N} = DenseDims{ntuple(_ -> true, Val{N}())}()
dense_dims(::Type{<:Tuple}) = DenseDims{(true,)}()
function dense_dims(::Type{<:Union{Transpose{T,A},Adjoint{T,A}}}) where {T,A<:AbstractMatrix{T}}
    dense = dense_dims(A)
    isnothing(dense) ? nothing : dense[Val{(2,1)}()]
end
function dense_dims(::Type{<:PermutedDimsArray{T,N,I1,I2,A}}) where {T,N,I1,I2,A<:AbstractArray{T,N}}
    dense = dense_dims(A)
    isnothing(dense) ? nothing : dense[Val{I1}()]
end
function dense_dims(::Type{S}) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    _dense_dims(S, dense_dims(A), stride_rank(A))
end
_dense_dims(::Any, ::Any) = nothing
@generated function _dense_dims(::Type{S}, ::DenseDims{D}, ::StrideRank{R}) where {D,R,N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}}
    still_dense = true
    sp = rank_to_sortperm(R)
    densev = Vector{Bool}(undef, NP)
    for np in 1:NP
        spₙ = sp[np]
        still_dense &= D[spₙ]
        densev[spₙ] = still_dense
        # a dim not being complete makes later dims not dense
        still_dense &= (I.parameters[spₙ] <: Base.Slice)::Bool
    end
    dense_tup = Expr(:tuple)
    for np in 1:NP
        spₙ = sp[np]
        if I.parameters[np] <: Base.Slice
            push!(dense_tup.args, densev[np])
        elseif I.parameters[np] <: AbstractUnitRange
            push!(dense_tup.args, densev[np])
        end
    end
    # If n != N, then an axis was indexed by something other than an integer or `AbstractUnitRange`, so we return `nothing`.
    length(dense_tup.args) == N ? Expr(:call, Expr(:curly, :DenseDims, dense_tup)) : nothing
end

permute(t::NTuple{N}, I::NTuple{N,Int}) where {N} = ntuple(n -> t[I[n]], Val{N}())
@generated function permute(t::Tuple{Vararg{Any,N}}, ::Val{I}) where {N,I}
    t = Expr(:tuple)
    foreach(i -> push!(t.args, Expr(:ref, :t, i)), I)
    Expr(:block, Expr(:meta, :inline), t)
end

"""
  size(A)

Returns the size of `A`. If the size of any axes are known at compile time,
these should be returned as `Static` numbers. For example:
```julia
julia> using StaticArrays, ArrayInterface

julia> A = @SMatrix rand(3,4);

julia> ArrayInterface.size(A)
(StaticInt{3}(), StaticInt{4}())
```
"""
size(A) = Base.size(A)
"""
  strides(A)

Returns the strides of array `A`. If any strides are known at compile time,
these should be returned as `Static` numbers. For example:
```julia
julia> A = rand(3,4);

julia> ArrayInterface.strides(A)
(StaticInt{1}(), 3)
```
"""
strides(A) = Base.strides(A)
"""
  offsets(A)

Returns offsets of indices with respect to 0. If values are known at compile time,
it should return them as `Static` numbers.
For example, if `A isa Base.Matrix`, `offsets(A) === (StaticInt(1), StaticInt(1))`.
"""
offsets(::Any) = (StaticInt{1}(),) # Assume arbitrary Julia data structures use 1-based indexing by default.
@inline strides(A::Vector{<:Any}) = (StaticInt(1),)
@inline strides(A::Array{<:Any,N}) where {N} = (StaticInt(1), Base.tail(Base.strides(A))...)
@inline strides(A::AbstractArray{<:Any,N}) where {N} = Base.strides(A)

@inline function offsets(x, i)
    inds = indices(x, i)
    start = known_first(inds)
    isnothing(start) ? first(inds) : StaticInt(start)
end
# @inline offsets(A::AbstractArray{<:Any,N}) where {N} = ntuple(n -> offsets(A, n), Val{N}())
# Explicit tuple needed for inference.
@generated function offsets(A::AbstractArray{<:Any,N}) where {N}
    quote
        $(Expr(:meta, :inline))
        Base.Cartesian.@ntuple $N n -> offsets(A, n)
    end
end

@inline size(B::Union{Transpose{T,A},Adjoint{T,A}}) where {T,A<:AbstractMatrix{T}} = permute(size(parent(B)), Val{(2,1)}())
@inline size(B::PermutedDimsArray{T,N,I1,I2,A}) where {T,N,I1,I2,A<:AbstractArray{T,N}} = permute(size(parent(B)), Val{I1}())
@inline size(A::AbstractArray, ::StaticInt{N}) where {N} = size(A)[N]
@inline size(A::AbstractArray, ::Val{N}) where {N} = size(A)[N]
@inline strides(B::Union{Transpose{T,A},Adjoint{T,A}}) where {T,A<:AbstractMatrix{T}} = permute(strides(parent(B)), Val{(2,1)}())
@inline strides(B::PermutedDimsArray{T,N,I1,I2,A}) where {T,N,I1,I2,A<:AbstractArray{T,N}} = permute(strides(parent(B)), Val{I1}())
@inline stride(A::AbstractArray, ::StaticInt{N}) where {N} = strides(A)[N]
@inline stride(A::AbstractArray, ::Val{N}) where {N} = strides(A)[N]
stride(A, i) = Base.stride(A, i)

size(B::S) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}} = _size(size(parent(B)), B.indices, map(static_length, B.indices))
strides(B::S) where {N,NP,T,A<:AbstractArray{T,NP},I,S <: SubArray{T,N,A,I}} = _strides(strides(parent(B)), B.indices)
@generated function _size(A::Tuple{Vararg{Any,N}}, inds::I, l::L) where {N, I<:Tuple, L}
    t = Expr(:tuple)
    for n in 1:N
        if (I.parameters[n] <: Base.Slice)
            push!(t.args, :(@inbounds(_try_static(A[$n], l[$n]))))
        elseif I.parameters[n] <: AbstractUnitRange
            push!(t.args, Expr(:ref, :l, n))
        end
    end
    Expr(:block, Expr(:meta, :inline), t)
end
@generated function _strides(A::Tuple{Vararg{Any,N}}, inds::I) where {N, I<:Tuple}
    t = Expr(:tuple)
    for n in 1:N
        if I.parameters[n] <: AbstractRange
            push!(t.args, Expr(:ref, :A, n))
        elseif !(I.parameters[n] <: Integer)
            return nothing
        end
    end
    Expr(:block, Expr(:meta, :inline), t)
end


text:
p = BigFloat(pi)
e = exp(BigFloat(1))
table = [
    ["PI", p],
    ["E", e],
    ["LOG2_E", log2(e)],
    ["LN_10", log(big"10")],
    ["LN_2", log(big"2")],
    ["LOG10_2", log10(big"2")],
    ["LOG10_E", log10(e)],
    ["LOG2_10", log2(big"10")],
    ["FRAC_1_PI", 1 / p],
    ["FRAC_1_SQRT_2", 1 / sqrt(big"2")],
    ["FRAC_2_PI", 2 / p],
    ["FRAC_2_SQRT_PI", 2 / sqrt(p)],
    ["FRAC_PI_2", p / 2],
    ["FRAC_PI_3", p / 3],
    ["FRAC_PI_4", p / 4],
    ["FRAC_PI_6", p / 6],
    ["FRAC_PI_8", p / 8],
    ["SQRT_2", sqrt(big"2")],
    ["TAU", p * 2],
]

for (name, value) in table
    println("#[inline]")
    hi = Float64(value)
    lo = Float64(value - hi)
    println("fn $name() -> Self { Self { hi: $hi, lo: $lo } }")
    println("")
end

max_stars_repo_name: roosephu/apc

max_issues_repo_name: roosephu/apc

max_forks_repo_name: roosephu/apc


text:
struct BlockIndexRange{N,R<:NTuple{N,AbstractUnitRange{Int}}}
    block::Block{N,Int}
    indices::R
end

"""
    BlockIndexRange(block, startind:stopind)

represents a cartesian range inside a block.
"""
BlockIndexRange

BlockIndexRange(block::Block{N}, inds::NTuple{N,AbstractUnitRange{Int}}) where {N} =
    BlockIndexRange{N,typeof(inds)}(inds)
BlockIndexRange(block::Block{N}, inds::Vararg{AbstractUnitRange{Int},N}) where {N} =
    BlockIndexRange(block,inds)

getindex(B::Block{N}, inds::Vararg{Int,N}) where N = BlockIndex(B,inds)
getindex(B::Block{N}, inds::Vararg{AbstractUnitRange{Int},N}) where N = BlockIndexRange(B,inds)

eltype(R::BlockIndexRange) = eltype(typeof(R))
eltype(::Type{BlockIndexRange{N}}) where {N} = BlockIndex{N}
eltype(::Type{BlockIndexRange{N,R}}) where {N,R} = BlockIndex{N}
IteratorSize(::Type{<:BlockIndexRange}) = Base.HasShape{1}()

first(iter::BlockIndexRange) = BlockIndex(iter.block.n, map(first, iter.indices))
last(iter::BlockIndexRange)  = BlockIndex(iter.block.n, map(last, iter.indices))

@inline function iterate(iter::BlockIndexRange)
    iterfirst, iterlast = first(iter), last(iter)
    if any(map(>, iterfirst.α, iterlast.α))
        return nothing
    end
    iterfirst, iterfirst
end
@inline function iterate(iter::BlockIndexRange, state)
    nextstate = BlockIndex(state.I, inc(state.α, first(iter).α, last(iter).α))
    nextstate.α[end] > last(iter.indices[end]) && return nothing
    nextstate, nextstate
end

size(iter::BlockIndexRange) = map(dimlength, first(iter).α, last(iter).α)
length(iter::BlockIndexRange) = prod(size(iter))

Block(bs::BlockIndexRange) = bs.block

### Views
Block(bs::BlockSlice{<:BlockIndexRange}) = Block(bs.block)

function _unblock(cum_sizes, I::Tuple{BlockIndexRange{1,R}, Vararg{Any}}) where {R}
    B = Block(first(I))
    range = cum_sizes[Int(B)]-1 .+ first(I).indices[1]

    BlockSlice(I[1], range)
end

to_index(::BlockIndexRange) = throw(ArgumentError("BlockIndexRange must be converted by to_indices(...)"))

@inline to_indices(A, inds, I::Tuple{BlockIndexRange{1,R}, Vararg{Any}}) where R =
    (unblock(A, inds, I), to_indices(A, _maybetail(inds), tail(I))...)

# splat out higher dimensional blocks
# this mimics view of a CartesianIndex
@inline to_indices(A, inds, I::Tuple{BlockIndexRange, Vararg{Any}}) =
    to_indices(A, inds, (BlockRange.(Block.(I[1].block.n), tuple.(I[1].indices))..., tail(I)...))

# In 0.7, we need to override to_indices to avoid calling linearindices
@inline to_indices(A, I::Tuple{BlockIndexRange, Vararg{Any}}) =
    to_indices(A, axes(A), I)

if VERSION >= v"1.2-"  # See also `reindex` definitions in views.jl
reindex(idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockIndexRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockIndexRange(Block(idxs[1].block.indices[1][Int(subidxs[1].block.block)]),
                                                            subidxs[1].block.indices),
                                            idxs[1].indices[subidxs[1].indices]),
                                 reindex(tail(idxs), tail(subidxs))...))
else  # if VERSION >= v"1.2-"
reindex(V, idxs::Tuple{BlockSlice{<:BlockRange}, Vararg{Any}},
        subidxs::Tuple{BlockSlice{<:BlockIndexRange}, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (BlockSlice(BlockIndexRange(Block(idxs[1].block.indices[1][Int(subidxs[1].block.block)]),
                                                            subidxs[1].block.indices),
                                            idxs[1].indices[subidxs[1].indices]),
                                    reindex(V, tail(idxs), tail(subidxs))...))
end  # if VERSION >= v"1.2-"

# De-reference blocks before creating a view to avoid taking `global2blockindex`
# path in `AbstractBlockStyle` broadcasting.
@inline function Base.unsafe_view(
        A::BlockArray{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    B = A[map(x -> x.block.block, I)...]
    return view(B, _splatmap(x -> x.block.indices, I)...)
end

@inline function Base.unsafe_view(
        A::PseudoBlockArray{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    return view(A.blocks, map(x -> x.indices, I)...)
end

@inline function Base.unsafe_view(
        A::ReshapedArray{<:Any, N, <:AbstractBlockArray{<:Any, M}},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N, M}
    @_propagate_inbounds_meta
    # Note: assuming that I[M+1:end] are verified to be singletons
    return reshape(view(A.parent, I[1:M]...), Val(N))
end

@inline function Base.unsafe_view(
        A::Array{<:Any, N},
        I::Vararg{BlockSlice{<:BlockIndexRange{1}}, N}) where {N}
    @_propagate_inbounds_meta
    return view(A, map(x -> x.indices, I)...)
end

# #################
# # support for pointers
# #################
#
# function unsafe_convert(::Type{Ptr{T}},
#                         V::SubArray{T, N, BlockArray{T, N, AT}, NTuple{N, BlockSlice{Block{1,Int}}}}) where AT <: AbstractArray{T, N} where {T,N}
#     unsafe_convert(Ptr{T}, parent(V).blocks[Int.(Block.(parentindices(V)))...])
# end


text:
#################### Monte Carlo Standard Errors ####################

function mcse{T<:Real}(x::Vector{T}, method::Symbol=:imse; args...)
  method == :bm ? mcse_bm(x; args...) :
  method == :imse ? mcse_imse(x) :
  method == :ipse ? mcse_ipse(x) :
    throw(ArgumentError("unsupported mcse method $method"))
end

function mcse_bm{T<:Real}(x::Vector{T}; size::Integer=100)
  n = length(x)
  m = div(n, size)
  m >= 2 ||
    throw(ArgumentError(
      "iterations are < $(2 * size) and batch size is > $(div(n, 2))"
    ))
  mbar = [mean(x[i * size + (1:size)]) for i in 0:(m - 1)]
  sem(mbar)
end

function mcse_imse{T<:Real}(x::Vector{T})
  n = length(x)
  m = div(n - 2, 2)
  ghat = autocov(x, [0, 1])
  Ghat = sum(ghat)
  value = -ghat[1] + 2 * Ghat
  for i in 1:m
    Ghat = min(Ghat, sum(autocov(x, [2 * i, 2 * i + 1])))
    Ghat > 0 || break
    value += 2 * Ghat
  end
  sqrt(value / n)
end

function mcse_ipse{T<:Real}(x::Vector{T})
  n = length(x)
  m = div(n - 2, 2)
  ghat = autocov(x, [0, 1])
  value = ghat[1] + 2 * ghat[2]
  for i in 1:m
    Ghat = sum(autocov(x, [2 * i, 2 * i + 1]))
    Ghat > 0 || break
    value += 2 * Ghat
  end
  sqrt(value / n)
end


text:
"""
VWN5 LDA correlation according to Vosko, Wilk, and Nusair, (DOI 10.1139/p80-159).
"""
function energy_per_particle(::Val{:lda_c_vwn}, ρ)
    # From https://math.nist.gov/DFTdata/atomdata/node5.html
    A = 0.0310907
    x0 = -0.10498
    b = 3.72744
    c = 12.9352
    rₛ = cbrt(3/(4π*ρ)) # τ in the above link
    x = sqrt(rₛ)
    Xx = x^2 + b*x + c
    Xx0 = x0^2 + b*x0 + c
    Q = sqrt(4c-b^2)
    A * (log(x^2 / Xx) + 2b/Q*atan(Q/(2x+b)) - b*x0/Xx0*(log((x-x0)^2/Xx) + 2*(b+2x0)/Q*atan(Q/(2x+b))))
end


text:
using BenchmarkTools
using DifferentialEquations

function main()
    x1 = 0.00001
    xf = 11.0
    eps = 1.0E-15

    y1 = [ 0.99998416139571, -1.58175227379914 ]

    u0 = [y1[1]; y1[2]]
    tspan = (x1, xf)
    prob = ODEProblem(equation!, u0, tspan)
    sol = solve(prob, abstol = eps, reltol = eps)
    y = [sol.u[size(sol.u)[1]][1], sol.u[size(sol.u)[1]][2]]

    return y
end

function equation!(du, u, p, t)
    du[1] = u[2]
    du[2] = u[1] * sqrt(u[1] / t)
end

@btime main()


text:
module FunctionalTables

using ArgCheck: @argcheck
using Base: issingletontype
using DocStringExtensions: SIGNATURES, TYPEDEF
using Parameters: @unpack
using IterTools: @ifsomething, imap
import Tables

include("keys.jl")
include("utilities.jl")
include("ordering.jl")
include("columns.jl")
include("functionaltables.jl")
include("tables-interface.jl")
include("sort.jl")
include("by.jl")

end # module


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

"""
    struct FloydWarshallState{T, U}

An [`AbstractPathState`](@ref) designed for Floyd-Warshall shortest-paths calculations.
"""
struct FloydWarshallState{T, U<:Integer}<:AbstractPathState
    dists::Matrix{T}
    parents::Matrix{U}
end

@doc_str """
floyd_warshall_shortest_paths(g, distmx=DefaultDistance())
Use the [Floyd-Warshall algorithm](http://en.wikipedia.org/wiki/Floyd–Warshall_algorithm)
to compute the shortest paths between all pairs of vertices in graph `g` using an
optional distance matrix `distmx`. Return a [`FloydWarshallState`](@ref) with relevant
traversal information.

### Performance
Space complexity is on the order of ``\\mathcal{O}(|V|^2)``.
"""
function floyd_warshall_shortest_paths{T}(
    g::AbstractGraph,
    distmx::AbstractMatrix{T} = DefaultDistance()
)
    U = eltype(g)
    n_v = nv(g)
    dists = fill(typemax(T), (Int(n_v),Int(n_v)))
    parents = zeros(U, (Int(n_v),Int(n_v)))

    # fws = FloydWarshallState(Matrix{T}(), Matrix{Int}())
    for v in 1:n_v
        dists[v,v] = zero(T)
    end
    undirected = !is_directed(g)
    for e in edges(g)
        u = src(e)
        v = dst(e)

        d = distmx[u,v]

        dists[u,v] = min(d, dists[u,v])
        parents[u,v] = u
        if undirected
            dists[v,u] = min(d, dists[v,u])
            parents[v,u] = v
        end
    end
    for w in vertices(g), u in vertices(g), v in vertices(g)
        if dists[u,w] == typemax(T) || dists[w,v] == typemax(T)
            ans = typemax(T)
        else
            ans = dists[u,w] + dists[w,v]
        end
        if dists[u,v] > ans
            dists[u,v] = dists[u,w] + dists[w,v]
            parents[u,v] = parents[w,v]
        end
    end
    fws = FloydWarshallState(dists, parents)
    return fws
end

function enumerate_paths(s::FloydWarshallState{T, U}, v::Integer) where T where U<:Integer
    pathinfo = s.parents[v,:]
    paths = Vector{Vector{U}}()
    for i in 1:length(pathinfo)
        if (i == v) || (s.dists[v,i] == typemax(T))
            push!(paths, Vector{U}())
        else
            path = Vector{U}()
            currpathindex = i
            while currpathindex != 0
                push!(path,currpathindex)
                currpathindex = pathinfo[currpathindex]
            end
            push!(paths, reverse(path))
        end
    end
    return paths
end

enumerate_paths(s::FloydWarshallState) = [enumerate_paths(s, v) for v in 1:size(s.parents,1)]
enumerate_paths(st::FloydWarshallState, s::Integer, d::Integer) = enumerate_paths(st, s)[d]

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
"""
    SMC(n_particles::Int)

Sequential Monte Carlo sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
SMC(1000)
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

sample(gdemo([1.5, 2]), SMC(1000))
```
"""
mutable struct SMC{T, F} <: InferenceAlgorithm
  n_particles           ::  Int
  resampler             ::  F
  resampler_threshold   ::  Float64
  space                 ::  Set{T}
  gid                   ::  Int
end
SMC(n) = SMC(n, resample_systematic, 0.5, Set(), 0)
function SMC(n_particles::Int, space...)
  _space = isa(space, Symbol) ? Set([space]) : Set(space)
  SMC(n_particles, resample_systematic, 0.5, _space, 0)
end
SMC(alg::SMC, new_gid::Int) = SMC(alg.n_particles, alg.resampler, alg.resampler_threshold, alg.space, new_gid)

Sampler(alg::SMC) = begin
  info = Dict{Symbol, Any}()
  info[:logevidence] = []
  Sampler(alg, info)
end

step(model, spl::Sampler{<:SMC}, vi::VarInfo) = begin
    particles = ParticleContainer{Trace}(model)
    vi.num_produce = 0;  # Reset num_produce before new sweep\.
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)

    push!(particles, spl.alg.n_particles, spl, vi)

    while consume(particles) != Val{:done}
      ess = effectiveSampleSize(particles)
      if ess <= spl.alg.resampler_threshold * length(particles)
        resample!(particles,spl.alg.resampler)
      end
    end

    ## pick a particle to be retained.
    Ws, _ = weights(particles)
    indx = randcat(Ws)
    push!(spl.info[:logevidence], particles.logE)

    particles[indx].vi
end

## wrapper for smc: run the sampler, collect results.
function sample(model::Model, alg::SMC)
  spl = Sampler(alg);

  particles = ParticleContainer{Trace}(model)
  push!(particles, spl.alg.n_particles, spl, VarInfo())

  while consume(particles) != Val{:done}
    ess = effectiveSampleSize(particles)
    if ess <= spl.alg.resampler_threshold * length(particles)
      resample!(particles,spl.alg.resampler)
    end
  end
  w, samples = getsample(particles)
  res = Chain(w, samples)

end


text:
function vfilter!(f::F, x::Vector{T}, y::AbstractArray{T}) where {F,T <: NativeTypes}
    W, Wshift = VectorizationBase.pick_vector_width_shift(T)
    N = length(y)
    Nrep = N >>> Wshift
    Nrem = N & (W - 1)
    j = 0
    st = VectorizationBase.static_sizeof(T)
    zero_index = MM(W, Static(0), st)
    GC.@preserve x y begin
        ptr_x = pointer(x)
        ptr_y = pointer(y)
        for _ ∈ 1:Nrep
            vy = vload(ptr_y, zero_index)
            mask = f(vy)
            VectorizationBase.compressstore!(gep(ptr_x, VectorizationBase.lazymul(st, j)), vy, mask)
            ptr_y = gep(ptr_y, register_size())
            j = vadd_fast(j, count_ones(mask))
        end
        rem_mask = VectorizationBase.mask(T, Nrem)
        vy = vload(ptr_y, zero_index, rem_mask)
        mask = rem_mask & f(vy)
        VectorizationBase.compressstore!(gep(ptr_x, VectorizationBase.lazymul(st, j)), vy, mask)
        j = vadd_fast(j, count_ones(mask))
        Base._deleteend!(x, N-j) # resize!(x, j)
    end
    x
end
vfilter!(f::F, x::Vector{T}) where {F, T<:NativeTypes} = vfilter!(f, x, x)
vfilter(f::F, y::AbstractArray{T}) where {F, T<:NativeTypes} = vfilter!(f, Vector{T}(undef, length(y)), y)
vfilter(f::F, y) where {F} = filter(f, y)
vfilter!(f::F, y) where {F} = filter!(f, y)

"""
    vfilter(f, a::AbstractArray)

SIMD-vectorized `filter`, returning an array containing the elements of `a` for which `f` return `true`.

This function requires AVX512 to be faster than `Base.filter`, as it adds compressstore instructions.
"""
vfilter

"""
    vfilter!(f, a::AbstractArray)

SIMD-vectorized `filter!`, removing the element of `a` for which `f` is false.
"""
vfilter!


text:
using SimpleTraits.BaseTraits, SimpleTraits

export Sum

struct Sum{N} <: CompositeBlock{N}
    list::Vector{AbstractBlock{N}}

    Sum{N}(list::Vector{AbstractBlock{N}}) where N = new{N}(list)
    Sum{N}(it::T) where {N, T} = Sum{N}(SimpleTraits.trait(IsIterator{T}), it)
    Sum{N}(::Type{<:IsIterator}, it) where N = new{N}(collect(AbstractBlock{N}, it))
end

Sum{N}(::Not, it) where N = error("expect an iterator/collection")

Sum{N}() where N = Sum(AbstractBlock{N}[])
Sum(blocks::Vector{<:AbstractBlock{N}}) where N = Sum{N}(blocks)
Sum(blocks::AbstractBlock{N}...) where N = Sum(collect(AbstractBlock{N}, blocks))

mat(::Type{T}, x::Sum) where T = mapreduce(x->mat(T, x), +, x.list)

chsubblocks(x::Sum{N}, it) where N = Sum{N}(it)

function apply!(r::AbstractRegister, x::Sum)
    isempty(x.list) && return r

    out = copy(r)
    apply!(out, first(x))
    for k in 2:length(x)
        out += apply!(copy(r), x[k])
    end
    copyto!(r, out)
    return r
end

export Sum

subblocks(x::Sum) = x.list
cache_key(x::Sum) = map(cache_key, x.list)

Base.length(x::Sum) = length(x.list)
Base.iterate(x::Sum) = iterate(x.list)
Base.iterate(x::Sum, st) = iterate(x.list, st)
Base.getindex(x::Sum, k) = getindex(x.list, k)

function Base.:(==)(lhs::Sum{N}, rhs::Sum{N}) where N
    for (a, b) in zip(lhs, rhs)
        a == b || return false
    end
    return true
end


text:
function solveSLAE(A)
    rangA = rang(A, isStep = true)
    rows, columns = size(A)
    if rangA == rows
        return zeros(rows)
    end
    k = 0
    for i in rows:-1:1
        if A[i, i] != 0
            k = i
            break
        end
    end

    A[k, k:columns] = A[k, k:columns]/A[k,k]
    for i in k-1:-1:firstindex(b)
        for j in columns:-1:i+1
            if isapproxzero(A[i, j])
                A[i,j] = 0
                continue
            end
            A[i, j:columns] = A[i, j:columns] - A[j, j:columns]/A[j,j]* A[i,j]
        end
        A[i, i:columns] = A[i, i:columns] / A[i,i]
    end

    c = []
    for i in 1:columns-k
        push!(c, "c"*"$i")
    end
    FRS = []
    for i in k:-1:1
        newline = "x$i = "
        for j in k+1:columns
            b = A[i,j]
            if b > 0
                newline = newline * " - $b *" *c[j-k]
            elseif  b < 0
                b *= -1
                newline = newline * " + $b * "*c[j-k]
            end
        end
        if newline == "x$i = "
            newline = newline * "0"
        end
        push!(FRS, newline)
    end

    for i in k+1:columns
        pushfirst!(FRS, "x$i = " * c[i - k])
    end

    return FRS
end


text:
include("sorted.jl")
function sortkey(key_values, a)
    indperm = sortperm(key_values)
    return a[indperm]
end

function sortkey(f::Function, a)
    key_values = []
    for i in a
        append!(key_values, f(i))
    end
    sortkey(key_values, a)
end

function sortkeymatrix(A::Matrix)
    vec_of_colons = [@view A[:,j] for j in 1:size(A,2)]
    sorted1_vec_of_colons = sortkey(a -> sum(a), vec_of_colons)
    sorted2_vec_of_colons = sortkey(a -> sum(map(iszero, a)), vec_of_colons)
    hcat(sorted1_vec_of_colons...)
    hcat(sorted2_vec_of_colons...)
end

function calcsort(values, array)
    pom = zeros(values[length(values)])
    for i in 1:length(array)
        pom[array[i]] += 1
    end
    pos = 1
    for i in 1:length(values)
        for _ in 1:pom[i]
            array[pos] = values[i]
            pos += 1
        end
    end
    return array
end

function calcsort(head, tail, array)
    values = []
    for i in head:tail
        append!(values, i)
    end
    return calcsort(values, array)
end

function insertsortperm(A)
    indecies = []
    for i in 1:length(A)
        append!(indecies, i)
    end
    for j = 2:length(A)
        key = A[j]
        ikey = indecies[j]
        i = j-1
        while (i > 0 && A[i] > key)
            A[i + 1] = A[i]
            indecies[i + 1] = indecies[i]
            i = i - 1
        end
        A[i+1] = key
        indecies[i+1] = ikey
    end
    println(A)
end

function insertsort(A)
    for j = 2:length(A)
        key = A[j]
        i = j-1
        while (i > 0 && A[i] > key)
            A[i + 1] = A[i]
            i = i - 1
        end
        A[i+1] = key
    end
    println(A)
end

function insertsortbinary(A)
    for j = 2:length(A)
        key = A[j]
        i = searchsortedfirst(A[1:j-1], key)
        for k in j-1:-1:i
            A[k + 1] = A[k]
        end
        A[i] = key
    end
    println(A)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
using aoc2019.wiring: manhattan_distance, signal_delay
using DelimitedFiles

input = joinpath(@__DIR__, "input")
ops = readdlm(input, String)

left = split(ops[1], ',')
right = split(ops[2], ',')

p1 = manhattan_distance(left, right)
@assert p1 == 1626

p2 = signal_delay(left, right)
@assert p2 == 27330

print("-----------------------------------------------------------------------\n")
print("crossed wires -- part one\n    distance : $p1\n")
print("crossed wires -- part two\n    delay    : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
## Makes the model associated with a fully adapted auxiliary particle filter for
## a linear Gaussian state space model. This was introduced in:
## M. K. Pitt and N. Shephard. Filtering via simulation: Auxiliary particle
## filters. J. Amer. Statist. Assoc., 94(446):590--599, 1999.

function makeLGAPFModel(theta::LGTheta, ys::Vector{Float64})
  n::Int64 = length(ys)

  RC2v0::Float64 = theta.R + theta.C * theta.C * theta.v0
  v0RoverRC2v0::Float64 = theta.v0*theta.R/RC2v0
  sqrtv0RoverRC2v0::Float64 = sqrt(v0RoverRC2v0)

  mu1::Float64 = v0RoverRC2v0*(theta.x0/theta.v0 + theta.C*ys[1]/theta.R)
  tmp1::Float64 = theta.C * theta.x0 - ys[1]
  lG1::Float64 = -0.5 * log(2 * π * RC2v0) - tmp1 * 0.5/RC2v0 * tmp1

  RC2Q::Float64 = theta.R + theta.C * theta.C * theta.Q
  invRC2Qover2::Float64 = 0.5/RC2Q
  QRoverRC2Q::Float64 = theta.Q * theta.R / RC2Q
  sqrtQRoverRC2Q::Float64 = sqrt(QRoverRC2Q)

  invRover2::Float64 = 0.5/theta.R
  logncG::Float64 = -0.5 * log(2 * π * RC2Q)

  @inline function lG(p::Int64, particle::Float64Particle, ::Nothing)
    if p == n
      return p == 1 ? lG1 : 0.0
    end
    v::Float64 = theta.C * theta.A * particle.x - ys[p+1]
    v = logncG - v * invRC2Qover2 * v
    p == 1 && (v += lG1)
    return v
  end
  @inline function M!(newParticle::Float64Particle, rng::RNG, p::Int64,
    particle::Float64Particle, ::Nothing)
    if p == 1
      newParticle.x = mu1 + sqrtv0RoverRC2v0*randn(rng)
    else
      mu::Float64 = QRoverRC2Q*(theta.A*particle.x/theta.Q + theta.C*ys[p]/theta.R)
      newParticle.x = mu + sqrtQRoverRC2Q*randn(rng)
    end
  end
  return SMCModel(M!, lG, n, Float64Particle, Nothing)
end


text:
"""
    ConfigTypes

Module containing ClimateMachine configuration types.
"""
module ConfigTypes

export ClimateMachineConfigType,
    AtmosLESConfigType,
    AtmosGCMConfigType,
    OceanBoxGCMConfigType,
    OceanSplitExplicitConfigType,
    SingleStackConfigType

abstract type ClimateMachineConfigType end
struct AtmosLESConfigType <: ClimateMachineConfigType end
struct AtmosGCMConfigType <: ClimateMachineConfigType end
struct OceanBoxGCMConfigType <: ClimateMachineConfigType end
struct OceanSplitExplicitConfigType <: ClimateMachineConfigType end
struct SingleStackConfigType <: ClimateMachineConfigType end

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
using Colors, ColorTypes, FixedPointNumbers

import Base: convert

function make_bitvector(v::Vector{UInt8})
    siz = sizeof(v)
    bv = falses(siz<<3)
    unsafe_copyto!(reinterpret(Ptr{UInt8}, pointer(bv.chunks)), pointer(v), siz)
    bv
end

"""
    convertY10BtoU16(imbuff::Vector{UInt8})
Convert a packed 10 bit grayscale buffer [`::Vector{UInt8}`] to 16 bit padded buffer [`::Vector{UInt16}`]
[x x x x x x b9 b8 b7 b6 b5 b4 b3 b2 b1 b0]
"""
function convertY10BtoU16(imbuff::Vector{UInt8} )

    L = length(imbuff)
    bv = reverse(reshape(make_bitvector(imbuff[(L):-1:1]),10,:)',dims=1)

    #convert bit array back to Array{UInt16}
    values = zeros(UInt16, size(bv,1))
    for i = 1:size(bv,1)
        val = 0x0000
        for j = 1:10
            val += bv[i,j] << (j-1)
        end
        values[i] = val
    end
    return values
end

"""
    convertUYVYtoYCbCr(imbuff::Vector{UInt8}, width::Int, height::Int)

Convert a UYVY 4:2:2 encoded YCbCr buffer [`::Vector{UInt8}`] to full YCbCr
"""
function convertUYVYtoYCbCr(imbuff::Vector{UInt8}, width::Int, height::Int)

    halfwidth = div(width,2)

    imY = reshape(imbuff[2:2:end],(width,height))'

    #recover colors
    imCbCr = imbuff[1:2:end]

    imCb = kron(reshape(imCbCr[1:2:end], (halfwidth,height))', [1 1])

    imCr = kron(reshape(imCbCr[2:2:end], (halfwidth,height))', [1 1])

    #create YCbCr array
    return YCbCr.(imY, imCb, imCr)

end

"""
    convertYUYVtoYCbCr(imbuff::Vector{UInt8}, width::Int, height::Int)

Convert a UYVY 4:2:2 encoded YCbCr buffer [`::Vector{UInt8}`] to full YCbCr
"""
function convertYUYVtoYCbCr(imbuff::Vector{UInt8}, width::Int, height::Int)

    halfwidth = div(width,2)

    imY = reshape(imbuff[1:2:end],(width,height))'

    #recover colors
    imCbCr = imbuff[2:2:end]

    imCb = kron(reshape(imCbCr[1:2:end], (halfwidth,height))', [1 1])

    imCr = kron(reshape(imCbCr[2:2:end], (halfwidth,height))', [1 1])

    #create YCbCr array
    return YCbCr.(imY, imCb, imCr)

end

"""
    convertUYVYtoArray(imbuff::Vector{UInt8}, width::Int, height::Int)

Convert a 4:2:2 encoded YCbCr buffer [`::Vector{UInt8}`] to full 4:4:4 buffer [`::Matrix{UInt8}[3]`]
"""
function convertUYVYtoArray(imbuff::Vector{UInt8}, width::Int, height::Int)

    halfwidth = div(width,2)

    imY = reshape(imbuff[2:2:end],(width,height))'

    #recover colors
    imCbCr = imbuff[1:2:end]

    imCb = kron(reshape(imCbCr[1:2:end], (halfwidth,height))', [1 1])

    imCr = kron(reshape(imCbCr[2:2:end], (halfwidth,height))', [1 1])

    #create YCbCr array
    return cat(3, imY, imCb, imCr)

end

"""
    convertYUYVtoArray(imbuff::Vector{UInt8}, width::Int, height::Int)

Convert a 4:2:2 encoded YCbCr buffer [`::Vector{UInt8}`] to full 4:4:4 buffer [`::Matrix{UInt8}[3]`]
"""
function convertYUYVtoArray(imbuff::Vector{UInt8}, width::Int, height::Int)

    halfwidth = div(width,2)

    imY = reshape(imbuff[1:2:end],(width,height))'

    #recover colors
    imCbCr = imbuff[2:2:end]

    imCb = kron(reshape(imCbCr[1:2:end], (halfwidth,height))', [1 1])

    imCr = kron(reshape(imCbCr[2:2:end], (halfwidth,height))', [1 1])

    #create YCbCr array
    return cat(imY, imCb, imCr, dims=3)

end

"""
    UYVYextractY(imbuff::Vector{UInt8}, width::Int, height::Int)

Extract the Y component out of a 4:2:2 YCbCr buffer
"""
function UYVYextractY(imbuff::Vector{UInt8}, width::Int, height::Int)
    reshape(imbuff[2:2:end],(width,height))'
end

"""
    UYVYextractY(imbuff::Vector{UInt8}, width::Int, height::Int)

Extract the Y component out of a 4:2:2 YCbCr buffer
"""
function YUYVextractY(imbuff::Vector{UInt8}, width::Int, height::Int)
    reshape(imbuff[1:2:end],(width,height))'
end

##

abstract type AbstractEncodings end

struct UYVYonlyY <: AbstractEncodings
    width::Int
    height::Int
    depth::Int
    datatype::DataType
end
UYVYonlyY(width::Int, height::Int) = UYVYonlyY(width, height, 1, UInt8)

function (decoder::UYVYonlyY)(imy::Array{UInt8,2})
    imbuff = copy_buffer_bytes(decoder.width * decoder.height * 2)
    imy[:,:] = reshape(imbuff[2:2:end],(decoder.width, decoder.height))'
    return nothing
end

struct UYVY <: AbstractEncodings
    width::Int
    height::Int
    depth::Int
    datatype::DataType
end
UYVY(width::Int, height::Int) = UYVY(width, height, 3, UInt8)

function (decoder::UYVY)(im444::Array{UInt8,3})
    imbuff = copy_buffer_bytes(decoder.width * decoder.height * 2)
    im444[:] = convertUYVYtoArray(imbuff, decoder.width, decoder.height)
    return nothing
end

struct YUYVonlyY <: AbstractEncodings
    width::Int
    height::Int
    depth::Int
    datatype::DataType
end
YUYVonlyY(width::Int, height::Int) = YUYVonlyY(width, height, 1, UInt8)

function (decoder::YUYVonlyY)(imy::Array{UInt8,2})
    imbuff = copy_buffer_bytes(decoder.width * decoder.height * 2)
    imy[:,:] = reshape(imbuff[1:2:end],(decoder.width, decoder.height))'
    return nothing
end

struct YUYV <: AbstractEncodings
    width::Int
    height::Int
    depth::Int
    datatype::DataType
end
YUYV(width::Int, height::Int) = YUYV(width, height, 3, UInt8)

function (decoder::YUYV)(im444::Array{UInt8,3})
    imbuff = copy_buffer_bytes(decoder.width * decoder.height * 2)
    im444[:] = convertYUYVtoArray(imbuff, decoder.width, decoder.height)
    return nothing
end

struct Y10B <: AbstractEncodings
    width::Int
    height::Int
    depth::Int
    datatype::DataType
end
Y10B(width::Int, height::Int) = Y10B(width, height, 1, UInt16)

function (decoder::Y10B)(imy10b::Array{UInt16,2})
    imbuff = copy_buffer_bytes(div(decoder.width * decoder.height * 10,8))
    imy10b[:] = reshape(convertY10BtoU16(imbuff), (decoder.width, decoder.height))'
    return nothing
end

#Convertions
function convert(::Type{Gray}, rf::RawFrame)
    if rf.pixelformat == :UYVY
        return reinterpret(Gray{N0f8}, UYVYextractY(rf.data, rf.width, rf.height))
    elseif rf.pixelformat == :YUYV
        return reinterpret(Gray{N0f8}, YUYVextractY(rf.data, rf.width, rf.height))
    elseif rf.pixelformat == :Y10B
        buflenght = div(rf.width * rf.height * 10, 8) #10 bits per pixel / 8 bits per byte
        return reinterpret(Gray{N0f16}, reshape(  convertY10BtoU16(rf.data[1:buflenght])  ,(rf.width,rf.height))')
    else
        error("Frame format $(rf.pixelformat) possibly not implemented yet")
    end
end

function convert(::Type{YCbCr}, rf::RawFrame)
    if rf.pixelformat == :UYVY
        return convertUYVYtoYCbCr(rf.data, rf.width, rf.height)
    elseif rf.pixelformat == :YUYV
        return convertYUYVtoYCbCr(rf.data, rf.width, rf.height)
    else
        error("Frame format $(rf.pixelformat) possibly not implemented yet, or not supported.")
    end
end


text:
# function to generate a Kitaev colorcode
function generateBondColorcodeKitaev(
        lattice :: L
    ) :: Dict{LB,Vector{Int64}} where {LS,D,LB,N,S<:AbstractSite{LS,D},B<:AbstractBond{LB,N},U,L<:AbstractLattice{S,B,U}}

    # put the default kitaev labels with default colors into a dictonary
    return Dict{LB,Vector{Int64}}(
        getDefaultLabelX(LB) => [255,0,0],
        getDefaultLabelY(LB) => [0,255,0],
        getDefaultLabelZ(LB) => [0,0,255]
    )
end


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using Printf
using Glob
using PyPlot
using Oceananigans

using ConvergenceTests
using ConvergenceTests.ForcedFlowFixedSlip: u

arch = CUDA.has_cuda() ? GPU() : CPU()

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

xy_filenames = glob("forced_fixed_slip_xy*", joinpath(@__DIR__, "data"))
xz_filenames = glob("forced_fixed_slip_xz*", joinpath(@__DIR__, "data"))

labels = ["(x, y)", "(x, z)"]

filenameses = [ xy_filenames,
                xz_filenames ]

errorses = [ ConvergenceTests.compute_errors((x, y, z, t) -> u(x, y, t), xy_filenames...),
             ConvergenceTests.compute_errors((x, y, z, t) -> u(x, z, t), xz_filenames...) ]

sizeses = [ ConvergenceTests.extract_sizes(xy_filenames...),
            ConvergenceTests.extract_sizes(xz_filenames...) ]

close("all")
fig, ax = subplots()

for j = 1:2
    filenames = filenameses[j]
    errors = errorses[j]
    sizes = sizeses[j]

    Nx = map(sz -> sz[1], sizes)

    L₁ = map(err -> err.L₁, errors)
    L∞ = map(err -> err.L∞, errors)

    @show size(L₁) size(Nx)

    ax.loglog(Nx, L₁, color=defaultcolors[j], alpha=0.6, mfc="None",
            linestyle="None", marker="o", label="\$L_1\$-norm, $(labels[j])")

    ax.loglog(Nx, L∞, color=defaultcolors[j], alpha=0.6, mfc="None",
            linestyle="None", marker="^", label="\$L_\\infty\$-norm, $(labels[j])")

    if j == 2
        L₁ = map(err -> err.L₁, errors)
        ii = sortperm(Nx)
        Nx = Nx[ii]
        L₁ = L₁[ii]
        ax.loglog(Nx, L₁[end] * (Nx[end] ./ Nx).^2, "k-", linewidth=1, alpha=0.6, label=L"\sim N_x^{-2}")
    end

    ax.set_xscale("log", base=2)
    ax.set_yscale("log", base=10)
end

filenames = filenameses[1]
errors = errorses[1]
sizes = sizeses[1]
Nx = map(sz -> sz[1], sizes)

legend()

xlabel(L"N_x")
ylabel("Norms of the absolute error, \$ | u_{\\mathrm{sim}} - u_{\\mathrm{exact}} | \$")

removespines("top", "right")
title("Convergence for forced fixed slip")

filename = "forced_fixed_slip_convergence_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480)

p = sortperm(Nx)
for (label, error) in zip(labels, errorses)
    L₁ = map(e -> e.L₁, error)
    L∞ = map(e -> e.L∞, error)
    name = "Forced flow fixed slip " * label
    test_rate_of_convergence(L₁[p], Nx[p], expected=-2.0, atol=0.05, name=name * " L₁")
    test_rate_of_convergence(L∞[p], Nx[p], expected=-2.0, atol=0.10, name=name * " L∞")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

__precompile__(true)

module Unicode

using Base.Unicode: normalize, graphemes, isassigned, textwidth, isvalid,
                    islower, isupper, isalpha, isdigit, isxdigit, isnumeric, isalnum,
                    iscntrl, ispunct, isspace, isprint, isgraph,
                    lowercase, uppercase, titlecase, lcfirst, ucfirst

export graphemes, textwidth, isvalid,
       islower, isupper, isalpha, isdigit, isxdigit, isnumeric, isalnum,
       iscntrl, ispunct, isspace, isprint, isgraph,
       lowercase, uppercase, titlecase, lcfirst, ucfirst

# BEGIN 0.7 deprecations

@deprecate isnumber(c::Char) Unicode.isnumeric(c)
@deprecate is_assigned_char(c::Char) Unicode.isassigned(c)
@deprecate normalize_string(s::AbstractString, nf::Symbol; kwargs...) Unicode.normalize(s, nf; kwargs...)
@deprecate normalize_string(s::AbstractString; kwargs...) Unicode.normalize(s; kwargs...)

# END 0.7 deprecations

end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

export CARGILLE, EYE

# all other glasses should follow the format below, new glasses must be added to OTHER_GLASSES and OTHER_GLASS_NAMES where the index in the array matches the numeric part of the GlassID

module CARGILLE
using ..AGFFileReader: Glass, GlassID, OTHER

const OG0608 = Glass(GlassID(OTHER, 1), -2, 1.4451400, 0.0043176, -1.80659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 800.0, -1.0, 0, -1.0, [(0.32, 0.03, 10.0), (0.365, 0.16, 100.0), (0.4047, 0.40, 100.0), (0.480, 0.71, 100.0), (0.4861, 0.72, 100.0), (0.5461, 0.80, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.95, 100.0), (0.6563, 0.96, 100.0), (0.6943, 0.99, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.74, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.16, 100.0)], 1.457518, -1.0, -1.0, 0, 57.18978, 0, 0.878, -1)
const OG0607 = Glass(GlassID(OTHER, 2), -2, 1.44503, 0.0044096, -2.85878e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.15, 10.0), (0.365, 0.12, 100.0), (0.4047, 0.42, 100.0), (0.480, 0.78, 100.0), (0.4861, 0.79, 100.0), (0.5461, 0.86, 100.0), (0.5893, 0.90, 100.0), (0.6328, 0.92, 100.0), (0.6439, 0.90, 100.0), (0.6563, 0.92, 100.0), (0.6943, 0.98, 100.0), (0.840, 0.99, 100.0), (0.10648, 0.61, 100.0), (0.1300, 0.39, 100.0), (0.1550, 0.11, 100.0)], 1.457587, -1.0, -1.0, 0, 57.19833, 0, 0.878, -1)
const OG081160 = Glass(GlassID(OTHER, 3), -2, 1.49614, 0.00692199, -8.07052e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.000885983052189022, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.014, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [(0.32, 0.04, 100.0), (0.365, 0.13, 100.0), (0.4047, 0.26, 100.0), (0.480, 0.48, 100.0), (0.4861, 0.49, 100.0), (0.5461, 0.60, 100.0), (0.5893, 0.68, 100.0), (0.6328, 0.71, 100.0), (0.6439, 0.73, 100.0), (0.6563, 0.74, 100.0), (0.6943, 0.76, 100.0), (0.840, 0.83, 100.0), (0.10648, 0.86, 100.0), (0.1300, 0.89, 100.0), (0.1550, 0.90, 100.0)], 1.515549, -1.0, -1.0, 0, 36.82493, 0, 1.11, -1)

end

module EYE
using ..AGFFileReader: Glass, GlassID, OTHER

"""    EYE.AQUEOUS
```
ID:                      OTHER:4
RI @ 587nm:              1.336981
Abbe Number:             52.658991
ΔPgF:                    0.0
TCE (÷1e-6):             0.0
Density:                 1.0g/m³
Valid wavelengths:       0.38μm to 0.78μm
Reference Temp:          20.0°C
```
"""
const AQUEOUS = Glass(GlassID(OTHER, 4), 5, 1.32107278, 0.00847113739, 0.000231825063, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.38, 0.78, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.32, 0.16, 25.0), (0.334, 0.42, 25.0), (0.35, 0.7, 25.0), (0.365, 0.85, 25.0), (0.37, 0.88, 25.0), (0.38, 0.92, 25.0), (0.39, 0.95, 25.0), (0.4, 0.963, 25.0), (0.42, 0.977, 25.0), (0.46, 0.988, 25.0), (0.5, 0.993, 25.0), (0.66, 0.996, 25.0), (1.06, 0.996, 25.0), (1.529, 0.975, 25.0), (1.9701, 0.93, 25.0), (2.325, 0.66, 25.0)], 1.336981, -1.0, -1.0, 0, 52.658991, 0, 1.0, 0)
"""    EYE.LENS
```
ID:                      OTHER:5
RI @ 587nm:              1.419976
Abbe Number:             51.226142
ΔPgF:                    0.0
TCE (÷1e-6):             0.0
Density:                 1.0g/m³
Valid wavelengths:       0.38μm to 0.78μm
Reference Temp:          20.0°C
```
"""
const LENS = Glass(GlassID(OTHER, 5), 5, 1.4014679, 0.00938901135, 0.000393175776, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.38, 0.78, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.32, 0.16, 25.0), (0.334, 0.42, 25.0), (0.35, 0.7, 25.0), (0.365, 0.85, 25.0), (0.37, 0.88, 25.0), (0.38, 0.92, 25.0), (0.39, 0.95, 25.0), (0.4, 0.963, 25.0), (0.42, 0.977, 25.0), (0.46, 0.988, 25.0), (0.5, 0.993, 25.0), (0.66, 0.996, 25.0), (1.06, 0.996, 25.0), (1.529, 0.975, 25.0), (1.9701, 0.93, 25.0), (2.325, 0.66, 25.0)], 1.419976, -1.0, -1.0, 0, 51.226142, 0, 1.0, 0)
"""    EYE.CORNEA
```
ID:                      OTHER:6
RI @ 587nm:              1.376981
Abbe Number:             56.279936
ΔPgF:                    0.0
TCE (÷1e-6):             0.0
Density:                 1.0g/m³
Valid wavelengths:       0.38μm to 0.78μm
Reference Temp:          20.0°C
```
"""
const CORNEA = Glass(GlassID(OTHER, 6), 5, 1.36313817, 0.00667127181, 0.000386916734, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.38, 0.78, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.32, 0.16, 25.0), (0.334, 0.42, 25.0), (0.35, 0.7, 25.0), (0.365, 0.85, 25.0), (0.37, 0.88, 25.0), (0.38, 0.92, 25.0), (0.39, 0.95, 25.0), (0.4, 0.963, 25.0), (0.42, 0.977, 25.0), (0.46, 0.988, 25.0), (0.5, 0.993, 25.0), (0.66, 0.996, 25.0), (1.06, 0.996, 25.0), (1.529, 0.975, 25.0), (1.9701, 0.93, 25.0), (2.325, 0.66, 25.0)], 1.376981, -1.0, -1.0, 0, 56.279936, 0, 1.0, 0)
"""    EYE.VITREOUS
```
ID:                      OTHER:7
RI @ 587nm:              1.335982
Abbe Number:             53.342173
ΔPgF:                    0.0
TCE (÷1e-6):             0.0
Density:                 1.0g/m³
Valid wavelengths:       0.38μm to 0.78μm
Reference Temp:          20.0°C
```
"""
const VITREOUS = Glass(GlassID(OTHER, 7), 5, 1.32238376, 0.00672767909, 0.000333967702, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.38, 0.78, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20.0, 0.0, -1.0, -1.0, 0.0, -1.0, 0, -1.0, [(0.32, 0.16, 25.0), (0.334, 0.42, 25.0), (0.35, 0.7, 25.0), (0.365, 0.85, 25.0), (0.37, 0.88, 25.0), (0.38, 0.92, 25.0), (0.39, 0.95, 25.0), (0.4, 0.963, 25.0), (0.42, 0.977, 25.0), (0.46, 0.988, 25.0), (0.5, 0.993, 25.0), (0.66, 0.996, 25.0), (1.06, 0.996, 25.0), (1.529, 0.975, 25.0), (1.9701, 0.93, 25.0), (2.325, 0.66, 25.0)], 1.335982, -1.0, -1.0, 0, 53.342173, 0, 1.0, 0)

end

const OTHER_GLASSES = [CARGILLE.OG0607, CARGILLE.OG0608, CARGILLE.OG081160, EYE.AQUEOUS, EYE.LENS, EYE.CORNEA, EYE.VITREOUS]
const OTHER_GLASS_NAMES = ["CARGILLE.OG0607", "CARGILLE.OG0608", "CARGILLE.OG081160", "EYE.AQUEOUS", "EYE.LENS", "EYE.CORNEA", "EYE.VITREOUS"]

max_stars_repo_name: rambunctiousapple/GlassCat

max_issues_repo_name: rambunctiousapple/GlassCat

max_forks_repo_name: rambunctiousapple/GlassCat


text:
using Literate

@show sources_dir = joinpath(@__DIR__, "..", "deps")
notebooks_dir = joinpath(@__DIR__, "..", "notebooks")

if !isdir(notebooks_dir)
    mkdir(notebooks_dir)
end

for t in readdir(sources_dir)
        @show t
    if occursin(r"tutorial.*.jl", t)
        # Literate.markdown(t, "."; documenter=false);
        Literate.notebook(joinpath(sources_dir, t), notebooks_dir; documenter=false, execute=false)
        cp(t,  joinpath(notebooks_dir, t), force=true)
    end
end

for insf in ["trunc_cyl_shell_0.nas", "block-w-hole.inp", "samplet4.jl", "t3_i1.png", "t4_i1.png", "t6_i1.png"]
    cp(insf,  joinpath(notebooks_dir, insf), force=true)
end


text:
using DroneRL
using Test

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


text:
type Spline
  d
  U #knots
  C #coefficients
  function Spline(d,U,C)
    p = sortperm(U)
    U,C = U[p], C[p]
    U = [ fill(U[1],d); U; fill(U[end],d) ]
    C = [ fill(C[1],d); C; fill(C[end],d) ]
    return new(d,U,C)
  end
end

test = Spline(1,linspace(0,1,10),rand(10))

function fastfind(A,x)
  i,j = 1,length(A)
  A[i] > x && return 0
  A[j] < x && return j
  while abs(i-j) > 1
    k = (i+j) ÷ 2
    if A[k] < x; i = k
    else j = k end
  end
  return i
end

function splinematrix(S::Spline,i,k,x)
  d,U = S.d,S.U
  M = zeros(k,k+1)
  for a in 1:k
    M[a,a] = (U[i+a] - x) / (U[i+a] - U[i+a-k])
    M[a,a+1] = (x - U[i+a-k]) / (U[i+a] - U[i+a-k])
  end
  return M
end

splinematrix(test,2,1,.1)

@noinline function call(S::Spline, x)
  core(S::Spline, x)
end

@noinline function core(S::Spline, x)
  d,U,C = S.d,S.U,S.C
  i = fastfind(U,x)
  # i = clamp(i,d+1,length(U)-d-1)
  i = max(i,d+1)
  println("x")
  c = C[i-d:i]
  for k in d:-1:1
    c = splinematrix(S,i,k,x) * c
  end
  return c[1]
end

using PyPlot
plot(linspace(0,1,length(test.C)), test.C)
plot(linspace(0,1,100), map(x->core(test, x), linspace(test.U[1],test.U[end],100)))

splinematrix(test, 10, 1, 1.0)

x = sort(rand(10))

fastfind(x,.0)


text:
using Documenter, ControlBenchmarks

DocMeta.setdocmeta!(ControlBenchmarks, :DocTestSetup, :(using ControlBenchmarks); recursive=true)

makedocs(;
    modules=[ControlBenchmarks],
    format=Documenter.HTML(prettyurls = false),
    repo="https://github.com/imciner2/ControlBenchmarks.jl/blob/{commit}{path}#L{line}",
    sitename="ControlBenchmarks.jl",
    authors="Ian McInerney, Imperial College London",
    pages = Any[
        "Home" => "index.md",
        "Usage" => "usage.md",
        "Problems" => [
            "ifac.md",
        ],
        "Extending" => "extending.md",
        "Internals" => "internals.md"
    ],
    assets=String[],
)

#deploydocs(;
#    repo="github.com/imciner2/ControlBenchmarks.jl",
#)


text: avg(x) = sum(x)/length(x)


text:
using ModelingToolkit
using OrdinaryDiffEq
using Plots

function duffing()
    @parameters t m c k μ F₀ ω
    @variables x₁(t) x₂(t)
    D = Differential(t)

    # mx″ + cx′ + kx + μx³ = F₀cos(ωt)
    eqns = [
        D(x₁) ~ x₂,
        D(x₂) ~ (F₀*cos(ω*t) - c*x₂ - k*x₁ - μ*x₁^3)/m
    ]
    sys = ODESystem(eqns)

    # Compile the ODEFunction with the states and parameters in the specified order
    return ODEFunction(sys, [x₁, x₂], [m, c, k, μ, F₀, ω])
end

function integrate(ode, u, p, t)
    # Explicitly convert times to Float64 in case integer values are used
    prob = ODEProblem(ode, u, Float64.(t), p)
    return solve(prob, Tsit5())
end

function example()
    ode = duffing()
    m = 1.0
    c = 0.01
    k = 2.0
    μ = 1
    F₀ = 0.05
    ω = 1.0
    sol = integrate(ode, [0, 0], [m, c, k, μ, F₀, ω], (0, 100))
    plot(sol.t, sol[1,:])
end

# Code to benchmark the timings

# using BenchmarkTools
# ode = duffing()
# m = 1.0
# c = 0.01
# k = 2.0
# μ = 1
# F₀ = 0.05
# ω = 1.0
# @benchmark integrate($ode, $([0, 0]), $([m, c, k, μ, F₀, ω]), (0, 100))

# Output on my laptop

# BenchmarkTools.Trial:
#   memory estimate:  168.48 KiB
#   allocs estimate:  1606
#   --------------
#   minimum time:     101.200 μs (0.00% GC)
#   median time:      107.600 μs (0.00% GC)
#   mean time:        141.544 μs (9.58% GC)
#   maximum time:     8.184 ms (97.48% GC)
#   --------------
#   samples:          10000
#   evals/sample:     1

max_stars_repo_name: dawbarton/julia-de-examples

max_issues_repo_name: dawbarton/julia-de-examples

max_forks_repo_name: dawbarton/julia-de-examples


text:
mutable struct Continuous_MountainCarEnv
    min_action::Float32
    max_action::Float32
    min_position::Float32
    max_position::Float32
    max_speed::Float32
    goal_position::Float32
    power::Float32

    state
    action_space::Box
    observation_space::Box
end

function Continuous_MountainCarEnv()
    min_action    =  -1f0
    max_action    =   1f0
    min_position  = -12f-1
    max_position  =   6f-1
    max_speed     =   7f-2
    goal_position =  45f-2 # was 0.5 in gym, 0.45 in Arnaud de Broissia's version
    power         =  15f-4
    action_space = Box(min_action, max_action, (1,), Float32)
    observation_space = Box([min_position, -max_speed], [max_position, max_speed], Float32)
    Continuous_MountainCarEnv(min_action, max_action, min_position, max_position,
                              max_speed, goal_position, power, nothing, action_space, observation_space)
end

function step!(env::Continuous_MountainCarEnv, action)
    @assert action ∈ env.action_space "Invalid action selected ($(action))"
    position, velocity = env.state[1:1], env.state[2:2]
    force = clamp.(action, env.min_action, env.max_action)

    v          = velocity .+ force * env.power .- 25f-4cos.(3f0position)
    velocity_  = clamp.(v, -env.max_speed, env.max_speed)
    x          = position .+ velocity_
    position_  = clamp.(x, env.min_position, env.max_position)
    if all(position_ .== env.min_position) && all(velocity_ .< 0)
        velocity_ = 0f0velocity_
    end

    done = all(position_ .≥ env.goal_position)

    r = [0f0]
    if done
        r = [1f2]
    end
    reward = r .- 1f-1action .^ 2

    env.state = vcat(position_, velocity_)
    return env.state, reward, done, Dict()
end

function reset!(env::Continuous_MountainCarEnv)
    env.state = param([2f-1rand(Float32) - 6f-1, 0f0])
end

show(io::IO, env::Continuous_MountainCarEnv) = print(io, "Continuous-MountainCarEnv")


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.6.3
#     language: julia
#     name: julia-1.6
# ---

# %%
using AbstractAlgebra
AbstractAlgebra.PrettyPrinting.set_html_as_latex(true)

R, (a, b, c, l, m, n) = ZZ["a", "b", "c", "l", "m", "n"]
Rx, x = R["x"]

# %%
f = a*x^2 + b*x + c

# %%
g = l*x^2 + m*x + n

# %%
resultant(f, g)

# %%

max_stars_repo_path: 0023/Sylvester's resultant.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0023/Sylvester's resultant.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0023/Sylvester's resultant.jl

max_forks_repo_name: genkuroki/public


text:
using ArgParse

include("functions_align.jl")
include("functions_pairsnp.jl")

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--infile1"
            help = "alignment 1"
            required = true
        "--infile2"
            help = "alignment 2"
        "--outfile", "-o"
            help = "the name of the file to write"
            default = "pairsnp.csv"
        end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()

    nucleotide_array_1, fasta_IDs_1 = populate_byte_array_get_names(parsed_args["infile1"])

    consensus_seq = consensus(nucleotide_array_1)

    if parsed_args["infile2"] != nothing
        nucleotide_array_2, fasta_IDs_2 = populate_byte_array_get_names(parsed_args["infile2"])
        snpmat = pairsnp_rectangle(nucleotide_array_1, nucleotide_array_2, consensus_seq)
        write_pairsnp_rectangle(parsed_args["outfile"], snpmat, fasta_IDs_1, fasta_IDs_2)
    else
        snpmat = pairsnp_square(nucleotide_array_1, consensus_seq)
        write_pairsnp_square(parsed_args["outfile"], snpmat, fasta_IDs_1)
    end
end

main()

max_stars_repo_name: snake-flu/julialign

max_issues_repo_name: snake-flu/julialign

max_forks_repo_name: snake-flu/julialign


text:
using Documenter, Hangul

makedocs(
    modules = [Hangul],
    clean = false,
    format = :html,
    sitename = "Hangul.jl",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
        "Letter" => "Letter.md",
        "YetJamo" => "YetJamo.md",
    ],
    html_prettyurls = !("local" in ARGS),
)


text:
using NCDatasets
if VERSION >= v"0.7.0-beta.0"
    using Test
else
    using Base.Test
end

fname = tempname()

ds = Dataset(fname,"c")
# Dimensions

ds.dim["time"] = 3

# Declare variables

ncLAT = defVar(ds,"LAT", Float64, ("time",))
ncLAT.attrib["standard_name"] = "latitude"
ncLAT.attrib["long_name"] = "Latitude coordiante"
ncLAT.attrib["units"] = "degree_north"
ncLAT.attrib["ancillary_variables"] = "QC_LAT"
ncLAT.attrib["axis"] = "Y"

ncQC_LAT = defVar(ds,"QC_LAT", Int8, ("time",))
ncQC_LAT.attrib["standard_name"] = "latitude status_flag"
ncQC_LAT.attrib["long_name"] = "Quality flag for latitude"
ncQC_LAT.attrib["_FillValue"] = Int8(10)
ncQC_LAT.attrib["flag_values"] = Int8[0, 1, 2, 3, 4, 6, 9]
ncQC_LAT.attrib["flag_meanings"] = "no_qc_performed good_data probably_good_data probably_bad_data bad_data spike missing_value"

ncDEPTH = defVar(ds,"DEPTH", Float64, ())
ncDEPTH.attrib["standard_name"] = "depth"
ncDEPTH.attrib["long_name"] = "Depth coordinate"
ncDEPTH.attrib["units"] = "m"
ncDEPTH.attrib["positive"] = "down"
ncDEPTH.attrib["axis"] = "Z"
ncDEPTH.attrib["reference_datum"] = "geographical coordinates, WGS84 projection"

ncLAT[:] = [1.,2.,3.]
ncQC_LAT[:] = [1,1,4]

close(ds)

ds = Dataset(fname,"r")

@test name(NCDatasets.ancillaryvariables(ds["LAT"],"status_flag")) == "QC_LAT"
@test isequal(NCDatasets.filter(ds["LAT"],:,accepted_status_flags = ["good_data","probably_good_data"]),
    [1.,2.,missing])


text:
# # ACOPF with [PowerSimulations.jl](https://github.com/NREL-SIIP/PowerSimulations.jl) using [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl)

# **Originally Contributed by**: Clayton Barrows

# ## Introduction

# PowerSimulations.jl supports non-linear AC optimal power flow through a deep integration
# with [PowerModels.jl](https://github.com/lanl-ansi/PowerModels.jl). This example shows a
# single multi-period optimization of economic dispatch with a full representation of
# AC optimal power flow.

# ## Dependencies
# We can use the a TAMU synthetic ERCOT dataset that is included in the PowerSystemsTestData.
using SIIPExamples
using PowerSystems
using PowerSimulations
using Dates

pkgpath = pkgdir(SIIPExamples)
PowerSystems.download(PowerSystems.TestData; branch = "master") # *note* add `force=true` to get a fresh copy
base_dir = pkgdir(PowerSystems);

# The TAMU data format relies on a folder containing `.m` or `.raw` files and `.csv`
# files for the time series data. We have provided a parser for the TAMU data format with
# the `TamuSystem()` function.

TAMU_DIR = joinpath(base_dir, "data", "ACTIVSg2000");
sys = TamuSystem(TAMU_DIR)
transform_single_time_series!(sys, 2, Hour(1))

# Since we'll be doing non-linear optimization, we need a solver that supports non-linear
# problems. Ipopt is quite good.
using Ipopt
solver = optimizer_with_attributes(Ipopt.Optimizer)

# In the [OperationsProblem example](../../notebook/3_PowerSimulations_examples/1_operations_problems.ipynb)
# we defined a unit-commitment problem with a copper plate representation of the network.
# Here, we want do define an economic dispatch (linear generation decisions) with an ACOPF
# network representation.
# So, starting with the network, we can select from _almost_ any of the endpoints on this
# tree:
#nb TypeTree(PSI.PM.AbstractPowerModel,  init_expand = 10, scopesep="\n")

# For now, let's just choose a standard ACOPF formulation.
devices = Dict(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
        :QLoads => DeviceModel(FixedAdmittance, StaticPowerLoad)
    )
ed_template = template_economic_dispatch(network = ACPPowerModel, devices = devices)

# Now we can build a 4-hour economic dispatch / ACOPF problem with the TAMU data.
problem = OperationsProblem(
    EconomicDispatchProblem,
    ed_template,
    sys,
    horizon = 1,
    optimizer = solver,
    balance_slack_variables = true,
)

# And solve it ...
solve!(problem)


text:
# Run this file to add packages to Project.toml

using Pkg

"""
    addpackage()

Install required packages
"""
function addpackage()
    Pkg.add("BenchmarkTools") # adds DataFrames and DelimitedFiles
    Pkg.add("CSV")
    Pkg.add("XLSX")
    Pkg.add("JLD")
    Pkg.add("NPZ")
    Pkg.add("RData")
    Pkg.add("MAT")
    Pkg.add("Images")
    Pkg.add("ImageMagick")
    Pkg.add("DataFrames")
    # Statistics
    Pkg.add("Statistics")
    Pkg.add("StatsBase")
    Pkg.add("RDatasets")
    Pkg.add("Plots")
    Pkg.add("Plotly")
    Pkg.add("StatsPlots")
    Pkg.add("KernelDensity")
    Pkg.add("Distributions")
    Pkg.add("HypothesisTests")
    Pkg.add("PyCall")
    Pkg.add("MLBase")
    # Dimensionality Reduction
    Pkg.add("UMAP")
    Pkg.add("Makie")
    Pkg.add("VegaDatasets")
    Pkg.add("MultivariateStats")
    Pkg.add("ScikitLearn")
    Pkg.add("Distances")
end
addpackage() # Do this only once


text:
@testset "HashIndex" begin
   a = [3, 1, 2, 1]
   b = accelerate(a, HashIndex)

   @test 2 ∈ b
   @test 0 ∉ b

   @test count(isequal(1), b) == 2
   @test count(isequal(2), b) == 1
   @test count(isequal(4), b) == 0

   @test findall(isequal(1), b) == [2, 4]
   @test findall(isequal(4), b) == []

   @test filter(isequal(1), b) == [1, 1]
   @test filter(isequal(4), b) == []

   @test issetequal(unique(b), [1,2,3])

   @test group(identity, b) == group(identity, a)
   @test groupfind(identity, b) == groupfind(identity, a)
   @test groupreduce(identity, +, b) == groupreduce(identity, +, a)

   @test issetequal(innerjoin(identity, identity, tuple, isequal, b, [0, 1, 2]),
                    innerjoin(identity, identity, tuple, isequal, a, [0, 1, 2]))
end


text:
using ControlSystems
using LinearAlgebra

"""
    ntf, stf = calculateTF(ABCD, k=1)

Calculate the NTF and STF of a delta-sigma modulator whose loop filter
is described by the `ABCD` matrix, assuming a quantizer gain of `k`.

The `ntf` and `stf` are TransferFunction objects.
"""
function calculateTF(ABCD, k=[1.0])
    nq = length(k)
    n = size(ABCD, 1) - nq
    nu = size(ABCD, 2) - size(ABCD, 1)
    A, B, C, D = partitionABCD(ABCD, nu+nq)
    B1 = B[:,1:nu]
    B2 = B[:,nu+1:end]
    D1 = D[:,1:nu]
    if any(x -> x != 0, D[:,nu+1:end])
        error("D must be zero")
    end
    K = diagm(k)

    # Find the noise transfer function by forming the closed-loop
    # system (sys_cl) in state-space form.
    Acl = A + B2*K*C
    Bcl = [(B1 + B2*K*D1) B2]
    Ccl = K*C
    Dcl = [K*D1 diagm(ones(nq))]

    sys_cl = ss(Acl, Bcl, Ccl, Dcl, 1)
    tfs = zpk(sys_cl)
    tol = min(1e-3, max(1e-6, eps(Float64)^(1/(size(ABCD,1)))))
    mtfs = minreal(tfs, tol)
    stf = TransferFunction(mtfs.matrix[:,1:nu], mtfs.Ts)
    ntf = TransferFunction(mtfs.matrix[:,nu+1:nu+nq], mtfs.Ts)

    return ntf, stf
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma