text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

using FEMCoupling: get_C
using Base.Test

@testset "Plain strain kinematic Coupling" begin
# plane strain problem
# Square shaped plane strain element(8x8) + one single node(2x2)
K=zeros(10,10)
K[1:8,1:8]=
  [500 225 -350 75 -250 -225 100 -75;
   225 500 -75 100 -225 -250 75 -350;
  -350 -75 500 -225  100 75 -250 225;
   75 100 -225  500 -75 -350 225 -250;
  -250 -225 100 -75   500 225 -350 75;
  -225 -250 75  -350  225 500 -75 100;
   100 75 -250 225 -350 -75 500 -225;
  -75 -350 225 -250 75 100 -225 500]

# Removing fixed DOFs (1 and 2) from nodes 1 and 4
K=K[1:end .!=1, 1:end .!=1]
K=K[1:end .!=1, 1:end .!=1]
K=K[1:end .!=5, 1:end .!=5]
K=K[1:end .!=5, 1:end .!=5]

# Force vector
f=zeros(10,1)
f[9,1]=1200

# Removing fixed DOFs
f=f[1:end .!= 1]
f=f[1:end .!= 1]
f=f[1:end .!= 5]
f=f[1:end .!= 5]

############### Calculations by hand

# Making the C matrix by hand
A=zeros(4,10)
A[1,3]-=1;  A[1,9]=1; A[2,4]=-1; A[2,10]=1;
A[3,5]-=1;  A[3,9]=1; A[4,6]=-1; A[4,10]=1

# Eliminating fixed dofs (1,2) from nodes 1 and 4
A=A[1:end , 1:end .!=1 ]
A=A[1:end , 1:end .!=1 ]
A=A[1:end , 1:end .!=5 ]
A=A[1:end , 1:end .!=5 ]

# Renaming variables to match with get_C.jl
C_expected=A
g_expected=zeros(4,1)
D_expected=zeros(4,4)

# Assembly for solving
K_expected= [K            C_expected';
             C_expected   D_expected]

f_expected= [f;
             g_expected]

u_expected = K_expected\f_expected

############### Calculating C,D and g with get_C.jl
#                  get_C(refnode,slaves,dofs,ndofs,K_size)
K_size=size(K,1)
C,D,g= FEMCoupling.get_C(5,[2,3],[1,2],2,10)

# Removing fixed DOFs
C=C[1:end , 1:end .!=1 ]
C=C[1:end , 1:end .!=1 ]
C=C[1:end , 1:end .!=5 ]
C=C[1:end , 1:end .!=5 ]

KK=[K C';
    C D]
ff=[f;
    g]
u = lufact(KK) \ full(ff)

@test isapprox(C,C_expected,rtol=0.0001)
@test isapprox(D,D_expected,rtol=0.0001)
@test isapprox(g,g_expected,rtol=0.0001)
@test isapprox(ff,f_expected,rtol=0.0001)
@test isapprox(KK,K_expected,rtol=0.0001)
@test isapprox(u,u_expected,rtol=0.0001)
# If the last test passes, all other tests will pass too.
# Other tests are made to help tracing why the last test doesn't pass.
end


text:
#=
Interpolation Methods

=#

## ------ Akima Interpolation  ---------

# struct used internally
struct Akima{TF, T1<:AbstractVector{TF}, T2<:AbstractVector{TF}}

    xdata::T1
    ydata::T2
    p0::Vector{TF}
    p1::Vector{TF}
    p2::Vector{TF}
    p3::Vector{TF}

end

"""
    akima_setup(xdata, ydata, delta_x=0.1)

Setup akima spline at node points: xdata, ydata.  This is a 1D spline that avoids
overshooting issues common with many other polynomial splines resulting in a more
natural curve.  It also only depends on local points (i-2...i+2) allow for more efficient
computation.  delta_x is the half width of a smoothing interval used for the absolute
value function.  Set delta_x=0 to recover the original akima spline.  The smoothing
is only useful if you want to differentiate xdata and ydata.  In many case the nodal points
are fixed so this is not needed.  Returns an akima spline object (Akima struct).
This function, in connection with akima_interp separates the construction from evaluation.
This is useful if you want to evaluate the same mesh at multiple different conditions.
"""
function akima_setup(xdata, ydata, delta_x=0.0)

    # setup
    eps = 1e-30
    n = length(xdata)

    # compute segment slopes
    m = OffsetVector(zeros(n+3), -1:n+1)
    for i = 1:n-1
        m[i] = (ydata[i+1] - ydata[i]) / (xdata[i+1] - xdata[i])
    end

    # estimation for end points
    m[0] = 2.0*m[1] - m[2]
    m[-1] = 2.0*m[0] - m[1]
    m[n] = 2.0*m[n-1] - m[n-2]
    m[n+1] = 2.0*m[n] - m[n-1]

    # slope at points
    t = zeros(n)
    for i = 1:n
        m1 = m[i-2]
        m2 = m[i-1]
        m3 = m[i]
        m4 = m[i+1]
        w1 = abs_smooth(m4 - m3, delta_x)
        w2 = abs_smooth(m2 - m1, delta_x)
        if (w1 < eps && w2 < eps)
            t[i] = 0.5*(m2 + m3)  # special case to avoid divide by zero
        else
            t[i] = (w1*m2 + w2*m3) / (w1 + w2)
        end
    end

    # polynomial cofficients
    p0 = zeros(n-1)
    p1 = zeros(n-1)
    p2 = zeros(n-1)
    p3 = zeros(n-1)
    for i = 1:n-1
        dx = xdata[i+1] - xdata[i]
        t1 = t[i]
        t2 = t[i+1]
        p0[i] = ydata[i]
        p1[i] = t1
        p2[i] = (3.0*m[i] - 2.0*t1 - t2)/dx
        p3[i] = (t1 + t2 - 2.0*m[i])/dx^2
    end

    return Akima(xdata, ydata, p0, p1, p2, p3)
end

"""
    akima_interp(x, spline)

Evaluate the spline setup by akima_setup.  x are the locations to evaluate the spline at
and spline is the object returned by akima_setup.
"""
function akima_interp(x, spline)

    # interpolate at each point
    n = length(x)
    y = zeros(n)

    for i = 1:n

        j = findlast(x[i] .>= spline.xdata[1:end-1])
        if isnothing(j)
            j = 1
        end

        # evaluate polynomial
        dx = x[i] - spline.xdata[j]
        y[i] = spline.p0[j] + spline.p1[j]*dx + spline.p2[j]*dx^2 + spline.p3[j]*dx^3

    end

    return y
end

"""
    akima(xdata, ydata, xpt)

A convenience method to perform both the setup and evaluation of the akima spline in one go.
xpt may be an array.
"""
function akima(xdata, ydata, xpt)
    spline = akima_setup(xdata, ydata)
    return akima_interp(xpt, spline)
end

"""
   interp2d(interp1d, xdata, ydata, fdata, xpt, ypt)

2D interpolation using recursive 1D interpolation.  This approach is likely less efficient than a more direct
2D interpolation method, especially one you can create separate creation from evaluation,
but it is generalizable to any spline approach and any dimension.

**Arguments**
- `interp1d`: any spline function of form: ypt = interp1d(xdata, ydata, xpt) where data are the known
    data(node) points and pt are the points where you want to evaluate the spline at.
- `xdata::Vector{Float}`, `ydata::Vector{Float}`: Define the 2D grid
- `fdata::Matrix{Float}`: where fdata[i, j] is the function value at xdata[i], ydata[j]
- `xpt::Vector{Float}`, `ypt::Vector{Float}`: the locations where you want to evaluate the spline

**Returns**
- `fhat::Matrix{Float}`: where fhat[i, j] is the estimate function value at xpt[i], ypt[j]
"""
function interp2d(interp1d, xdata, ydata, fdata, xpt, ypt)

    ny = length(ydata)
    nxpt = length(xpt)
    nypt = length(ypt)
    yinterp = zeros(ny, nxpt)
    output = zeros(nxpt, nypt)

    for i = 1:ny
        yinterp[i, :] = interp1d(xdata, fdata[:, i], xpt)
    end
    for i = 1:nxpt
        output[i, :] = interp1d(ydata, yinterp[:, i], ypt)
    end

    return output
end

"""
    interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Same as interp2d, except in three dimension.
"""
function interp3d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

    nz = length(zdata)
    nxpt = length(xpt)
    nypt = length(ypt)
    nzpt = length(zpt)
    zinterp = zeros(nz, nxpt, nypt)
    output = zeros(nxpt, nypt, nzpt)

    for i = 1:nz
        zinterp[i, :, :] = interp2d(interp1d, xdata, ydata, fdata[:, :, i], xpt, ypt)
    end
    for j = 1:nypt
        for i = 1:nxpt
            output[i, j, :] = interp1d(zdata, zinterp[:, i, j], zpt)
        end
    end

    return output
end

"""
    interp4d(interp1d, xdata, ydata, zdata, fdata, xpt, ypt, zpt)

Same as interp3d, except in four dimensions.
"""
function interp4d(interp1d, xdata, ydata, zdata, tdata, fdata, xpt, ypt, zpt, tpt)

    nt = length(tdata)
    nxpt = length(xpt)
    nypt = length(ypt)
    nzpt = length(zpt)
    ntpt = length(tpt)
    tinterp = zeros(nt, nxpt, nypt, nzpt)
    output = zeros(nxpt, nypt, nzpt, ntpt)

    for i = 1:nt
        tinterp[i, :, :, :] = interp3d(interp1d, xdata, ydata, zdata, fdata[:, :, :, i], xpt, ypt, zpt)
    end
    for k = 1:nzpt
        for j = 1:nypt
            for i = 1:nxpt
                output[i, j, k, :] = interp1d(tdata, tinterp[:, i, j, k], tpt)
            end
        end
    end

    return output
end

ext: jl


text:
using Distributions

𝕀(b) = b ? 1 : 0 # indicator function

function sierra_mixture(; 𝛍=[0,0], σ=3, 𝐬=[+σ,-σ], δ=2, η=6, decay=true)
    𝚺 = [σ 0.0; 0.0 σ]
    decay = 𝕀(decay)
    origin = [0.0, 0.0]
    𝐒 = MvNormal[MvNormal(𝛍, 𝚺/(σ*η))]

    for g in [[+δ, +δ], [+δ, -δ], [-δ, +δ], [-δ, -δ]]
        for (i,𝐩) in enumerate([origin, [+1, +1], [+2, 0], [+3, +1], [0, +2], [+1, +3]])
            for s in 𝐬
                push!(𝐒, MvNormal(g + s*𝐩 + 𝛍, i^decay/η * 𝚺))
            end
        end
    end

    return MixtureModel(𝐒)
end

global MMCACHE = Dict() # hyperparameters => MixtureModel

function sierra(x; minimize=true, kwargs...)
    global MMCACHE

    hpk = kwargs
    if haskey(MMCACHE, hpk)
        𝐒 = MMCACHE[hpk]
    else
        𝐒 = sierra_mixture(; kwargs...)
    end

    # mixture model with equal weights
    sgn = minimize ? -1 : +1
    return sgn*pdf(𝐒, x)
end


text:
using ClassicalOrthogonalPolynomials, FastGaussQuadrature

P = Legendre()
x = axes(P,1)
x .* P


text:
using BioSymbols
using BioSequences
using IterTools
using Base.Test

using Rifraf

include("test_utils.jl")

import Rifraf.sample_from_template,
       Rifraf.random_seq,
       Rifraf.ScoredProposal,
       Rifraf.Proposal,
       Rifraf.Substitution,
       Rifraf.Insertion,
       Rifraf.Deletion,
       Rifraf.apply_proposals,
       Rifraf.rbase,
       Rifraf.smart_forward_moves!,
       Rifraf.backward,
       Rifraf.phred_to_log_p,
       Rifraf.initial_state,
       Rifraf.estimate_probs,
       Rifraf.has_single_indels,
       Rifraf.score_proposal,
       Rifraf.single_indel_proposals,
       Rifraf.resample!,
       Rifraf.realign_rescore!,
       Rifraf.get_candidates,
       Rifraf.alignment_proposals,
       Rifraf.align_moves,
       Rifraf.moves_to_aligned_seqs,
       Rifraf.alignment_error_probs,
       Rifraf.Prob,
       Rifraf.LogProb

srand(1234)

@testset "scoring proposals" begin
    const errors = ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0)
    const scores = Scores(errors)

    function test_random_proposal(proposal, template_len)
        template = random_seq(template_len)

        template_error_p = Prob[0.1 for _=1:length(template)]
        codon_moves = rand([true, false])
        if codon_moves
            local_errors = ErrorModel(2.0, 0.1, 0.1, 3.0, 3.0)
        else
            local_errors = ErrorModel(2.0, 4.0, 4.0, 0.0, 0.0)
        end
        local_scores = Scores(local_errors)
        phred_scale = 3.0
        log_actual_error_std = 0.5
        log_reported_error_std = 0.5

        (seq, actual, phreds, sbools,
         tbools) = sample_from_template(template,
                                        template_error_p,
                                        errors,
                                        phred_scale,
                                        log_actual_error_std,
                                        log_reported_error_std)
        log_p = LogProb[float(q) / (-10.0) for q in phreds]
        bandwidth = max(5 * abs(length(template) - length(seq)), 30)
        pseq = RifrafSequence(seq, log_p, bandwidth, local_scores)

        new_template = apply_proposals(template, Proposal[proposal])
        new_template_error_p = Prob[0.1 for _=1:length(new_template)]

        Anew = Rifraf.forward(new_template, pseq)
        Bnew = backward(new_template, pseq)
        check_all_cols(Anew, Bnew, codon_moves)

        A = Rifraf.forward(template, pseq)
        B = backward(template, pseq)
        check_all_cols(A, B, codon_moves)

        newcols = zeros(size(A)[1], Rifraf.CODON_LENGTH + 1)
        score = score_proposal(proposal, A, B, template, pseq, newcols)

        @test score ≈ Anew[end, end]
    end

    @testset "random substitutions" begin
        for i = 1:1000
            template_len = rand(30:50)
            pos = rand(1:template_len)
            proposal = Substitution(pos, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random insertions" begin
        for i = 1:1000
            template_len = rand(30:50)
            pos = rand(0:template_len)
            proposal = Insertion(pos, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random deletions" begin
        for i = 1:1000
            template_len = rand(30:50)
            pos = rand(1:template_len)
            proposal = Deletion(pos)
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random deletions at beginning" begin
        for i = 1:10
            template_len = rand(30:50)
            proposal = Deletion(1)
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random deletions at end" begin
        for i = 1:10
            template_len = rand(30:50)
            proposal = Deletion(template_len)
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random substitutions at beginning" begin
        # subs at the beginning
        for i = 1:10
            template_len = rand(30:50)
            proposal = Substitution(1, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random substitutions at end" begin
        # subs at the end
        for i = 1:10
            template_len = rand(30:50)
            proposal = Substitution(template_len, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random insertions at beginning" begin
        for i = 1:10
            template_len = rand(30:50)
            proposal = Insertion(0, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
    @testset "random insertions at end" begin
        for i = 1:10
            template_len = rand(30:50)
            proposal = Insertion(template_len, rbase())
            test_random_proposal(proposal, template_len)
        end
    end
end

@testset "no single indels" begin
    reference = DNASeq("AAAGGGTTT")
    ref_log_p = fill(log10(0.01), length(reference))
    bandwidth = 6
    local_errors = Rifraf.normalize(ErrorModel(2.0, 0.5, 0.5, 1.0, 1.0))
    local_scores = Scores(local_errors)
    rseq = RifrafSequence(reference, ref_log_p, bandwidth, local_scores)

    template = DNASeq("AAACCCGGGTTT")
    @test !has_single_indels(template, rseq)

    template = DNASeq("AAACCCGGGTTTT")
    @test has_single_indels(template, rseq)

    template = DNASeq("AAA")
    @test !has_single_indels(template, rseq)
end

@testset "single indel proposals" begin
    ref_errors = ErrorModel(10.0, 1e-10, 1e-10, 1.0, 1.0)
    ref_scores = Scores(ref_errors)

    refseq = DNASeq("CGGCGATTT")
    ref_lps = fill(-1.0, length(refseq))
    ref = RifrafSequence(refseq, ref_lps, 10, ref_scores)

    consensus = DNASeq("CTGCCGA")

    proposals = single_indel_proposals(consensus, ref)
    expected = [Deletion(2), Deletion(4), Deletion(5)]
    @test length(proposals) == 1
    @test proposals[1] in expected
end

@testset "alignment proposals" begin
    function _test_alignment_proposals(consensus, seqs, lps, expected;
                                  do_indels::Bool=true)
        params = RifrafParams(bandwidth=6, batch_fixed=false,
                              batch_size=length(seqs))

        errors = ErrorModel(1.0, 5.0, 5.0, 0.0, 0.0)
        scores = Scores(errors)
        if length(lps) == 0
            lps =  Vector{LogProb}[fill(-9.0, length(s)) for s in seqs]
        end
        pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                               for (s, p) in zip(seqs, lps)]
        state = initial_state(consensus, pseqs, DNASeq(), params)
        resample!(state, params)
        realign_rescore!(state, RifrafParams())

        proposals = alignment_proposals(state.Amoves, state.consensus,
                                        pseqs, do_indels)
        @test length(symdiff(Set(proposals), Set(expected))) == 0
    end

    @testset "fast proposals 1" begin
        consensus = DNASeq("ACGAG")
        seqs = [DNASeq("CGTAC"),
                DNASeq("CGAC"),
                DNASeq("CGTAG")]
        expected = [Deletion(1),
                    Insertion(3, DNA_T),
                    Substitution(5, DNA_C)]
        _test_alignment_proposals(consensus, seqs, [], expected)
    end

    @testset "fast proposals 2" begin
        consensus = DNASeq("AA")
        seqs = [DNASeq("AAG"),
                DNASeq("AA"),
                DNASeq("AAG")]
        expected = [Insertion(2, DNA_G)]
        _test_alignment_proposals(consensus, seqs, [], expected)
    end

    @testset "fast proposals 3" begin
        consensus = DNASeq("AA")
        seqs = [DNASeq("GAA"),
                DNASeq("AA"),
                DNASeq("GAA")]
        expected = [Insertion(0, DNA_G)]
        _test_alignment_proposals(consensus, seqs, [], expected)
    end

    @testset "fast proposals 4" begin
        consensus = DNASeq("AA")
        seqs = [DNASeq("AGA"),
                DNASeq("AA"),
                DNASeq("AGA")]
        expected = [Insertion(1, DNA_G)]
        _test_alignment_proposals(consensus, seqs, [], expected)
    end

    @testset "fast proposals - highly confident base" begin
        consensus = DNASeq("AA")
        seqs = [DNASeq("GAA"),
                DNASeq("AA"),
                DNASeq("AA")]
        lps = Vector{LogProb}[[-10.0, -5.0, -5.0],
                                   [-3.0, -5.0],
                                   [-3.0, -50]]
        expected = [Insertion(0, DNA_G)]
        _test_alignment_proposals(consensus, seqs, lps, expected)
    end
end

@testset "candidate scores" begin
    params = RifrafParams(bandwidth=9,
                          do_alignment_proposals=true)

    function _test_candidate_scores(consensus, pseqs, expected)
        state = initial_state(consensus, pseqs, DNASeq(), params)
        resample!(state, params)
        realign_rescore!(state, RifrafParams())
        cands = get_candidates(state, params)
        @test length(cands) == length(expected)
        if length(cands) == length(expected)
            cand_scores = sort([c.score for c in cands])
            exp_scores = sort([c.score for c in expected])
            @test all([a ≈ b for (a, b) in zip(cand_scores, exp_scores)])
        end
    end

    @testset "substitutions" begin
        consensus = DNASeq("TTT")
        seqs = [DNASeq("TAT")]
        lps =  Vector{LogProb}[fill(-1.0, length(s)) for s in seqs]

        scores = Scores(ErrorModel(1.0, 2.0, 2.0, 0.0, 0.0))
        pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                               for (s, p) in zip(seqs, lps)]

        expected = [ScoredProposal(Substitution(2, DNA_A),
                                   sum(pseqs[1].match_scores))]
        _test_candidate_scores(consensus, pseqs, expected)
    end

    @testset "deletion" begin
        consensus = DNASeq("TTT")
        seqs = [DNASeq("TT")]
        lps =  Vector{LogProb}[fill(-1.0, length(s)) for s in seqs]

        scores = Scores(ErrorModel(1.0, 2.0, 2.0, 0.0, 0.0))
        pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                               for (s, p) in zip(seqs, lps)]

        expected = [ScoredProposal(Deletion(3),
                                   sum(pseqs[1].match_scores))]
        _test_candidate_scores(consensus, pseqs, expected)
    end

    @testset "insertion" begin
        # test insertion
        consensus = DNASeq("TT")
        seqs = [DNASeq("TAT")]
        lps =  Vector{LogProb}[fill(-1.0, length(s)) for s in seqs]

        scores = Scores(ErrorModel(1.0, 2.0, 2.0, 0.0, 0.0))
        pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                               for (s, p) in zip(seqs, lps)]

        expected = [ScoredProposal(Insertion(2, DNA_A),
                                   sum(pseqs[1].match_scores))]
        _test_candidate_scores(consensus, pseqs, expected)
    end
end

@testset "full model" begin
    # can't guarantee this test actually passes, since it is random
    n_seqs = 5
    len = 30

    ref_sample_errors = ErrorModel(8.0, 0.0, 0.0, 1.0, 1.0)
    ref_errors = ErrorModel(8.0, 0.1, 0.1, 1.0, 1.0)
    ref_scores = Scores(ref_errors)

    seq_errors = ErrorModel(1.0, 2.0, 2.0, 0.0, 0.0)
    seq_scores = Scores(seq_errors)

    sample_params = Dict(:ref_error_rate => 0.1,
                         :ref_errors => ref_sample_errors,
                         :error_rate => 0.005,
                         :alpha => 1.0,
                         :phred_scale => 1.5,
                         :actual_std => 3.0,
                         :reported_std => 0.3,
                         :seq_errors => seq_errors)

    @testset "full model" begin
        for (use_ref,
             do_alignment_proposals,
             seed_indels,
             indel_correction_only,
             batch_size) in product([true, false],
                                    [true, false],
                                    [true, false],
                                    [true, false],
                                    [3, 6])

            (reference, template, template_error, reads, actual, phreds, sbools,
             tbools) = sample_sequences(n_seqs, len; sample_params...)
            if !use_ref
                reference = DNASeq()
            end

            params = RifrafParams(scores=seq_scores,
                                  ref_scores=ref_scores,
                                  do_alignment_proposals=do_alignment_proposals,
                                  seed_indels=seed_indels,
                                  indel_correction_only=indel_correction_only,
                                  batch_size=batch_size)

            result = rifraf(reads, phreds;
                            reference=reference, params=params)
            @test result.consensus == template
        end
    end
end

@testset "base_probs" begin
    consensus = DNASeq("CGTAC")
    seqs = [DNASeq("CGAC"),
            DNASeq("CGAC"),
            DNASeq("CGAC")]
    lps = Vector{LogProb}[[-9.0, -9.0, -9.0, -9.0],
                               [-9.0, -9.0, -9.0, -9.0],
                               [-9.0, -9.0, -9.0, -9.0]]

    params = RifrafParams(bandwidth=6)
    errors = Rifraf.normalize(ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0))
    scores = Scores(errors)
    pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                           for (s, p) in zip(seqs, lps)]
    state = initial_state(consensus, pseqs, DNASeq(), params)
    resample!(state, params)
    realign_rescore!(state, RifrafParams())
    probs = estimate_probs(state, false)
    @test probs.sub[1, 2] > 0.9
    @test probs.del[1] < 1e-9
    @test probs.del[3] > 0.9
end

@testset "ins_probs" begin
    consensus = DNASeq("CGAT")
    seqs = [DNASeq("CGTAT"),
            DNASeq("CGTAT"),
            DNASeq("CGTAT")]
    params = RifrafParams(bandwidth=6)
    errors = Rifraf.normalize(ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0))
    scores = Scores(errors)

    lps = Vector{LogProb}[[-9.0, -9.0, -9.0, -9.0, -9.0],
                          [-9.0, -9.0, -9.0, -9.0, -9.0],
                          [-9.0, -9.0, -9.0, -9.0, -9.0]]
    pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                           for (s, p) in zip(seqs, lps)]
    state = initial_state(consensus, pseqs, DNASeq(), params)
    resample!(state, params)
    realign_rescore!(state, RifrafParams())
    probs = estimate_probs(state, false)
    @test maximum(probs.ins[1, :]) < 1e-9
    @test probs.ins[3, 4] > 0.9
end

@testset "alignment_probs" begin
    consensus = DNASeq("ACGT")
    seqs = [DNASeq("ACGT"),
            DNASeq("CGT"),
            DNASeq("CCGT")]

    params = RifrafParams(bandwidth=6)
    errors = Rifraf.normalize(ErrorModel(1.0, 1.0, 1.0, 0.0, 0.0))
    scores = Scores(errors)

    ps = [[0.1, 0.1, 0.1, 0.1],
          [0.2, 0.1, 0.1],
          [0.2, 0.1, 0.1, 0.1]]
    lps = [log10.(p) for p in ps]
    pseqs = RifrafSequence[RifrafSequence(s, p, params.bandwidth, scores)
                           for (s, p) in zip(seqs, lps)]
    state = initial_state(consensus, pseqs, DNASeq(), params)
    resample!(state, params)
    realign_rescore!(state, RifrafParams())

    result = alignment_error_probs(length(consensus),
                                   pseqs, state.Amoves)
    indices = sortperm(result)
    expected = [4, 3, 2, 1]
    @test indices == expected
end

@testset "smart_forward_moves" begin
    consensus = DNASeq("AAACCCGGGTTT")

    # make sequence with long insertion
    seq = DNASeq("AAAGGGTTTCCC")
    errors = fill(0.3, length(seq))
    errors[end-3:end] = 0.45
    err_log_p = log10.(errors)
    bandwidth = 1
    scores = Scores(ErrorModel(1.0, 10.0, 10.0, 0.0, 0.0))
    rseq = RifrafSequence(seq, err_log_p, bandwidth, scores)

    # make arrays
    shape = (length(seq) + 1, length(consensus) + 1)
    A = BandedArray(Score, shape, bandwidth, default=-Inf)
    B = BandedArray(Score, shape, bandwidth, default=-Inf)
    Amoves = BandedArray(Rifraf.Trace, shape, bandwidth)
    pvalue = 0.1
    smart_forward_moves!(consensus, rseq, A, B, Amoves, pvalue)
    @test rseq.bandwidth > bandwidth
end


text:
# ---
# title: 41. First Missing Positive
# id: problem41
# author: Indigo
# date: 2021-03-10
# difficulty: Hard
# categories: Array
# link: <https://leetcode.com/problems/first-missing-positive/description/>
# hidden: true
# ---
#
# Given an unsorted integer array `nums`, find the smallest missing positive
# integer.
#
# **Follow up:**  Could you implement an algorithm that runs in `O(n)` time and
# uses constant extra space.?
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,0]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,4,-1,1]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: nums = [7,8,9,11,12]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 300`
#   * `-231 <= nums[i] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function first_missing_positive(nums::Vector{Int})
    len = length(nums)
    for idx in 1:len
        while nums[idx] <= len && nums[idx] > 0 && nums[nums[idx]] != nums[idx]
            num = nums[idx]
            nums[idx], nums[num] = nums[num], nums[idx]
        end
    end
    for i in 1:len
        (nums[i] != i) && return i
    end
end
## @lc code=end


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

using Markdown
using InteractiveUtils

# ╔═╡ 4133e4a2-ed4d-11ea-11cd-0da48c09c2e8
using Pkg

# ╔═╡ 54bd3938-ed4d-11ea-29e4-c1c7020c1400
begin
	using StatisticalRethinking
end

# ╔═╡ 0c5fe59e-ed4f-11ea-35a8-45c3aff83362
md"## Broadcasting.jl"

# ╔═╡ e5a410a0-e61f-4590-9ad6-1228d6060ecd
Pkg.status()

# ╔═╡ 54ccfc74-ed4d-11ea-366d-b738d54ff79a
p_grid = range(0, 1, length = 5)

# ╔═╡ 54d9bce8-ed4d-11ea-279e-3bc3bf465e28
prior = ones(5)

# ╔═╡ 54da91d4-ed4d-11ea-04f5-e9a4f8c659a3
likelihood = pdf.(Binomial.(9, p_grid), 6)

# ╔═╡ 54e45ff4-ed4d-11ea-0288-e3475c709007
posterior = likelihood .* prior

# ╔═╡ 54e8f906-ed4d-11ea-2fee-316c2df1d8f6
posterior ./= sum(posterior)

# ╔═╡ 54f1592c-ed4d-11ea-0210-0be7741d2604
md"##### Broadcasting is another important difference between Julia and R.

You don't have automatic
broadcasting (making functions work the same whether you run it on a single
number or on an array of numbers). However, fixing this is very easy,
you only have to annotate your function call with a dot:
  `f(single_number)` -> `f.(array_of_numbers)`.

That means the translation of the line `likelihood = pdf.(Binomial.(9, p_grid), 6)` goes something like:
1. First or every value in `p_grid` make a binomial distribution with that p value.
2. For every distribution then take the pdf of 6.

The same in the next line, this is elementwise multiplication.

For more infos read [this](https://julialang.org/blog/2017/01/moredots)."

# ╔═╡ e822d97e-ed4f-11ea-2c18-17ecfc81fdbd
md"## End of broadcasting.jl"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
StatisticalRethinking = "2d09df54-9d0f-5258-8220-54c2a3d4fbee"

[compat]
Distributions = "~0.25.16"
StatisticalRethinking = "~4.0.9"
"""

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

julia_version = "1.7.0-rc1"
manifest_format = "2.0"

[[deps.ANSIColoredPrinters]]
git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c"
uuid = "a4c015fc-c6ff-483c-b24f-f7ea428134e9"
version = "0.0.1"

[[deps.AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

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

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "b8d49c34c3da35f220e7295659cd0bab8e739fed"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.33"

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

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "a4d07a1c313392a77042855df46c5f534076fab9"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.0"

[[deps.AxisKeys]]
deps = ["AbstractFFTs", "CovarianceEstimation", "IntervalSets", "InvertedIndices", "LazyStack", "LinearAlgebra", "NamedDims", "OffsetArrays", "Statistics", "StatsBase", "Tables"]
git-tree-sha1 = "8b382307c6195762a5473ba3522a2830c3014620"
uuid = "94b1ba4f-4ee9-5380-92f1-94cde586c3c5"
version = "0.1.19"

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

[[deps.BitTwiddlingConvenienceFunctions]]
deps = ["Static"]
git-tree-sha1 = "652aab0fc0d6d4db4cc726425cadf700e9f473f1"
uuid = "62783981-4cbd-42fc-bca8-16325de8dc4b"
version = "0.1.0"

[[deps.CPUSummary]]
deps = ["Hwloc", "IfElse", "Static"]
git-tree-sha1 = "38d0941d2ce6dd96427fd033d45471e1f26c3865"
uuid = "2a0fbf3d-bb9c-48f3-b0a9-814d99fd7ab9"
version = "0.1.5"

[[deps.CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "3a877c2fc5c9b88ed7259fd0bdb7691aad6b50dc"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.4"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "bd4afa1fdeec0c8b89dad3c6e92bc6e3b0fec9ce"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.6.0"

[[deps.CloseOpenIntervals]]
deps = ["ArrayInterface", "Static"]
git-tree-sha1 = "ce9c0d07ed6e1a4fecd2df6ace144cbd29ba6f37"
uuid = "fb6a15b2-703c-40df-9091-08a04967cfa9"
version = "0.1.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[deps.Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

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

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "bc3930158d2be029e90b7c40d1371c4f54fa04db"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.6"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

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

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "7220bc21c33e990c14f4a9a319b1d242ebc5b269"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.3.1"

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

[[deps.Distributions]]
deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "f4efaa4b5157e0cdb8283ae0b5428bc9208436ed"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.16"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[deps.Documenter]]
deps = ["ANSIColoredPrinters", "Base64", "Dates", "DocStringExtensions", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
git-tree-sha1 = "8b43e37cfb4f4edc2b6180409acc0cebce7fede8"
uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
version = "0.27.7"

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

[[deps.EllipsisNotation]]
deps = ["ArrayInterface"]
git-tree-sha1 = "8041575f021cba5a099a456b4163c9a08b566a02"
uuid = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
version = "1.1.0"

[[deps.FFTW]]
deps = ["AbstractFFTs", "FFTW_jll", "LinearAlgebra", "MKL_jll", "Preferences", "Reexport"]
git-tree-sha1 = "463cb335fa22c4ebacfd1faba5fde14edb80d96c"
uuid = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
version = "1.4.5"

[[deps.FFTW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c6033cc3892d0ef5bb9cd29b7f2f0331ea5184ea"
uuid = "f5851436-0d7a-5f13-b9de-f02708fd171a"
version = "3.3.10+0"

[[deps.FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "6d4b609786127030d09e6b1ee0e2044ec20eb403"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.11"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "29890dfbc427afa59598b8cfcc10034719bd7744"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.6"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.HostCPUFeatures]]
deps = ["BitTwiddlingConvenienceFunctions", "IfElse", "Libdl", "Static"]
git-tree-sha1 = "3169c8b31863f9a409be1d17693751314241e3eb"
uuid = "3e5b6fbb-0976-4d2c-9146-d79de83f2fb0"
version = "0.1.4"

[[deps.Hwloc]]
deps = ["Hwloc_jll"]
git-tree-sha1 = "92d99146066c5c6888d5a3abc871e6a214388b91"
uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
version = "2.0.0"

[[deps.Hwloc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3395d4d4aeb3c9d31f5929d32760d8baeee88aaf"
uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
version = "2.5.0+0"

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

[[deps.IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[deps.IntelOpenMP_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d979e54b71da82f3a65b62553da4fc3d18c9004c"
uuid = "1d5cc7b8-4909-519e-a0f8-d0f5ad9712d0"
version = "2018.0.3+2"

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

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "61aa005707ea2cebf47c8d780da8dc9bc4e0c512"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.4"

[[deps.IntervalSets]]
deps = ["Dates", "EllipsisNotation", "Statistics"]
git-tree-sha1 = "3cc368af3f110a767ac786560045dceddfc16758"
uuid = "8197267c-284f-5f27-9208-e0e47529a953"
version = "0.5.3"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

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

[[deps.KernelDensity]]
deps = ["Distributions", "DocStringExtensions", "FFTW", "Interpolations", "StatsBase"]
git-tree-sha1 = "591e8dc09ad18386189610acafb970032c519707"
uuid = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
version = "0.6.3"

[[deps.LayoutPointers]]
deps = ["ArrayInterface", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static"]
git-tree-sha1 = "d2bda6aa0b03ce6f141a2dc73d0bcb7070131adc"
uuid = "10f19ff3-798f-405d-979b-55457f8fc047"
version = "0.1.3"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LazyStack]]
deps = ["LinearAlgebra", "NamedDims", "OffsetArrays", "Test", "ZygoteRules"]
git-tree-sha1 = "a8bf67afad3f1ee59d367267adb7c44ccac7fdee"
uuid = "1fad7336-0346-5a1a-a56f-a06ba010965b"
version = "0.0.7"

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

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

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

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

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

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "34dc30f868e368f8a17b728a1238f3fcda43931a"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.3"

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

[[deps.LoopVectorization]]
deps = ["ArrayInterface", "CPUSummary", "CloseOpenIntervals", "DocStringExtensions", "HostCPUFeatures", "IfElse", "LayoutPointers", "LinearAlgebra", "OffsetArrays", "PolyesterWeave", "Requires", "SLEEFPirates", "Static", "ThreadingUtilities", "UnPack", "VectorizationBase"]
git-tree-sha1 = "23b22559423dbf0217f28ea58ac05242a7770194"
uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
version = "0.12.79"

[[deps.MCMCDiagnosticTools]]
deps = ["AbstractFFTs", "DataAPI", "Distributions", "LinearAlgebra", "MLJModelInterface", "Random", "SpecialFunctions", "Statistics", "StatsBase", "Tables"]
git-tree-sha1 = "f3f0c23f0ebe11db62ff1e81412919cf7739053d"
uuid = "be115224-59cd-429b-ad48-344e309966f0"
version = "0.1.1"

[[deps.MKL_jll]]
deps = ["Artifacts", "IntelOpenMP_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "Pkg"]
git-tree-sha1 = "5455aef09b40e5020e1520f551fa3135040d4ed0"
uuid = "856f044c-d86e-5d09-b602-aeab76dc8ba7"
version = "2021.1.1+2"

[[deps.MLJModelInterface]]
deps = ["Random", "ScientificTypesBase", "StatisticalTraits"]
git-tree-sha1 = "91ef121a2c458806973c8aaeb502c57b2f0d74b3"
uuid = "e80e1ace-859a-464e-9ed9-23947d8ae3ea"
version = "1.3.2"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "5a5bc6bf062f0f95e62d0fe0a2d99699fed82dd9"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.8"

[[deps.ManualMemory]]
git-tree-sha1 = "9cb207b18148b2199db259adfa923b45593fe08e"
uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
version = "0.1.6"

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

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

[[deps.MonteCarloMeasurements]]
deps = ["Distributed", "Distributions", "LinearAlgebra", "MacroTools", "Random", "RecipesBase", "Requires", "SLEEFPirates", "StaticArrays", "Statistics", "StatsBase", "Test"]
git-tree-sha1 = "fcefbe707597ce0628f19cb355757f281f37f15e"
uuid = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
version = "1.0.2"

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

[[deps.NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[deps.NamedArrays]]
deps = ["Combinatorics", "DataStructures", "DelimitedFiles", "InvertedIndices", "LinearAlgebra", "Random", "Requires", "SparseArrays", "Statistics"]
git-tree-sha1 = "2fd5787125d1a93fbe30961bd841707b8a80d75b"
uuid = "86f7a689-2022-50b4-a561-43c23ac3c673"
version = "0.9.6"

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "fb4530603a1e62aa5ed7569f283d4b47c2a92f61"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.38"

[[deps.NamedTupleTools]]
git-tree-sha1 = "63831dcea5e11db1c0925efe5ef5fc01d528c522"
uuid = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
version = "0.13.7"

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

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "c0e9e582987d36d5a61e650e6e543b9e44d9914b"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.7"

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

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.ParetoSmooth]]
deps = ["AxisKeys", "Distributions", "DocStringExtensions", "FFTW", "InteractiveUtils", "LinearAlgebra", "LoopVectorization", "MCMCDiagnosticTools", "NamedDims", "Polyester", "PrettyTables", "Printf", "Random", "Requires", "Statistics", "StatsFuns", "Tullio"]
git-tree-sha1 = "ee4bb422cf79f696da410cb033e774c03c1c3023"
uuid = "a68b5a21-f429-434e-8bfa-46b447300aac"
version = "0.6.6"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "9d8c00ef7a8d110787ff6f170579846f776133a9"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.4"

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

[[deps.Polyester]]
deps = ["ArrayInterface", "BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "ManualMemory", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities"]
git-tree-sha1 = "21d8a7163d0f3972ade36ca2b5a0e8a27ac96842"
uuid = "f517fe37-dbe3-4b94-8317-1923a5111588"
version = "0.4.4"

[[deps.PolyesterWeave]]
deps = ["BitTwiddlingConvenienceFunctions", "CPUSummary", "IfElse", "Static", "ThreadingUtilities"]
git-tree-sha1 = "371a19bb801c1b420b29141750f3a34d6c6634b9"
uuid = "1d0040c9-8b98-4ee7-8388-3f51789ca0ad"
version = "0.1.0"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "0d1245a357cc61c8cd61934c07447aa569ff22e6"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.1.0"

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "01d341f502250e81f6fec0afe662aa861392a3aa"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.2"

[[deps.RecipesBase]]
git-tree-sha1 = "44a75aa7a527910ee3d1751d1f0e4148698add9e"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.2"

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

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

[[deps.SIMDTypes]]
git-tree-sha1 = "330289636fb8107c5f32088d2741e9fd7a061a5c"
uuid = "94e857df-77ce-4151-89e5-788b33177be4"
version = "0.1.0"

[[deps.SLEEFPirates]]
deps = ["IfElse", "Static", "VectorizationBase"]
git-tree-sha1 = "2e8150c7d2a14ac68537c7aac25faa6577aff046"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
version = "0.6.27"

[[deps.ScientificTypesBase]]
git-tree-sha1 = "9c1a0dea3b442024c54ca6a318e8acf842eab06f"
uuid = "30f210dd-8aff-4c5f-94ba-8e64358c1161"
version = "2.2.0"

[[deps.SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "54f37736d8934a12a200edea2f9206b03bdf3159"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.7"

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

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

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

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "ad42c30a6204c74d264692e633133dcea0e8b14e"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.2"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "a8f30abc7c64a39d389680b74e749cf33f872a70"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.3.3"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3240808c6d463ac46f1c1cd7638375cd22abbccb"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.12"

[[deps.StatisticalRethinking]]
deps = ["AxisKeys", "DataFrames", "Distributions", "DocStringExtensions", "Documenter", "Formatting", "KernelDensity", "LinearAlgebra", "MonteCarloMeasurements", "NamedArrays", "NamedTupleTools", "OrderedCollections", "Parameters", "ParetoSmooth", "PrettyTables", "Random", "Requires", "Statistics", "StatsBase", "StatsFuns", "StructuralCausalModels", "Tables", "Test", "Unicode"]
git-tree-sha1 = "f0280939946d2af83a175202f5e12b80724a761f"
uuid = "2d09df54-9d0f-5258-8220-54c2a3d4fbee"
version = "4.0.9"

[[deps.StatisticalTraits]]
deps = ["ScientificTypesBase"]
git-tree-sha1 = "730732cae4d3135e2f2182bd47f8d8b795ea4439"
uuid = "64bff920-2084-43da-a3e6-9bb72801c0c9"
version = "2.1.0"

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

[[deps.StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[deps.StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

[[deps.StrideArraysCore]]
deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "ManualMemory", "Requires", "SIMDTypes", "Static", "ThreadingUtilities"]
git-tree-sha1 = "1258e25e171aec339866f283a11e7d75867e77d7"
uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
version = "0.2.4"

[[deps.StructuralCausalModels]]
deps = ["CSV", "Combinatorics", "DataFrames", "DataStructures", "Distributions", "DocStringExtensions", "LinearAlgebra", "NamedArrays", "Reexport", "Statistics"]
git-tree-sha1 = "2ff1d30658057a17f3fca1ef501d8911ae66566e"
uuid = "a41e6734-49ce-4065-8b83-aff084c01dfd"
version = "1.1.0"

[[deps.SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "1162ce4a6c4b7e31e0e6b14486a6986951c73be9"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.5.2"

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

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

[[deps.ThreadingUtilities]]
deps = ["ManualMemory"]
git-tree-sha1 = "03013c6ae7f1824131b2ae2fc1d49793b51e8394"
uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
version = "0.4.6"

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[deps.Tullio]]
deps = ["ChainRulesCore", "DiffRules", "LinearAlgebra", "Requires"]
git-tree-sha1 = "0288b7a395fc412952baf756fac94e4f28bfec65"
uuid = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
version = "0.3.2"

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

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

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

[[deps.VectorizationBase]]
deps = ["ArrayInterface", "CPUSummary", "HostCPUFeatures", "Hwloc", "IfElse", "LayoutPointers", "Libdl", "LinearAlgebra", "SIMDTypes", "Static"]
git-tree-sha1 = "3e2385f4ec895e694c24a1d5aba58cb6d27cf8b6"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.21.10"

[[deps.WeakRefStrings]]
deps = ["DataAPI", "Parsers"]
git-tree-sha1 = "4a4cfb1ae5f26202db4f0320ac9344b3372136b0"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.3.0"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "59e2ad8fd1591ea019a5259bd012d7aee15f995c"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.3"

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

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

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

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

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

# ╔═╡ Cell order:
# ╟─0c5fe59e-ed4f-11ea-35a8-45c3aff83362
# ╠═4133e4a2-ed4d-11ea-11cd-0da48c09c2e8
# ╠═e5a410a0-e61f-4590-9ad6-1228d6060ecd
# ╠═54bd3938-ed4d-11ea-29e4-c1c7020c1400
# ╠═54ccfc74-ed4d-11ea-366d-b738d54ff79a
# ╠═54d9bce8-ed4d-11ea-279e-3bc3bf465e28
# ╠═54da91d4-ed4d-11ea-04f5-e9a4f8c659a3
# ╠═54e45ff4-ed4d-11ea-0288-e3475c709007
# ╠═54e8f906-ed4d-11ea-2fee-316c2df1d8f6
# ╟─54f1592c-ed4d-11ea-0210-0be7741d2604
# ╟─e822d97e-ed4f-11ea-2c18-17ecfc81fdbd
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module GitREPL

import Git
import ReplMaker

export gitrepl

const GIT_REPL_MODE_NAME        = "GitREPL.jl Git REPL mode"
const GIT_REPL_MODE_PROMPT_TEXT = "git> "
const GIT_REPL_MODE_START_KEY   = ','

function _gitrepl_parser(repl_input::AbstractString)
    return quote
        repl_input = $(Expr(:quote, repl_input))
        run(`$(GitREPL.Git.git()) $(split(repl_input))`)
        return nothing
    end
end

"""
    gitrepl()
    gitrepl(; kwargs...)

Set up the Git REPL mode.

## Optional keyword arguments

| Name          | Type             | Default Value                        |
|:------------- | :--------------- |:-------------------------------------|
| `mode_name`   | `AbstractString` | `$(repr(GIT_REPL_MODE_NAME))`        |
| `prompt_text` | `AbstractString` | `$(repr(GIT_REPL_MODE_PROMPT_TEXT))` |
| `start_key`   | `AbstractChar`   | `$(repr(GIT_REPL_MODE_START_KEY))`   |

## Descriptions of optional keyword arguments:

| Name          | Description                                         |
|:------------- | :-------------------------------------------------- |
| `mode_name`   | Name of the REPL mode                               |
| `prompt_text` | Prompt text to display when the REPL mode is active |
| `start_key`   | Key to press to enter the REPL mode                 |

"""
function gitrepl(; mode_name::AbstractString   = GIT_REPL_MODE_NAME,
                   prompt_text::AbstractString = GIT_REPL_MODE_PROMPT_TEXT,
                   start_key::AbstractChar     = GIT_REPL_MODE_START_KEY,
                   kwargs...)
    ReplMaker.initrepl(
        _gitrepl_parser;
        mode_name,
        prompt_text,
        start_key,
        kwargs...
    )
    return nothing
end

function __init__()
    try
        gitrepl()
    catch ex
        msg = "Unable to automatically initialize the Git REPL mode"
        @debug(msg, exception=(ex, catch_backtrace()))
        @warn(msg)
    end
    return nothing
end

end # module


text:
using Test
using TopOpt
using TopOpt.TopOptProblems: getE
using TopOpt.TrussTopOptProblems: load_truss_json, load_truss_geo
using Base.Iterators
# if get(ENV, "CI", nothing) != "true"
#     import Makie
#     @eval using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
# end

ins_dir = joinpath(@__DIR__, "instances", "ground_meshes");

@testset "Test parsing $file_format" for file_format in [".geo", ".json"]
    file_name = "tim_2d" * file_format
    problem_file = joinpath(ins_dir, file_name)
    mat = TrussFEAMaterial(1.0, 0.3)
    crossec = TrussFEACrossSec(800.0)
    if file_format == ".geo"
        node_points, elements, fixities, load_cases = load_truss_geo(problem_file)
        loads = load_cases[1]
    else
        node_points, elements, _, _, fixities, load_cases = load_truss_json(problem_file)
        loads = load_cases["0"]
    end
    problem = TrussProblem(
        Val{:Linear}, node_points, elements, loads, fixities, mat, crossec
    )

    solver = FEASolver(Direct, problem)
    solver()
    @test !all(@. isnan(solver.u))
end

@testset "Tim problem to solve" for (problem_dim, lc_ind) in product(["2d", "3d"], [0, 1])
    problem_dim = "2d"
    file_name = "tim_2d.json"
    lc_ind = 1

    file_name = "tim_$(problem_dim).json"
    problem_file = joinpath(ins_dir, file_name)

    node_points, elements, mats, crosssecs, fixities, load_cases = load_truss_json(
        problem_file
    )
    loads = load_cases[string(lc_ind)]

    problem = TrussProblem(
        Val{:Linear}, node_points, elements, loads, fixities, mats, crosssecs
    )

    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    @test getE(problem) == [m.E for m in mats]

    V = 0.3 # volume fraction
    xmin = 0.001 # minimum density
    rmin = 4.0 # density filter radius

    penalty = TopOpt.PowerPenalty(1.0) # 1
    solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)
    ## call solver to trigger assemble!
    solver()

    # * Compliance
    comp = TopOpt.Compliance(problem, solver)
    obj = comp
    volfrac = TopOpt.Volume(problem, solver)
    constr = x -> volfrac(x) - V

    options = MMAOptions(; maxiter=3000, tol=Nonconvex.Tolerance(; kkt=0.001))
    convcriteria = Nonconvex.KKTCriteria()
    x0 = fill(V, length(solver.vars))
    nelem = length(x0)

    m = Model(obj)
    addvar!(m, zeros(nelem), ones(nelem))
    add_ineq_constraint!(m, constr)

    TopOpt.setpenalty!(solver, penalty.p)
    result = Nonconvex.optimize(m, MMA87(), x0; options=options)

    println("="^10)
    println(
        "tim-$(problem_dim) - LC $(lc_ind) - #elements $(ncells), #dof: $(ncells*ndim): opt iter $(result.iter)",
    )
    println("$(result.convstate)")

    # if get(ENV, "CI", nothing) != "true"
    #     solver()
    #     fig = visualize(
    #         problem, solver.u; topology = result.topology, vector_arrowsize = 0.1,
    #         vector_linewidth=0.8, default_exagg_scale=ndim == 3 ? 1.0 : 0.01,
    #         exagg_range = ndim == 3 ? 10.0 : 0.1,
    #     )
    #     Makie.display(fig)
    # end
end # end testset

@testset "PointLoadCantileverTruss" for dim in [2, 3]
    nels = dim == 2 ? (10, 4) : (10, 4, 4)
    cell_size = Tuple(ones(Float32, dim))

    problem = PointLoadCantileverTruss(nels, cell_size; k_connect=1)

    V = 0.1 # volume fraction
    xmin = 0.001 # minimum density
    rmin = 4.0 # density filter radius

    penalty = TopOpt.PowerPenalty(1.0) # 1
    solver = FEASolver(Direct, problem; xmin=xmin, penalty=penalty)
    ## call solver to trigger assemble!
    solver()

    # * Compliance
    comp = TopOpt.Compliance(problem, solver)
    obj = comp
    volfrac = TopOpt.Volume(problem, solver)
    constr = x -> volfrac(x) - V

    options = MMAOptions(; maxiter=3000, tol=Nonconvex.Tolerance(; kkt=0.001))
    convcriteria = Nonconvex.KKTCriteria()
    x0 = fill(V, length(solver.vars))
    nelem = length(x0)

    m = Model(obj)
    addvar!(m, zeros(nelem), ones(nelem))
    add_ineq_constraint!(m, constr)

    TopOpt.setpenalty!(solver, penalty.p)
    result = Nonconvex.optimize(m, MMA87(), x0; options=options)

    # if get(ENV, "CI", nothing) != "true"
    #     fig = visualize(
    #         problem; topology = result.minimizer)
    #     Makie.display(fig)
    # end
end # end testset


text:
"""
    fivePointStencil(f,x,n::Integer,h)

Evaluates function f and its derivatives up to order n ∈ [0,4] at x:
``f(x),f'(x),...,f^{(n)}(x)``
The result is an array of length n+1.
Derivatives are numerically computed using the 5-point stencil method
with h≠0 being the grid spacing:
[https://en.wikipedia.org/wiki/Five-point_stencil](https://en.wikipedia.org/wiki/Five-point_stencil)

"""
function fivePointStencil(f,x,n::Integer,h)
    if ((n<0) || (n>4))
        error("In fivePointStencil: Invalid order $n")
    end
    if (h==0.0)
        error("In fivePointStencil: Invalid grid spacing $h")
    end
    fm2 = f(x-2h)
    fm1 = f(x-h)
    fn0 = f(x)
    fp1 = f(x+h)
    fp2 = f(x+2h)

    res = Vector{typeof(fn0)}(undef,n+1)

    res[1] = fn0

    (n==0) && return res

    res[2] = (-fp2+8fp1-8fm1+fm2)/(12*h)

    (n==1) && return res

    h2 = h * h
    res[3] = (-fp2+16fp1-30fn0+16fm1-fm2)/(12*h2)

    (n==2) && return res

    h3 = h2 * h
    res[4] = (fp2-2fp1+2fm1-fm2)/(2*h3)

    (n==3) && return res

    h4 = h3 * h
    res[5] = (fp2-4fp1+6fn0-4fm1+fm2)/h4

    return res
end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "src")
AUX_DIR = joinpath(SRC_DIR, "auxiliary")
OUT_DIR=joinpath(Base.source_dir(), "..", "output")
mkpath(OUT_DIR)

include(joinpath(SRC_DIR, "fitzHughNagumo.jl"))
include(joinpath(SRC_DIR, "fitzHughNagumo_conjugateUpdt.jl"))

include(joinpath(SRC_DIR, "types.jl"))
include(joinpath(SRC_DIR, "vern7.jl"))
include(joinpath(SRC_DIR, "tsit5.jl"))
include(joinpath(SRC_DIR, "rk4.jl"))
include(joinpath(SRC_DIR, "ralston3.jl"))
include(joinpath(SRC_DIR, "priors.jl"))
include(joinpath(SRC_DIR, "guid_prop_bridge.jl"))
include(joinpath(SRC_DIR, "random_walk.jl"))
include(joinpath(SRC_DIR, "blocking_schedule.jl"))
include(joinpath(SRC_DIR, "starting_pt.jl"))
include(joinpath(SRC_DIR, "mcmc.jl"))
include(joinpath(SRC_DIR, "path_to_wiener.jl"))

include(joinpath(AUX_DIR, "read_and_write_data.jl"))
include(joinpath(AUX_DIR, "transforms.jl"))


text:
__precompile__(false)

module Hydra
  using MacroTools
  using MacroTools: @forward
  using IRTools: meta, varargs!, argnames!, spliceargs!, update!, IR, @code_ir, iscontrol
  using IRTools: blocks, map, xcall, isgoto, block, blockidx, blocks, successors, label, Block
  using IRTools: block!
  using IRTools.Wrap: CFG
  import Core.Compiler: userefs

  export @spmd, lanewidth, lane, lanesum

  include("vec.jl")

  include("interface.jl")
  include("lib/general.jl")
  include("lib/numeric.jl")
  include("lib/array.jl")
  include("lib/lane.jl")
  include("uniform.jl")
  include("pass.jl")
end


text:
function rel_err(sol, ref)
    maximum(abs.((sol.d[:,end] .- ref) ./ ref))
end

function l2_err(sol, ref)
    sqrt(sum((sol.d[:,end] .- ref).^2) / sum(ref.^2))
end


text:
# All the VAR.jl, SVAR.jl, SVEC.jl (will) have the following components
# - estimateXXXX function
# - XXXXEstimate type
# - λ function meaning the log-likelihood
# - Phi to get the Wold representation
# - Psi to get the orthogonalized Wold representation
# - show function that gives the output

# Additional machinery such as tests should be properly referenced to other files

# In the long run, the following should be supported:
# - Analytical plots using the PyPlot library
# - Automatic selection using information criteria (or at least giving table with the numbers)
# - Macro language for the restricted VAR using the GLS.

module VARmodels
	using Formatting
	using Distributions

	abstract VARRepresentation

	include("dataManipulations.jl")
	include("VAR.jl")
	include("restrictions.jl")
	include("bindings.jl")
	include("fevd.jl")
	include("plots.jl")
	include("cointegrationJohansen.jl")

	# Main VAR machinery
	export
		varEstimate,
		varEstimateShrink,
		Psi,
		restrictVAR,
		restrictVAR2,
		getR,
		Phi

	# Main VECM machinery
	export
		cointegrationJohansen,
		johansenTeststat,
		CointegrationEstimate,
		getVARrepresentation,
		varRepresentationVECM

	# FEVD functions
	export
		fevd,
		genFEVD,
		fftFEVD,
		fftGenFEVD

	# Simulation functions
	export
		simulateVAR

	# Plotting
	export
		plotIRF

end # module


text:
__precompile__()

module ForwardDiff

using DiffBase
using DiffBase: DiffResult
using StaticArrays

import NaNMath
import SpecialFunctions
import RealInterface
import CommonSubexpressions

#############################
# types/functions/constants #
#############################

const NANSAFE_MODE_ENABLED = false

const REAL_TYPES = (AbstractFloat, Irrational, Integer, Rational, Real)

const DEFAULT_CHUNK_THRESHOLD = 10

struct Chunk{N} end

function Chunk(input_length::Integer, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    N = pickchunksize(input_length, threshold)
    return Chunk{N}()
end

function Chunk(x::AbstractArray, threshold::Integer = DEFAULT_CHUNK_THRESHOLD)
    return Chunk(length(x), threshold)
end

# Constrained to `N <= threshold`, minimize (in order of priority):
#   1. the number of chunks that need to be computed
#   2. the number of "left over" perturbations in the final chunk
function pickchunksize(input_length, threshold = DEFAULT_CHUNK_THRESHOLD)
    if input_length <= threshold
        return input_length
    else
        nchunks = round(Int, input_length / DEFAULT_CHUNK_THRESHOLD, RoundUp)
        return round(Int, input_length / nchunks, RoundUp)
    end
end

############
# includes #
############

include("partials.jl")
include("dual.jl")
include("config.jl")
include("utils.jl")
include("derivative.jl")
include("gradient.jl")
include("jacobian.jl")
include("hessian.jl")
include("deprecated.jl")

export DiffBase

end # module


text:
module m4thsymmcheck1
using FinEtools
using FinEtoolsDeforLinear: tens4checksymmetry, tens4symmt6x6tot!, tens4symmtto6x6t!
using LinearAlgebra
using Test
function test()
	C = rand(6, 6)
	C = C + C'
	Co = deepcopy(C)
	t = fill(0.0, 3, 3, 3, 3)
	tens4symmt6x6tot!(t, C)
	@test tens4checksymmetry(t)
	tens4symmtto6x6t!(C, t)
	@assert norm(C - C') <= eps(1.0)
	@assert norm(C - Co) <= eps(1.0)
end
end
using .m4thsymmcheck1
m4thsymmcheck1.test()

module mlumpedmass1
using FinEtools
using FinEtoolsDeforLinear
using Test
import Arpack: eigs
function test()
    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 10;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;

    MR = DeforModelRed3D
    fens,fes  = H20block(a,b,h, na,nb,nh)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material=MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,2)), material)
    K = stiffness(femm, geom, u)

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    @test abs(fs[7]-0.26259869196259) < 1.0e-5

    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3,3)), material)
    M = mass(femm, SysmatAssemblerSparseHRZLumpingSymm(), geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    @test abs(fs[7]- 0.2598164590380608) < 1.0e-5
    # println("Eigenvalues: $fs [Hz]")

    # mode = 17
    # scattersysvec!(u, v[:,mode])
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])

    true
end
end
using .mlumpedmass1
mlumpedmass1.test()


text: my_f(x, y) = 2x + y


text:
struct Uploader <: System end

function Overseer.requested_components(::Uploader)
    return (Mesh, BufferColor, DefaultVao, DefaultProgram, LineVao, LineProgram,
            PeelingProgram, PeelingVao, LineGeometry, Alpha, Visible)
end

shaders(::Type{DefaultProgram}) = default_shaders()
shaders(::Type{PeelingProgram}) = peeling_shaders()
shaders(::Type{InstancedDefaultProgram}) = instanced_default_shaders()
shaders(::Type{InstancedPeelingProgram}) = instanced_peeling_shaders()
shaders(::Type{PeelingCompositingProgram}) = peeling_compositing_shaders()
shaders(::Type{CompositingProgram}) = compositing_shaders()
shaders(::Type{BlendProgram}) = blending_shaders()
shaders(::Type{FXAAProgram}) = fxaa_shaders()
#TODO there is some preparation doubling here!

#TODO cleanup: really not the place to do this
function Overseer.prepare(::Uploader, dio::Diorama)
    for prog in components(dio, RenderProgram)
        if isempty(prog)
            ProgType = eltype(prog)
            dio[Entity(1)] = ProgType(Program(shaders(ProgType)))
        end
    end
end

function Overseer.update(::Uploader, m::AbstractLedger)
    mesh, bcolor, ucolor = m[Mesh], m[BufferColor], m[UniformColor]
    default_vao = m[DefaultVao]
    peeling_vao = m[PeelingVao]
    line_vao = m[LineVao]
    peeling_prog = singleton(m, PeelingProgram)
    default_prog = singleton(m, DefaultProgram)
    line_prog = singleton(m, LineProgram)
    uc = singleton(m, UpdatedComponents)
    alpha = m[Alpha]

    #Buffer color entities are always not instanced
    for e in @entities_in(mesh && bcolor && !default_vao && !peeling_vao)
        gen_vao = prog -> begin
            buffers = [generate_buffers(prog, e.mesh);
                       generate_buffers(prog, GEOMETRY_DIVISOR; color = e.color)]
            return VertexArray(buffers, map(x -> x .- GLint(1), faces(e.mesh)))
        end

        if e in alpha && alpha[e].α < 1
            peeling_vao[e] = PeelingVao(gen_vao(peeling_prog))
        else
            default_vao[e] = DefaultVao(gen_vao(default_prog))
        end
    end
    line_geom = m[LineGeometry]
    #Line Entities
    for e in @entities_in(line_geom)
        vert_loc = attribute_location(line_prog.program, :vertices)
        color_loc = attribute_location(line_prog.program, :color)

        if e in ucolor
            color_vec = fill(ucolor[e].color, length(e.points))
        elseif e in bcolor
            color_vec = bcolor[e].color
        else
            continue
        end

        if !(e in line_vao)
            color_attach = BufferAttachmentInfo(:color, color_loc, Buffer(color_vec),
                                                GEOMETRY_DIVISOR)
            points_attach = BufferAttachmentInfo(:vertices, vert_loc, Buffer(e.points),
                                                 GEOMETRY_DIVISOR)
            line_vao[e] = LineVao(VertexArray([points_attach, color_attach],
                                              GL_LINE_STRIP_ADJACENCY))
        else
            GLA.upload_data!(GLA.bufferinfo(line_vao[e].vertexarray, :vertices).buffer,
                             e.points)
            GLA.upload_data!(GLA.bufferinfo(line_vao[e].vertexarray, :color).buffer,
                             color_vec)
        end
    end

    # Add rendered entities to the Visible component if it's not existing yet.
    # I.e. if not set to be invisible they will be visible
    vis = m[Visible]
    for comp in (default_vao, peeling_vao, line_vao)
        for e in @entities_in(comp && !vis)
            vis[e] = Visible()
        end
    end
end

struct InstancedUploader <: System end

function Overseer.update(::InstancedUploader, m::AbstractLedger)
    default_prog = m[InstancedDefaultProgram][1].program
    peeling_prog = m[InstancedPeelingProgram][1].program
    default_vao  = m[InstancedDefaultVao]
    peeling_vao  = m[InstancedPeelingVao]
    idc          = m[Selectable]
    mesh         = m[Mesh]
    if isempty(mesh)
        return
    end
    ucolor = m[UniformColor]
    modelmat = m[ModelMat]
    material = m[Material]
    alpha = m[Alpha]
    vis = m[Visible]

    default_entities = @entities_in(mesh && ucolor && modelmat && material && !alpha)
    peeling_entities = @entities_in(mesh && ucolor && modelmat && material && alpha)
    n_default = length(default_entities)
    n_peeling = length(peeling_entities)

    ndef = isempty(default_vao) ? 0 : sum(x -> x.vertexarray.ninst, default_vao)
    npeel = isempty(peeling_vao) ? 0 : sum(x -> x.vertexarray.ninst, peeling_vao)

    if n_default == ndef && n_peeling == npeel
        return
    end

    max_entities = maximum(mesh.group_size)
    modelmats    = Vector{ModelMat}(undef, max_entities)
    materials    = Vector{Material}(undef, max_entities)
    ids          = Vector{Entity}(undef, max_entities)
    colors       = Vector{UniformColor}(undef, max_entities)
    idcolors     = Vector{RGBf0}(undef, max_entities)
    alphas       = Vector{Alpha}(undef, max_entities)

    for (i, m) in enumerate(mesh)
        default_it = @entities_in(entity_group(mesh, i) &&
                                  ucolor &&
                                  modelmat &&
                                  material &&
                                  !alpha)
        peeling_it = @entities_in(entity_group(mesh, i) &&
                                  ucolor &&
                                  modelmat &&
                                  material &&
                                  alpha)
        for (it, vao, prog, vao_T) in
            zip((default_it, peeling_it), (default_vao, peeling_vao),
                (default_prog, peeling_prog), (InstancedDefaultVao, InstancedPeelingVao))
            tot = length(it)
            if i <= length(vao.group_size) && tot == vao.data[i].vertexarray.ninst
                # Nothing to be done for this meshgroup
                continue
            end

            for (ie, e) in enumerate(it)
                modelmats[ie] = e[ModelMat]
                materials[ie] = e[Material]
                ids[ie] = e.e
                colors[ie] = e[UniformColor]
                if e in idc
                    idcolors[ie] = idc[e].color
                else
                    idcolors[ie] = RGBf0(1, 1, 1)
                end
                if !(e in vis)
                    vis[e] = Visible()
                end
                if e in alpha
                    alphas[ie] = vis[e].visible ? alpha[e] : Alpha(0)
                else
                    alphas[ie] = vis[e].visible ? Alpha(1) : Alpha(0)
                end
            end

            if (tot > 0 && length(vao.group_size) < i)
                # Mesh needs to be uploaded, only time when a vao gets created
                buffers = [generate_buffers(prog, m.mesh);
                           generate_buffers(prog, GLA.UNIFORM_DIVISOR;
                                            color = view(colors, 1:tot),
                                            modelmat = view(modelmats, 1:tot),
                                            material = view(materials, 1:tot),
                                            object_id_color = view(idcolors, 1:tot),
                                            alpha = view(alphas, 1:tot))]
                vao[ids[1]] = vao_T(VertexArray(buffers,
                                                map(x -> x .- GLint(1), m.mesh.faces), tot))
                for e in ids[2:tot]
                    vao[e] = ids[1]
                end

            elseif i <= length(vao.group_size) && tot != vao.data[i].vertexarray.ninst
                # buffers need to be reuploaded because entities were added
                upload_buffer!(vao.data[i], view(colors, 1:tot))
                upload_buffer!(vao.data[i], view(modelmats, 1:tot))
                upload_buffer!(vao.data[i], view(materials, 1:tot))
                upload_buffer!(vao.data[i], view(idcolors, 1:tot), :object_id_color)
                upload_buffer!(vao.data[i], view(alphas, 1:tot))
                vao.data[i].vertexarray.ninst = tot
                to_remove = Entity[]
                for e in entity_group(vao, i)
                    if !in(e, it)
                        push!(to_remove, e)
                    end
                end
                delete!(vao, to_remove)
                p = parent(vao, i)
                for e in ids
                    if e == p
                        continue
                    end
                    vao[e] = p
                end
            end
        end
    end
end

struct UniformUploader <: System end

function Overseer.requested_components(::UniformUploader)
    return (InstancedDefaultVao, InstancedPeelingVao, ModelMat, Selectable, UniformColor,
            UpdatedComponents)
end

# function find_contiguous_bounds(indices)
#     ranges = UnitRange[]
#     i = 1
#     cur_start = indices[1]
#     while i <= length(indices) - 1
#         id = indices[i]
#         id_1 = indices[i + 1]
#         if id_1 - id != 1
#             push!(ranges, cur_start:id)
#             cur_start = id_1
#         end
#         i += 1
#     end
#     push!(ranges, cur_start:indices[end])
#     return ranges
# end

"Uploads buffer using glBufferSubData."
function upload_buffer!(vao, vec::AbstractVector{T}, sym = shader_symbol(T)) where {T}
    binfo = GLA.bufferinfo(vao.vertexarray, sym)
    if binfo !== nothing
        datsize = sizeof(T)
        b = binfo.buffer
        GLA.bind(b)
        n = length(vec)
        s = n * datsize
        if n <= length(b)
            glBufferSubData(b.buffertype, 0, s, pointer(vec, 1))
        else
            glBufferData(b.buffertype, s, pointer(vec, 1), b.usage)
        end
        b.len = n
        GLA.unbind(binfo.buffer)
    end
end

function upload_to_vao!(f::Function, vao_comp, comp,
                        buffer = Vector{eltype(comp)}(undef, maximum(vao_comp.group_size)))
    for (i, v) in enumerate(vao_comp)
        for (ie, e) in enumerate(entity_group(vao_comp, i))
            buffer[ie] = f(comp, e)
        end
        upload_buffer!(v, view(buffer, 1:vao_comp.group_size[i]))
    end
end
function upload_to_vao!(vao_comp::Overseer.AbstractComponent, args...)
    return upload_to_vao!((c, e) -> (@inbounds c[e]), vao_comp, args...)
end

function Overseer.update(::UniformUploader, m::AbstractLedger)
    uc = m[UpdatedComponents][1]
    vis = m[Visible]
    alpha = m[Alpha]
    peel = m[InstancedPeelingVao]
    default = m[InstancedDefaultVao]
    maxsize = max(isempty(peel) ? 0 : maximum(peel.group_size),
                  isempty(default) ? 0 : maximum(default.group_size))
    if maxsize == 0
        return
    end

    # always reupload visible
    push!(uc.components, Alpha)
    for c in filter(x -> shader_symbol(x) != :none, uc.components)
        comp = m[c]
        vec  = Vector{eltype(comp)}(undef, maxsize)
        for vao in (default, peel)
            if isempty(vao)
                continue
            end
            if c == Alpha
                upload_to_vao!(vao, comp, vec) do a, e
                    if e in a
                        return @inbounds vis[e].visible ? a[e] : Alpha(0.0)
                    else
                        return @inbounds vis[e].visible ? Alpha(1.0) : Alpha(0.0)
                    end
                end
            else
                upload_to_vao!(vao, comp, vec)
            end
        end
    end
end