text:
#################### Hamiltonian Monte Carlo ####################

#################### Types and Constructors ####################

mutable struct HMCTune <: SamplerTune
    logf::Union{Function,Missing}
    epsilon::Float64
    L::Int
    SigmaL::Union{UniformScaling{Bool},LowerTriangular{Float64}}

    HMCTune() = new()

    HMCTune(x, epsilon::Real, L::Integer) = new(missing, epsilon, L, I)

    HMCTune(x, epsilon::Real, L::Integer, logfgrad::Function) = new(logfgrad, epsilon, L, I)

    function HMCTune(x, epsilon::Real, L::Integer, Sigma::Matrix{T}) where {T<:Real}
        new(missing, epsilon, L, cholesky(Sigma).L)
    end

    function HMCTune(
        x,
        epsilon::Real,
        L::Integer,
        Sigma::Matrix{T},
        logfgrad::Function,
    ) where {T<:Real}
        new(logfgrad, epsilon, L, cholesky(Sigma).L)
    end

    function HMCTune(
        x,
        epsilon::Real,
        L::Integer,
        Sigma::UniformScaling{Bool},
        logfgrad::Function,
    )
        new(logfgrad, epsilon, L, Sigma)
    end
end

const HMCVariate = Sampler{HMCTune,T} where {T}

validate(v::HMCVariate) = validate(v, v.tune.SigmaL)

validate(v::HMCVariate, SigmaL::UniformScaling) = v

function validate(v::HMCVariate, SigmaL::LowerTriangular)
    n = length(v)
    size(SigmaL, 1) == n ||
        throw(ArgumentError("Sigma dimension differs from variate length $n"))
    v
end

#################### Sampler Constructor ####################
"""
    HMC(params::ElementOrVector{Symbol}, epsilon::Real, L::Integer; args...)

Construct a `Sampler` object for HMC sampling. Parameters are assumed to be
continuous, but may be constrained or unconstrained.

Returns a `Sampler{HMCTune}` type object.

* `params`: stochastic node(s) to be updated with the sampler. Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic `unlist()` function.

* `epsilon`: step size.

* `L`: number of steps to take in the Leapfrog algorithm.

* `Sigma`: covariance matrix for the multivariate normal proposal distribution. The covariance matrix is relative to the unconstrained parameter space, where candidate draws are generated. If omitted, the identity matrix is assumed.

"""
function HMC(params::ElementOrVector{Symbol}, epsilon::Real, L::Int)
    tune = HMCTune(Float64[], epsilon, L, I, logpdfgrad!)
    Sampler(Float64[], params, tune, Symbol[], true)
end

function HMC(params::ElementOrVector{Symbol}, epsilon::Real, L::Int, Sigma::Matrix{<:Real})
    tune = HMCTune(Float64[], epsilon, L, Sigma, logpdfgrad!)
    Sampler(Float64[], params, tune, Symbol[], true)
end

#################### Sampling Functions ####################

"""
    sample!(v::HMCVariate, logfgrad::Function)

Draw one sample from a target distribution using the HMC sampler. Parameters
are assumed to be continuous and unconstrained.

Returns `v` updated with simulated values and associated tuning parameters.
"""
function sample!(v::HMCVariate, logfgrad::Function; kwargs...)
    tune = v.tune

    x1 = v.value[:]
    logf0, grad0 = logf1, grad1 = logfgrad(x1)

    ## Momentum variables
    p0 = p1 = tune.SigmaL * randn(length(v))

    ## Make a half step for a momentum at the beginning
    p1 += 0.5 * tune.epsilon * grad0

    ## Alternate full steps for position and momentum
    for i = 1:tune.L
        ## Make a full step for the position
        x1 += tune.epsilon * p1

        logf1, grad1 = logfgrad(x1)

        ## Make a full step for the momentum
        p1 += tune.epsilon * grad1
    end

    ## Make a half step for momentum at the end
    p1 -= 0.5 * tune.epsilon * grad1

    ## Negate momentum at end of trajectory to make the proposal symmetric
    p1 *= -1.0

    ## Evaluate potential and kinetic energies at start and end of trajectory
    SigmaLinv = inv(tune.SigmaL)
    Kp0 = 0.5 * sum(abs2, SigmaLinv * p0)
    Kp1 = 0.5 * sum(abs2, SigmaLinv * p1)

    if rand() < exp((logf1 - Kp1) - (logf0 - Kp0))
        v[:] = x1
    end

    v
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
module KeyOutputs

using CSV
using DataFrames

function getKeyoutputs(file)
    df = CSV.read(file,
                  delim="\t",
                  datarow=1,
                  quotechar="\\",
                  nullable=false,
                  header=["pathway_id",
                          "pathway_name",
                          "node_id",
                          "node_name"],
                  types=[String,
                         String,
                         String,
                         String])

    keyoutputList = Dict()
    for row in eachrow(df)
        if !isnothing(row[:pathway_name])
            pathwayName = row[:pathway_name]
            if haskey(keyoutputList, pathwayName) == false
               keyoutputList[pathwayName] = Dict()
            end
            keyoutputList[pathwayName][row[:node_id]] = row[:node_name]
        end
    end

    return keyoutputList
end

end

ext: jl

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
module TestFidelityUMI

using Base.Test
NEW_PATH=join(split(Base.source_path(),"/")[1:end-3],"/")
push!(LOAD_PATH, NEW_PATH)
using BioAlignments
using AutoHashEquals
using DataStructures
using DataFrames
using BioSequences
using BioAlignments
using FidelityUMI

println(LOAD_PATH)

bam_file1 = NEW_PATH*"/FidelityUMI.jl/tests/1_duplicated_subSortbyName.bam"
bam_file2 = NEW_PATH*"/FidelityUMI.jl/tests/2_duplicated_subSortbyName.bam"
reference_file = NEW_PATH*"/FidelityUMI.jl/tests/umi_fidelity.fasta"
gene_csv = NEW_PATH*"/FidelityUMI.jl/tests/genes.csv"

@testset "FidelityUMI" begin
    referenceSruct = ReferenceSequences(reference_file)
    cvCoverage = Coverage(referenceSruct)
    coverage_counter!(cvCoverage, bam_file1; mappingQuality=1)
    clust_coverages = Coverage(referenceSruct)

    @test String["lambda"] == cvCoverage.chromosomes
    @test Dict("lambda"=>379) == cvCoverage.chrom_lengths
    @test Dict("lambda"=>[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 8.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 17.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 18.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 16.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
              ) == cvCoverage.coverage

    all_variants = Dict(Any[])
    reader = open(BAM.Reader, bam_file2)
    record = BAM.Record()
    while !eof(reader)
        read!(reader, record)
        rname = BAM.tempname(record)
        flag1 = BAM.flag(record)
        variants = call_variants(record; quality=0, verbose=false)

        if flag1&64 == 64
            rname = "R1:"*rname
        elseif flag1&128 == 128
            rname = "R2:"*rname
        end

        if length(variants) > 0
            if rname in keys(all_variants)
                push!(all_variants[rname], variants)
            else
                all_variants[rname] = [variants]
            end
        end
    end
    close(reader)

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('D', "test", 8, 8, "G", "*"),
           FidelityUMI.Mutation('M', "test", 9, 9, "C", "A"),
           FidelityUMI.Mutation('I', "test", 3, 4, "-", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:insert:r1:4A:deletion:8G:mistmatch:r1:C9A_CCCCCCCCCC"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "A"),
           FidelityUMI.Mutation('I', "test", 3, 4, "-", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:insert:r1:4A:mistmatch:r1:G8A_CCCCCCCCCC"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "C")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:mistmatch:G8C_GGGGGGGGGG"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 4, 4, "T", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:mistmatch:r1:T4A_CCCCCCCCCC"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 46, 46, "T", "A")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:1:1:1:mistmatch:r2:T5A_CCCCCCCCCC"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 4, 4, "T", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:soft:clip:r1:1G:mistmatch:r1:T4A_CCCCCCCCCC"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "C")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:mistmatch:G8C_GGGGGGGGGG"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "C")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:2:mistmatch:G8C_GGGGGGGGGG"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "C")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:3:mistmatch:G8C_GGGGGGGGGG"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "C")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:4:mistmatch:G8C_GGGGGGGGGG"]

    @test [
           FidelityUMI.Mutation[FidelityUMI.Mutation('D', "test", 4, 4, "T", "*"),
           FidelityUMI.Mutation('D', "test", 8, 8, "G", "*"),
           FidelityUMI.Mutation('M', "test", 14, 14, "T", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:deletion:r1:4T:8G:mistmatch:r1:T14A_CCCCCCCCCC"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "test", 8, 9, "GC", "*"),
            FidelityUMI.Mutation('M', "test", 14, 14, "T", "A"),
            FidelityUMI.Mutation('I', "test", 3, 4, "-", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:insert:r1:4A:deletion:8GC:mistmatch:r1:T14A_CCCCCCCCCC"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 14, 14, "T", "A"),
            FidelityUMI.Mutation('I', "test", 7, 8, "-", "AAAA")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:insert:r1:4AAAA:mismatch:T14A_CCCCCCCCCC"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 8, 8, "G", "A"),
            FidelityUMI.Mutation('I', "test", 3, 4, "-", "TA")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:insert:r1:4TA:mistmatch:r1:G8A_CCCCCCCCCC"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 4, 4, "T", "A"),
            FidelityUMI.Mutation('I', "test", 10, 11, "-", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:mismatch:r1:T4A:insert:r1:10G_CCCCCCCCCC"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "test", 4, 4, "T", "A"),
            FidelityUMI.Mutation('M', "test", 7, 7, "T", "A")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:1:1:1:mistmatch:r1:T4A:T7A_CCCCCCCCCC"]

    all_variants = Dict(Any[])
    reader = open(BAM.Reader, bam_file1)
    record = BAM.Record()
    while !eof(reader)
        read!(reader, record)
        rname = BAM.tempname(record)
        variants = call_variants(record; quality=0, verbose=false)
        flag1 = BAM.flag(record)

        if flag1&64 == 64
            rname = "R1:"*rname
        elseif flag1&128 == 128
            rname = "R2:"*rname
        end

        if length(variants) > 0
            if rname in keys(all_variants)
                push!(all_variants[rname], variants)
            else
                all_variants[rname] = [variants]
            end
        end
    end
    close(reader)

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2101:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2101:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2102:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2102:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2103:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2103:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2104:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2104:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2105:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*"),
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2105:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*"),
            FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2106:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*"),
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 143, 144, "AA", "*")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2106:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*"),
            FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*"),
            FidelityUMI.Mutation('I', "lambda", 45, 46, "-", "AA")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2107:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*"),
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R2:M00827:86:000000000-BPT2L:1:2107:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2108:17692:1620_ATTACAGGGGACCATGGCAA"]

    all_variants = Dict(Any[])
    reader = open(BAM.Reader, bam_file1)
    record = BAM.Record()
    while !eof(reader)
        read!(reader, record)
        rname = BAM.tempname(record)
        variants = call_variants(record; quality=1, verbose=false)
        flag1 = BAM.flag(record)

        if flag1&64 == 64
            rname = "R1:"*rname
        elseif flag1&128 == 128
            rname = "R2:"*rname
        end

        if length(variants) > 0
            if rname in keys(all_variants)
                push!(all_variants[rname], variants)
            else
                all_variants[rname] = [variants]
            end
        end
    end
    close(reader)

    @test [
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")]
    ] == all_variants["R1:M00827:86:000000000-BPT2L:1:2109:17692:1620_ATTACAGGGGACCATGGCAA"]

    @test true == check_operations([BioAlignments.Operation(OP_MATCH)])
    @test true == check_operations([BioAlignments.Operation(OP_INSERT)])
    @test true == check_operations([BioAlignments.Operation(OP_DELETE)])
    @test true == check_operations([BioAlignments.Operation(OP_SOFT_CLIP)])
    @test true == check_operations([BioAlignments.Operation(OP_HARD_CLIP)])
    @test false == check_operations([BioAlignments.Operation(OP_SKIP)])
    @test false == check_operations([BioAlignments.Operation(OP_PAD)])
    @test false == check_operations([BioAlignments.Operation(OP_SEQ_MATCH)])
    @test false == check_operations([BioAlignments.Operation(OP_SEQ_MISMATCH)])
    @test false == check_operations([BioAlignments.Operation(OP_BACK)])
    @test false == check_operations([BioAlignments.Operation(OP_START)])

    (errors,
    mut_coverages,
    total_bases1,
    total_reads,
    filtered_reads,
    filtered_R1,
    filtered_R2,
    clusts_with_vars
    ) = collect_variants(bam_file1, clust_coverages, 0, 10, 0, "_", 2, false, "BX"; verbose=false)

    @test String["lambda"] == clust_coverages.chromosomes
    @test Dict("lambda"=>379) == clust_coverages.chrom_lengths
    @test Dict("lambda"=>[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
                ) == clust_coverages.coverage

    (true_errors,
    err_distr,
    clust_size_distr,
    found_mutations,
    found_insertions,
    found_deletions,
    found_ns,
    total_clusters,
    conf_bases,
    clusts_with_vars,
    passing_vars,
    cluster_total_bp,
    conf_coverage
    ) = filter_variants(errors, mut_coverages, cvCoverage, referenceSruct,
                        2, 0.9; verbose=false)

    test_conf_coverage = Dict("lambda"=>UInt32[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
    )

    @test test_conf_coverage == conf_coverage.coverage

    test_errors = Dict(
                "lambda:33:164:ATTACAGGGGACCATGGCAA"=>FidelityUMI.MutationsCluster(
                    Dict(
                        FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*")=>3,
                        FidelityUMI.Mutation('D', "lambda", 143, 144, "AA", "*")=>1,
                        FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*")=>2,
                        FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")=>3,
                        FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*")=>3,
                        FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")=>3,
                        FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>3,
                        FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>3,
                        FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*")=>3,
                        FidelityUMI.Mutation('I', "lambda", 45, 46, "-", "AA")=>1
                    ),
                    3,
                    132,
                    92,
                    10,
                    [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0],
                    [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0],
                    33,
                    164
                ),
                "lambda:63:164:ATTACAGGGGACCATGGCAA"=>FidelityUMI.MutationsCluster(
                    Dict(
                        FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*")=>1,
                        FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")=>1,
                        FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*")=>1,
                        FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")=>1,
                        FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*")=>1,
                        FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>1
                    ),
                    1,
                    80,
                    80,
                    6,
                    [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
                    [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0],
                    63,
                    164
                ),
                "lambda:33:154:ATTACAGGGGACCATGGCAA"=>FidelityUMI.MutationsCluster(
                    Dict(
                        FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")=>5,
                        FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>5,
                        FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>5
                    ),
                    5,
                    122,
                    92,
                    3,
                    [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0],
                    [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0],
                    33,
                    154
                )
            )

    @test test_errors == errors
    @test 18 == total_reads
    @test 1890 == total_bases1
    @test 426 == conf_bases
    @test 1706 == cluster_total_bp
    @test Dict("lambda"=>Dict(
        FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*")=>4,
        FidelityUMI.Mutation('D', "lambda", 143, 144, "AA", "*")=>1,
        FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*")=>2,
        FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")=>18,
        FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*")=>4,
        FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")=>4,
        FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>8,
        FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>9,
        FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*")=>4,
        FidelityUMI.Mutation('I', "lambda", 45, 46, "-", "AA")=>1)
        ) == mut_coverages

    @test 8 == true_errors
    @test Float32[1.0, 0.33, 0.67, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.33, 1.0, 1.0, 1.0] == err_distr
    @test [3, 1, 5] == clust_size_distr
    @test 4 == found_mutations
    @test 0 == found_insertions
    @test 4 == found_deletions
    @test 0 == found_ns
    @test 3 == total_clusters

    test_passing_vars = Dict(
        "lambda:33:164:ATTACAGGGGACCATGGCAA"=>FidelityUMI.MutationsCluster(
            Dict(
                FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*")=>1,
                FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*")=>1,
                FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")=>1,
                FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>1,
                FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>1,
                FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*")=>1
            ),
            3,
            132,
            92,
            6,
            [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0],
            [3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0],
            33,
            164
        ),
        "lambda:33:154:ATTACAGGGGACCATGGCAA"=>FidelityUMI.MutationsCluster(
            Dict(
                FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>1,
                FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>1
            ),
            5,
            122,
            92,
            2,
            [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0],
            [5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0],
            33,
            154
        )
    )

    @test test_passing_vars == passing_vars
    aa_subs, frame_shifts = translate_mutations(referenceSruct, gene_csv, passing_vars; verbose=false)
    @test Dict("lambda"=>Array{String,1}[String["Q27P", "C38W"]]) == aa_subs
    @test 1 == frame_shifts

    test_stringify = [
        "G71*",
        "AA143*",
        "GA57*",
        "T100G",
        "TA114*",
        "C142*",
        "A85C",
        "C119G",
        "AT98*",
        "-45AA",
        "G71*",
        "T100G",
        "TA114*",
        "C142*",
        "AT98*",
        "C119G",
        "T100G",
        "A85C",
        "C119G",
    ]

    ct = 0
    for id in keys(test_errors)
        for var in keys(test_errors[id].mutation_count)
            ct += 1
            @test test_stringify[ct] == stringify(var)
        end
    end

    test_mapp_mutations = Dict(
            FidelityUMI.Mutation[FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]=>1,
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*"),
            FidelityUMI.Mutation('D', "lambda", 143, 144, "AA", "*"),
            FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*"),
            FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G"),
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*"),
            FidelityUMI.Mutation('I', "lambda", 45, 46, "-", "AA")]=>1,
            FidelityUMI.Mutation[FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*"),
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G"),
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*"),
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*"),
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*"),
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")]=>1
    )
    #TODO: cluster=false case
    @test test_mapp_mutations == mapp_mutations(test_errors)

    @test String["P27Q", "W38C"] == compare_seq(
                                        "MQFAHKVSGKYRGVAKLEGNTKAKVLPVLATFAYADYWRSA",
                                        "MQFAHKVSGKYRGVAKLEGNTKAKVLQVLATFAYADYCRSA"
                                    )

    cvCoverage = Coverage(referenceSruct)
    test_conf_coverage = Dict("lambda"=>
        UInt32[
            0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        ]
    )

    conf_coverage!(cvCoverage, test_errors)
    @test test_conf_coverage == cvCoverage.coverage

    clust_coverages = Coverage(referenceSruct)
    (errors_paired,
    mut_coverages_paired,
    total_bases1_paired,
    total_reads_paired,
    filtered_reads_paired,
    filtered_R1_paired,
    filtered_R2_paired,
    clusts_with_vars_paired
    ) = collect_variants(bam_file1, clust_coverages, 0, 10, 0, "_", 2, false, "BX"; verbose=false, algo="paired")

    test_errors_paired = Dict(
        "lambda"=>Dict(
            FidelityUMI.Mutation('D', "lambda", 71, 71, "G", "*")=>4,
            FidelityUMI.Mutation('D', "lambda", 143, 144, "AA", "*")=>1,
            FidelityUMI.Mutation('D', "lambda", 57, 58, "GA", "*")=>2,
            FidelityUMI.Mutation('M', "lambda", 100, 100, "T", "G")=>18,
            FidelityUMI.Mutation('D', "lambda", 114, 115, "TA", "*")=>4,
            FidelityUMI.Mutation('D', "lambda", 142, 142, "C", "*")=>4,
            FidelityUMI.Mutation('M', "lambda", 85, 85, "A", "C")=>8,
            FidelityUMI.Mutation('M', "lambda", 119, 119, "C", "G")=>9,
            FidelityUMI.Mutation('D', "lambda", 98, 99, "AT", "*")=>4,
            FidelityUMI.Mutation('I', "lambda", 45, 46, "-", "AA")=>1
            )
        )

end
end

max_stars_repo_name: thermofisher/rnd-it-molbio-fidelity-umi

max_issues_repo_name: thermofisher/rnd-it-molbio-fidelity-umi

max_forks_repo_name: thermofisher/rnd-it-molbio-fidelity-umi


text:
using Plots, LaTeXStrings; pyplot()

xGrid = 0:0.01:10
uGrid = 0:0.01:1
busy = 0.8

F(t)= t<=0 ? 0 : 1 - busy*exp(-(1-busy)t)

infimum(B) = isempty(B) ? Inf : minimum(B)
invF(u) = infimum(filter((x) -> (F(x) >= u),xGrid))

p1 = plot(xGrid,F.(xGrid), c=:blue, xlims=(-0.1,10), ylims=(0,1),
	xlabel=L"x", ylabel=L"F(x)")

p2 = plot(uGrid,invF.(uGrid), c=:blue, xlims=(0,0.95), ylims=(0,maximum(xGrid)),
	xlabel=L"u", ylabel=L"F^{-1}(u)")

plot(p1, p2, legend=:none, size=(800, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
# Most constants are defined here, except:
#
# BUF           src/Types/SeisIOBuf.jl
# KW            src/Types/KWDefs.jl
# PhaseCat      src/Types/Quake/PhaseCat.jl
# flat_resp     src/Types/InstResp.jl
# RespStage     src/Types/InstResp.jl
# type_codes    src/Types/Methods/0_type_codes.jl
#
# ...and submodule-specific constants, which are in the submodule declaration files (e.g. SEED.jl)

# Type aliases
const FloatArray  = Union{Array{Float64,1}, Array{Float32,1}}
const TimeSpec    = Union{Real,DateTime,String}
const bad_chars = Dict{String, Any}(
  "File" => (0x22, 0x24, 0x2a, 0x2f, 0x3a, 0x3c, 0x3e, 0x3f, 0x40, 0x5c, 0x5e, 0x7c, 0x7e, 0x7f),
  "HTML" => (0x22, 0x26, 0x27, 0x3b, 0x3c, 0x3e, 0xa9, 0x7f),
  "Julia" => (0x24, 0x5c, 0x7f),
  "Markdown" => (0x21, 0x23, 0x28, 0x29, 0x2a, 0x2b, 0x2d, 0x2e, 0x5b, 0x5c, 0x5d, 0x5f, 0x60, 0x7b, 0x7d),
  "SEED" => (0x2e, 0x7f),
  "Strict" => (0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0x28, 0x29, 0x2a,
               0x2b, 0x2c, 0x2d, 0x2e, 0x2f, 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f,
               0x40, 0x5b, 0x5c, 0x5d, 0x5e, 0x60, 0x7b, 0x7c, 0x7d, 0x7e, 0x7f) )
const datafields = (:id, :name, :loc, :fs, :gain, :resp, :units, :src, :notes, :misc, :t, :x)
const days_per_month = Int32[31,28,31,30,31,30,31,31,30,31,30,31]
const dtconst = 62135683200000000
const regex_chars = String[Sys.iswindows() ? "/" : "\\", "\$", "(", ")", "+", "?", "[", "\\0",
"\\A", "\\B", "\\D", "\\E", "\\G", "\\N", "\\P", "\\Q", "\\S", "\\U", "\\U",
"\\W", "\\X", "\\Z", "\\a", "\\b", "\\c", "\\d", "\\e", "\\f", "\\n", "\\n",
"\\p", "\\r", "\\s", "\\t", "\\w", "\\x", "\\x", "\\z", "]", "^", "{", "|", "}"]
const show_os = 8
const sac_keys = (  String[ "delta", "depmin", "depmax", "scale", "odelta",
                            "b", "e", "o", "a", "internal1",
                            "t0", "t1", "t2", "t3", "t4",
                            "t5", "t6", "t7", "t8", "t9",
                            "f", "resp0", "resp1", "resp2", "resp3",
                            "resp4", "resp5", "resp6", "resp7", "resp8",
                            "resp9", "stla", "stlo", "stel", "stdp",
                            "evla", "evlo", "evel", "evdp", "mag",
                            "user0", "user1", "user2", "user3", "user4",
                            "user5", "user6", "user7", "user8", "user9",
                            "dist", "az", "baz", "gcarc", "internal2",
                            "internal3", "depmen", "cmpaz", "cmpinc", "xminimum",
                            "xmaximum", "yminimum", "ymaximum", "unused1", "unused2",
                            "unused3", "unused4", "unused5", "unused6", "unused7" ],
                    String[ "nzyear", "nzjday", "nzhour", "nzmin", "nzsec",
                            "nzmsec", "nvhdr", "norid", "nevid", "npts",
                            "internal4", "nwfid", "nxsize", "nysize", "unused8",
                            "iftype", "idep", "iztype", "unused9", "iinst",
                            "istreg", "ievreg", "ievtyp", "iqual", "isynth",
                            "imagtyp", "imagsrc", "unused10", "unused11", "unused12",
                            "unused13", "unused14", "unused15", "unused16", "unused17",
                            "leven", "lpspol", "lovrok", "lcalda", "unused18" ],
                    String[ "kstnm", "kevnm", "khole", "ko", "ka", "kt0", "kt1", "kt2",
                            "kt3", "kt4", "kt5", "kt6", "kt7", "kt8", "kt9", "kf", "kuser0",
                            "kuser1", "kuser2", "kcmpnm", "knetwk", "kdatrd", "kinst" ] )
const sac_nul_c = UInt8[0x2d, 0x31, 0x32, 0x33, 0x34, 0x35, 0x20, 0x20]
const sac_nul_f = -12345.0f0
const sac_nul_i = Int32(-12345)
const sac_nul_start = 0x2d
const sac_nul_Int8 = UInt8[0x31, 0x32, 0x33, 0x34, 0x35]
const segy_ftypes  = Array{DataType, 1}([UInt32, Int32, Int16, Any, Float32, Any, Any, Int8]) # Note: type 1 is IBM Float32
const seis_inst_codes = ('H', 'J', 'L', 'M', 'N', 'P', 'Z')
const seisio_file_begin = UInt8[0x53, 0x45, 0x49, 0x53, 0x49, 0x4f]
const sμ = 1000000.0
const vSeisIO = Float32(0.53)
const unindexed_fields = (:c, :n)
const webhdr = Dict("User-Agent" => "Mozilla/5.0 (X11; Linux x86_64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/78.0.3904.70 Safari/537.36") # lol
const xml_endtime = 19880899199000000
const μs = 1.0e-6


text:
"The simplest viable forward mode a AD, only supports `Float64`"
module ForwardDiffZero
using ChainRulesCore
using ChainRulesOverloadGeneration
# resolve conflicts while this code exists in both.
const on_new_rule = ChainRulesOverloadGeneration.on_new_rule
const refresh_rules = ChainRulesOverloadGeneration.refresh_rules

using Test

#########################################
# Initial rule setup
@scalar_rule x + y (1, 1)
@scalar_rule x - y (1, -1)
##########################
# Define the AD

# Note that we never directly define Dual Number Arithmetic on Dual numbers
# instead it is automatically defined from the `frules`
struct Dual <: Real
    primal::Float64
    partial::Float64
end

primal(d::Dual) = d.primal
partial(d::Dual) = d.partial

primal(d::Real) = d
partial(d::Real) = 0.0

# needed for `^` to work from having `*` defined
Base.to_power_type(x::Dual) = x

function define_dual_overload(sig)
    sig = Base.unwrap_unionall(sig)  # Not really handling most UnionAlls
    opT, argTs = Iterators.peel(sig.parameters)
    opT isa Type{<:Type} && return  # not handling constructors
    fieldcount(opT) == 0 || return  # not handling functors
    all(argT isa Type && Float64 <: argT for argT in argTs) || return  # only handling purely Float64 ops.

    N = length(sig.parameters) - 1  # skip the op
    fdef = quote
        # we use the function call overloading form as it lets us avoid namespacing issues
        # as we can directly interpolate the function type into to the AST.
        function (op::$opT)(dual_args::Vararg{Union{Dual, Float64}, $N}; kwargs...)
            ȧrgs = (NoTangent(),  partial.(dual_args)...)
            args = (op, primal.(dual_args)...)
            y, ẏ = frule(ȧrgs, args...; kwargs...)
            return Dual(y, ẏ)  # if y, ẏ are not `Float64` this will error.
        end
    end
    eval(fdef)
end

# !Important!: Attach the define function to the `on_new_rule` hook
on_new_rule(define_dual_overload, frule)

"Do a calculus. `f` should have a single input."
function derv(f, arg)
    duals = Dual(arg, one(arg))
    return partial(f(duals...))
end

# End AD definition
################################

# add a rule later also
function ChainRulesCore.frule((_, Δx, Δy), ::typeof(*), x::Number, y::Number)
    return (x * y, Δx * y + x * Δy)
end

# Manual refresh needed as new rule added in same file as AD after the `on_new_rule` call
refresh_rules();

@testset "ForwardDiffZero" begin
    foo(x) = x + x
    @test derv(foo, 1.6) == 2

    bar(x) = x + 2.1 * x
    @test derv(bar, 1.2) == 3.1

    baz(x) = 2.0 * x^2 + 3.0*x + 1.2
    @test derv(baz, 1.7) == 2*2.0*1.7 + 3.0

    qux(x) = foo(x) + bar(x) + baz(x)
    @test derv(qux, 1.7) == (2*2.0*1.7 + 3.0) + 3.1 + 2

    function quux(x)
        y = 2.0*x + 3.0*x
        return 4.0*y + 5.0*y
    end
    @test derv(quux, 11.1) == 4*(2+3) + 5*(2+3)
end

end  # module


text:
import ObjFileBase: Section, DebugSections, endianness

immutable DIETreeRef
    oh
    strtab
    tree
end

function show(io::IO, ref::DIETreeRef)
    show(io::IO, ref.tree; strtab = ref.strtab)
end

function dies(oh::ObjectHandle)
    dbgs = debugsections(oh)
    if dbgs.debug_str !== nothing
        return DIETreeRef(oh, load_strtab(dbgs.debug_str), read(oh, DIETree; dbgs = dbgs))
    else
        return read(oh, DIETree; dbgs = dbgs)
    end
end

function read{T<:ObjectHandle}(oh::T,::Type{DWARF.DIETree}; dbgs = debugsections(oh))
    seek(dbgs.oh, ObjFileBase.sectionoffset(dbgs.debug_info))
    s = read(dbgs.oh, DWARF.DWARFCUHeader)
    pos = position(dbgs.oh)
    abbrev = read(dbgs.oh, ObjFileBase.deref(dbgs.debug_abbrev), s, DWARF.AbbrevTableSet)
    seek(dbgs.oh, pos)
    DIETree = read(oh, s, abbrev, zero(DWARF.DIETree), DWARF.DIETreeNode, endianness(oh));
    DIETree
end

function read{T<:ObjectHandle,S}(oh::T,sec::Section{S},::Type{DWARF.ARTable})
    @assert T <: S
    seek(oh, sec)
    ret = DWARF.ARTable(Array(DWARF.ARTableSet,0))
    while position(oh) < sectionoffset(sec) + sectionsize(sec)
        push!(ret.sets, read(oh, DWARF.ARTableSet, endianness(oh)))
    end
    ret
end

function read{T<:ObjectHandle,S}(oh::T,sec::Section{S},::Type{DWARF.PUBTable})
    @assert T <: S
    seek(oh, sec)
    ret = DWARF.PUBTable(Array(DWARF.PUBTableSet,0))
    while position(oh) < sectionoffset(sec) + sectionsize(sec)
        push!(ret.sets,read(oh,DWARF.PUBTableSet, endianness(oh)))
    end
    ret
end

function read{T<:ObjectHandle,S}(oh::T, sec::Section{S},
        ::Type{DWARF.AbbrevTableSet})
    @assert T <: S
    seek(oh, sec)
    read(oh, AbbrevTableSet, endianness(oh))
end

function read{T<:ObjectHandle,S}(oh::T,debug_info::Section{S},
        s::DWARF.PUBTableSet,::Type{DWARF.DWARFCUHeader})
    @assert T <: S
    seek(oh,sectionoffset(debug_info)+s.header.debug_info_offset)
    read(oh,DWARF.DWARFCUHeader, endianness(oh))
end

function read{T<:ObjectHandle,S}(oh::T, debug_info::Section{S}, s::DWARF.DWARFCUHeader,
        ::Type{DWARF.AbbrevTableSet})
    @assert T <: S
    seek(oh,sectionoffset(debug_info)+s.debug_abbrev_offset)
    read(oh,DWARF.AbbrevTableSet, endianness(oh))
end

function read{T<:ObjectHandle,S}(oh::T,
    debug_info::Section{S}, debug_abbrev::Section{S},
    s::DWARF.PUBTableSet, e::DWARF.PUBTableEntry,
    header::DWARF.DWARFCUHeader, ::Type{DWARF.DIETree})

    @assert T <: S
    ats = read(oh,debug_abbrev,header,DWARF.AbbrevTableSet)
    seek(oh,sectionoffset(debug_info)+s.header.debug_info_offset+e.offset)
    ret = DWARF.DIETree(Array(DWARF.DIETreeNode,0))
    read(oh,header,ats,ret,DWARF.DIETreeNode,:LittleEndian)
    ret
end

function findindexbyname(x::DebugSections,name;
        pubtable=read(x.oh, deref(x.debug_pubnames), DWARF.PUBTable))
    symbols = map(y->map(x->x.name,y.entries),pubtable.sets)
    for i = 1:length(symbols)
        ind = findfirst(symbols[i],name)
        if ind != 0
            return (i,ind)
        end
    end
    return (0,0)
end

function findcubyname(x::DebugSections, name;
    pubtable = read(x.oh, deref(x.debug_pubnames), DWARF.PUBTable))
    (si,ei) = findindexbyname(x, name; pubtable = pubtable)
    if si == ei == 0
        error("Not Found")
    end
    read(x.oh,deref(x.debug_info),pubtable.sets[si],DWARF.DWARFCUHeader)
end

function finddietreebyname(x::DebugSections, name;
    pubtable = read(x.oh, deref(x.debug_pubnames), DWARF.PUBTable))
    (si,ei) = findindexbyname(x, name; pubtable = pubtable)
    if si == ei == 0
        error("Not Found")
    end
    pubset = pubtable.sets[si]
    pubentry = pubset.entries[ei]
    cu = read(x.oh,deref(x.debug_info),pubset,DWARF.DWARFCUHeader)
    d = read(x.oh,deref(x.debug_info),deref(x.debug_abbrev),pubset,pubentry,
        cu,DWARF.DIETree)
end

function read(x::DebugSections, ::Type{DWARF.ARTableSet})
    read(x.oh, deref(x.debug_aranges), DWARF.ARTableSet)
end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
#Very simple visualization output.
#Uses template from http://www.paraview.org/Wiki/ParaView/Data_formats#Displaying_data_as_structured_grid
function tocsv(xs,ys,zs,fn,filename)
    f=open(filename,"w");
    print(f,"x coord,y coord,z coord,field\n");
    for x in xs
        for y in ys
            for z in zs
                of=map(abs,fn(x,y,z));
                out=hypot(hypot(of[1],of[2]),of[3]);
                print(f,"$(x), $(y), $(z), $(out)\n");
            end
        end
    end
    close(f);
end

max_stars_repo_name: ReidAtcheson/incns

max_issues_repo_name: ReidAtcheson/incns

max_forks_repo_name: ReidAtcheson/incns


text:
null_theme = Theme(
              minor_label_font_size=0pt, #Size of ticks
              major_label_font_size=0pt, #Size of axis text
              line_width=0pt, #Plot line
              panel_fill=Colors.RGB(.95,.95,.95),
              background_color=Colors.RGB(.95,.95,.95)
              )
default_theme = Theme(
              minor_label_font_size=16pt, #Size of ticks
              major_label_font_size=26pt, #Size of axis text
              line_width=5pt, #Plot line
              panel_fill=Colors.RGB(.99,.99,.99),
              background_color=Colors.RGB(.99,.99,.99)
              )

#Almost empty plot
na_plot = Immerse.plot(x=[0.0], y=[0.0], Geom.line, null_theme)

function Base.push!(obj::GtkInteract.CairoGraphic, p::Void)
    push!(obj, na_plot)
end

function Base.push!(obj::GtkInteract.CairoGraphic, p::Gadfly.Plot)
    if obj.obj != nothing
      push!(obj.obj, p)
    end
end

function get_log_ticks(x)
    min = ceil(log10(minimum(x)))
    max = floor(log10(maximum(x)))
    major = 10.^collect(min:max)
    majorText = ["\$10^{$(round(Int64,i))}\$" for i=min:max]
    if length(major) < 7
      minor = [j*10^i for i=(min-1):(max+1) for j=2:9]
      minor = minor[find(minimum(x) .<= minor .<= maximum(x))]
    else
      minor = Float64[]
    end
    if length(major) < 2
      minorText = ["\$j*10^\{$i\}\$" for i=(min-1):(max+1) for j=2:9]
    else
      minorText = fill("", length(minor))
    end
    major, minor, majorText, minorText
end

function labelf(x)
  if abs(round(Int64, x) - x) < 1e-5
    @sprintf("10<sup>%1.0f</sup>", x)
  else
    ""
  end
end

function major_lines(x, y, xticks, yticks)
  layer(
    x=x, y=y,
    yintercept=yticks, xintercept=xticks, Geom.nil,
    Geom.hline, Geom.vline, Theme(line_width=0.2pt, default_color=Colors.RGB(.1,.1,.1))
  )
end


text:
using BrainAnnotationMapping
using Test

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


text:
abstract type AbstractListofStates end

struct VoidListofStates <: AbstractListofStates end

"""
Type: list of States

Constructor:

`ListofStates(:: AbstractState)`

`ListofStates(n :: Int, :: Val{AbstractState})`

`ListofStates(n :: Int, list :: Array{AbstractState,1})`

`ListofStates(state :: S; n :: Int = -1, kwargs...)`

Note:
- If `n != -1`, then it stores at most n `AbstractState`.
- `ListofStates` recursively handles sub-list of states as the attribute list is
an array of pair whose first component is a, `AbstractState` and the second
component is a `ListofStates` (or `VoidListofStates`).

Examples:
`ListofStates(state)`
`ListofStates(state, n = 2)`
`ListofStates(-1, Val{NLPAtX}())`
`ListofStates(-1, [(state1, VoidListofStates), (state2, VoidListofStates)], 2)`
"""
mutable struct ListofStates{S <: AbstractState, T} <: AbstractListofStates

  n     :: Int #If length of the list is knwon, -1 if unknown
  i     :: Int #current index in the list/length
  list  :: Array{Tuple{S, T},1}
                 #Tanj: \TODO Tuple instead of an Array would be better, I think

end

state_type(:: ListofStates{S, T}) where {S <: AbstractState, T} = S

function ListofStates(n :: T, :: Val{S}) where {T <: Int, S <: AbstractState}
  list = Array{Tuple{S, VoidListofStates},1}(undef, 0)
  i = 0
  return ListofStates(n, i, list)
end

function ListofStates(n :: Ti, list :: Array{S,1}) where {S <: AbstractState, Ti <: Int}
  i = length(list)
  tuple_list = Array{Tuple{S, VoidListofStates},1}(undef, 0)
  for j=1:i
    push!(tuple_list, (list[j], VoidListofStates()))
  end
  return ListofStates(n, i, tuple_list)
end

function ListofStates(n :: Ti, list :: Array{Tuple{S, T},1}) where {S <: AbstractState, T, Ti <: Int}
  i   = length(list)
  return ListofStates(n, i, list)
end

function ListofStates(state :: S; n :: Int = -1, kwargs...) where S <: AbstractState
  i =  1
  list = [(copy_compress_state(state; kwargs...), VoidListofStates())]
  return ListofStates(n, i, list)
end

"""
add\\_to\\_list!: add a State to the list of maximal size n.
If a n+1-th State is added, the first one in the list is removed.
The given is State is compressed before being added in the list (via State.copy\\_compress\\_state).

`add_to_list!(:: AbstractListofStates, :: AbstractState; kwargs...)`

Note:
 -  kwargs are passed to the compress_state call.
 -  does nothing for `VoidListofStates`

see also: ListofStates, State.compress\\_state, State.copy\\_compress\\_state
"""
function add_to_list!(list :: AbstractListofStates, state :: AbstractState; kwargs...)

  if typeof(list.n) <: Int && list.n > 0 #If n is a natural number
    if list.i + 1 > list.n
      popfirst!(list.list) #remove the first item
      list.i = list.n
    else
      list.i += 1
    end
    cstate = copy_compress_state(state; kwargs...)
    push!(list.list, (cstate, VoidListofStates()))
  else
    push!(list.list, (copy_compress_state(state; kwargs...), VoidListofStates()))
    list.i += 1
  end

  return list
end

function add_to_list!(list :: VoidListofStates, state :: AbstractState; kwargs...)
  return list
end

import Base.length
"""
length: return the number of States in the list.

`length(:: ListofStates)`

see also: print, add_to_list!, ListofStates
"""
function length(list :: AbstractListofStates)
  return list.i
end

import Base.print
"""
print: output formatting. return a DataFrame.

`print(:: ListofStates; verbose :: Bool = true, print_sym :: Union{Nothing,Array{Symbol,1}})`

Note:
- set `verbose` to false to avoid printing.
- if `print_sym` is an Array of Symbol, only those symbols are printed. Note that
the returned DataFrame still contains all the columns.
- More information about DataFrame: http://juliadata.github.io/DataFrames.jl

see also: add\\_to\\_list!, length, ListofStates
"""
function print(list :: AbstractListofStates;
               verbose :: Bool = true,
               print_sym :: Union{Nothing,Array{Symbol,1}} = nothing)

  tab = zeros(0, length(list.list))#Array{Any,2}(undef, length(fieldnames(typeof(list.list[1,1]))))
  for k in fieldnames(typeof(list.list[1,1]))
    tab = vcat(tab, [getfield(i[1], k) for i in list.list]');
  end
  df = DataFrame(tab, :auto)

  if isnothing(print_sym)
    verbose && print(df)
  else
    verbose && print(df[!, print_sym])
  end

  return df
end

import Base.getindex
"""
`getindex(:: ListofStates, :: Int)`
`getindex(:: ListofStates, :: Int, :: Int)`

Example:
stop_lstt.listofstates.list[3]
stop_lstt.listofstates.list[3,1]
"""
function getindex(list :: AbstractListofStates, i :: Int)
  return list.list[i]
end

function getindex(list :: AbstractListofStates, i :: Int, j :: Int)
  return list.list[i][j]
end


text:
function diagonal_elements(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  A = randn(n,n)
  B = randn(n,n)

  Benchmarker.add_data(csv, "diag_mult_nai", Benchmarker.measure(reps, implicit_diagonal_mult, A, B))
  Benchmarker.add_data(csv, "diag_add_nai", Benchmarker.measure(reps, implicit_diagonal, A, B))
  Benchmarker.add_data(csv, "diag_add_rec", Benchmarker.measure(reps, explicit_diagonal, A, B))

  @test isapprox(implicit_diagonal(A, B)[2], explicit_diagonal(A, B)[2])
end

function implicit_diagonal_mult(A, B)
  result = zeros(1, size(A,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = diag(A*B)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function implicit_diagonal(A, B)
  result = zeros(1, size(A,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = diag(A+B)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function explicit_diagonal(A, B)
  result = zeros(1, size(A,1))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = diag(A) + diag(B)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end


text:
function test6()
  filename = "test6.asc"
  exc = ""
  # exectuablepath = null string will not run LTspice.exe.  Test parsing only.
  sim = LTspiceSimulation(filename,executablepath="")
  @test LTspice.does_circuitfilearray_file_match(sim)
  show(IOBuffer(),sim)

  @test measurementnames(sim) == ()
  @test stepnames(sim) == ()
  @test logpath(sim) != ""

  @static if Sys.iswindows()
      @test circuitpath(sim) == filename
  end

  @test parametervalues(sim) == []
  @test measurementvalues(sim) == []
  @test executablepath(sim) == ""
  @test length(sim) == 0

  pli = perlineiterator(sim,header=true)
  for line in pli
    nothing
  end

  show(IOBuffer(),sim)
end
test6()


text:
# BigWig Writer
# =============

mutable struct WriterState
    # section info
    datatype::UInt8
    count::UInt64
    chromid::UInt32
    chromstart::UInt32
    chromend::UInt32
    itemstep::UInt32
    itemspan::UInt32

    # previous record info
    chromstart_prev::UInt32
    chromend_prev::UInt32

    # section state
    started::Bool
    buffer::IOBuffer
    compressed::Vector{UInt8}

    # global state and stats
    blocks::Vector{BBI.Block}
    cov::UInt32
    min::Float32
    max::Float32
    sum::Float32
    ssq::Float32

    function WriterState(datatype::Symbol)
        return new(
            # section info
            encode_datatype(datatype), 0, 0, 0, 0, 0, 0,
            # last record info
            0, 0,
            # section data
            false, IOBuffer(), UInt8[],
            # global info
            BBI.Block[],
            0, +Inf32, -Inf32, 0.0f0, 0.0f0)
    end
end

struct Writer <: BioCore.IO.AbstractWriter
    # output stream
    stream::IO

    # number of zoom levels
    zoomlevels::Int

    # maximum size of uncompressed buffer
    max_buffer_size::UInt64

    # file offsets
    summary_offset::UInt64
    chrom_tree_offset::UInt64
    data_offset::UInt64

    # chrom name => (ID, length)
    chroms::Dict{String,Tuple{UInt32,UInt32}}

    # chrom ID => length
    chromlens::Dict{UInt32,UInt32}

    # mutable state
    state::WriterState

    # zoom buffer
    zoombuffer::BBI.ZoomBuffer
end

const ZOOM_SCALE_FACTOR = 4

"""
    BigWig.Writer(output::IO, chromlist; binsize=64, datatype=:bedgraph)

Create a data writer of the bigWig file format.

Arguments
---------

* `output`: data sink
* `chromlist`: chromosome list with length
* `binsize=64`: size of a zoom with the highest resolution
* `datatype=:bedgraph`: encoding of values (`:bedgraph`, `:varstep` or `:fixedstep`)

Examples
--------

```julia
output = open("data.bw", "w")
writer = BigWig.Writer(output, [("chr1", 12345), ("chr2", 9100)])
write(writer, ("chr1", 501, 600, 1.0))
write(writer, ("chr2", 301, 450, 3.0))
close(writer)
```
"""
function Writer(output::IO, chromlist::Union{AbstractVector,AbstractDict}; binsize::Integer=64, datatype::Symbol=:bedgraph)
    # write dummy header (filled later)
    write_zeros(output, BBI.HEADER_SIZE)

    # write dummy zoom headers (filled later)
    chromlist_with_id = BBI.add_chrom_ids(chromlist)
    maxlen = Base.maximum(x[3] for x in chromlist_with_id)
    zoomlevels = BBI.determine_zoomlevels(maxlen, binsize, ZOOM_SCALE_FACTOR)
    write_zeros(output, BBI.ZOOM_HEADER_SIZE * zoomlevels)

    # write dummy total summary (filled later)
    summary_offset = position(output)
    write_zeros(output, BBI.SUMMARY_SIZE)

    # write chromosome B+-tree
    chrom_tree_offset = position(output)
    BBI.write_btree(output, chromlist_with_id)

    # write dummy data count (filled later)
    data_offset = position(output)
    write_zeros(output, sizeof(UInt64))

    # initialize zoom buffer (use temporary file?)
    max_block_size = 16 * 2^10
    chromlens = Dict(id => len for (name, id, len) in chromlist_with_id)
    zoombuffer = BBI.ZoomBuffer(chromlens, binsize, max_block_size)

    return Writer(
        output,
        zoomlevels,
        max_block_size,
        summary_offset,
        chrom_tree_offset,
        data_offset,
        Dict(name => (id, len) for (name, id, len) in chromlist_with_id),
        chromlens,
        WriterState(datatype),
        zoombuffer)
end

function Base.write(writer::Writer, record::Tuple{String,Integer,Integer,Real})
    chromname, chromstart, chromend, value = record
    chromid = writer.chroms[chromname][1]
    return write_impl(writer, chromid, UInt32(chromstart - 1), UInt32(chromend), Float32(value))
end

function Base.write(writer::Writer, interval::GenomicFeatures.Interval{T}) where T<:Real
    chromid = writer.chroms[seqname(interval)][1]
    return write_impl(writer, chromid, UInt32(leftposition(interval) - 1), UInt32(rightposition(interval)), Float32(metadata(interval)))
end

function finalize_file(writer::Writer)
    state = writer.state
    if state.started
        finish_section!(writer)
    end

    # write R-tree
    stream = writer.stream
    data_index_offset = position(stream)
    BBI.write_rtree(writer.stream, state.blocks)

    # fill header
    header = BBI.Header(
        BBI.WIG_MAGIC,
        3,  # version
        writer.zoomlevels,
        writer.chrom_tree_offset,
        writer.data_offset,
        data_index_offset,
        0,  # field count (0 for bigWig)
        0,  # predefined field count (0 for bigWig?)
        0,  # autoSQL offset (0 for bigWig?)
        writer.summary_offset,
        writer.max_buffer_size,
        0   # reserved
    )
    seekstart(stream)
    write(stream, header)

    # fill summary
    seek(stream, writer.summary_offset)
    write(stream, BBI.Summary(state.cov, state.min, state.max, state.sum, state.ssq))

    # fill data count
    seek(stream, writer.data_offset)
    write(stream, UInt64(length(state.blocks)))

    # write zoom
    seekend(stream)
    zoomheaders = BBI.write_zoom(stream, writer.zoombuffer, writer.zoomlevels, ZOOM_SCALE_FACTOR)
    @assert length(zoomheaders) == writer.zoomlevels

    # fill zoom headers
    seek(stream, BBI.HEADER_SIZE)
    for zheader in zoomheaders
        write(stream, zheader)
    end
end

function Base.close(writer::Writer)
    finalize_file(writer)
    close(writer.stream)
    return
end

function BioCore.IO.stream(writer::Writer)
    return writer.stream
end

function write_impl(writer::Writer, chromid::UInt32, chromstart::UInt32, chromend::UInt32, value::Float32)
    n = 0
    state = writer.state
    if state.started && (chromid != state.chromid || position(state.buffer) ≥ writer.max_buffer_size - sizeof(UInt32) * 3)
        finish_section!(writer)
    end
    if !state.started
        n += start_section!(writer, chromid, chromstart, chromend)
    end

    # check new record
    check_interval(state, chromid, chromstart, chromend)

    # infer the step size from the first and second records
    if isfixedstep(state.datatype) && state.count == 1
        state.itemstep = chromstart - state.chromstart_prev
    end

    # write data to the buffer
    if isbedgraph(state.datatype)
        n += write(state.buffer, chromstart, chromend, value)
    elseif isvarstep(state.datatype)
        n += write(state.buffer, chromstart, value)
    elseif isfixedstep(state.datatype)
        n += write(state.buffer, value)
    else
        assert(false)
    end

    # update state
    state.count += 1
    state.chromend = max(state.chromend, chromend)
    state.chromstart_prev = chromstart
    state.chromend_prev = chromend

    size = chromend - chromstart
    state.cov += size
    state.min = min(state.min, value)
    state.max = max(state.max, value)
    state.sum += value   * size
    state.ssq += value^2 * size

    BBI.add_value!(writer.zoombuffer, chromid, chromstart, chromend, value)

    return n
end

function check_interval(state::WriterState, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    if chromstart ≥ chromend
        throw(ArgumentError("empty interval"))
    end
    if BBI.compare_intervals((chromid, chromstart, chromend), (state.chromid, state.chromstart_prev, state.chromend_prev)) != 1
        throw(ArgumentError("disordered intervals"))
    end
    if isvarstep(state.datatype) && chromend - chromstart != state.itemspan
        throw(ArgumentError("inconsistent interval span"))
    end
    if isfixedstep(state.datatype)
        if chromend - chromstart != state.itemspan
            throw(ArgumentError("inconsistent interval span"))
        end
        if state.count > 1 && chromstart - state.chromstart_prev != state.itemstep
            throw(ArgumentError("inconsistent intreval step"))
        end
    end
end

function start_section!(writer::Writer, chromid::UInt32, chromstart::UInt32, chromend::UInt32)
    state = writer.state
    @assert !state.started

    # initialize section info
    state.chromid = chromid
    state.chromstart = chromstart
    state.chromend = chromend
    if isfixedstep(state.datatype) || isvarstep(state.datatype)
        # inter the span from the first record
        state.itemspan = chromend - chromstart
    end
    state.count = 0

    # initialize last record info
    state.chromstart_prev = 0
    state.chromend_prev = 0

    # write dummy section header (filled later)
    seekstart(state.buffer)
    truncate(state.buffer, SECTION_HEADER_SIZE)
    seek(state.buffer, SECTION_HEADER_SIZE)
    resize!(state.compressed, div(writer.max_buffer_size * 11, 10))
    state.started = true
    return SECTION_HEADER_SIZE
end

function finish_section!(writer::Writer)
    state = writer.state

    # fill section header
    seekstart(state.buffer)
    write(
        state.buffer,
        SectionHeader(
            state.chromid,
            state.chromstart,
            state.chromend,
            state.itemstep,
            state.itemspan,
            state.datatype,
            0x00,  # reserved
            state.count))

    # write compressed section
    size = BBI.compress!(state.compressed, take!(state.buffer))
    offset = position(writer.stream)
    unsafe_write(writer.stream, pointer(state.compressed), size)

    # record block
    push!(
        state.blocks,
        BBI.Block((state.chromid, state.chromstart), (state.chromid, state.chromend), offset, size))
    state.started = false
    return
end

function write_zeros(stream::IO, n::Integer)
    for _ in 1:n
        write(stream, 0x00)
    end
    return n
end


text:
using SSAGES
using Test

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


text:
using Plots
using Oceananigans

using Statistics

"""
    make_vertical_slice_movie(model::Model, nc_writer::NetCDFOutputWriter,
                              var_name, Nt, Δt, var_offset=0, slice_idx=1)

Make a movie of a vertical slice produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀). A `slice_idx` can be specified to select
the index of the y-slice to be plotted (useful when plotting vertical slices
from a 3D model, it should be set to 1 for 2D xz-slice models).
"""
function make_vertical_slice_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0, slice_idx=1)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.heatmap(model.grid.xC, model.grid.zC, rotl90(var[:, slice_idx, :]) .- var_offset,
                      color=:balance, clims=(-0.01, 0.01), title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

"""
    make_horizontal_slice_movie(model::Model, nc_writer::NetCDFOutputWriter,
                                var_name, Nt, Δt, var_offset=0)

Make a movie of a horizontal slice produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀).
"""
function make_horizontal_slice_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.heatmap(model.grid.xC, model.grid.yC, var[:, :, 1] .- var_offset,
                      color=:balance, clims=(-0.01, 0.01),
                      title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

"""
    make_vertical_profile_movie(model::Model, nc_writer::NetCDFOutputWriter,
                                var_name, Nt, Δt, var_offset=0)

Make a movie of a vertical profile produced by `model` with output being saved by
`nc_writer`. The variable name `var_name` can be either of "u", "v", "w", "T",
or "S". `Nt` is the number of model iterations (or time steps) taken and ``Δt`
is the time step. A plotting offset `var_offset` can be specified to be
subtracted from the data before plotting (useful for plotting e.g. small
temperature perturbations around T₀).
"""
function make_vertical_profile_movie(model::Model, nc_writer::NetCDFOutputWriter, var_name, Nt, Δt, var_offset=0)
    freq = nc_writer.output_frequency
    N_frames = Int(Nt/freq)

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")
        var = read_output(nc_writer, var_name, freq*n)
        Plots.plot(var[1, 1, :] .- var_offset, model.grid.zC,
                   title="t=$(freq*n*Δt) s ($(round(freq*n*Δt/86400; digits=2)) days)")
    end

    mp4(animation, nc_writer.filename_prefix * "$(round(Int, time())).mp4", fps=30)
end

using NetCDF

function make_avg_temperature_profile_movie()
    Nt, dt = 86400, 0.5
    freq = 3600
    N_frames = Int(Nt/freq)
    filename_prefix = "convection"
    var_offset = 273.15

    Nz, Lz = 128, 100
    dz = Lz/Nz
    zC = -dz/2:-dz:-Lz

    print("Producing movie... ($N_frames frames)\n")
    Plots.gr(dpi=150)

    animation = @animate for n in 0:N_frames
        print("\rframe = $n / $N_frames   ")

        filepath = filename_prefix * lpad(freq*n, 9, "0") * ".nc"
        field_data = ncread(filepath, "T")
        ncclose(filepath)

        T_profile = mean(field_data; dims=[1,2])

        Plots.plot(reshape(T_profile, Nz) .- var_offset, zC,
                   title="t=$(freq*n*dt) s ($(round(freq*n*dt/86400; digits=2)) days)")
    end

    mp4(animation, filename_prefix * "$(round(Int, time())).mp4", fps=30)
end


text:
module ReinterpretTests

using Test
using InplaceArrays.TensorValues

v = VectorValue(1.3,2.1)

V = fill(v,5)

A = reinterpret(V)

R = [1.3 1.3 1.3 1.3 1.3; 2.1 2.1 2.1 2.1 2.1]

@test A == R

end # module ReinterpretTests


text:
using Test
using PyBaMM
using SparseArrays, LinearAlgebra
using Sundials
using PyCall
using OrderedCollections

pybamm = pyimport("pybamm")

@testset "SPMe L2 loss voltage" begin
    # load model
    model = pybamm.lithium_ion.SPMe(name="SPMe")
    parameter_values = model.default_parameter_values
    parameter_values.update(
        Dict(
            "Cation transference number" => pybamm.InputParameter("t_plus"),
            "Electrolyte conductivity [S.m-1]" => pybamm.InputParameter("kappa_e"),
        )
    )
    sim = pybamm.Simulation(model, parameter_values=parameter_values)

    inputs = OrderedDict{String,Float64}(
        "t_plus" => 0.4,
        "kappa_e" => 0.7,
    )

    # Build ODE problem
    prob = get_ode_problem(sim, inputs)
    t = prob.tspan[2] / 100
    # Generate data
    sol = solve(prob, CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1), reltol=1e-6, abstol=1e-6, saveat=t);
    V_data = get_variable(sim, sol, "Terminal voltage [V]", inputs)

    # Get loss function
    loss = get_l2loss_function(sim, "Terminal voltage [V]", inputs, V_data)

    # Loss of solution used to generate data should be zero
    @test loss(sol) == 0.0

    # Loss should get bigger as parameters get further away from truth
    prob2 = remake(prob; p=[0.5,0.7])
    sol2 = solve(prob2, CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1), reltol=1e-6, abstol=1e-6, saveat=t);

    prob3 = remake(prob; p=[0.6,0.7])
    sol3 = solve(prob3, CVODE_BDF(linear_solver=:Band,jac_upper=1,jac_lower=1), reltol=1e-6, abstol=1e-6, saveat=t);
    @test loss(sol2) < loss(sol3)
end


text:
using FastHistograms
using FastHistograms: BinType, FixedWidth, VariableWidth, UnboundedWidth
using FastHistograms: BinSearchAlgorithm, Arithmetic, BinarySearch, HashFunction
using FastHistograms: HistogramParallelization, NoParallelization, SIMD, PrivateThreads
using Test
using Random
import StatsBase
using BenchmarkTools
using Downloads

const run_1D_tests = true
const run_2D_tests = true

function get_test_resources()
    path1 = joinpath(@__DIR__, "pride_and_prejudice.txt")
    if !isfile(path1)
        Downloads.download("https://www.gutenberg.org/files/1342/1342-0.txt", path1)
    end
end

get_test_resources()

include("test_fixed_width.jl")
include("test_variable_width.jl")
include("test_text_hist.jl")

invalid_combination(::BinType, ::BinSearchAlgorithm, ::HistogramParallelization) = false
# SIMD binary search needs a special implementation that I don't have yet
invalid_combination(::BinType, ::BinarySearch, ::SIMD) = true
# Arithmetic only works for fixed width
invalid_combination(::Union{VariableWidth,UnboundedWidth}, ::Arithmetic, ::HistogramParallelization) = true
# Text histograms can only use UnboundedWidth, HashFunction, and NoParallelization or PrivateThreads. SIMD not implemented yet.
invalid_combination(::UnboundedWidth, ::HashFunction, ::SIMD) = true

@testset "FastHistograms.jl" begin
    bin_types = [FixedWidth(), VariableWidth()]
    search_algorithms = [Arithmetic(), BinarySearch()]
    parallelizations = [NoParallelization(), SIMD(), PrivateThreads()]

    @testset "histogram computations BinType=$(bin_type)" for bin_type in bin_types
        @testset "BinSearchAlgorithm=$(search_algorithm)" for search_algorithm in search_algorithms
            @testset "HistogramParallelization=$(parallelization)" for parallelization in parallelizations
                if invalid_combination(bin_type, search_algorithm, parallelization)
                    continue
                end

                test_parameterized_hist(bin_type, search_algorithm, parallelization)
            end
        end
    end

    @testset "zero bins" begin
        h = create_fast_histogram(FixedWidth(), Arithmetic(), NoParallelization(), [(0x00, 0xff, 16), (0x00, 0xff, 16)])

        increment_bins!(h, rand(UInt8, 10, 10), rand(UInt8, 10, 10))

        @test any(counts(h) .!= 0)

        zero!(h)
        @test all(counts(h) .== 0)
    end

    @testset "default hist bin type" begin
        h = create_fast_histogram(FixedWidth(), Arithmetic(), NoParallelization(), [(0x00, 0xff, 16), (0x00, 0xff, 16)])

        @test eltype(h) == UInt8
    end

    @testset "custom hist bin type" begin
        h = create_fast_histogram(
            FixedWidth(),
            Arithmetic(),
            NoParallelization(),
            [(Int32(0), Int32(16), 16), (Int32(0), Int32(16), 16)],
        )

        @test eltype(h) == Int32
    end

    @testset "AHTL fixed data" begin
        # Generate the data in the same way the AHTL test does
        data = rand(Float32, 100000000) .* (512 - Float32(0.01))

        img_vec = vec(data)
        sb_counts = StatsBase.fit(
            StatsBase.Histogram,
            img_vec,
            range(Float32(0.0); stop = Float32(512.0), length = 513),
        ).weights

        h = create_fast_histogram(
            FastHistograms.FixedWidth(),
            FastHistograms.Arithmetic(),
            FastHistograms.PrivateThreads(),
            [(Float32(0.0), Float32(512.0), 512)],
        )
        increment_bins!(h, data)
        fh_counts = counts(h)
        @test fh_counts == sb_counts
    end

    let bin_type = UnboundedWidth(), search_algorithm = HashFunction()
        @testset "HistogramParallelization=$(parallelization)" for parallelization in parallelizations
            if invalid_combination(bin_type, search_algorithm, parallelization)
                continue
            end

            test_text_hist(parallelization)
        end
    end
end


text:
module Guide

using Color
using Compose
using Gadfly
using JSON

import Gadfly.render, Gadfly.escape_id

# Where the guide should be placed in relation to the plot.
abstract GuidePosition
type TopGuidePosition    <: GuidePosition end
type RightGuidePosition  <: GuidePosition end
type BottomGuidePosition <: GuidePosition end
type LeftGuidePosition   <: GuidePosition end
type UnderGuidePosition  <: GuidePosition end

const top_guide_position    = TopGuidePosition()
const right_guide_position  = RightGuidePosition()
const bottom_guide_position = BottomGuidePosition()
const left_guide_position   = LeftGuidePosition()
const under_guide_position  = UnderGuidePosition()

type PanelBackground <: Gadfly.GuideElement
end

const background = PanelBackground()

function render(guide::PanelBackground, theme::Gadfly.Theme,
                aess::Vector{Gadfly.Aesthetics})
    back = compose(canvas(order=-1),
                   rectangle(),
                   svgclass("guide background"),
                   stroke(theme.panel_stroke),
                   fill(theme.panel_fill))

    {(back, under_guide_position)}
end

type ColorKey <: Gadfly.GuideElement
end

const colorkey = ColorKey()

# A helper for render(::ColorKey) for rendering guides for discrete color
# scales.
function render_discrete_color_key(colors::Vector{ColorValue},
                                   labels::Dict{ColorValue, String},
                                   title_canvas::Canvas,
                                   title_width::Measure,
                                   theme::Gadfly.Theme)
    # Key entries
    n = length(colors)

    entry_width, entry_height = text_extents(theme.minor_label_font,
                                             theme.minor_label_font_size,
                                             values(labels)...)
    entry_width += entry_height # make space for the color swatch

    # Rewrite to put toggleable things in a group.
    swatch_padding = 1mm
    swatch_size = 1cy - swatch_padding
    swatch_canvas = canvas(0w, 0h + title_canvas.box.height,
                           1w, n * (entry_height + swatch_padding),
                           unit_box=Units(0, 0, 1, n))
    for (i, c) in enumerate(colors)
        swatch_square = compose(rectangle(0, i - 1, swatch_size, swatch_size),
                                fill(c),
                                stroke(theme.highlight_color(c)),
                                linewidth(theme.highlight_width))

        label = labels[c]
        swatch_label = compose(text(1cy, (i - 1)cy + entry_height/2,
                                    label, hleft, vcenter),
                               stroke(nothing),
                               fill(theme.minor_label_color))
        swatch = swatch_square | swatch_label

        color_class = @sprintf("color_%s", escape_id(label))
        swatch <<= svgclass(@sprintf("guide %s", color_class))

        swatch <<= d3embed(
            @sprintf(".on(\"click\", guide_toggle_color(parent_id, \"%s\"))",
                     color_class))

        swatch_canvas <<= swatch
    end
    swatch_canvas <<= font(theme.minor_label_font) |
                      fontsize(theme.minor_label_font_size)
    c = canvas(0, 0, max(title_width, entry_width) + 3swatch_padding,
               swatch_canvas.box.height + title_canvas.box.height) <<
        pad(canvas() << swatch_canvas << title_canvas, 2mm)
    c
end

# A helper for render(::ColorKey) for rendering guides for continuous color
# scales.
function render_continuous_color_key(colors::Vector{ColorValue},
                                   labels::Dict{ColorValue, String},
                                   title_canvas::Canvas,
                                   title_width::Measure,
                                   theme::Gadfly.Theme)

    # Key entries
    entry_width, entry_height = text_extents(theme.minor_label_font,
                                             theme.minor_label_font_size,
                                             values(labels)...)
    entry_width += entry_height # make space for the color swatch

    unlabeled_swatches = 0
    for c in colors
        if labels[c] == ""
            unlabeled_swatches += 1
        end
    end

    unlabeled_swatch_height = 1.0mm
    swatch_padding = 1mm

    swatch_canvas = canvas(0w, 0h + title_canvas.box.height, 1w,
                           unlabeled_swatches * unlabeled_swatch_height +
                           (length(colors) - unlabeled_swatches) * entry_height)

    # Nudge things to overlap slightly avoiding any gaps.
    nudge = 0.1mm

    y = 0cy
    for (i, c) in enumerate(colors)
        if labels[c] == ""
            swatch_square = compose(rectangle(0, y,
                                              entry_height,
                                              unlabeled_swatch_height + nudge),
                                    fill(c),
                                    linewidth(theme.highlight_width))

            swatch_canvas <<= swatch_square

            y += unlabeled_swatch_height
        else
            swatch_square = compose(rectangle(0, y,
                                              entry_height,
                                              entry_height + nudge),
                                    fill(c),
                                    linewidth(theme.highlight_width))
            swatch_label = compose(text(entry_height + swatch_padding,
                                        y + entry_height / 2,
                                        labels[c],
                                        hleft, vcenter),
                                   fill(theme.minor_label_color))

            swatch_canvas <<= swatch_square | swatch_label

            y += entry_height
        end
    end

    swatch_canvas <<= font(theme.minor_label_font) |
                      fontsize(theme.minor_label_font_size) |
                      stroke(nothing)

    c = canvas(0, 0, max(title_width, entry_width) + 3swatch_padding,
               swatch_canvas.box.height + title_canvas.box.height) <<
        pad(canvas() << swatch_canvas << title_canvas, 2mm)
    c
end

function render(guide::ColorKey, theme::Gadfly.Theme,
                aess::Vector{Gadfly.Aesthetics})
    used_colors = Set{ColorValue}()
    colors = Array(ColorValue, 0) # to preserve ordering
    labels = Dict{ColorValue, Set{String}}()

    continuous_guide = false
    guide_title = nothing

    for aes in aess
        if guide_title === nothing && !is(aes.color_key_title, nothing)
            guide_title = aes.color_key_title
        end

        if aes.color_key_colors === nothing
            continue
        end

        if !is(aes.color_key_continuous, nothing) && aes.color_key_continuous
            continuous_guide = true
        end

        for color in aes.color_key_colors
            label = aes.color_label(color)
            if !contains(used_colors, color)
                add!(used_colors, color)
                push!(colors, color)
                labels[color] = Set{String}(label)
            else
                add!(labels[color], label)
            end
        end
    end

    if guide_title === nothing
        guide_title = ""
    end

    pretty_labels = Dict{ColorValue, String}()
    for (color, label) in labels
        pretty_labels[color] = join(labels[color], ", ")
    end

    # Key title
    title_width, title_height = text_extents(theme.major_label_font,
                                             theme.major_label_font_size,
                                             guide_title)

    title_padding = 2mm
    title_canvas = compose(canvas(0w, 0h, 1w, title_height + title_padding),
                           text(0.5w, title_height, guide_title, hcenter, vbottom),
                           stroke(nothing),
                           font(theme.major_label_font),
                           fontsize(theme.major_label_font_size),
                           fill(theme.major_label_color))

    if continuous_guide
        c = render_continuous_color_key(colors, pretty_labels, title_canvas,
                                        title_width, theme)
    else
        c = render_discrete_color_key(colors, pretty_labels, title_canvas,
                                      title_width, theme)
    end

    {(c, right_guide_position)}
end

type XTicks <: Gadfly.GuideElement
end

const x_ticks = XTicks()

function render(guide::XTicks, theme::Gadfly.Theme,
                aess::Vector{Gadfly.Aesthetics})

    ticks = Dict{Float64, String}()
    grids = Set()
    for aes in aess
        if Gadfly.issomething(aes.xtick)
            for val in aes.xtick
                ticks[val] = aes.xtick_label(val)
            end
        end

        if Gadfly.issomething(aes.xgrid)
            for val in aes.xgrid
                add!(grids, val)
            end
        end
    end

    # grid lines
    grid_lines = compose(canvas(),
                         [lines((t, 0h), (t, 1h)) for t in grids]...,
                         stroke(theme.grid_color),
                         linewidth(theme.grid_line_width),
                         svgclass("guide xgridlines"))

    # tick labels

    (_, height) = text_extents(theme.minor_label_font,
                               theme.minor_label_font_size,
                               values(ticks)...)
    #padding = 1mm
    padding = 0mm

    tick_labels = compose(canvas(0, 0, 1w, height + 2padding, order=-1),
                          [text(tick, 1h - padding, label, hcenter, vbottom)
                           for (tick, label) in ticks]...,
                          stroke(nothing),
                          fill(theme.minor_label_color),
                          font(theme.minor_label_font),
                          fontsize(theme.minor_label_font_size),
                          svgclass("guide xlabels"))

    {(grid_lines, under_guide_position),
     (tick_labels, bottom_guide_position)}
end

type YTicks <: Gadfly.GuideElement
end

const y_ticks = YTicks()

function render(guide::YTicks, theme::Gadfly.Theme,
                aess::Vector{Gadfly.Aesthetics})
    ticks = Dict{Float64, String}()
    grids = Set()
    for aes in aess
        if Gadfly.issomething(aes.ytick)
            for val in aes.ytick
                ticks[val] = aes.ytick_label(val)
            end
        end

        if Gadfly.issomething(aes.ygrid)
            for val in aes.ygrid
                add!(grids, val)
            end
        end
    end

    # grid lines
    grid_lines = compose(canvas(),
                         [lines((0w, t), (1w, t)) for t in grids]...,
                         stroke(theme.grid_color),
                         linewidth(theme.grid_line_width),
                         svgclass("guide ygridlines"))

    # tick labels
    (width, _) = text_extents(theme.minor_label_font,
                              theme.minor_label_font_size,
                              values(ticks)...)
    padding = 1mm
    width += 2padding

    tick_labels = compose(canvas(0, 0, width, 1cy, order=-1),
                          [text(width - padding, t, label, hright, vcenter)
                           for (t, label) in ticks]...,
                          stroke(nothing),
                          fill(theme.minor_label_color),
                          font(theme.minor_label_font),
                          fontsize(theme.minor_label_font_size),
                          svgclass("guide ylabels"))

    {(grid_lines, under_guide_position),
     (tick_labels, left_guide_position)}
end

# X-axis label Guide
type XLabel <: Gadfly.GuideElement
    label::String
end

function render(guide::XLabel, theme::Gadfly.Theme,
                aess::Vector{Gadfly.Aesthetics})
    (_, text_height) = text_extents(theme.major_label_font,
                                    theme.major_label_font_size,
                                    guide.label)

    padding = 2mm
    c = compose(canvas(0, 0, 1w, text_height + 2padding),
                text(0.5w, 1h - padding, guide.label, hcenter, vbottom),
                stroke(nothing),
                fill(theme.major_label_color),
                font(theme.major_label_font),
                fontsize(theme.major_label_font_size))

    {(c, bottom_guide_position)}
end

# Y-axis label Guide
type YLabel <: Gadfly.GuideElement
    label::String
end

function render(guide::YLabel, theme::Gadfly.Theme, aess::Vector{Gadfly.Aesthetics})
    (text_width, text_height) = text_extents(theme.major_label_font,
                                             theme.major_label_font_size,
                                             guide.label)
    padding = 2mm
    c = compose(canvas(0, 0, text_height + 2padding, 1cy,
                       rotation=Rotation(-0.5pi, 0.5w, 0.5h)),
                text(0.5w, 0.5h, guide.label, hcenter, vcenter),
                stroke(nothing),
                fill(theme.major_label_color),
                font(theme.major_label_font),
                fontsize(theme.major_label_font_size))

    {(c, left_guide_position)}
end

# Arrange a plot with its guides
#
# Args:
#   plot_canvas: A canvas containing the plot graphics.
#   guides: Tuples of guide canvases each with a GuidePosition giving
#           where the guide should be placed relative to the plot.
#
# Returns:
#   A new canvas containing the plot with guides layed out in the specified
#   manner.
function layout_guides(plot_canvas::Canvas,
                       theme::Gadfly.Theme,
                       guides::(Canvas, GuidePosition)...)

    # Every guide is updated to use the plot's unit box.
    guides = [(set_unit_box(guide, plot_canvas.unit_box), pos)
              for (guide, pos) in guides]

    # Group by position
    top_guides    = Canvas[]
    right_guides  = Canvas[]
    bottom_guides = Canvas[]
    left_guides   = Canvas[]
    under_guides  = Canvas[]
    for (guide, pos) in guides
        if pos === top_guide_position
            push!(top_guides, guide)
        elseif pos === right_guide_position
            push!(right_guides, guide)
        elseif pos === bottom_guide_position
            push!(bottom_guides, guide)
        elseif pos === left_guide_position
            push!(left_guides, guide)
        elseif pos === under_guide_position
            push!(under_guides, guide)
        end
    end

    # Since top/right/bottom/left guides are drawn without overlaps, we use the
    # canvas z-order to determine ordering, with lowest ordered guides placed
    # nearest to the panel.
    canvas_order = canvas -> canvas.order
    rev_canvas_order = canvas -> -canvas.order
    sort!(by=canvas_order,     top_guides)
    sort!(by=canvas_order,     right_guides)
    sort!(by=canvas_order,     bottom_guides)
    sort!(by=rev_canvas_order, left_guides)

    # Stack the guides on edge edge of the plot
    top_guides    = vstack(0, 0, 1, [(g, hcenter) for g in top_guides]...)
    right_guides  = hstack(0, 0, 1, [(g, vcenter) for g in right_guides]...)
    bottom_guides = vstack(0, 0, 1, [(g, hcenter) for g in bottom_guides]...)
    left_guides   = hstack(0, 0, 1, [(g, vcenter) for g in left_guides]...)

    # Reposition each guide stack, now that we know the extents
    t = top_guides.box.height
    r = right_guides.box.width
    b = bottom_guides.box.height
    l = left_guides.box.width

    pw = 1cx - l - r # plot width
    ph = 1cy - t - b # plot height

    top_guides    = set_box(top_guides,    BoundingBox(l, 0, pw, t))
        # TODO: clip path

    right_guides  = set_box(right_guides,  BoundingBox(l + pw, t, r, ph))
        # TODO: clip path

    bottom_guides = set_box(bottom_guides, BoundingBox(l, t + ph, pw, b)) <<
                      clip((0cx, 1cy - b), (1cx, 1cy - b),
                           (1cx, 1cy), (0cx, 1cy))

    left_guides   = set_box(left_guides,   BoundingBox(0, t, l, ph)) <<
                      clip((0cx, 0cy), (0cx + l, 0cy),
                           (0cx + l, 1cy), (0cx, 1cy))

    compose(canvas(),
            (canvas(l, t, pw, ph),
                {canvas(units_inherited=true, order=-1, clip=true), under_guides...},
                (canvas(units_inherited=true, order=1, clip=true),  plot_canvas),
                d3embed(@sprintf(
                    ".on(\"mouseover\", guide_background_mouseover(parent_id, %s))",
                    json(theme.highlight_color(theme.grid_color)))),
                d3embed(@sprintf(
                    ".on(\"mouseout\", guide_background_mouseout(parent_id, %s))",
                    json(theme.grid_color))),
                d3embed(".call(zoom_behavior(parent_id, t))")),
            top_guides, right_guides, bottom_guides, left_guides)
end

end # module Guide