text:
"""
    kde(x,w)
Weighted kernel density estimate of the data `x` ∈ ℜN with weights `w` ∈ ℜN
`xi, densityw, density = kde(x,w)`
The number of grid points is chosen automatically and will approximately be equal to N/3
The bandwidth of the Gaussian kernel is chosen based on Silverman's rule of thumb
returns both weighted and non-weighted densities `densityw, density`
"""
function kde(x,w::AbstractVector)
    @assert all(x->x>0, w) "All weights must be non-negative"
    e   = StatsBase.histrange(x,ceil(Int,length(x)/3))
    nb  = length(e)-1
    np  = length(x)
    I   = sortperm(x)
    s   = (rand()/nb+0):1/nb:e[end]
    j   = zeros(Float64,nb)
    bo = 1
    ilast = 0
    for i = 1:np
        ii = I[i]
        for b = bo:nb
            if x[ii] <= e[b+1]
                j[b] += w[ii]
                bo = b
                ilast = i
                break
            end
        end
    end
    j[end] += sum(w[I[ilast+1:np]])
    @assert all(j .>= 0) "j"
    xi      = e[1:end-1] .+ 0.5*(e[2]-e[1])
    σ       = std(x)
    if σ == 0 # All mass on one point
        return x,ones(x),ones(x)
    end
    h        = 1.06σ*np^(-1/5)
    K(x)     = exp(-(x/h).^2/2) / √(2π)
    densityw = [1/h*sum(j.*K.(xi.-xi[i])) for i=1:nb] # This density is effectively normalized by nb due to sum(j) = 1
    density  = [1/(np*h)*sum(K.(x.-xi[i])) for i=1:nb]

    # @assert all(densityw .>= 0) "densityw"
    # @assert all(density  .>= 0) "density"

    return xi, densityw, density

end

@userplot DensityPlot
@recipe function f(dp::DensityPlot)
    x = dp.args[1]
    w = length(dp.args) > 1 ? dp.args[2] : ones(x)/length(x)
    xi, densityw, density = kde(x,w)
    if maximum(x)-minimum(x) > 0
        title --> "Kernel density estimate"
        seriestype := :path
        if length(dp.args) > 1
            @series begin
                label --> "Weighted density"
                linecolor --> :blue
                xi,densityw
            end
        end
        @series begin
            label --> "Non-weighted density"
            linecolor --> :red
            xi,density
        end
    end
end

"""
    densityplot(x,[w])
Plot (weighted) particles densities
"""
densityplot

@userplot HeatboxPlot
@recipe function f(p::HeatboxPlot; nbinsy=30)
    x,t = p.args[1:2]
    seriestype := :histogram2d
    nbins --> (1,nbinsy)
    if length(p.args) >= 3
        weights --> p.args[3]
    end
    fill(t, length(x)), x
end

"""
    pdata = pplot(x, w, y, yhat, a, t, pdata; kwargs...)
    pdata = pplot(pf, y, pdata; kwargs...)

To be called inside a particle filter, plots either particle density (`density=true`) or individual particles (`density=false`) \n
Will plot all the real states in `xindices` as well as the expected vs real measurements of `yindices`.
# Arguments:
- `x`: `Vector{Vector}(N)`. The states for each particle where `N` number of Particles
- `w`: `Vector(N)`. weight of each particle
- `y`: `Vector{Vector}(T)`. All true outputs. `T` is total number of time steps (will only use index `t`)
- `yhat`: `Vector{Vector}(N)` The expected output per particle. `N` number of Particles
- `a`, `Vector(N)`, reorderng of particles (e.g. `1:N`)
- `t`, Current time step
- `xreal`: `Vector{Vector}(T)`. All true states if available. `T` is total number of time steps (will only use index `t`)
- `xprev`: Same as `x`, but for previous time step, only used when `!density` to show states origins
- `pdata`: Persistant data for plotting. Set to `nothing` in first call and pdata on remaining \n
- `density = true` To only plot the particle trajectories, set (`leftonly=false`)\n
- `leftonly = true`: only plot the left column\n
- `xindices = 1:n_states`\n
- `yindices = 1:n_measurements`\n
Returns: `pdata`
"""
function pplot(pf, y, args...; kwargs...)
    s = state(pf)
    pplot(s.x, s.we, y, LowLevelParticleFilters.measurement(pf).(s.x, s.t[]), s.j, s.t[], args...; xprev=s.xprev, kwargs...)
end

function pplot(x, w, y, yhat, a, t, pdata; xreal=nothing, xprev=nothing,  density = true, leftonly = true, xindices = 1:length(x[1]), yindices = 1:length(y[1]), lowpass=0.9)

    (t == 1 || t % 35 == 0) &&  @printf "Time     Surviving    Effective nbr of particles\n--------------------------------------------------------------\n"

    T = length(y)
    N = length(x)
    D = length(x[1])
    x = reduce(hcat, x)
    y = reduce(hcat, y)
    yhat = reduce(hcat, yhat)

    cols = leftonly ? 1 : 2
    grd = (r,c) -> (r-1)*cols+c
    @printf("t: %5d %7.3f %9.1f\n", t, (N-length(setdiff(Set(1:N),Set(a))))/N, effective_particles(w))
    plotvals = [x;yhat]
    realVals = xreal === nothing ? [fill(Inf,D,T);y] : [reduce(hcat,xreal);y] # Inf turns plotting off
    if !density
        plotvalsOld = xprev === nothing ? [Inf*x;yhat] : [reduce(hcat, xprev);yhat]
    end

    plotindices = [xindices; size(x,1) .+ yindices]
    if pdata === nothing
        pdata = plot(layout=(length(plotindices),cols)), zeros(length(plotindices),2)
    end
    p, minmax = pdata
    dataMin = minimum(plotvals[plotindices,:])
    dataMax = maximum(plotvals[plotindices,:])
    minmax = [min.(minmax[:,1], dataMin)*lowpass .+ (1-lowpass)*dataMin max.(minmax[:,2], dataMax)*lowpass .+ (1-lowpass)*dataMax]

    for (i, pind) in enumerate(plotindices)
        #Plot the heatmap on the left plot
        if density
            heatboxplot!(plotvals[pind,:], t, w, ylims=minmax[i,:], subplot=grd(i,1), reuse=false)
            if !leftonly
                densityplot!( plotvals[pind,:], w , ylims = tuple(minmax[i,:]...), c=:blue, subplot=grd(i,2), reuse=true, legend=false)
            end
        else
            #Plot the line on the left plot
            plot!(repeat([t-0.5 t+0.5], N)', [plotvalsOld[pind,:] plotvalsOld[pind,:]]', legend=false, subplot=grd(i,1), l=(:black, 0.1))
            # plot!(repeat([t-1 t-0.5], N)', [plotvalsOld[pind,a] plotvals[pind,:]]', legend=false, subplot=grd(i,1), l=(:black, 0.1))
        end
        #Plot Real State Here
        scatter!( [t+0.5], [realVals[pind,t]], subplot=grd(i,1), legend=false, m=(:cyan,))
    end
    # gui(p)
    (p, minmax)
end

function commandplot(f)
    res = f(nothing)
    display(res[1])
    while true
        print("Waiting for command. q to Quit, s NN to skip NN steps:\n")
        line = readline()
        skip = 1
        if line[1] == 'q'
            return
        elseif occursin("s", line)
            ss = split(strip(line,'\n'))
            skip = parse(Int,ss[2])
        end
        foreach(1:skip) do i
            res = f(res)
        end
        display(res[1])
    end
end

"""
    commandplot(pf, u, y; kwargs...)

Produce a helpful plot. For customization options (`kwargs...`), see `?pplot`.
After each time step, a command from the user is requested.
- q: quit
- s n: step `n` steps
"""
function commandplot(pf, u, y; kwargs...)
    # pdata = nothing
    reset!(pf)
    pfp = pf isa AuxiliaryParticleFilter ? pf.pf : pf
    commandplot() do pdata
        pdata = pplot(pfp, y, pdata; kwargs...)
        t = index(pf)
        LowLevelParticleFilters.update!(pf,u[t],y[t])
        pdata
    end
end

"""
    debugplot(pf, u, y; runall=false, kwargs...)

Produce a helpful plot. For customization options (`kwargs...`), see `?pplot`.
- ` runall=false:` if true, runs all time steps befor displaying (faster), if false, displays the plot after each time step.

The generated plot becomes quite heavy. Initially, try limiting your input to 100 time steps to verify that it doesn't crash.
"""
function debugplot(pf, u, y; runall=false, kwargs...)
    pdata = nothing
    reset!(pf)
    pfp = pf isa AuxiliaryParticleFilter ? pf.pf : pf
    for i = 1:length(y)
        pdata = pplot(pfp, y, pdata; kwargs...)
        t = index(pf)
        LowLevelParticleFilters.update!(pf,u[t],y[t])
        runall || display(pdata[1])
    end
    display(pdata[1])
end

@userplot TrajectoryDensity
@recipe function f(p::TrajectoryDensity; nbinsy=30, xreal=nothing)
    pf,x,w,y = p.args[1:4]
    N,T = size(x)
    D = length(x[1])
    P = length(y[1])
    if sum(w) ≉ T
        w = exp.(w)
        w ./= sum(w, dims=1)
    end
    layout := D+P
    label := ""
    markercolor --> :cyan
    title --> reshape([["State $d" for d = 1:D];["Measurement $d" for d = 1:P]], 1, :)
    for d = 1:D
        subplot := d
        @series begin
            seriestype := :histogram2d
            bins --> (0.5:1:T+0.5,nbinsy)
            weights --> vec(w)
            repeat((1:T)' .-0.5,N)[:], vec(getindex.(x,d))
        end
        xreal === nothing || @series begin
            seriestype := :scatter
            1:T, getindex.(xreal,d)
        end
    end
    yhat = measurement(pf).(x,0)
    for d = 1:P
        subplot := d+D
        @series begin
            seriestype := :histogram2d
            bins --> (0.5:1:T+0.5,nbinsy)
            weights --> vec(w)
            repeat((1:T)' .-0.5,N)[:], vec(getindex.(yhat,d))
        end
        @series begin
            seriestype := :scatter
            1:T, getindex.(y,d)
        end
    end
end

@userplot DimensionDensity
@recipe function f(p::DimensionDensity; nbinsy=30, xreal=nothing)
    length(p.args) >=5 || throw(ArgumentError("Supply arguments: pf,x,w,y,dimension::Int"))
    pf,x,w,y,d = p.args[1:5]
    N,T = size(x)
    D = length(x[1])
    P = length(y[1])

    if sum(w) ≉ T
        w = exp.(w)
        w ./= sum(w, dims=1)
    end

    label := ""
    markercolor --> :cyan
    if d <= D
        @series begin
            seriestype := :histogram2d
            bins --> (0.5:1:T+0.5,nbinsy)
            weights --> vec(w)
            repeat((1:T)' .-0.5,N)[:], vec(getindex.(x,d))
        end
        xreal === nothing || @series begin
            seriestype := :scatter
            1:T, getindex.(xreal,d)
        end
    else
        d -= D
        yhat = measurement(pf).(x,0)
        @series begin
            seriestype := :histogram2d
            bins --> (0.5:1:T+0.5,nbinsy)
            weights --> vec(w)
            repeat((1:T)' .-0.5,N)[:], vec(getindex.(yhat,d))
        end
        @series begin
            seriestype := :scatter
            1:T, getindex.(y,d)
        end

    end

end

"""
    dimensiondensity(pf,x,we,y, dimension, nbinsy=30, xreal=nothing)
Same as trajectorydensity but only plots subplot `dimension`.
"""
dimensiondensity

"""
    trajectorydensity(pf,x,we,y, nbinsy=30, xreal=nothing)
Plots particle densities along trajectory. See the readme for an example.
"""
trajectorydensity


text:
function ner_datadep_register()
    register(DataDep("NER Model Weights",
        """
        The weights for NER Sequence Labelling Model.
        """,
        "https://github.com/Ayushk4/NER.jl/releases/download/0.0.0.1/ner_weights.tar.xz",
        "6eda5cd778af99f57a0a0b7eb4d5bc46a5a61c214e3e515e620b7db6b76ce3aa",
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "ner_weights"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))

    register(DataDep("NER Model Dicts",
        """
        The character and words dict for NER Sequence Labelling Model.
        """,
        "https://github.com/Ayushk4/NER.jl/releases/download/0.0.0.1/ner_dicts.tar.xz",
        "49619d793a5974dd41859e68d73eae68e58f8b264d49ba98489ab6ed74bf5f86",
        post_fetch_method = function(fn)
            unpack(fn)
            dir = "model_dicts"
            innerfiles = readdir(dir)
            mv.(joinpath.(dir, innerfiles), innerfiles)
            rm(dir)
        end
    ))
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 4ea44e3a-0730-11ec-1592-4dd06da99612
using Statistics

# ╔═╡ b8014e77-3a78-4bae-9dd5-11dc52bb5157
# Julian Statistics
# We will go through some basic functions and their uses from the Statistics package shipped with JuliaLang.

# ╔═╡ fe54d8f2-fe7d-4f81-a5df-dde04fad9a12
mean(30:60) # Intt64 can be used

# ╔═╡ f134f937-d2b6-443f-87fd-aca05b14ed8a
mean(30.5:45.6) #Float64 can be used

# ╔═╡ c259168a-7fc3-4f88-92c5-600ec733b207
mean([1, missing, 3, 4, 5, missing]) #returns missing

# ╔═╡ 60192011-3913-4257-a1f8-38af13f3beab
mean([1, NaN, 3, 4, 5, NaN, 7, 8]) #returns NaN

# ╔═╡ 87c16954-78cd-4b68-b234-9956c5ba374a
# Trivia Q'
# What do you mean will be the output of the following?

mean([1, missing, 3, 4, NaN, 6, 7, 8])#missing takes precedence if both NaN, missing are present

# ╔═╡ f9646321-6cb1-4693-972b-e39347235fe3
mean(skipmissing([1, missing, 3, 4, 5, 6])) #use skipmissing to remove "missing" and get mean of rest of the data

# ╔═╡ b20a6304-8dc5-4cd3-b5dc-87474dd2bb8e
mat = [3 4 6; 5 7 8] # a 2D matrix

# ╔═╡ 4dd9726f-1ce3-488f-a9ef-f52827ee8e0b
mean!([1., 1.], mat) # get mean along the rows

# ╔═╡ a806544f-24d8-49f1-9fa5-3e4465605716
# Trivia Q'
# What do you mean will be the output of the following?

mean!([7., 7.], mat) # this can be any numeric constant

# ╔═╡ 0b949f0a-08f5-4e02-8174-e630dd165156
mean!([1. 1. 1.], mat) # get mean along the columns

# ╔═╡ 3246fbe2-71f2-48e3-b207-1c49e5aa989c
mean(mat, dims=1) # get mean along one dimension

# ╔═╡ a1b703a6-ce76-4b49-9b23-7169acc90ed5
mean(mat, dims=2) # get mean along two dimensions

# ╔═╡ dfe50e29-1b0a-407d-b7b2-eedcb4d810eb
mean(-, [3, 4, 5]) # apply the function to the data, then compute mean

# ╔═╡ 58c1f4ab-6379-44d8-a047-5d2649e36c9a
median([1, 2, 3])

# ╔═╡ 5bf65349-2990-4f14-95eb-71c753bf2a63
median([1, 2, missing, 4])

# ╔═╡ 7ac95528-2672-4e8a-a51e-efc9287a39d9
median(mat, dims=1)

# ╔═╡ 8e251afb-74be-4f2c-a7f7-c4e03f7c29b7
middle(mat)

# ╔═╡ 19579d09-6d78-4236-925c-0fff03185dd4
middle([1.4, 2.76, 3, 4.27])

# ╔═╡ 1d616172-129c-4e4a-ba1b-83bfff25c40c
quantile!([0, 1, 5, 6, 3, 7, 8], 0.5)

# ╔═╡ ea2f09ad-91ab-4e22-a5bb-6f09ec47056e
quantile(0:20, [0.1, 0.3, 0.9, 0.5, 0.7])

# ╔═╡ 40deb75f-7f21-414e-98b0-606c2ad6bf68
quantile(skipmissing([7, 36, missing]), 0.5)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
"""

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═b8014e77-3a78-4bae-9dd5-11dc52bb5157
# ╠═4ea44e3a-0730-11ec-1592-4dd06da99612
# ╠═fe54d8f2-fe7d-4f81-a5df-dde04fad9a12
# ╠═f134f937-d2b6-443f-87fd-aca05b14ed8a
# ╠═c259168a-7fc3-4f88-92c5-600ec733b207
# ╠═60192011-3913-4257-a1f8-38af13f3beab
# ╠═87c16954-78cd-4b68-b234-9956c5ba374a
# ╠═f9646321-6cb1-4693-972b-e39347235fe3
# ╠═b20a6304-8dc5-4cd3-b5dc-87474dd2bb8e
# ╠═4dd9726f-1ce3-488f-a9ef-f52827ee8e0b
# ╠═a806544f-24d8-49f1-9fa5-3e4465605716
# ╠═0b949f0a-08f5-4e02-8174-e630dd165156
# ╠═3246fbe2-71f2-48e3-b207-1c49e5aa989c
# ╠═a1b703a6-ce76-4b49-9b23-7169acc90ed5
# ╠═dfe50e29-1b0a-407d-b7b2-eedcb4d810eb
# ╠═58c1f4ab-6379-44d8-a047-5d2649e36c9a
# ╠═5bf65349-2990-4f14-95eb-71c753bf2a63
# ╠═7ac95528-2672-4e8a-a51e-efc9287a39d9
# ╠═8e251afb-74be-4f2c-a7f7-c4e03f7c29b7
# ╠═19579d09-6d78-4236-925c-0fff03185dd4
# ╠═1d616172-129c-4e4a-ba1b-83bfff25c40c
# ╠═ea2f09ad-91ab-4e22-a5bb-6f09ec47056e
# ╠═40deb75f-7f21-414e-98b0-606c2ad6bf68
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    getPhenotype(tree::GEPGenotype)
documentation
"""
function getPhenotype(tree::GEPGenotype)

    representation = tree._representation
    lenRep = length(representation)

    phenotype = Array{String}(undef, 0)

    visitedFromRoot = zeros(Integer, lenRep, 4)
    visitedFromRoot[:, 1] = fill(-1, lenRep)
    visitedFromRoot[:, 4] = fill(-1, lenRep)
    visitedNodes = Array{Integer}(undef, 0)
    step = getArity(representation[1])
    currLevel = 0
    i = 1
    visitedFromRoot[i, 1] = currLevel

    while step != 0
        currLevel += 1
        i += 1
        currLevelNodes = collect(i:(i+step-1))
        visitedFromRoot[currLevelNodes, 1] = fill(currLevel, step)
        i += step - 1
        step = isempty(currLevelNodes) ? 0 : reduce(+, getArity.(representation[currLevelNodes]))
    end

    i = 1

    while i != 0

        if typeof(representation[i]) <: FunctionNode

            if findfirst(x->x==i, visitedNodes) == nothing
                push!(visitedNodes, i)
                push!(phenotype, getFunctionName(representation[i]) * "(")
            end

            if visitedFromRoot[i, 3] < getArity(representation[i])

                if visitedFromRoot[i, 4] == -1
                    prevSameLevel = representation[findprevIndexes(x -> visitedFromRoot[i, 1] == x, visitedFromRoot[:, 1], i-1)]
                    postSameLevel = findnextIndexes(x -> visitedFromRoot[i, 1] == x, visitedFromRoot[:, 1], i+1)
                    visitedFromRoot[i, 4] = isempty(prevSameLevel) ? length(postSameLevel) + 1 : reduce(+, getArity.(prevSameLevel)) + length(postSameLevel) + 1
                end

                parent = i
                i += visitedFromRoot[i, 4] + visitedFromRoot[i, 3]
                visitedFromRoot[i, 2] = parent
            else
                push!(phenotype, "),")
                i = visitedFromRoot[i, 2]
                if i != 0
                    visitedFromRoot[i, 3] += 1
                end
            end
        else

            if typeof(representation[i]) <: VariableNode
                push!(phenotype, getVariableName(representation[i]) * ",")

            elseif typeof(representation[i]) <: ConstantNode
                push!(phenotype, string(eval(representation[i])) * ",")

            elseif typeof(representation[i]) <: NoArgsFunctionNode
                push!(phenotype, getNoArgsFunctionName(representation[i]) * "(),")
            end

            i = visitedFromRoot[i, 2]
            if i != 0
                visitedFromRoot[i, 3] += 1
            end
        end
    end

    phenotype[end] = phenotype[end][1:end-1]

    if length(phenotype) == 1
        phenotype = Meta.parse(phenotype[1])
        return quote $phenotype end
    else
        return Meta.parse(join(phenotype))
    end
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
module MarkupGenerator

import JSON, Random
export PKG_ROOT_DIR

const PKG_ROOT_DIR = normpath(joinpath(@__DIR__, ".."))

include("./Utils.jl")
include("./Markup.jl")
include("./Abstract.jl")
include("./SVG.jl")

end # end of module


text:
using Test
using Random
using LinearAlgebra

using Celerite2

function get_term()
  term = Celerite2.SHOTerm(1.5, 0.8, 2.3)
  term += Celerite2.SHOTerm(1.5, 0.8, 0.3)
  return term
end

function get_coeffs()
  return Celerite2.get_coefficients(get_term())
end

function get_matrices()
  N, x, diag, Z = get_data()
  return Celerite2.get_matrices(get_term(), x, diag)..., Z
end

function get_data()
  Random.seed!(42)
  N = 50
  x = 100 * sort(Random.rand(N))
  diag = 1 .+ Random.rand(N)
  Z = Random.randn(50, 2)
  return N, x, diag, Z
end

function test_to_dense()
  N, x, diag, Z = get_data()
  ar, cr, ac, bc, cc, dc = get_coeffs()

  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  # Brute force
  tau = abs.(x .- transpose(x))
  K0 = zeros(N, N)
  for j in 1:size(ar, 1)
    K0 = K0 + ar[j] .* exp.(-cr[j] .* tau)
  end
  for j in 1:size(ac, 1)
    K0 = K0 + (ac[j] .* cos.(dc[j] .* tau) + bc[j] .* sin.(dc[j] .* tau)) .* exp.(-cc[j] .* tau)
  end
  K0 += Diagonal(diag)

  return maximum(abs.(K - K0)) < 1e-12
end

function test_matmul()
  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  Y0 = K * Z
  Y = Celerite2.matmul(a, U, V, P, Z)
  return maximum(abs.(Y - Y0)) < 1e-12
end

function test_factor()
  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  Celerite2.factor!(U, P, a, V)
  L0 = cholesky(K).L
  L = LowerTriangular(Celerite2.to_dense(ones(size(a, 1)), U, V, P))

  return (
    maximum(abs.(L * Diagonal(a) * transpose(L) - K)) < 1e-12 &&
    maximum(abs.(L * Diagonal(sqrt.(a)) - L0)) < 1e-12
  )
end

function test_dot_tril()
  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  Celerite2.factor!(U, P, a, V)
  L = cholesky(K).L

  Y0 = L * Z
  Celerite2.dot_tril!(U, P, a, V, Z)
  return maximum(abs.(Z - Y0)) < 1e-12
end

function test_solve()
  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  Celerite2.factor!(U, P, a, V)
  L = cholesky(K).L

  Y0 = K \ Z
  Celerite2.solve!(U, P, a, V, Z)
  return maximum(abs.(Z - Y0)) < 1e-12
end

@testset "core" begin

  @test test_to_dense()
  @test test_matmul()
  @test test_factor()
  @test test_dot_tril()
  @test test_solve()

end

function test_logp()
  N, x, diag, Z = get_data()
  y = Z[:, 1]
  a, U, V, P, Z = get_matrices()
  K = Celerite2.to_dense(a, U, V, P)

  term = get_term()
  logp = Celerite2.logp(term, x, y, sqrt.(diag))

  chol = cholesky(K)
  logp0 = -0.5 * transpose(y) * (chol \ y)
  logp0 -= 0.5 * logdet(chol)
  logp0 -= 0.5 * N * log(2 * pi)

  return abs(logp - logp0) < 1e-12
end

@testset "gp" begin

  @test test_logp()

end


text:
# # Modal analysis of a free-floating steel circle

# Source code: [`circle_modal_tut.jl`](circle_modal_tut.jl)

# ## Description

# Vibration analysis of a beam simply supported in one plane, and clamped in
# another. The results are compared with analytical expressions. This is a
# benchmark from the NAFEMS Selected Benchmarks for Natural Frequency Analysis,
# publication: Test VM09: Circular Ring --  In-plane and Out-of-plane
# Vibration.

# ## Reference frequencies

# There will be 6 rigid body modes (zero natural frequencies).

# The numerical results are due to the publication:
# NAFEMS Finite Element Methods & Standards, Abbassian, F., Dawswell, D. J., and
# Knowles, N. C. Selected Benchmarks for Natural Frequency Analysis, Test No.
# 6. Glasgow: NAFEMS, Nov., 1987.

# The reference values were analytically determined (Blevins, FORMULAS FOR
# DYNAMICS, ACOUSTICS AND VIBRATION, Table 4.16). Note that shear flexibility
# was neglected when computing the reference values.

# | Mode       |         Reference Value (Hz)  |  NAFEMS Target Value (Hz) |
# | -------   |     -------  |  ------- |
# | 7, 8 | (out of plane)   |        51.85          |         52.29  |
# | 9, 10 |  (in plane)       |       53.38         |          53.97  |
# | 11, 12 |  (out of plane)   |     148.8          |         149.7  |
# | 13, 14 |  (in plane)       |     151.0          |         152.4  |
# | 15, 16 |  (out of plane)   |     287.0          |         288.3  |
# |  17, 18 |  (in plane)      |      289.5         |          288.3  |

# ## Goals

# - Show convergence relative to reference values.
# - Demonstrate the optimization of eigenvalue accuracy by choosing mass type.

##
# ## Definition of the basic inputs

# The finite element code realize on the basic functionality implemented in this
# package.
using FinEtools

# The material parameters may be defined with the specification of the units.
# The elastic properties are:
E = 200.0 * phun("GPa")
nu = 0.3;

# The mass density is
rho = 8000 * phun("kg/m^3")
# Here are the cross-sectional dimensions and the length of the beam between
# supports.
radius = 1.0 * phun("m"); diameter = 0.1 * phun("m");

##
# ## Cross-section

# Cross-sectional properties are incorporated in the cross-section property. The
# three arguments supplied are functions. All are returning "constants". In
# particular the first two functions each return the dimension of the
# cross-section as a constant(the beam has a uniform cross-section); the third
# function defines the orientation of the cross-section in the global Cartesian
# coordinates. `[1.0, 0.0, 0.0]` is the vector that together with the tangent
# to the midline curve of the beam spans the $x_1x_2$ plane of the local
# coordinates for the beam.
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionCircle
cs = CrossSectionCircle(s -> diameter/2, s -> [1.0, 0.0, 0.0])
@show cs.parameters(0.0)

# For instance the the first out of plane mode is listed in the reference cited
# above as
R = radius
I = cs.parameters(0.0)[4]
m = rho * cs.parameters(0.0)[1]
J = cs.parameters(0.0)[2]
G = E/2/(1+nu)
i = 2 # the first non-rigid body mode
@show i*(i^2-1)/(2*pi*R^2)*sqrt(E*I/m/(i^2+E*I/G/J))

# The first "ovaling" (in-plane) mode is:
i = 2 # the first ovaling mode
@show i*(i^2-1)/(2*pi*R^2*(i^2+1)^(1/2))*sqrt(E*I/m)

# The purpose of the numerical model is to calculate approximation to the reference frequencies.

neigvs = 18;

# We will generate
n = 20
# beam elements along the member.
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member
tolerance = radius/n/1000;
fens, fes = frame_member([0 0 0; 2*pi 0 0], n, cs)
for i in 1:count(fens)
    a = fens.xyz[i, 1]
    fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
end
fens, fes = mergenodes(fens, fes, tolerance, [1, n+1])

##
# ## Material

# Material properties can be now used to create a material: isotropic elasticity model of the `FinEtoolsDeforLinear` package is instantiated.
using FinEtoolsDeforLinear
material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

##
# ## Fields

# Now we start constructing the discrete finite element model.
# We begin by constructing the requisite fields, geometry and displacement.
# These are the so-called "configuration variables", all initialized to 0.
# This is that geometry field.
geom0 = NodalField(fens.xyz)
# This is the displacement field, three unknown displacements per node.
u0 = NodalField(zeros(size(fens.xyz, 1), 3))
# This is the rotation field, three unknown rotations per node are represented
# with a rotation matrix, in total nine numbers. The utility function
# `initial_Rfield`
using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield
Rfield0 = initial_Rfield(fens)
# Here we verify the number of nodes and the number of degrees of freedom in the
# rotation field per node.
@show nents(Rfield0)
@show ndofs(Rfield0)

# Finally, this is the displacement and rotation field for incremental changes,
# incremental displacements and incremental rotations. In total, 6 unknowns per
# node.
dchi = NodalField(zeros(size(fens.xyz, 1), 6))

# There are no support conditions.
applyebc!(dchi)
# The  the number of free
# (unknown) degrees of freedom is equal to the total number of degrees of freedom in the system.
numberdofs!(dchi);

##
# ## Assemble the global discrete system

using FinEtoolsFlexBeams.FEMMCorotBeamModule: FEMMCorotBeam
femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material);

# For disambiguation we will refer to the stiffness and mass functions by qualifying them with the corotational-beam module, `FEMMCorotBeamModule`.
using FinEtoolsFlexBeams.FEMMCorotBeamModule
CB = FEMMCorotBeamModule
# Thus we can construct the stiffness and mass matrix as follows:
# Note that the finite element machine is the first argument. This provides
# access to the integration domain. The next argument is the geometry field,
# followed by the displacement, rotations, and incremental
# displacement/rotation fields.
K = CB.stiffness(femm, geom0, u0, Rfield0, dchi);
M = CB.mass(femm, geom0, u0, Rfield0, dchi);
# We can compare the size of the stiffness matrix with the number of degrees of
# freedom that are unknown (20).
@show size(K)

##
# ## Solve the free-vibration problem

oshift = (2*pi*15)^2

# The Arnoldi algorithm implemented in the well-known `Arpack` package is used
# to solve the generalized eigenvalue problem with the sparse matrices. As is
# common in structural dynamics, we request the smallest eigenvalues in
# absolute value (`:SM`).
using Arpack
evals, evecs, nconv = eigs(K + oshift * M, M; nev=neigvs, which=:SM);
# First  we should check that the requested eigenvalues actually converged:
@show nconv == neigvs

# The eigenvalues (i. e. the squares of the angular frequencies) are returned in
# the vector `evals`. The mode shapes constitute the columns of the matrix `evecs`.
@show size(evecs)
# The natural frequencies are obtained from the squares of the angular
# frequencies. We note the use of `sqrt.` which broadcast the square root over
# the array `evals`.
fs = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);

##
# ## Comparison of computed and analytical results

# The approximate and analytical frequencies are now reported.
sigdig(n) = round(n * 10000) / 10000
println("Approximate frequencies: $(sigdig.(fs)) [Hz]")

##
# ## Set up the visualization of the vibration modes

# The animation will show one of the vibration modes overlaid on the undeformed geometry. The configuration during the animation needs to reflect rotations. The function `update_rotation_field!` will update the rotation field given a vibration mode.
using FinEtoolsFlexBeams.RotUtilModule: update_rotation_field!

# The visualization utilities take advantage of the PlotlyJS library.
using PlotlyJS
using FinEtoolsFlexBeams.VisUtilModule: plot_space_box, plot_solid, render, react!, default_layout_3d, save_to_json

# The magnitude of the vibration modes (displacements  and rotations) will be amplified with this scale factor:
scale = 1.5

# In order to handle variables inside loops correctly, we create a local scope with the `let end` block.
vis(mode) = let
    # The extents of the box will be preserved during animation in order to eliminate changes in the viewing parameters.
    tbox = plot_space_box(reshape(inflatebox!(boundingbox(fens.xyz), 0.3 * radius), 2, 3))
    # This is the geometry of the structure without deformation (undeformed). It is displayed as gray, partially transparent.
    tenv0 = plot_solid(fens, fes; x=geom0.values, u=0.0 .* dchi.values[:, 1:3], R=Rfield0.values, facecolor="rgb(125, 155, 125)", opacity=0.3);
    # Initially the plot consists of the box and the undeformed geometry.
    plots = cat(tbox, tenv0; dims=1)
    # Create the layout for the plot. Set the size of the window.
    layout = default_layout_3d(;width=600, height=600)
    # Set the aspect mode to get the correct proportions.
    layout[:scene][:aspectmode] = "data"
    # Render the undeformed structure
    pl = render(plots; layout=layout, title="Mode $(mode)")
    sleep(2.115)
    # This is the animation loop.
    # 1. Distribute a fraction of the selected eigenvector into the incremental displacement/rotation field.
    # 2. Create the deformed configuration by defining displacement field `u1` and rotation field `Rfield1`.
    # 3. Create the plot for the deformed configuration, and add it to the list of plots.
    # 4. Call the `react!` function to update the display. Sleep for a brief period of time to give the display a chance to become current.
    for xscale in scale .* sin.(collect(0:1:89) .* (2 * pi / 21))
        scattersysvec!(dchi, xscale .* evecs[:, mode])
        u1 = deepcopy(u0)
        u1.values .= dchi.values[:, 1:3]
        Rfield1 = deepcopy(Rfield0)
        update_rotation_field!(Rfield1, dchi)
        tenv1 = plot_solid(fens, fes; x=geom0.values, u=dchi.values[:, 1:3], R=Rfield1.values, facecolor="rgb(50, 55, 125)");
        plots = cat(tbox, tenv0, tenv1; dims=1)
        react!(pl, plots, pl.plot.layout)
        sleep(0.115)
    end
end

##
# ## Visualize vibration mode

# Animate the harmonic motion of the mode given as argument:
# vis(7)

using FinEtoolsFlexBeams.FESetCorotBeamModule: MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA,
MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA,
MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA,
MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA

results = let
    results = Dict()
    for mtype in [
        MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA,
        MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA,
        MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA,
        MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA]
        M = CB.mass(femm, geom0, u0, Rfield0, dchi; mass_type = mtype);
        evals, evecs, nconv = eigs(K + oshift * M, M; nev=neigvs, which=:SM);
        results[mtype] = evals, evecs
    end
    results
end

colors = [
"rgb(125, 15, 15)",
"rgb(15, 155, 15)",
"rgb(15, 15, 155)",
"rgb(125, 115, 115)"
]

evals = results[MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA][1]
x = 1:length(evals); y = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);
tc0 = scatter(;x=x, y=y, mode="markers", name = "cons, wo", line_color = "rgb(215, 15, 15)", marker = attr(size = 9, symbol = "diamond-open"))
evals = results[MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA][1]
x = 1:length(evals); y = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);
tc1 = scatter(;x=x, y=y, mode="markers", name = "cons, w", line_color = "rgb(15, 215, 15)", marker = attr(size = 9, symbol = "triangle-down"))
evals = results[MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA][1]
x = 1:length(evals); y = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);
tc2 = scatter(;x=x, y=y, mode="markers", name = "lumped, wo", line_color = "rgb(15, 15, 215)", marker = attr(size = 9, symbol = "x-open"))
evals = results[MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA][1]
x = 1:length(evals); y = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);
tc3 = scatter(;x=x, y=y, mode="markers", name = "lumped, w", line_color = "rgb(165, 165, 15)", marker = attr(size = 9, symbol = "square"))

# 7, 8 (out of plane)         51.85                 52.29
# 9, 10 (in plane)            53.38                 53.97
# 11, 12 (out of plane)      148.8                 149.7
# 13, 14 (in plane)          151.0                 152.4
# 15, 16 (out of plane)      287.0                 288.3
# 17, 18 (in plane)          289.5                 288.3
rfs = vec(Float64[0 0 0 0 0 0 51.85 51.85 53.38 53.38 148.8 148.8 151.0 151.0 287.0 287.0 289.5 289.5])

rtc = scatter(;x=collect(1:length(rfs)), y=rfs, mode="lines", name = "ref", line_color = "rgb(15, 15, 15)")

# Set up the layout:
layout = Layout(;width=650, height=400, xaxis=attr(title="Mode", type = "linear"), yaxis=attr(title="Frequency [hertz]", type = "linear"), title = "Comparison of mass types")
# Plot the graphs:
pl = plot([rtc, tc0, tc1, tc2, tc3], layout; options = Dict(
        :showSendToCloud=>true,
        :plotlyServerURL=>"https://chart-studio.plotly.com"
        ))
display(pl)

##
# ## "Mixed" mass matrix

# As the graph suggests we can try mixing together mass matrices computed from
# different assumptions (lumped versus consistent, with or without rotation
# inertia). No particular justification can be provided for these numbers other
# than best low-frequency matching for linear bars suggests that 50% can be
# taken lumped and 50% consistent. Here it is obviously more complicated due to
# the presence of rotation inertia.
M = 0.5 .* CB.mass(femm, geom0, u0, Rfield0, dchi; mass_type = MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA) +
    0.5 .* CB.mass(femm, geom0, u0, Rfield0, dchi; mass_type = MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA);
# With this mixed mass matrix we solve the free vibration problem again.
evals, evecs, nconv = eigs(K + oshift * M, M; nev=neigvs, which=:SM);

# Plotting the newly obtained data on top of the previously presented data, we
# can observe sometimes substantial improvement of accuracy of the mixed-matrix
# formulation relative to the individual mass matrix types.
x = 1:length(evals); y = sqrt.([max(0, e - oshift) for e in evals]) / (2 * pi);
mtc = scatter(;x=x, y=y, mode="markers", name = "mixed", line_color = "rgb(215, 15, 215)", marker = attr(size = 9, symbol = "circle"))
pl = plot([rtc, tc0, tc1, tc2, tc3, mtc], layout; options = Dict(
        :showSendToCloud=>true,
        :plotlyServerURL=>"https://chart-studio.plotly.com"
        ))
display(pl)


text:
using ExecutableSpecifications
using Base.Test

include("runsmalltests.jl")


text:
using Comonicon
using Test
using Pkg

@testset "cast" begin
    include("cast.jl")
end

@testset "markdown" begin
    include("markdown.jl")
end

@testset "projects" begin
    include("projects.jl")
end

empty!(ARGS)
push!(ARGS, "arg", "--opt1=2", "--opt2", "3", "-f")
@test Base.include(Main, "scripts/hello.jl") == 0

empty!(ARGS)
push!(ARGS, "activate", "-h")
@test Base.include(Main, "scripts/pkg.jl") == 0

empty!(ARGS)
push!(ARGS, "activate", "path", "--shared")
@test Base.include(Main, "scripts/pkg.jl") == 0

empty!(ARGS)
push!(ARGS, "Author - Year.pdf")
@test Base.include(Main, "scripts/searchpdf.jl") == 0


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*cos(e+f*x))^2*(-1c+c*cos(e+f*x))^-1*sec(e+f*x)^2, x) == :(-1*a^2*c^-1*f^-1*tan(e+f*x)+-3*a^2*c^-1*f^-1*arctanh(sin(e+f*x))+4*a^2*c^-1*f^-1*(1+-1*cos(e+f*x))^-1*sin(e+f*x))

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.2.3 (g cos)^p (a+b cos)^m (c+d cos)^n.jl


text:
#https://doi.org/10.1128/AEM.00062-07
#The Classifier uses a feature space consisting of all possible 8-base
#subsequences (words).
# The psition of the word in a sequence is ignored.
# Only the words occurring in the query contribute to the score.

base_vals = Base.ImmutableDict(DNA_A=>0,DNA_C=>1,DNA_G=>2,DNA_T=>3)

word_prior(n,N) = (n+0.5)/(N+1)
function word_priors!(priors ::Vector,N)
    @turbo for i in 1:N
        priors[i] = word_prior(priors[i],N)
    end
end

conditional_prob(m,P,M) = (m+P)/(M+1)

function get_mer_idx(mer)
    idx = 0
    for base in mer
        idx = 4*idx + base_vals[base]
    end
    return idx +1
end

function count_ref_mers!(ref_seq,column_idx,mer_array,mer_vec,k = 8)
    for mer in each(DNAMer{k}, ref_seq)
        current_mer = mer.fw
        current_mer_idx = get_mer_idx(current_mer)
        mer_array[current_mer_idx,column_idx] = 1
        mer_vec[current_mer_idx] +=1
    end
end

function count_seq_mers(seq,k = 8)
    mask = zeros(Bool,4^k)
    for mer in each(DNAMer{k}, seq)
        current_mer = mer.fw
        current_mer_idx = get_mer_idx(current_mer)
        mask[current_mer_idx] = 1
    end
    return mask, (1:4^k)[mask]
end

function count_mers(refs, k = 8)
    n_refs = length(refs)
    mer_vec= zeros(Float64,4^k)
    mer_array= zeros(Bool,4^k,n_refs)
    for i in 1:n_refs
        count_ref_mers!(refs[i],i,mer_array,mer_vec,k)
    end
    return mer_vec, mer_array
end

function count_mers(refs, genus_vec ::Vector, k = 8)
    n_refs = length(refs)
    mer_vec= zeros(Float64,4^k)
    genus_array= zeros(Int64,4^k,maximum(genus_vec))
    mer_array= zeros(Bool,4^k,n_refs)
    for i in 1:n_refs
        count_ref_mers!(refs[i],i,mer_array,mer_vec,k)
        genus_array[:,genus_vec[i]] .+= mer_array[:,i]
    end
    return mer_vec, genus_array
end

function assign(seq_mask,log_probs)
    cond_probs = vec(sum(log_probs[seq_mask,:], dims = 1))
    println(maximum(cond_probs))
    return findmax(cond_probs)[2]
end

function bootstrap(seq_inds,log_probs,assignment, sample_size)
    hits = 0
    for i in 1:100
        inds = rand(seq_inds,sample_size)
        cond_probs = vec(sum(log_probs[inds,:], dims = 1))
        if findmax(cond_probs)[2] == assignment
            hits +=1
        end
    end
    return hits/100
end

function naieve_bayes(seqs,refs)
    N = length(refs)
    assignments = []
    confs = []
    priors, a =count_mers(refs)
    word_priors!(priors,N)
    log_probs = log.(conditional_prob.(a,priors,1))
    for seq in seqs
        seq_mask, seq_inds = count_seq_mers(seq)
        sample_size = sum(seq_mask) ÷8
        assignment = assign(seq_mask,log_probs)
        push!(assignments,assignment)
        push!(confs, bootstrap(seq_inds,log_probs,assignment,sample_size))
    end
    return assignments, confs
end

function naieve_bayes(seqs,refs,genus_vec)
    N = length(refs)
    assignments = []
    confs = []
    M = reverse([sum(genus_vec .==i) for i in 1:maximum(genus_vec)]) # reverse gives right answer. but not sure why
    priors, a =count_mers(refs,genus_vec)
    word_priors!(priors,N)
    log_probs = log.(conditional_prob.(a,priors,M')) #M` gave wrong results. need to figure out right way...
    for seq in seqs
        seq_mask, seq_inds = count_seq_mers(seq)
        sample_size = sum(seq_mask)  ÷8
        assignment = assign(seq_mask,log_probs)
        push!(assignments,assignment)
        push!(confs, bootstrap(seq_inds,log_probs,assignment,sample_size))
    end
    return assignments, confs
end

function get_genera(fasta)
    genera =[]
    open(FASTA.Reader,fasta) do reader
        for record in reader
            push!(genera,split(identifier(record),";")[5])
        end
    end
    unq = unique(genera)
    n = length(unq)
    mapping = Base.ImmutableDict([unq[i] => i for i in 1:n]...)
    reverse_mapping = Base.ImmutableDict([i => unq[i] for i in 1:n]...)
    return [mapping[genus] for genus in genera], reverse_mapping
end

function get_genera(genera ::Vector)
    unq = unique(genera)
    n = length(unq)
    mapping = Base.ImmutableDict([unq[i] => i for i in 1:n]...)
    reverse_mapping = Base.ImmutableDict([i => unq[i] for i in 1:n]...)
    return [mapping[genus] for genus in genera], reverse_mapping
end

function assignTaxonomy(seqs,ref_fasta)
    genera = []
    refs = []
    open(FASTA.Reader, ref_fasta) do reader
        for record in reader
            if length(findall(";",identifier(record))) >4
                push!(refs,sequence(record))
                push!(genera,split(identifier(record),";")[6])
            end
        end
    end
    genus_vec, genus_mapping = get_genera(genera)

    assignments, confs = naieve_bayes(seqs,refs,genus_vec)
    return hcat(assignments,[genus_mapping[a] for a in assignments],confs)

end

using FASTX


text:
if basename(pwd()) == "aoc"
    cd("2016/6")
end
using StatsBase: countmap

function loadcs(filename)
    reshape(
        [c for c in readchomp(filename) if !isspace(c)],
        length(readline(filename)),
        :
    )
end

function part1(filename::AbstractString)
    join(findmax(countmap(row))[2] for row in eachrow(loadcs(filename)))
end

@assert part1("example.txt") == "easter"
part1("input.txt")

function part2(filename::AbstractString)
    join(findmin(countmap(row))[2] for row in eachrow(loadcs(filename)))
end

@assert part2("example.txt") == "advent"
part2("input.txt")

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
function train_network2d(datatyp, n_dims, η, β, number_epochs, Sb, L,
                        hidden_units_1, hidden_units_2, hidden_units_3,
                        hidden_units_4, hidden_units_5)

    #Load data
    x_train = h5open("datasets/$(n_dims)d/traindata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "X")
    end
    y_train = h5open("datasets/$(n_dims)d/traindata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "Y")
    end
    x_valid = h5open("datasets/$(n_dims)d/validdata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "X")
    end
    y_valid = h5open("datasets/$(n_dims)d/validdata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "Y")
    end

    @info("Building model...")
    input_size = size(x_train)[1]
    model2d = Flux.Chain(
                    Flux.Dense(input_size, hidden_units_1, leakyrelu, initW=Flux.glorot_normal),
                    Flux.Dense(hidden_units_1, hidden_units_2, leakyrelu, initW=Flux.glorot_normal),
                    Flux.Dense(hidden_units_2, hidden_units_3, leakyrelu, initW=Flux.glorot_normal),
                    Flux.Dense(hidden_units_3, hidden_units_4, leakyrelu, initW=Flux.glorot_normal),
                    Flux.Dense(hidden_units_4, hidden_units_5, leakyrelu, initW=Flux.glorot_normal),
                    Flux.Dense(hidden_units_5, 2, initW=Flux.glorot_normal),
                    Flux.softmax)

    opt = Flux.ADAM(η)
    ps = Flux.params(model2d)[1:2:11]
    accuracy(ŷ, y) = mean(onecold(ŷ) .== onecold(y))
    sqnorm(x) = sum(abs2, x)
    loss(x, y) = crossentropy(model2d(x), y) + β * sum(sqnorm, ps)

    @info("Beginning training loop...")
    best_acc = 0.0
    last_improvement = 0
    for epoch_idx in 1:number_epochs
        # Create the full dataset
        train_data = DataLoader((x_train, y_train) ; batchsize=Sb, shuffle=true)
        Flux.train!(loss, params(model2d), train_data, opt)
        acc = accuracy(model2d(x_valid), y_valid)
        acc2 = accuracy(model2d(x_train), y_train)
        loss_idx = loss(x_train, y_train)
        parasum = sum(sqnorm, params(model2d))
        @show epoch_idx, acc, loss_idx, parasum, acc2
        #@save "networks/$(n_dims)d/model$(n_dims)d$(datatyp)-$(acc)-$(β).bson" model2d

        if acc > best_acc
            best_acc = acc
            last_improvement = epoch_idx
            @save "networks/$(n_dims)d/model$(n_dims)d$(datatyp)-$(acc)-$(β).bson" model2d
        end

        if epoch_idx - last_improvement >= L
            @warn("Early stopping")
            @show best_acc
            break
        end
    end
end

function train_network2dcnn(datatyp, n_dims, η, β, number_epochs, Sb, L)

    #Load data
    x_train = h5open("datasets/$(n_dims)d/traindata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "X")
    end
    y_train = h5open("datasets/$(n_dims)d/traindata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "Y")
    end
    x_valid = h5open("datasets/$(n_dims)d/validdata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "X")
    end
    y_valid = h5open("datasets/$(n_dims)d/validdata$(n_dims)d$(datatyp).h5", "r") do file
        read(file, "Y")
    end

    x_train = Flux.unsqueeze(x_train, 3)
    x_valid = Flux.unsqueeze(x_valid, 3)

    @info("Building model...")
    model2dcnn = Chain(
                        Conv((3, 3), 1=>32, pad=(1,1), leakyrelu),
                        Conv((3, 3), 32=>16, pad=(1,1),leakyrelu),
                        Conv((3, 3), 16=>8, pad=(1,1), leakyrelu),
                        Conv((3, 3), 8=>4, pad=(1,1), leakyrelu),
                        Conv((3, 3), 4=>1, pad=(1,1), leakyrelu),
                        flatten,
                        Dense(16, 8, leakyrelu),
                        Dense(8, 2),
                        softmax
                    )

    opt = ADAM(η)
    accuracy(ŷ, y) = mean(onecold(ŷ) .== onecold(y))
    sqnorm(x) = sum(abs2, x)
    loss(x, y) = Flux.crossentropy(model2dcnn(x), y)

    @info("Beginning training loop...")
    best_acc = 0.0
    last_improvement = 0
    for epoch_idx in 1:number_epochs
        # Create the full dataset
        train_data = DataLoader((x_train, y_train) ; batchsize=Sb, shuffle=true)
        Flux.train!(loss, params(model2dcnn), train_data, opt)
        acc = accuracy(model2dcnn(x_valid), y_valid)
        acc2 = accuracy(model2dcnn(x_train), y_train)
        loss_idx = loss(x_train, y_train)
        parasum = sum(sqnorm, params(model2dcnn))
        @show epoch_idx, acc, loss_idx, parasum, acc2
        @save "networks/$(n_dims)d/model$(n_dims)d$(datatyp)-$(acc2)-$(β).bson" model2dcnn

        if acc > best_acc
            best_acc = acc
            last_improvement = epoch_idx
            @save "networks/$(n_dims)d/model$(n_dims)d$(datatyp)-$(acc)-$(β).bson" model2dcnn
        end

        if epoch_idx - last_improvement >= L
            @warn("Early stopping")
            @show best_acc
            break
        end
    end
end


text:
# Dog <-> World
function move( dog::Dog, world::World, dt, motion::Motion )
    new_position, new_velocity = step( motion, dog.position, dog.velocity, dt )
    new_position = pbc( new_position, world )
    return Dog( dog; position = new_position, velocity = new_velocity )
end # function

function move_dogs!( world::World, dt )
    dogs = world.dogs
    dogs.member .= move.( dogs.member, Ref(world), dt, Ref(dogs.motion) )
    return nothing
end # function


text:
module ONSUD

using ZipFile
using CSV
using Serialization
using Index1024
using DataFrames
using Proj4

export check_for_update

export generate

export uprn_data, create_index1024, index_by_postcode, open_pcode_index, pcode_info, save, load
export postcode_to_UInt64, UInt64_to_postcode

const mask = 0xffff000000000000 # 64k files should be enough for anybody
const shift = 48

const LaLoDim = typeof((la=0.0, lo=0.0, dim=zero(UInt64)))
const LaLo = typeof((la=0.0, lo=0.0))

tag(n, v) = UInt64(v) | (UInt64(n) << shift)
tag(v::UInt64) = UInt32((UInt64(v) & mask) >> shift)
key(n::Index1024.NodeInfo) = key(n.tagged_key)
key(v::UInt) = UInt64(v) & ~mask

const Grid = NamedTuple{(:la, :lo, :pc), Tuple{Float64, Float64, UInt64}}

const GEODIR = "/home/matt/wren/UkGeoData"
const DATADIR = joinpath(GEODIR, "ONSUD_NOV_2021/Data")

const Field = Symbol

struct UPRNDB
    grid::Dict{Int64, Grid}
    field2uprn::Dict{Field, Dict{String, Vector{Int64}}}
    uprn2dimension::Dict{Int64, UInt64}
    dimensions::Dict
    UPRNDB() = new(Dict{Int64, Grid}(), Dict{Field, Dict{String, Vector{Int64}}}(), Dict{Int64, UInt64}(), Dict())
end

str(v) = ismissing(v) ? "missing" : String(v)

function dimension(db::UPRNDB, row)
    dim = Dict{Field, String}()
    for k in keys(db.field2uprn)
        dim[Field(k)] = str(row[k])
    end
    (;dim...)
end

function dimension_index(db::UPRNDB, row)
    dim = dimension(db, row)
    hsh = hash(dim)
    if ! (hsh in keys(db.dimensions))
        db.dimensions[hsh] = dim
    end
    hsh
end

function add!(db::UPRNDB, row, en2lola)
    uprn = parse(Int, row.uprn)
    e = parse(UInt64, row.gridgb1e)
    n = parse(UInt64, row.gridgb1n)
    lo, la = en2lola((e, n))
    db.grid[uprn] = (la=la, lo=lo, pc=postcode_to_UInt64(row.pcds))
    db.uprn2dimension[uprn] = dimension_index(db, row)
    for f in keys(db.field2uprn)
        v = str(row[f])
        if ! (v in keys(db.field2uprn[f]))
            db.field2uprn[f][v] = [uprn]
        else
            push!(db.field2uprn[f][v], uprn)
        end
    end
end

function uprninfo(db::UPRNDB, uprn)
    if uprn in keys(db.grid)
        (grid=db.grid[uprn], db.dimensions[db.uprn2dimension[uprn]]...)
    end
end

check_for_update() = println("Visit https://geoportal.statistics.gov.uk/search?sort=-created&tags=onsud")

zipped_row_readers(zipfile) = map(file->(file.name, ()->CSV.Rows(read(file))), filter(f->startswith(f.name, "Data/") && endswith(f.name, ".csv"), ZipFile.Reader(zipfile).files))
row_readers(datadir) = map(n->(n, ()->CSV.Rows(read(n))), readdir(datadir, join=true))

# @pipe ONSUD.row_readers(ONSUD.DATADIR) |> generate |> ONSUD.save(joinpath(ONSUD.GEODIR, "nov_2021.uprndb"), _)

function generate(reader::Tuple)
    println(reader[1])
    en2lola = Proj4.Transformation("EPSG:27700", "EPSG:4326", always_xy=true)
    db = UPRNDB()
    rows = reader[2]()
    foreach(f->db.field2uprn[f] = Dict{String, Vector{Int64}}(), filter(f->!(f in [:uprn, :gridgb1e, :gridgb1n, :pcds]), rows.names))
    foreach(row->add!(db, row, en2lola), rows)
    db
end

function generate(readers)
    dbs = Vector{UPRNDB}(undef, length(readers))

    Threads.@threads for i in 1:length(readers)
        dbs[i] = generate(readers[i])
    end
    println(stderr, "generated")
    uprndb = UPRNDB()
    for db in dbs
        merge!(uprndb.grid, db.grid)
        merge!(uprndb.field2uprn, db.field2uprn)
        merge!(uprndb.uprn2dimension, db.uprn2dimension)
        merge!(uprndb.dimensions, db.dimensions)
    end
    uprndb
end

function save(io::IO, object)
    serialize(io, object)
    object
end

function save(memofile::AbstractString, object)
    try
        open(memofile, "w+") do io
           return save(io, object)
        end
    catch e
        println(stderr, e)
    end
    object
end

function load(memofile)
    if filesize(memofile) > 0
        open(memofile, "r") do io
            return deserialize(io)
        end
    end
end

############### Index1024 stuff

function postcode_to_UInt64(pc)
    if ismissing(pc)
        return 0
    end
    m = match(r"([A-Z]+)([0-9]+([A-Z]+)?) ?([0-9]+)([A-Z]+)", replace(pc, " "=>""))
    if m == nothing || m[1] === nothing || m[2] === nothing || m[4] === nothing || m[5] === nothing
        return 0
    end
    reduce((a,c) -> UInt64(a) << 8 + UInt8(c), collect(lpad(m[1], 2) * lpad(m[2], 2) * m[4] * m[5]), init=0)
end

function UInt64_to_postcode(u)
    if u == 0
        return ""
    end
    cs = Char[]
    while u > 0
        push!(cs, Char(u & 0xff))
        u >>= 8
    end
    part(cs) = replace(String(reverse(cs)), " "=>"")
    "$(part(cs[4:end])) $(part(cs[1:3]))"
end

function pc_index(db::UPRNDB, dim_positions)
    pcds = Dict{UInt64, Set{LaLoDim}}()
    for uprn in keys(db.grid)
        (;la,lo,pc) =  db.grid[uprn]
        dim = dim_positions[db.uprn2dimension[uprn]]
        if !haskey(pcds, pc)
            pcds[pc] = Set{LaLoDim}([(;la,lo,dim)])
        else
            push!(pcds[pc], (;la,lo,dim))
        end
    end
    pcds
end

Base.write(io::IO, lalodim::LaLoDim) = write(io, lalodim.la) + write(io, lalodim.lo) + write(io, lalodim.dim)
Base.read(io::IO, ::Type{LaLoDim}) = (la=read(io, Float64), lo=read(io, Float64), dim=read(io, UInt64))
Base.write(io::IO, lalodims::Set{LaLoDim}) = reduce((a,lalodim)->a+write(io, lalodim), lalodims, init=write(io, UInt64(length(lalodims))))
Base.read(io::IO, ::Type{Set{LaLoDim}}) = Set{LaLoDim}([read(io, LaLoDim) for i in 1:read(io, UInt64)])

Base.write(io::IO, lalo::LaLo) = write(io, lalo.la) + write(io, lalo.lo)
Base.read(io::IO, ::Type{LaLo}) = (la=read(io, Float64), lo=read(io, Float64))
Base.write(io::IO, lalos::Set{LaLo}) = reduce((a,lalo)->a+write(io, lalo), lalos, init=write(io, UInt64(length(lalos))))
Base.read(io::IO, ::Type{Set{LaLo}}) = Set{LaLo}([read(io, LaLo) for i in 1:read(io, UInt64)])

function write_and_index_pcode_data(io::IO, pcds)
    kvs = Dict{UInt64, Index1024.DataAux}()
    for pc in keys(pcds)
        kvs[pc] = (data=position(io), aux=length(pcds[pc])) # may as well put something in aux
        write(io, pcds[pc])
    end
    kvs
end

function write_dimensions(io::IO, uprndb::UPRNDB)
    fields = sort(collect(keys(uprndb.field2uprn)))
    write(io, length(fields))
    dim_positions = Dict{UInt64, UInt64}()
    for (hsh, dim) in uprndb.dimensions
        dim_positions[hsh] = position(io)
        foreach(f->println(io, dim[f]), fields)
    end
    dim_positions
end

#==
    index_by_postcode("nov_2021.uprndb")
    pdb = ONSUD.open_pcode_index(joinpath(ONSUD.GEODIR, "pcode.db.index"))
    pcode_info("S17 3BB", pdb)

    index_by_postcode("test.uprndb"; pcindexfile="pctest.db.index", datadir="ONSUD_NOV_2021/Test")
    pdb = ONSUD.open_pcode_index(joinpath(ONSUD.GEODIR, "pctest.db.index"))
    pcode_info("S17 3BB", pdb)

    index_by_postcode("bbtest.uprndb"; pcindexfile="bbtest.db.index", datadir="ONSUD_NOV_2021/BB")
    pdb = ONSUD.open_pcode_index(joinpath(ONSUD.GEODIR, "bbtest.db.index"))
    pcode_info("S17 3BB", pdb)
==#

function index_by_postcode(uprndbfile::AbstractString; pcindexfile="pcode.db.index", geodir="", datadir="")
    if geodir == ""
        geodir = ONSUD.GEODIR
    end
    if datadir == ""
        datadir = ONSUD.DATADIR
    end
    uprndb = load(joinpath(geodir, uprndbfile))
    if uprndb === nothing
        uprndb = save(joinpath(geodir, uprndbfile), generate(row_readers(joinpath(geodir, datadir))))
    end
    index_by_postcode(uprndb; pcindexfile, geodir)
end

function index_by_postcode(uprndb::UPRNDB; pcindexfile="pcode.db.index", geodir="")
    if geodir == ""
        geodir = ONSUD.GEODIR
    end
    open(joinpath(geodir, pcindexfile), "w+") do io
        write(io, zero(Int64)) # placeholder for index_pos
        write(io, zero(Int64)) # placeholder for dimension_pos
        dimension_pos = position(io)
        @assert dimension_pos > 0
        dim_positions = write_dimensions(io, uprndb)
        kvs = write_and_index_pcode_data(io, pc_index(uprndb, dim_positions)) # pcode=>(data=offset, aux=dimcount)
        index_pos = position(io)
        @assert index_pos > 0
        build_index_file(io, kvs; meta=map(String, sort(collect(keys(uprndb.field2uprn)))))
        seekstart(io)
        write(io, index_pos)
        write(io, dimension_pos)
    end
end

function open_pcode_index(pcodefn)
    io = open(pcodefn, "r")
    index_pos = read(io, Int64)
    seek(io, index_pos)
    open_index(io)
end

pcode_info(pcode, pcodefn::AbstractString) = pcode_info(pcode, open_pcode_index(pcodefn))

function pcode_info(pcode, pcodedb::Index)
    node = search(pcodedb, postcode_to_UInt64(pcode))
    if node === nothing
        return nothing
    end
    seek(pcodedb.io, node.data)
    lalodims = read(pcodedb.io, Set{LaLoDim})
    dim0 = (;Dict([Symbol(f)=>"" for f in pcodedb.meta])...)
    dims = Dict{Int64, typeof(dim0)}()
    function get_dim(dp)
        if !haskey(dims, dp)
            seek(pcodedb.io, dp)
            dims[dp] = (;Dict{Symbol, String}([Symbol(f)=>readline(pcodedb.io) for f in pcodedb.meta])...)
        end
        dims[dp]
    end
    [(la=lalodim.la, lo=lalodim.lo, get_dim(lalodim.dim)...) for lalodim in lalodims]
end

function pcode_inbound_info(pcodedb::Index, pcode)
    min_p = postcode_to_UInt64("$pcode 1AA")
    max_p = postcode_to_UInt64("$pcode 9ZZ")
    if min_p === nothing && max_p === nothing
        println(stderr, "Inbound postcode not found $pcode")
        return
    end
    nodes = Index1024.node_range(pcodedb, min_p, max_p)
end

function by_uprn!(io, n, lk, kvs)
    readline(io) # header
    while ! eof(io)
        pos = position(io)
        uprn = parse(UInt64, readuntil(io, ","))
        lock(lk)
        try
            kvs[uprn] = (data=pos, aux=n)
        finally
            unlock(lk)
        end
        readline(io)
    end
end

function index_datadir(datadir, index_by!)
    kvs = Dict{UInt64, Index1024.DataAux}()
    lk = ReentrantLock()
    files = readdir(datadir)
    @sync for n in 1:length(files)
        Threads.@spawn open(joinpath(datadir, files[n]), "r") do io
            index_by!(io, n, lk, kvs)
        end
    end
    files, kvs
end

#==
@time ONSUD.create_index1024(ONSUD.DATADIR, "/home/matt/wren/UkGeoData/onsud_nov_2021.uprn.index", ONSUD.by_uprn!)
using BenchmarkTools
@benchmark ONSUD.uprn_data(ONSUD.DATADIR, "/home/matt/wren/UkGeoData/onsud_nov_2021.uprn.index", 10015278860)
==#

function create_index1024(datadir, indexfile, index_by!)
    meta, kvs = index_datadir(datadir, index_by!)
    build_index_file(indexfile, kvs; meta)
end

uprn_data(datadir::AbstractString, indexfile::AbstractString, uprn) = uprn_data(datadir, open_index(indexfile), uprn)

function csv(io::IO, offset)
    buff = IOBuffer()
    write(buff, readline(io; keep=true))
    seek(io, offset)
    write(buff, readline(io; keep=true))
    seekstart(buff)
    CSV.File(buff)
end

function uprn_data(datadir, idx::Index, uprn)
    node = search(idx, uprn)
    if node === nothing
        return nothing
    end
    open(joinpath(datadir, idx.meta[node.aux])) do io
        csv(io, node.data)
    end
end

###
end


text:
using UUIDs
using Pkg

const JULIA_UUID = UUID("1222c4b2-2114-5bfd-aeef-88e4692bbb3e")

const stdlib_to_uuid = Dict(
    "Dates" => UUID("ade2ca70-3891-5945-98fb-dc099432e06a"),
    "FileWatching" => UUID("7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"),
    "SHA" => UUID("ea8e919c-243c-51af-8825-aaa63cd721ce"),
    "Sockets" => UUID("6462fe0b-24de-5631-8697-dd941f90decc"),
    "DelimitedFiles" => UUID("8bb1440f-4735-579b-a4ab-409b98df4dab"),
    "Random" => UUID("9a3f8284-a2c9-5f02-9a11-845980a1fd5c"),
    "Distributed" => UUID("8ba89e20-285c-5b6f-9357-94700520ee1b"),
    "LinearAlgebra" => UUID("37e2e46d-f89d-539d-b4ee-838fcccc9c8e"),
    "Mmap" => UUID("a63ad114-7e13-5084-954f-fe012c677804"),
    "SparseArrays" => UUID("2f01184e-e22b-5df5-ae63-d93ebab69eaf"),
    "UUIDs" => UUID("cf7118a7-6976-5b1a-9a39-7adc72f591a4"),
    "Logging" => UUID("56ddb016-857b-54e1-b83d-db4d58db5568"),
    "Statistics" => UUID("10745b16-79ce-11e8-11f9-7d13ad32a3b2"),
    "CRC32c" => UUID("8bf52ea8-c179-5cab-976a-9e18b702a9bc"),
    "Pkg" => UUID("44cfe95a-1eb2-52ea-b672-e2afdf69b78f"),
    "Serialization" => UUID("9e88b42a-f829-5b0c-bbe9-9e923198166b"),
    "Markdown" => UUID("d6f4376e-aef5-505a-96c1-9c027394607a"),
    "SuiteSparse" => UUID("4607b0f0-06f3-5cda-b6b1-a6196a1729e9"),
    "REPL" => UUID("3fa0cd96-eef1-5676-8a61-b3b8758bbffb"),
    "Libdl" => UUID("8f399da3-3557-5675-b5ff-fb832c97cbdb"),
    "Base64" => UUID("2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"),
    "Profile" => UUID("9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"),
    "InteractiveUtils" => UUID("b77e0a4c-d291-57a0-90e8-8db25a27a240"),
    "SharedArrays" => UUID("1a1011a3-84de-559e-8e89-a11a2f7dc383"),
    "Printf" => UUID("de0858da-6303-5e67-8744-51eddeeeb8d7"),
    "Future" => UUID("9fa8497b-333b-5362-9e8d-4d0656e87820"),
    "LibGit2" => UUID("76f85450-5226-5b5a-8eaa-529ad045b433"),
    "Test" => UUID("8dfed614-e22c-5e08-85e1-65c5234f0b40"),
    "Unicode" => UUID("4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"),
    "julia" => JULIA_UUID
)

const stdlib_uuids = Set([
    "ade2ca70-3891-5945-98fb-dc099432e06a",
    "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee",
    "ea8e919c-243c-51af-8825-aaa63cd721ce",
    "6462fe0b-24de-5631-8697-dd941f90decc",
    "8bb1440f-4735-579b-a4ab-409b98df4dab",
    "9a3f8284-a2c9-5f02-9a11-845980a1fd5c",
    "8ba89e20-285c-5b6f-9357-94700520ee1b",
    "37e2e46d-f89d-539d-b4ee-838fcccc9c8e",
    "a63ad114-7e13-5084-954f-fe012c677804",
    "2f01184e-e22b-5df5-ae63-d93ebab69eaf",
    "cf7118a7-6976-5b1a-9a39-7adc72f591a4",
    "56ddb016-857b-54e1-b83d-db4d58db5568",
    "10745b16-79ce-11e8-11f9-7d13ad32a3b2",
    "8bf52ea8-c179-5cab-976a-9e18b702a9bc",
    "44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
    "9e88b42a-f829-5b0c-bbe9-9e923198166b",
    "d6f4376e-aef5-505a-96c1-9c027394607a",
    "4607b0f0-06f3-5cda-b6b1-a6196a1729e9",
    "3fa0cd96-eef1-5676-8a61-b3b8758bbffb",
    "8f399da3-3557-5675-b5ff-fb832c97cbdb",
    "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f",
    "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79",
    "b77e0a4c-d291-57a0-90e8-8db25a27a240",
    "1a1011a3-84de-559e-8e89-a11a2f7dc383",
    "de0858da-6303-5e67-8744-51eddeeeb8d7",
    "9fa8497b-333b-5362-9e8d-4d0656e87820",
    "76f85450-5226-5b5a-8eaa-529ad045b433",
    "8dfed614-e22c-5e08-85e1-65c5234f0b40",
    "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
])

is_stdlib(uuid) = uuid in stdlib_uuids
is_jll(name) = endswith(name, "_jll")

function build_uuid_name_map(version = VERSION; registry=joinpath(homedir(), ".julia/registries/General"))
    allpkgs = installable_on_version(version, registry=registry)
    name_to_uuid = Dict()
    for (uuid, pkg) in allpkgs
        name_to_uuid[pkg.name] = UUID(uuid)
    end
    merge!(name_to_uuid, stdlib_to_uuid)
    name_to_uuid
end

"""
    dependencies_per_package(registry=joinpath(homedir(), ".julia/registries/General"))

Find all declared (direct) dependencies for each package in `registry`.
"""
function dependencies_per_package(registry=joinpath(homedir(), ".julia/registries/General"), depmap = Dict())
    name_uuid_map = build_uuid_name_map(registry = registry)

    for dir in readdir(registry)
        startswith(dir, ".") && continue

        dirpath = joinpath(registry, dir)

        isdir(dirpath) || continue

        for pkg in readdir(dirpath)
            startswith(dir, ".") && continue

            pkgpath = joinpath(dirpath, pkg)

            isdir(pkgpath) || continue

            pkgtoml = Pkg.TOML.parsefile(joinpath(pkgpath, "Package.toml"))
            versionstoml = Pkg.TOML.parsefile(joinpath(pkgpath, "Versions.toml"))

            depsperversion = Dict{String, Dict}()

            if isfile(joinpath(pkgpath, "Deps.toml"))
                depstoml = Pkg.TOML.parsefile(joinpath(pkgpath, "Deps.toml"))
            else
                depstoml = Dict()
            end

            if isfile(joinpath(pkgpath, "Compat.toml"))
                compattoml = Pkg.TOML.parsefile(joinpath(pkgpath, "Compat.toml"))
            else
                compattoml = Dict()
            end

            for version in keys(versionstoml)
                deps = Dict{String, Dict}()

                for depsver in keys(depstoml)
                    if VersionNumber(version) in Pkg.Types.VersionRange(depsver)
                        for (dep, uuid) in depstoml[depsver]
                            deps[dep] = Dict{String, Any}(
                                "name" => dep,
                                "uuid" => uuid,
                                "is_stdlib" => is_stdlib(uuid),
                                "is_jll" => is_jll(dep),
                                "slug" => Base.package_slug(UUID(uuid), 5),
                                "versions" => "*"
                            )

                        end
                    end
                end

                for compatver in keys(compattoml)
                    if VersionNumber(version) in Pkg.Types.VersionRange(compatver)
                        for (dep, vers) in compattoml[compatver]
                            if !haskey(name_uuid_map, dep)
                                @error "UUID not found for" dep
                                continue
                            end
                            depdict = get!(deps, dep) do
                                uuid = string(name_uuid_map[dep])
                                Dict{String, Any}(
                                    "name" => dep,
                                    "uuid" => uuid,
                                    "is_stdlib" => is_stdlib(uuid),
                                    "is_jll" => is_jll(dep),
                                    "slug" => Base.package_slug(UUID(uuid), 5),
                                )
                            end

                            depdict["versions"] = vers
                        end
                    end
                end

                depsperversion[version] = deps
            end

            name = pkgtoml["name"]
            uuid = pkgtoml["uuid"]

            depmap[uuid] = Dict(
                "uuid" => uuid,
                "name" => name,
                "slug" => Base.package_slug(UUID(uuid), 5),
                "deps" => depsperversion
            )
        end
    end
    return depmap
end

function dependencies_per_package(registries::Vector)
    deps = Dict()
    for reg in registries
        dependencies_per_package(reg, deps)
    end
    return deps
end

function reverse_dependencies_per_package(deps_per_pkg)
    reversedeps = Dict()

    for (uuid, d) in deps_per_pkg
        for (ver, deps) in d["deps"]
            for (dep, depdict) in deps
                depuuid = depdict["uuid"]

                rdeps = get!(reversedeps, depuuid, Dict())

                push!(get!(rdeps, depdict["versions"], Set([])), Dict(
                    "uuid" => uuid,
                    "name" => d["name"],
                    "slug" => Base.package_slug(UUID(uuid), 5),
                    "version" => ver
                ))
            end
        end
    end

    return reversedeps
end

function _alldeps(uuid, version, deps_per_pkg, deps, seen = Set([]), isdirect=true; directonly = false)
    uuid in seen && return deps
    haskey(deps_per_pkg, uuid) || return deps

    push!(seen, uuid)

    for (dep, depdict) in get(get(deps_per_pkg[uuid], "deps", Dict()), version, [])
        depuuid = depdict["uuid"]

        depentry = get!(deps, (depuuid, isdirect), Dict(
            "uuid" => depuuid,
            "name" => depdict["name"],
            "slug" => depdict["slug"],
            "direct" => isdirect,
            "versions" => vcat(depdict["versions"])
        ))

        sort!(unique!(append!(depentry["versions"], vcat(depdict["versions"]))))

        if !directonly && haskey(deps_per_pkg, depuuid)
            versions = (filter(collect(keys(deps_per_pkg[depuuid]["deps"]))) do ver
                VersionNumber(ver) in Pkg.Types.VersionSpec(depdict["versions"])
            end)
            isempty(versions) && continue
            ver = last(versions)
            _alldeps(depuuid, ver, deps_per_pkg, deps, seen, false)
        end
    end
end

function alldeps(uuid, version, deps_per_pkg; directonly = false)
    deps = Dict()
    _alldeps(uuid, version, deps_per_pkg, deps; directonly = directonly)

    sort!(collect(values(deps)), by = x -> x["name"])
end

directdeps(uuid, version, deps_per_pkg) = alldeps(uuid, version, deps_per_pkg; directonly = true)

function _allreversedeps(uuid, version, reversedeps, rdeps, seen = Set([]), isdirect=true; directonly = false)
    uuid in seen && return rdeps
    haskey(reversedeps, uuid) || return rdeps

    push!(seen, uuid)

    for (versionrange, deps) in reversedeps[uuid]
        if VersionNumber(version) in Pkg.Types.VersionSpec(versionrange)
            for dep in deps
                depentry = get!(rdeps, (dep["uuid"], isdirect), merge(dep, Dict("direct" => isdirect)))
                depentry["version"] = unique!(vcat(depentry["version"], dep["version"]))

                if !directonly
                    _allreversedeps(dep["uuid"], dep["version"], reversedeps, rdeps, seen, false)
                end
            end
        end
    end
end

function allreversedeps(uuid, version, reversedeps; directonly = false)
    rdeps = Dict()
    _allreversedeps(uuid, version, reversedeps, rdeps; directonly = directonly)
    rdeps = sort!(collect(values(rdeps)), by = x -> x["name"])

    foreach(rdeps) do dep
        sort!(dep["version"], by=VersionNumber)
    end
    rdeps
end

directreversedeps(uuid, version, reversedeps) = allreversedeps(uuid, version, reversedeps; directonly = true)

"""
    installable_on_version(version = VERSION; registry=joinpath(homedir(), ".julia/registries/General"))

Returns a Dict mapping from `uuid` to named tuples containing `(name, url, uuid, path, versions)`
of packages in `registry` compatible with Julia version `version`.
"""
function installable_on_version(version = VERSION; registry=joinpath(homedir(), ".julia/registries/General"))
    allpkgs = Dict()
    for initial in filter!(isdir, joinpath.(registry, readdir(registry)))
        for pkg in filter!(isdir, joinpath.(registry, initial, readdir(initial)))
            "Compat.toml" in readdir(pkg) || continue
            pkgtoml = Pkg.TOML.parsefile(joinpath(pkg, "Package.toml"))
            versions = Pkg.TOML.parsefile(joinpath(pkg, "Versions.toml"))
            compat = Pkg.TOML.parsefile(joinpath(pkg, "Compat.toml"))
            compatible_versions = VersionNumber[]
            for pkgver in keys(compat)
                try
                    if haskey(compat[pkgver], "julia")
                        if any(in.(version, Pkg.Types.VersionRange.(compat[pkgver]["julia"])))
                            append!(compatible_versions, [
                                                VersionNumber(v) for v in keys(versions) if
                                                    VersionNumber(v) in Pkg.Types.VersionRange(pkgver)
                                               ])
                        end
                    else
                        append!(compatible_versions, [
                                            VersionNumber(v) for v in keys(versions) if
                                                VersionNumber(v) in Pkg.Types.VersionRange(pkgver)
                                           ])
                    end
                catch err
                    @error err
                end
            end
            allpkgs[pkgtoml["uuid"]] = (
                name = pkgtoml["name"],
                url = pkgtoml["repo"],
                uuid = pkgtoml["uuid"],
                path = pkg,
                versions = compatible_versions
            )
        end
    end
    allpkgs
end


text:
@reexport module ProgressiveHedging

# Standard library
using LinearAlgebra
using SparseArrays
using Distributed
using Printf

# External libraries
using Parameters
using JuMP
using MathOptInterface
using StochasticPrograms
using StochasticPrograms: AcceptableTermination
using StochasticPrograms: UnspecifiedInstantiation, HorizontalStructure, AbstractScenarioProblems, ScenarioProblems, DistributedScenarioProblems
using StochasticPrograms: AbstractExecution, ExecutionParameter, Serial, Synchronous, Asynchronous
using StochasticPrograms: AbstractStructuredOptimizer, RelativeTolerance, SubproblemOptimizer
using StochasticPrograms: get_decisions, set_known_decision!, SingleKnownSet
using StochasticPrograms: add_subscript
using StochasticPrograms: AbstractPenaltyterm, Quadratic, InfNorm, ManhattanNorm, initialize_penaltyterm!, update_penaltyterm!, remove_penalty!
using ProgressMeter

import Base: show, put!, wait, isready, take!, fetch
import StochasticPrograms: supports_structure, default_structure, load_structure!, restore_structure!, optimize!, optimizer_name, master_optimizer, subproblem_optimizer, num_subproblems

const MOI = MathOptInterface
const MOIU = MOI.Utilities

export
    AbstractProgressiveHedgingAttribute,
    ProgressiveHedgingAlgorithm,
    PrimalTolerance,
    DualTolerance,
    Penalizer,
    PenalizationParameter,
    RawPenalizationParameter,
    set_penalization_attribute,
    set_penalization_attributes,
    Penaltyterm,
    Fixed,
    Adaptive

# Include files
include("types/types.jl")
include("penalties/penalization.jl")
include("execution/execution.jl")
include("solver.jl")
include("MOI_wrapper.jl")

end # module


text:
module LocalSearchSolvers

using Dictionaries: include
using CompositionalNetworks: include
using Constraints: include
using Base.Threads
using CompositionalNetworks
using ConstraintDomains
using Constraints
using Dictionaries
using Distributed
using JSON
using Lazy

# Exports internal
export constraint!, variable!, objective!, add!, add_var_to_cons!, add_value!
export delete_value!, delete_var_from_cons!, domain, variable, constraint, objective
export length_var, length_cons, constriction, draw, describe, get_variable
export get_variables, get_constraint, get_constraints, get_objective, get_objectives
export get_cons_from_var, get_vars_from_cons, get_domain, get_name, solution

# export for CBLS.jl
export _set_domain!, is_sat, max_domains_size, get_value, update_domain!, has_solution
export set_option!, get_option, sense, sense!

# Exports Model
export model

# Exports error/predicate/objective functions
export o_dist_extrema, o_mincut

# Exports Solver
export solver, solve!, specialize, specialize!, Options, get_values, best_values
export best_value, time_info, status

# Include utils
include("utils.jl")

# Include model related files
include("variable.jl")
include("constraint.jl")
include("objective.jl")
include("model.jl")

# Include solver state and pool of configurations related files
include("configuration.jl")
include("pool.jl")
include("fluct.jl")
include("state.jl")

# Include strategies
include("strategies/move.jl")
include("strategies/neighbor.jl")
include("strategies/objective.jl")
include("strategies/parallel.jl")
include("strategies/perturbation.jl")
include("strategies/portfolio.jl")
include("strategies/tabu.jl") # preceed restart.jl
include("strategies/restart.jl")
include("strategies/selection.jl")
include("strategies/solution.jl")
include("strategies/termination.jl")
include("strategy.jl") # meta strategy methods and structures

# Include solvers
include("options.jl")
include("time_stamps.jl")
include("solver.jl")
include("solvers/sub.jl")
include("solvers/meta.jl")
include("solvers/lead.jl")
include("solvers/main.jl")

# Include usual objectives
include("objectives/extrema.jl")
include("objectives/cut.jl")

end


text:
# Copyright Robin Harper 2019-2020

module PEEL

using LsqFit, Hadamard, DelimitedFiles, LinearAlgebra
using PyPlot, Statistics, SparseArrays
import Base.show

# Defines a helper function in jupyter notebooks its /oplus<tab>
⊕ = (x,y)->mod.(x+y,2)

export peelBack,checkAndExtractSingleton,constructTheKeyIndexMUB,getIndexOf
export generateFromPVecSamples4N,twoPattern,fourPattern,simplexUp,generatePauliFromClifford
export getχ,bitArrayToNumber

function getχ(a,b)
    return count_ones(a&b) %2 == 0 ? +1 : -1
end

function bitArrayToNumber(toS)
    return foldl((x,y)->x<<1+y,toS)
end

# The way I ended up doing it in the FULL MUB is probably better, but this is easy to follow

"""
fourPattern(p)

## Arguments
    p::Array{Array{Array{Int64,1},1},1} represents a series of commuting two qubit paulis.

Given a MUB works out what binary string will go in what bin
This takes a 4 bit (two Pauli) set of commuting Paulis, eg.

```
1-element Array{Array{Array{Int64,1},1},1}:
 [[0, 0, 0, 0], [1, 0, 1, 1], [0, 1, 1, 0], [1, 1, 0, 1]]
```

And then returns the Array that shows what bin each of the 16 two qubit Paulis go into e.g.

```
[[["0000", "0111", "1001", "1110"], ["0001", "0110", "1000", "1111"], ["0011", "0100", "1010", "1101"], ["0010", "0101", "1011", "1100"]]]
```

Which tells us e.g. that "1001" goes into bin 0 and "1101" goes into bin 2.

## Returns
    an Array, containing a single array of 4 vectors each of size 4 --- splitting all 16 Paulis (rep as binary strings) into 4 bins.

If you want to think of this as a hashing function, it shows how to take 4 bits -> 2 bits using the supplied pattern.
"""
function fourPattern(p)
    allPatterns=[]
    for x in p
        patterns =[Hadamard.hadamard(16)[x+1,:] for x in [foldl((x,y)->x<<1+y,t) for t in x]]
        push!(allPatterns,
            map(y->reverse.(string.(findall(x->x!=0,y).-1,base=2,pad=4)),Hadamard.hadamard(4)*patterns))
    end
    return allPatterns
end

"""
twoPattern(p)

### Arguments
     p Array{Array{Int64,1},1} represents a series of commuting two qubit paulis.

Given a MUB works out what binary string will go in what bin
This takes a 2 bit (one Pauli) set of commuting Paulis, eg.

```
1-element Array{Array{Array{Int64,1},1},1}:
 [0, 0], [1, 1]
```

And then returns a two element Array that shows what bin each of the 16 two qubit Paulis go into e.g.

```
["00", "11"]
["01", "10"]
```

Which tells us e.g. that "11" goes into bin 0 and "01" goes into bin 2.

## Returns
    an Array, containing an array of 2 vectors each of size 2 --- splitting all 4 Paulis (rep as binary strings) into 2 bins.

If you want to think of this as a hashing function, it shows how to take 2 bits -> 1 bits using the supplied pattern.
"""
function twoPattern(ps)
    patterns =[Hadamard.hadamard(4)[x+1,:] for x in [foldl((x,y)->x<<1+y,t) for t in ps]]
    return map(y->reverse.(string.(findall(x->x!=0,y).-1,base=2,pad=2)),Hadamard.hadamard(2)*patterns)
end

"""
getIndexOf(hadamardMap,no,dictMap=Dict())
### Arguments
	hadamardMap: An array of hashing bins (see later)
	no: Integer  the index of the Pauli we want to know which bin it gotes into
	dictMap: Dictionary optional mapping eg qubit 2->qubit 4.

 Given a map (constructed from fourPattern and twoPattern), this will
 split the qubits up into the pairs/singles identified by the map and
 then determine the hash destingation of a particular Pauli, i.e. which bin it goes into.
 The optional dictMapping, allows us to map qubits to different qubits. This might be
 relevant for different topologies. E.g. we might want the pair to span qubits 1 and 14, rather
 than 1 and 2, so we would map 2 onto 14. There should be examples later.

 **NOTE** because we use the least significant (rh) qubit - there is a subte reversal of bits. See example.

### Example

The hadamard Map for a six qubit system might look like this
```
4-element Array{Any,1}:
 [["00", "10"], ["01", "11"]]
 [["0000", "0110", "1101", "1011"], ["1100", "1010", "0001", "0111"], ["1000", "1110", "0101", "0011"], ["0100", "0010", "1001", "1111"]]
 [["0000", "1110", "1001", "0111"], ["1100", "0010", "0101", "1011"], ["0100", "1010", "1101", "0011"], ["1000", "0110", "0001", "1111"]]
 [["00", "01"], ["10", "11"]]
```

 Where we have used twoPattern for qubits 0 and 5 and fourPattern for qubits 1&2 and 3&4

 getIndexOf(patternAbove, 23), would then return 6.

 Because 23 = "000000010111"
 Then we map the bins as
 00-> reversed to 00 -> bin 0 -> 0
 0000 -> reversed to 0000 ->  bin 0 -> 00
 0101 -> reversed to 1010 ->  bin 2 -> 10
 11 ->  reversed to 11 -> bin 1 -> 01

 Therefore the index is 0 00 10 01 = 5, as we one index we return 6.

 ## Returns

 Index of bin (1 indexed).

"""
function getIndexOf(hadamardMap,no,dictMapping=Dict())
    bitsize = sum([length(x)>1 ? length(x) : floor(Int,log2(x[1])) for x in hadamardMap]) # 2* number of qubits
    bitst = (string(no,base=2,pad=bitsize))
    if length(dictMapping)!=0
        currentValues = [parse(Int,x) for x in bitst]
        newB = copy(currentValues)
        for i = 1:length(dictMapping)
               currentValues[dictMapping[i]*2-1] = newB[i*2-1]
               currentValues[dictMapping[i]*2] = newB[i*2]
        end
        bitst = string(foldl((x,y)->x<<1+y,currentValues),base=2,pad=bitsize)
    end
    bitsize = sum([length(x)>1 ? length(x) : floor(Int,log2(x[1])) for x in hadamardMap]) # 2* number of qubits
    #lshift = round(Int,bitsize/2)-1
    occupancy = [length(x)>1 ? length(x) : x[1] for x in hadamardMap]
    cpos = 1
    indexArray = []
    for hm in hadamardMap
        if (length(hm)>1) # We assume that we have a list of MUBS
            bits = bitst[cpos:cpos+length(hm)-1]
            mv = findfirst(x->x==true,in.(reverse(bits),hm))-1
            push!(indexArray,mv)
            cpos = cpos+length(hm)
        else # we have a binary watchamacallit
            numberOfBits = floor(Int,log2(hm[1]))
            bits = bitst[cpos:cpos+numberOfBits-1]
            push!(indexArray,parse(Int,bits,base=2))
            cpos = cpos+numberOfBits
        end

    end

    #print(indexArray)

    val = indexArray[1]
    hmLengths = [floor(Int,log2(t)) for t in occupancy]
    st = ""
    for ix in 1:length(hmLengths)
        st = st*string(indexArray[ix],base=2,pad=hmLengths[ix])
    end
    return parse(Int,st,base=2)+1
end

"""
constructTheKeyIndexMUB(indx,listOfPs,map=[])

## Arguments
    indx The index of a Pauli, we want to ascertain which bins that Pauli was hashed into.
	listOfPs The used to sample/bin the Paulis
	map A map of qubit substitutions e.g mapping qubit 2 to 14

Given a Paulis to test, returns the index into all the bins that Paulis is hashed into
(it is assumed we have a number a sub-sampling groups).
Basically just iterates over them. If we supply a map, it should be an
array of maps, one for each different set of n-qubit MUBS supplied.

It uses getIndexOf to do all the heavy lifting, the documentation for that function
provides an example of the type of map it requires, this is an array of these maps.

Returns and array of offsets representing the index the Pauli gets hashed to.

"""
function constructTheKeyIndexMUB(indx,listOfPs;map=[])
    if map == []
        return [getIndexOf(x,indx) for x in listOfPs]
    else
        return [getIndexOf(x,indx,map[ix]) for (ix,x) in enumerate(listOfPs)]
    end
end

"""
generateFromPVecSamples4N(pvec,d=[],dictMapping)

## Arguments
	pvec The list of Paulis which are used in our experiment. Given as an array of stabilisers (see example)
	d The offset, this is the binary string 'hashed' into the Paulis given by pvec.
	map an optional dictionary mapping qubits to other qubits.

Generates the fidelities we are going to sample from for a given offset
vector (d).

By way of example, if we supplied as a pvec

```
[[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 0, 1], [1, 1, 1, 0]]

```

 Then we would get 0,7,9 and 14.

 If we also suplied a d, [1 0 0 0], then we get 8, 15, 1 6

 A pvec as

```
 [[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 0, 1], [1, 1, 1, 0]]
 [[0, 0, 0, 0], [0, 1, 1, 1], [1, 0, 0, 1], [1, 1, 1, 0]]
```

 Would give us 0, 7, 9, 14, 112, 119, 121, 126, 144, 151, 153, 158, 224, 231, 233, 238

 ### Returns
 	Array{Int64,1} with the relevant fidelity indexes. (zero indexed).

"""
function generateFromPVecSamples4N(pvec,d=[],dictMapping=Dict())
    subSamples = []
    if d == []
        d = zeros(Int64,sum([length(x[1]) for x in pvec]))
    end

    maxValue = foldl(*,[length(x) for x in pvec])
    occupancy = vcat(map(x->length(x),pvec),1)
    # Tuple of what we need to divide by and what we need to mod by
    # So for example if the lengths given were [4,8,4]
    # Then we need to to range from 0:(4*8*4=256)-1
    # We need to divide our 'c' by (8*4)%4 and (4)%8 and (1)%4 to get the three indexes needed.
    divs = [(foldl(*,occupancy[ix:end]),occupancy[ix-1]) for ix=length(occupancy):-1:2]
    for c = 0:maxValue-1
        # express c as a four  array of length 1
        iv = reverse(map(x->floor(Int,c/x[1])%x[2]+1,divs))
#           print("$iv\n")
        xs_o = foldl(vcat,map(x->pvec[x[1]][x[2]],enumerate(iv)))
        xs = copy(xs_o)
        if length(dictMapping) != 0
           for i = 1:length(dictMapping)
               xs[dictMapping[i]*2-1] = xs_o[i*2-1]
               xs[dictMapping[i]*2] = xs_o[i*2]
           end
        else
                xs = xs_o
        end
            #println(xs)
        nxs = xs ⊕ d
            #print("$(itoP(foldl((x,y)->x<<1+y,nxs)))\n")
        push!(subSamples,foldl((x,y)->x<<1+y,nxs))
    end
    return [x for x in subSamples]
end

"""
closeToZero(bin,BITS,cutoff=1e-15)

## Arguments
	bin An array of values in the bin and its various offsets.
	bits The number of bits (typically 2*#ofqubits)
	cutoff the entropy in the bin, below which we will assume it is zero.

Checks if the values in the bin have an absolute value close to zero, where we sum the
squares of the entries in the bin, divide by (BITS+1) and compare to the cutoff.

## Returns
	true or false.
"""
function closeToZero(indexes,BITS;cutoff=1e-15)
    # checks to see if the supplied [U[k̂] for randomDs] are close to zero. i.e. no weights
    # Just now , this is effectively ||U||^2
    U = [ix^2 for ix in indexes]
    if (sum(U)/(BITS+1) < cutoff)
        return true
    else
        return false
    end
end

"""
checkAndExtractSingleton(indexes,BITS;cutoff=0.000002)

## Arguments
	indexes: A list of arrays, being the bins and their offsets.
	Bits: the number of bits (being 2*qubits)
	cutoff: what the entropy has to be less than
Checks the 'entropy' in the bin (ie sum of the square of the fluctations around the mean of the abs value)/BITS
If it is less than the cutoff admits it as a singleton.
Extracts the Pauli (using majority vote if given multiple samples)
-  Note if you wanted to implement the "Error Correcting Code in the offsets" as described in the appendix of the paper then this funciton would be different.

## Returns
    A tuple of: (found,index, value), where found is boolean, the index is an integer and value is float (between 0 and 1).

"""
function checkAndExtractSingleton(indexes,BITS;cutoff=0.000002)
    stringsFound = zeros(BITS)
    threshold = floor((length(indexes)-1)/2)
    for (ii,i) in enumerate(indexes)
        #print("Check and Extract $ii : $i\n")
        pcheck = i[1]
        toAdd = ""
        for (ix,xx) in enumerate(i[2:end])
            if pcheck >=0
                if xx >=0
                    stringsFound[ix]+=1
                end
            else
                if xx < 0
                    stringsFound[ix]+=1
                end
            end
        end
    end
    #print("$stringsFound\n")
    kstr =  [x > threshold ? "0" : "1" for x in stringsFound]
    kii = foldl(*,kstr)
    kval = parse(Int,kii,base=2)
    # so our estimate of this particular string should be indexes[1][1]
    Uk̂ = indexes[1][1]
    Uk̂bar = mean(mean([abs.(x) for x in indexes]))
    entropy = sum([sum([(abs.(x)-Uk̂)^2 for x in ix])/(BITS+1) for ix in indexes])
    #print(entropy)
    if entropy < cutoff
       # print("Singleton $(Uk̂) and $(mean(mean([abs.(x) for x in extracted])))")
        return (true,kii,Uk̂bar)
    else
        return (false,0,0)
    end

end

"""
peelBack(listOfX,listOfPs,singletonBits,singletonValue,found,ds,mappings)

## Arguments
	listOfX - the bins and their offsets <- this gets modified
	listOfPs - the stabiliser groups used to create the bins (see generateFromPVecSamples4N)
	singletonBits - the bit index of the Pauli to peel back
	singletonVlaue - the value of the Pauli to peel back
	found - a Dict of Pauli bitstrings that have been found and their value
	ds - the offsets used.
	mappings - an optional mapping of qubit->qubit.

Peel back algorithm tailored for noise version.

Check if we have already found the supplied Pauli (is it in found?)
-   If so check if the value we think we have found is greater than the previous sample
    -   If not, then just return - we had already found it, and the noise has just confused us!
    -   If yes, then the one we found was probably wrong and a result of the noise, just update the value.

*[Query this logic, it seems to work but perhaps, we may want to update the bins with the difference.]*

- Otherwise (i.e. we hadn't found it before.)
    -    Get the index of that Pauli into each of the other bin sets
    -    And add or subtract (depending on relevant offset) the supplied value of the Pauli
    -    Update found by adding in this Pauli and the value.
"""
function peelBack(listOfX,listOfPs,singletonBits,singletonValue,found,ds,mappings)
    if (singletonBits in keys(found))
        #print("We already done that one!\n")
        vval = parse(Int,singletonBits,base=2)
        vstring = singletonBits
        if singletonValue > mean(found[vstring])
                found[vstring] = [singletonValue]
        end
        return
        #@assert false
    end
    vval = parse(Int,singletonBits,base=2)
    vstring = singletonBits
    found[vstring] = vcat(get(found,vstring,[]),[singletonValue])
    # Get the keys/bin number for all the other hash/MUBs we have.
    theKeys = constructTheKeyIndexMUB(vval,listOfPs,map=mappings)
    for dictIndex in 1:length(theKeys)
        pv = theKeys[dictIndex]
        #print("\t****** Stepping through left to right nodes (so to speak) for identified singleton\n")
        extractValue = pv#parse(Int,pv,base=2)+1
        extracted = []
            push!(extracted,[x[extractValue] for x in listOfX[dictIndex]])
        newIndex = []
        for (sidx,bd) in enumerate(ds)
               push!(newIndex,extracted[1][sidx] + -1*PEEL.getχ(PEEL.bitArrayToNumber(bd),vval)*singletonValue)
        end
         theX = listOfX[dictIndex]
         for (ix1) in 1:length(theX)
             theX[ix1][extractValue] = newIndex[ix1]
         end
    end
end

# ======================================================================================
# The following are just convience functions for giving us labels we can use
# ======================================================================================

"""
probabilityLabels(index; qubits=2)

## Arguments
    index (Int) The index (0 index) for which we want the label.
    qubits (Int) The number of qubits in the label, defaults to 2
Because we use the Walsh-Hadmard transform based on bits (i.e. the standard 'natural' transform rather than a Pauli-commuting transform)
The transform changes the way we need to label our qubits. In the probability regime we have the following bit labelling:
```
	00-> I
	01-> Y
	10-> X
	11-> Z
```
The first qubit is the least signifcant number (i.e. in binary would appear on the right hand side).
So the number 1001 is XY, where qubit 0 = Y.

### Returns
The probability label for the index.

"""
function probabilityLabels(x;qubits=2)
    str = string(x,base=4,pad=qubits)
    paulis = ['I','Y','X','Z']
    return map(x->paulis[parse(Int,x)+1],str)
end

"""
fidelityLabels(index; qubits=2)

## Arguments
    index (Int) The index (0 index) for which we want the label.
    qubits (Int) The number of qubits in the label, defaults to 2
Because we use the Walsh-Hadmard transform based on bits (i.e. the standard 'natural' transform rather than a Pauli-commuting transform)
The transform changes the way we need to label our qubits. In the probability regime we have the following bit labelling:
```
	00-> I
	01-> X
	10-> Y
	11-> Z
```
The first qubit is the least signifcant number (i.e. in binary would appear on the right hand side).
So the number 1001 is YX, where qubit 0 = X.

### Returns
The fideltiy label for the index.

"""
function fidelityLabels(x;qubits=2)
    str = string(x,base=4,pad=qubits)
    paulis = ['I','X','Y','Z']
    return map(x->paulis[parse(Int,x)+1],str)
end

# ======================================================================================
# The following functions are to generate a Pauli probability from a "Clifford" probability
# ======================================================================================

"""
simplexUp(n)

## Arguments

   `n (Int)`  = number to project onto

Projects onto the n-simplex. i.e. returns n numbers such that the sum of them is 1.

## Returns

vector of n-numbers that add up to 1.

"""
function simplexUp(n)
    return diff(vcat(0,sort(rand(n-1)),1))
end

"""
noOfOnesIn(ix)

## Arguments

	`ix (Int)` A number representing an index in Julia

Return the number of '1' in the binary form of the number -1.
Note that we subtract 1 from the supplied index to take into account Juila's 1 indexing.
The 'index' 1 actually represents the number 0, which has 0 '1's in its binary form.
The indexes 2 and 3, both have 1 index in their binary form (e.g 01 and 10), index 4
(which is binary 11) has 2 and so on.

	Return the number of 1s in the binary representation of ix-1.
"""
function noOfOnesIn(ix)
    return length(filter(x->x=='1',string(ix-1,base=2)))
end

"""
gen(n)

## Arguments
     `n (Int) ` Number of qubits

Calculates and returns all possible (non-identity) Paulis  for a certain system size. The Paulis
are represented by the numbers 1,2 and 3.

### Returns

A vector of the Paulis.
"""
function gen(n)
    if n == 1
        return [1,2,3]
    else
        x = gen(n-1)
        return vcat(map(y->vcat([1],y),x),map(y->vcat([2],y),x),map(y->vcat([3],y),x))
    end
end

"""
splitUp(n)

## Arguments
    `n Integer`: An 'index' into the Cliffo

Splits the "Clifford average" into the Paulis that were averaged into it.
The starting point here is that if we assume we had single-Clifford twirls on an n-qubit machine
We have our noise supoeroperator, diagonalised into a superoperator matrix, Pauli basis, with 2^n distinct
entries (4^n entries in total).  If we imagine we labelled these distinct entries sequentially then the first entry would represent the 'all I' eigenvalue.
(These distinct entries are in fact the eigenvalues we recover from a Hadamard transform of the measurement outcomes, where we measure in the computational basis).
Next there would then be three averaged eigenvalues representing the I...IX, I...IY, I....IZ eigenvalues and so on. What we are interested in is
splitting that up to retrieve the eigenvalues that were averaged to give that number.

For example, the no-error starts IIIII.IIII = index 1

-	splitUp(1) = 1
-	splitUp(2) (Z) = 2,3,4 i.e. IX, IY, IZ
-	splitUp(3)  (ZI) = 5,9,13 ie. XI, YI, ZI
-	splitUp(4) (ZZ) = 6,7,8,10,11,12,14,15,16 i.e. XX,XY, XZ, YX,YY,YZ, ZX, ZY, ZZ

## Returns
A vector of the Paulis that would have been averaged together to give that specific Clifford.
"""
function splitUp(n)
    if n == 1
        return [1]
    end
    vals = []
    for subs in gen(noOfOnesIn(n))
        original = string(n-1,base=2)
        rec = zeros(Int64,length(original))
        count = 1
        for x in subs
            while original[count] != '1'
                rec[count]= 0
                count += 1
            end
            rec[count] = x
            count+=1
        end
        push!(vals,1+foldl((x,y)->x*4+y,rec))
    end
    return vals
end

"""
generatePauliFromClifford(n)

## Arguments
    cliffordPs: Array{Float64,1}

Splits the "Clifford average" which we recover from the protocol in Efficient Learning of quantum Noise arXiv:1907.13022. into
a fake Pauli distribution, that would average to the same Clifford average. We are going to use this to test the algorithm.

We do this by realising that each member of our ``2^n`` distribution was in fact an average of certain Paulis. We work out how many
Paulis went into the distribution (```splitUp(n)```) and then project the number onto an appropriate simplex (```simplexUp```).

This allows us to recreate a ``4^n`` distribution that would average down to the supplied distribution. Of course there is a random
element with the projection up stage, but it is *experimentally* inspired if the original averaged protocol came from an experiment.

## Returns
A Pauli probability distribution that result in the supplied probability distribution if they were the errors caused by the
avarged noise in the machine and that machine was subject to single qubit Clifford twirls.
"""
function generatePauliFromClifford(cliffordPs)
    s = Dict{Int64,Float64}()
    nos=0
    totalToSub = 0
    for (ix,x) in enumerate(cliffordPs)
        if x>0
            split  = splitUp(ix)
            # Number of of paulis without a Z
            nsplit = 2^round(Int,(log(3,length(split))))
            increaseRatio = (length(split)+nsplit/2)/length(split)
            totalToSub += x*(increaseRatio-1)
            vals = simplexUp(length(splitUp(ix))).*(x*increaseRatio)
            for (jx,j) in enumerate(splitUp(ix))
                s[j]=vals[jx]
            end
        end
    end
    s[1] -= totalToSub
    lastValue = length(cliffordPs).^2
    s[lastValue] = get(s,lastValue,0)
    sp = sparsevec(s)
    return Vector(sp)
end

end # of module.