text:
type mask_Weight <: masks_Forward_Backward
  Dim1 :: Tuple{Int64,Vararg{Int64}}
  Weight :: Tuple{Array{Float64,1},Vararg{Array{Float64,1}}}
end

function mask_Weight(Dim1 :: Int64, Weight :: Array{Float64,1} )
  return mask_Weight(tuple(Dim1), tuple(Weight))
end

function W_Forward(xarr :: Array{Float64}, M :: mask_Weight)
  ya = copy(xarr)
  for (d, w) in zip(M.Dim1, M.Weight)
    ds = fill(1, ndims(ya))
    ds[d] = length(w)
    we = reshape(w, ds...)
    ya = broadcast(*, we, ya)
  end
  return (ya, M)
end
function W_Forward(xarr :: Array{Float64}, Dim1, Weight)
  M = mask_Weight(Dim1, Weight)
  return W_Forward(xarr,M)
end

export mask_Weight, W_Forward

function BackwardD2(xarr :: Array{Float64}, M :: mask_Weight;
                    Dim :: Tuple{Int64,Vararg{Int64}} = M.Dim1)
  ya = copy(xarr)
  for (d,w) in zip(Dim, M.Weight)
    ds = fill(1, ndims(ya))
    ds[d] = length(w)
    we = reshape(w, ds...)
    we = we
    ya = broadcast(/, ya, we)
  end
  return ya
end
export BackwardD2

function Wevaluate(SS :: _SVD_Functor,
                     xarr :: Array{Float64}, yarr :: Array{Float64},
                     Mx :: mask_Weight ,
                     My :: mask_Weight ,
                     NanDim1x :: Array{Int64,1};
                     NanDim1y :: Array{Int64,1} = NanDim1x,nsv :: Int64 = 3)

  M2x = mask_Weight(Mx.Dim1, map(sqrt, Mx.Weight))
  M2y = mask_Weight(My.Dim1, map(sqrt, My.Weight))
  (xa,) = W_Forward(xarr, M2x)
  (ya,) = W_Forward(yarr, M2y)
  (S1x, S1y, ratio, V, Pctx, Pcty) = DFevaluate(SS, xarr, yarr, NanDim1x, NanDim1y; nsv = nsv)
  Sx = BackwardD2(S1x,M2x)
  Sy = BackwardD2(S1y,M2y)
  return (Sx,Sy,ratio,V,Pctx,Pcty)
end

function Wevaluate(EE :: _EOF_Functor,
                     xarr :: Array{Float64},
                     Mx :: mask_Weight ,
                     NanDim1x :: Array{Int64,1};
                     nsv :: Int64 = 3)
  M2 = mask_Weight(Mx.Dim1, map(sqrt,Mx.Weight))
  (xa,) = W_Forward(xarr, M2)
  (S1x, V,D,ratio) = DFevaluate(EE, xarr, NanDim1x; nsv = nsv)
  Sx = BackwardD2(S1x, M2)
  return (Sx,V,D, ratio)
end

function Wevaluate( MM :: _Mean_Functor, xarr :: Array{Float64},
                     Mx :: mask_Weight, Dim1 :: Array{Int64,1})
  ref = fill(1.0, size(xarr)); ref[isnan(xarr)] = NaN;
  (refw,) = W_Forward(ref , Mx)
  (xarw,) = W_Forward(xarr, Mx)
  return DFevaluate(MM, xarw,  Dim1) ./ DFevaluate(MM, refw,  Dim1)
end

function WMTIndex( pattern :: Array{Float64}, xarr :: Array{Float64},
                  Mx :: mask_Weight, PDim :: Array{Int64,1} = [1:ndims(pattern);] )
  M2 = mask_Weight(Mx.Dim1, map(sqrt,Mx.Weight))
  (x1,) = W_Forward(xarr, M2)
  (p1,) = W_Forward(pattern,M2)
  indc = broadcast(/,MTIndex(p1,x1,PDim), MTIndex(p1,p1, PDim))
end

export Wevaluate, WMTIndex


text:
# by Lazaro Alonso
using CairoMakie, LaTeXStrings, SpecialFunctions
CairoMakie.activate!() # HIDE
let
    fig = Figure(resolution = (600,400), font = "CMU Serif") # probably you need to install this font in your system
    ax = Axis(fig, xlabel = L"x", ylabel = L"J_{\nu}(x)", ylabelsize = 22,
        xlabelsize= 22, xgridstyle=:dash, ygridstyle=:dash, xtickalign = 1,
        xticksize=10, ytickalign=1, yticksize=10,  xlabelpadding = -10)
    x = 0:0.1:15
    for ν in 0:6
        lines!(ax, x, besselj.(ν, x), label = latexstring("J_{$(ν)}(x)"))
    end
    axislegend(; nbanks = 3, framecolor = (:grey,0.5))
    fig[1,1] = ax
    #save("filename.png", fig, px_per_unit = 2.0)
    save(joinpath(@__DIR__, "output", "besselFunctions.png"), fig, px_per_unit = 2.0) # HIDE
    display(fig)
end
using Pkg # HIDE
Pkg.status() # HIDE

ext: jl

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using CheapThreads
using Documenter

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

makedocs(;
    modules=[CheapThreads],
    authors="Chris Elrod <elrodc@gmail.com> and contributors",
    repo="https://github.com/JuliaSIMD/CheapThreads.jl/blob/{commit}{path}#{line}",
    sitename="CheapThreads.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaSIMD.github.io/CheapThreads.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaSIMD/CheapThreads.jl",
)


text:
using JuMP, JuMPChance
using Ipopt
using Gurobi
using Distributions

function Solve_CCLinDist(buses, lines, generators, error_variances)
info("Starting Chance Constraint LinDist Model")

bus_set = collect(keys(buses))
line_set = collect(keys(lines))
gen_set = collect(keys(generators))
gen_bus = []
for g in keys(generators)
    push!(gen_bus, generators[g].bus_idx)
end
root_bus = 0
for b in keys(buses)
    if buses[b].is_root
        root_bus = b
    end
end

lines_to = Dict()
for l in keys(lines)
    lines_to[lines[l].to_node] = lines[l]
end

v_root = 1

error_var_sum = sum(error_variances)

# m = ChanceModel(solver=IpoptSolver())
m = ChanceModel(solver=GurobiSolver())

@variable(m, v[bus_set] >= 0) #variable for voltage square
@variable(m, fp[bus_set]) #variable for active power flow
@variable(m, fq[bus_set]) # variable for reactive power flow
@variable(m, gp[bus_set]) #variable for active power generation
@variable(m, gq[bus_set]) #variable for reactive power generation
@variable(m, α[bus_set] >= 0) # variable to distribute load deviations (affine frequency control)

@indepnormal(m, e[1:length(bus_set)], mean=zeros(length(bus_set)), var=error_variances)
# @indepnormal(m, e2, mean=0, var=1)

@objective(m, Min, sum((gp[g]^2 + α[g]^2*error_var_sum) * buses[g].generator.cost for g in gen_bus))
# Expected Costs

# Control variable
@constraint(m, sum(α) == 1)

#Voltage and flow constraints for root node
@constraint(m, v[root_bus] == v_root)
@constraint(m, fp[root_bus] == 0)
@constraint(m, fq[root_bus] == 0)

for b in setdiff(bus_set, gen_bus)
# All buses without generation
    @constraint(m, α[b] == 0)
    @constraint(m, gp[b] == 0)
    @constraint(m, gq[b] == 0)
end

for b in bus_set
# All buses
    @constraint(m, buses[b].d_P - gp[b] + sum(fp[k] for k in buses[b].children) == fp[b])
    @constraint(m, buses[b].d_Q - gq[b] + sum(fq[k] for k in buses[b].children) == fq[b])
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    b_ancestor = buses[b].ancestor[1]
    @constraint(m, v[b] == v[b_ancestor] - 2*(lines_to[b].r * fp[b] + lines_to[b].x * fq[b]))
    @constraint(m, fp[b]^2 + fq[b]^2 <= lines_to[b].s_max^2)
end

# CHANCE CONSTRAINTS:
for b in gen_bus
# All buses with generation
    @constraint(m, gp[b] + α[b]*sum(e[bb] for bb in bus_set) <= buses[b].generator.g_P_max, with_probability = (1 - η_g))
    @constraint(m, gp[b] + α[b]*sum(e[bb] for bb in bus_set) >= 0, with_probability = (1 - η_g))

    @constraint(m, gq[b] + α[b]*sum(e[bb] * buses[bb].tanphi for bb in bus_set) <=  buses[b].generator.g_Q_max, with_probability = (1 - η_g))
    @constraint(m, gq[b] + α[b]*sum(e[bb] * buses[bb].tanphi for bb in bus_set) >= -buses[b].generator.g_Q_max, with_probability = (1 - η_g))
end

for b in setdiff(bus_set, [root_bus])
# All buses without root
    l = lines_to[b]
    line_arr = [lines[l] for l in line_set]
    @constraint(m, v[b] - 2 * (sum(rPTDF[l.index, b] * l.r * sum(rPTDF[l.index, j] * (e[j] - α[j]*sum(e[z] for z in bus_set)) for j in bus_set) for l in line_arr) +
                    sum(rPTDF[l.index, b] * l.x * sum(rPTDF[l.index, j] * (e[j]*buses[j].tanphi - α[j]*sum(e[z]*buses[j].tanphi for z in bus_set)) for j in bus_set) for l in line_arr))
                    <= buses[b].v_max, with_probability = (1 - η_v))
    @constraint(m, v[b] - 2 * (sum(rPTDF[l.index, b] * l.r * sum(rPTDF[l.index, j] * (e[j] - α[j]*sum(e[z] for z in bus_set)) for j in bus_set) for l in line_arr) +
                    sum(rPTDF[l.index, b] * l.x * sum(rPTDF[l.index, j] * (e[j]*buses[j].tanphi - α[j]*sum(e[z]*buses[j].tanphi for z in bus_set)) for j in bus_set) for l in line_arr))
                    >= buses[b].v_min, with_probability = (1 - η_v))
end

method = :Cuts
# method = :Reformulate
status = solve(m, method=method)
objective = getobjectivevalue(m)
if status == :Optimal
    info("CC solved to optimality with method $(method)")
    info("Objective: $(objective)")
else
    warn("CC not optimal, termination status $(status)")
end

# Prepare Results
bus_results = DataFrame(bus = Any[], d_P = Any[], d_Q = Any[], cosphi=Any[], tanphi = [],  g_P = Any[], g_Q = Any[], g_cost=Any[], v_squared = Any[], α=Any[])
line_results = DataFrame(line = Any[], from = Any[], to = Any[], f_P = Any[], f_Q = Any[], a = Any[])
for b in bus_set
    cost = 0
    if in(b, gen_bus)
        cost = buses[b].generator.cost
    end

    # Numerical smoothing
    α_res = getvalue(α[b])
    α_res < 1e-10 ? α_res = 0 : nothing

    res = [b, buses[b].d_P, buses[b].d_Q, buses[b].cosphi, buses[b].tanphi,  getvalue(gp[b]), getvalue(gq[b]), cost, getvalue(v[b]), α_res]
    push!(bus_results, res)

    # Calc square current on line
    a_res = 0
    fp_res = getvalue(fp[b])
    fq_res = getvalue(fq[b])
    if b != root_bus
        v_res = getvalue(v[buses[b].ancestor[1]])
        a_res = (fp_res^2 + fq_res^2)/v_res^2
        lres = [lines_to[b].index, buses[b].ancestor[1], b, fp_res, fq_res, a_res]
        push!(line_results, lres)
    end
end

sort!(bus_results, cols=:bus)
sort!(line_results, cols=:to)

return objective, bus_results, line_results
end #Solve_CCLinDist


text:
using Luxor

scale = 2
ratio = (√5 -1)/2

w, h = scale * 128 * 2, scale * 128 * 2

r1 = w/2
r2 = r1 * ratio
r3 = r2 * ratio
r4 = r3 * ratio

c1 = Point(0, 0)
c2 = c1 + Point(r1-r2, r1-r2) * √2/2
c3 = c1 + Point(r1-r3, r1-r3) * √2/2
c4 = c1 + Point(r1-r4, r1-r4) * √2/2

Drawing(w, h, "logo.svg")
background(1, 1, 1, 0)
Luxor.origin()

setcolor(1,1,1)
circle(c1, r1, :fill)
setcolor(0.251, 0.388, 0.847)  # dark blue
circle(c1, r1-4*scale, :fill)

setcolor(1,1,1)
circle(c2, r2, :fill)
setcolor(0.796, 0.235, 0.2)  # dark red
circle(c2, r2-4*scale, :fill)

setcolor(1,1,1)
circle(c3, r3, :fill)
setcolor(0.22, 0.596, 0.149) # dark green
circle(c3, r3-4*scale, :fill)

setcolor(1,1,1)
circle(c4, r4, :fill)
setcolor(0.584, 0.345, 0.698) # dark purple
circle(c4, r4-4*scale, :fill)

finish()


text:
using Distributions
using KernelDensityEstimate, KernelDensityEstimatePlotting
using IncrementalInference
using Gadfly, DataFrames

# requried for overloading with new factors
import IncrementalInference: getSample

include(joinpath(@__DIR__, "SquareRootTypes.jl"))

## FIXED POINT ILLUSTRATION

function extractCycleProjections(FGl::Vector{FactorGraph}, FR; N::Int=1000)
  itersl = length(FGl)
  numfr = length(FR)
  ALLXX = zeros(numfr, itersl, 2)
  ALLXY = zeros(numfr, itersl, 2)
  PL = Array{Any,2}(itersl,3)
  p = Dict{Symbol, ManifoldKernelDensity}()
  PP = Dict{Symbol, Vector{ManifoldKernelDensity}}()
  DIV = Dict{Symbol, Vector{Float64}}()
  DIVREF = Dict{Symbol, Vector{Float64}}()
  PP[:x] = ManifoldKernelDensity[]
  PP[:xy] = ManifoldKernelDensity[]
  DIV[:x] = Float64[]
  DIV[:xy] = Float64[]
  DIVREF[:x] = Float64[]
  DIVREF[:xy] = Float64[]

  for i in 1:itersl
    p[:x] = getBelief(FGl[i], :x)
    p[:xy] = getBelief(FGl[i], :xy)

    push!(PP[:x], deepcopy(p[:x]))
    push!(PP[:xy], deepcopy(p[:xy]))

    for j in 1:numfr
      ALLXX[j,i,1] = approxHilbertInnerProd(p[:x], (x) -> phic(x, f=FR[j]), N=N)
      ALLXX[j,i,2] = approxHilbertInnerProd(p[:x], (x) -> phis(x, f=FR[j]), N=N)
      ALLXY[j,i,1] = approxHilbertInnerProd(p[:xy], (x) -> phic(x, f=FR[j]), N=N)
      ALLXY[j,i,2] = approxHilbertInnerProd(p[:xy], (x) -> phis(x, f=FR[j]), N=N)
    end

    PL[i,1] = plotKDE(p[:x],N=N, extend=0.2)
    PL[i,2] = plotKDE(p[:xy],N=N, extend=0.2)
    PL[i,3] = plotKDE([p[:x],p[:xy]],c=["red","green"],N=N, extend=0.2)

    if i > 1
      # KL-divergence
      mkl = minimum([kld(p[:x],getBelief(FGl[i-1],:x))[1]; kld(getBelief(FGl[i-1],:x),p[:x])[1]])
      push!(DIV[:x], abs(mkl))
      mklxy = minimum([kld(p[:xy],getBelief(FGl[i-1],:xy))[1]; kld(getBelief(FGl[i-1],:xy),p[:xy])[1]])
      push!(DIV[:xy], abs(mklxy))
    end

    # reference KL-divergence
    mkl = minimum([kld(p[:x],getBelief(FGl[end],:x))[1]; kld(getBelief(FGl[end],:x),p[:x])[1]])
    push!(DIVREF[:x], abs(mkl))
    mklxy = minimum([kld(p[:xy],getBelief(FGl[end],:xy))[1]; kld(getBelief(FGl[end],:xy),p[:xy])[1]])
    push!(DIVREF[:xy], abs(mklxy))
  end

  cycle = Dict{Symbol, Any}()
  # save everything
  cycle[:N] = N
  cycle[:iters] = itersl-1
  cycle[:FR] = FR
  cycle[:numfr] = numfr
  cycle[:PP] = PP
  cycle[:ALLXX] = ALLXX
  cycle[:ALLXY] = ALLXY
  cycle[:DIV] = DIV
  cycle[:DIVREF] = DIVREF
  return cycle
end

function runFullBatchIterations(;N=100, iters=50)
  fg = initfg()

  x0 = 0.5 .- rand(1,N)  #[1.0+0.1*randn(N);10+0.1*randn(N)]'
  addVariable!(fg, :x, ContinuousScalar)
  # addVariable!(fg, :y, x0,N=N)

  # TODO make Mixture instead
  pts = rand(Distributions.Normal(4.0,0.05),N)   #;rand(Distributions.Normal(144.0,0.05),N)]
  md = kde!(pts)

  npx = NumbersPrior(md)
  addVariable!(fg, :xy, ContinuousScalar)

  addFactor!(fg, [:xy;], npx)
  #
  # xey = AreEqual(Distributions.Normal(0.0,0.01))
  # addFactor!(fg, [getVariable(fg, :x);getVariable(fg, :y)], xey)

  # xty = ProductNumbers(Distributions.Normal(0.0,0.01))
  # addFactor!(fg, [getVariable(fg, :x);getVariable(fg, :y);getVariable(fg, :xy)], xty)

  xty = Square(Distributions.Normal(0.0,0.01))
  addFactor!(fg, [:x,:xy], xty)

  FG = deepcopy(fg)
  for i in 1:iters
    tree = solveTree!(fg)
    push!(FG, deepcopy(fg))
  end
  return FG
end

runFullBatchIterations(iters=1)

# Do all the runs
# frequencies of interest

FR = range(0.5/(2pi),stop=3.0/(2pi), length=8)
# FR = range(0.5/(2pi),stop=3.0/(2pi), length=8)

mc = 3

# data containers
FG = []
CYCLE = Vector{Dict}(undef, mc)

for i in 1:mc
  push!(FG, runFullBatchIterations(;N=100, iters=50))
  CYCLE[i] = extractCycleProjections(FG[i], FR, N=2000)
end

# Analyse the data

function plotXandXYFreq(CYCLE, FR, fridx)
  XXfrs = DataFrame[]
  for i in 1:mc
    push!(XXfrs, DataFrame(
      x=CYCLE[i][:ALLXX][fridx,:,1],
      y=CYCLE[i][:ALLXX][fridx,:,2],
      MC="$(i)"
    ))
  end

  XXrepeat = vcat(XXfrs...)

  plxrep = Gadfly.plot(XXrepeat,
    Geom.path(),
    Guide.title("X"),
    Guide.title("μ=0, frequency=$(round(FR[fridx],3))"),
    Guide.ylabel("p ⋅ ϕs"),
    Guide.xlabel("p ⋅ ϕc"),
    x=:x,
    y=:y,
    color=:MC
  )

  XYfrs = DataFrame[]
  for i in 1:mc
    push!(XYfrs, DataFrame(
      x=CYCLE[i][:ALLXY][fridx,:,1],
      y=CYCLE[i][:ALLXY][fridx,:,2],
      MC="$(i)"
    ))
  end

  XYrepeat = vcat(XYfrs...)

  plxyrep = Gadfly.plot(XYrepeat,
    Geom.path(),
    Guide.title("X^2"),
    Guide.title("μ=0, frequency=$(round(FR[fridx],3))"),
    Guide.ylabel("p ⋅ ϕs"),
    Guide.xlabel("p ⋅ ϕc"),
    x=:x,
    y=:y,
    color=:MC
  )
  return plxyrep,plxrep
end

for i in 1:length(FR)
  plxyrep,plxrep = plotXandXYFreq(CYCLE, FR, i)
  Gadfly.draw(PDF("sqrtexamplexxtrajFR$(i).pdf",15cm,7cm),hstack(plxyrep,plxrep))
end
# @async run(`evince sqrtexamplexxtrajFR2.pdf`)

0

#
#
# plotKDE(
# [CYCLE[1][:PP][:x][1];
# CYCLE[1][:PP][:x][2];
# CYCLE[1][:PP][:x][3]],
# c=["black","red","green"]
# )

DFs = DataFrame[]

for i in [1,2,3,5,13]
  p = CYCLE[1][:PP][:x][i]
  mxmx = getKDERange(p)
  x = [range(mxmx[1], stop=mxmx[2], length=2000);]
  push!(DFs, DataFrame(
    x = x,
    y = clamp(evaluateDualTree(p,x),0,4),
    Iteration="$(i-1)"
  ))
end

plx = Gadfly.plot(vcat(DFs...) , x=:x, y=:y, color=:Iteration,
  Geom.line,
  Guide.xlabel("X"),
  Guide.ylabel("pdf")
)

DFs = DataFrame[]

for i in [1,2,3,5,13]
  p = CYCLE[1][:PP][:xy][i]
  mxmx = getKDERange(p, extend=0.4)
  x = [range(mxmx[1], stop=mxmx[2], length=2000);]
  push!(DFs, DataFrame(
    x = x,
    y = clamp(evaluateDualTree(p,x),0,6),
    Iteration="$(i-1)"
  ))
end

plxy = Gadfly.plot(vcat(DFs...) , x=:x, y=:y, color=:Iteration,
  Geom.line,
  Guide.xlabel("X^2"),
  Guide.ylabel("pdf")
)

plh = vstack(plxy, plx)

Gadfly.draw(PDF("sqrtexamplebeliefs.pdf",12cm,9cm),plh)

@async run(`evince sqrtexamplebeliefs.pdf`)

0

# Also plot the KL divergences

DFs = DataFrame[]

for i in 1:mc
  push!(DFs, DataFrame(
    x = 1:length(CYCLE[i][:DIV][:xy]),
    y = CYCLE[i][:DIV][:xy],
    MC="$(i), X^2"
  ))
end

for i in 1:mc
  push!(DFs, DataFrame(
    x = 1:length(CYCLE[i][:DIV][:x]),
    y = CYCLE[i][:DIV][:x],
    MC="$(i), X"
  ))
end

pldiv = Gadfly.plot(vcat(DFs...) , x=:x, y=:y, color=:MC,
  Geom.line,
  Guide.title("Relative between iterations"),
  Guide.xlabel("Iterations"),
  Guide.ylabel("KL-Divergence")
)

Gadfly.draw(PDF("sqrtexamplekldrelative.pdf",10cm,6cm),pldiv)
# @async run(`evince sqrtexamplekldrelative.pdf`)

DFs = DataFrame[]
# x0 = [1.0+0.1*randn(100);10+0.1*randn(100)]'
# addVariable!(fg, :x, ContinuousScalar, N=N)
#
# pts = [rand(Distributions.Normal(4.0,0.05),100);rand(Distributions.Normal(144.0,0.05),100)]
# md = kde!(pts)
# npx = NumbersPrior(md)
# pts0 = getSample(npx,N)[1]
# addVariable!(fg, :xy, ContinuousScalar, N=N)

for i in 1:mc
  push!(DFs, DataFrame(
    x = 1:length(CYCLE[i][:DIVREF][:xy]),
    y = CYCLE[i][:DIVREF][:xy],
    MC="$(i), X^2"
  ))
end

for i in 1:mc
  push!(DFs, DataFrame(
    x = 1:length(CYCLE[i][:DIVREF][:x]),
    y = CYCLE[i][:DIVREF][:x],
    MC="$(i), X"
  ))
end

pldivref = Gadfly.plot(vcat(DFs...) , x=:x, y=:y, color=:MC,
  Geom.line,
  Guide.title("Referenced to final belief"),
  Guide.xlabel("Iterations"),
  Guide.ylabel("KL-Divergence")
)

Gadfly.draw(PDF("sqrtexamplekldreferenced.pdf",10cm,6cm),pldivref)
@async run(`evince sqrtexamplekldreferenced.pdf`)

0

#
# Gadfly.plot(x=cycle1[:ALLXX][1,:,1],y=cycle1[:ALLXX][1,:,2],
# Geom.path(), Guide.title("μ=0, frequency=$(round(FR[1],3))"))
# Gadfly.plot(x=ALLXY[1,:,1],y=ALLXY[1,:,2],
# Geom.path(), Guide.title("μ=0, frequency=$(round(FR[1],3))"))
#
# Gadfly.plot(x=ALLXX[2,:,1],y=ALLXX[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[2],3))"))
# Gadfly.plot(x=ALLXY[2,:,1],y=ALLXY[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[2],3))"))
#
# Gadfly.plot(x=ALLXX[3,:,1],y=ALLXX[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[3],3))"))
# Gadfly.plot(x=ALLXY[3,:,1],y=ALLXY[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[3],3))"))
#
# Gadfly.plot(x=ALLXX[4,:,1],y=ALLXX[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[4],3))"))
# Gadfly.plot(x=ALLXY[4,:,1],y=ALLXY[2,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[4],3))"))
#
# Gadfly.plot(x=ALLXX[end,:,1],y=ALLXX[3,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[5],3))"))
# Gadfly.plot(x=ALLXY[end,:,1],y=ALLXY[3,:,2], Geom.path(), Guide.title("μ=0, frequency=$(round(FR[5],3))"))
#
#
#
# # plotKDE(getBelief(fg,:y))
#
# PL[1,3]
# PL[2,3]
# PL[3,3]
# PL[4,3]
# PL[6,3]
#
# PL[8,3]
#
# PL[10,3]
#
# PL[12,3]
#
# PL[20,3]
#

# Now evaluate the value function for studing the Bellman equation

using KernelDensityEstimatePlotting

plot(getBelief(fg,:xy),N=2000)

plot(getBelief(fg,:x),N=2000)

# plot(getBelief(fg,:y))

plot(md)

#


text:
export iauAtic13
"""
Transform star RA,Dec from geocentric CIRS to ICRS astrometric.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    ri,di  double  CIRS geocentric RA,Dec (radians)
    date1  double  TDB as a 2-part...
    date2  double  ...Julian Date (Note 1)

Returned:
    rc,dc  double  ICRS astrometric RA,Dec (radians)
    eo     double  equation of the origins (ERA-GST, Note 4)

Notes:

    1. The TDB date date1+date2 is a Julian Date, apportioned in any
    convenient way between the two arguments.  For example,
    JD(TDB)=2450123.7 could be expressed in any of these ways, among
    others:

        date1          date2

        2450123.7           0.0       (JD method)
        2451545.0       -1421.3       (J2000 method)
        2400000.5       50123.2       (MJD method)
        2450123.5           0.2       (date & time method)

    The JD method is the most natural and convenient to use in cases
    where the loss of several decimal digits of resolution is
    acceptable.  The J2000 method is best matched to the way the
    argument is handled internally and will deliver the optimum
    resolution.  The MJD method and the date & time methods are both
    good compromises between resolution and convenience.  For most
    applications of this function the choice will not be at all
    critical.

    TT can be used instead of TDB without any significant impact on
    accuracy.

    2. Iterative techniques are used for the aberration and light
    deflection corrections so that the functions iauAtic13 (or
    iauAticq) and iauAtci13 (or iauAtciq) are accurate inverses;
    even at the edge of the Sun's disk the discrepancy is only about
    1 nanoarcsecond.

    3. The available accuracy is better than 1 milliarcsecond, limited
    mainly by the precession-nutation model that is used, namely
    IAU 2000A/2006.  Very close to solar system bodies, additional
    errors of up to several milliarcseconds can occur because of
    unmodeled light deflection;  however, the Sun's contribution is
    taken into account, to first order.  The accuracy limitations of
    the SOFA function iauEpv00 (used to compute Earth position and
    velocity) can contribute aberration errors of up to
    5 microarcseconds.  Light deflection at the Sun's limb is
    uncertain at the 0.4 mas level.

    4. Should the transformation to (equinox based) J2000.0 mean place
    be required rather than (CIO based) ICRS coordinates, subtract the
    equation of the origins from the returned right ascension:
    RA = RI - EO.  (The iauAnp function can then be applied, as
    required, to keep the result in the conventional 0-2pi range.)

Called:
    iauApci13    astrometry parameters, ICRS-CIRS, 2013
    iauAticq     quick CIRS to ICRS astrometric

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauAtic13(double ri, double di, double date1, double date2,
#                double *rc, double *dc, double *eo)

function iauAtic13(ri::Real, di::Real, date1::Real, date2::Real)
   # Allocate return values
   ref_rc = Ref{Float64}(0.0)
   ref_dc = Ref{Float64}(0.0)
   ref_eo  = Ref{Float64}(0.0)

   status = ccall((:iauAtic13, libsofa_c), Cint,
            (Cdouble, Cdouble, Cdouble, Cdouble,
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, ri), convert(Float64, di),
            convert(Float64, date1), convert(Float64, date2),
            ref_rc, ref_dc, ref_eo)

    return ref_rc[], ref_dc[], ref_eo[]
end


text:
abstract type AbstractMovingWindowSolver <: AbstractInitialValueSolver end
abstract type AbstractMovingWindowSolution <: AbstractInitialValueSolution end
abstract type AbstractMovingWindowCache <: AbstractInitialValueCache end
abstract type AbstractMovingWindowParameters <: AbstractInitialValueParameters end
abstract type AbstractBudget <: AbstractMovingWindowParameters end
abstract type AbstractAdaptiveParameters <: AbstractMovingWindowParameters end


text:
# This file defines the CompoundNode type. This type defines a node which is a
#   collection of several base nodes.

export CompoundNode
"""
The `CompoundNode` type defines a compound node which groups a number of base nodes in a manpower simulation. This is used to generate a hierarchy of levels that personnel members can occupy. For example: assume that there are compound nodes for each combination of service branch (Land, Air, Marine, Medical) and grade category (Volunteer, Non-com, Officer). The compound nodes can then group all personnel members by service brance, so each compound node corresponds to a service branch.

The type contains the following fields:
* `name::String`: the name of the compound node.
* `baseNodeList::Vector{String}`: the names of the (base) nodes the compound node consists of. Remark that the simulation must ensure that all the component base nodes actually exist.
* `nodeTarget::Int`: the target number of personnel members in the compound node, where no target is represented by -1. Default = -1

Constructor:
```
CompoundNode( name::String )
```
This constructor creates a `CompoundNode` object with name `name`.
"""
mutable struct CompoundNode

    name::String
    baseNodeList::Vector{String}
    nodeTarget::Int  # Necessary? It's not used in the simulation (yet).

    function CompoundNode( name::String )

        newCompNode = new()
        newCompNode.name = name
        newCompNode.baseNodeList = Vector{String}()
        newCompNode.nodeTarget = -1
        return newCompNode

    end  # CompoundNode( name )

end  # mutable struct CompoundNode


text:
## datastore
using Distributed

struct DRef
    owner::Int
    id::Int
    size::UInt
end

mutable struct RefState
    size::UInt64
    data::Union{Some{Any}, Nothing}
    file::Union{String, Nothing}
    destroyonevict::Bool # if true, will be removed from memory
end

isinmemory(x::RefState) = x.data !== nothing
isondisk(x::RefState) = x.file !== nothing

const datastore = Dict{Int,RefState}()
const id_counter = Ref(0)

function poolset(@nospecialize(x), pid=myid(); size=approx_size(x), destroyonevict=false, file=nothing)
    if pid == myid()
        id = id_counter[] += 1
        #lru_free(size)
        datastore[id] = RefState(size,
                                 Some{Any}(x),
                                 file,
                                 destroyonevict)
        #lru_touch(id, size)
        DRef(myid(), id, size)
    else
        # use our serialization
        remotecall_fetch(pid, MMWrap(x)) do wx
            # todo: try to evict here before deserializing...
            poolset(unwrap_payload(wx), pid)
        end
    end
end

function forwardkeyerror(f)
    try
        f()
    catch err
        if isa(err, RemoteException) && isa(err.captured.ex, KeyError)
            rethrow(err.captured.ex)
        end
        rethrow(err)
    end
end

const file_to_dref = Dict{String, DRef}()
const who_has_read = Dict{String, Vector{DRef}}() # updated only on master process
const enable_who_has_read = Ref(true)
const enable_random_fref_serve = Ref(true)
const wrkrips = Dict{IPv4,Vector{Int}}() # cached IP-pid map to lookup FileRef locations

is_my_ip(ip::String) = getipaddr() == IPv4(ip)
is_my_ip(ip::IPv4) = getipaddr() == ip

function get_wrkrips()
    d = Dict{IPv4,Vector{Int}}()
    for w in Distributed.PGRP.workers
        wip = IPv4(isa(w, Distributed.Worker) ? w.config.bind_addr : w.bind_addr)
        if wip in keys(d)
            if enable_random_fref_serve[]
                push!(d[wip], w.id)
            else
                d[wip] = [min(d[wip][1], w.id)]
            end
        else
            d[wip] = [w.id]
        end
    end

    loopback = ip"127.0.0.1"
    if (loopback in keys(d)) && (length(d) > 1)
        # there is a chance that workers on loopback are actually on same ip as master
        realip = remotecall_fetch(getipaddr, first(d[loopback]))
        if realip in keys(d)
            append!(d[realip], d[loopback])
        else
            d[realip] = d[loopback]
        end
        delete!(d, loopback)
    end
    d
end

get_worker_at(ip::String) = get_worker_at(IPv4(ip))
get_worker_at(ip::IPv4) = rand(get_workers_at(ip))

function get_workers_at(ip::IPv4)
    isempty(wrkrips) && merge!(wrkrips, remotecall_fetch(get_wrkrips, 1))
    wrkrips[ip]
end

function poolget(r::FileRef)
    # since loading a file is expensive, and often
    # requested in quick succession
    # we add the data to the pool
    if haskey(file_to_dref, r.file)
        ref = file_to_dref[r.file]
        if ref.owner == myid()
            return poolget(ref)
        end
    end

    if is_my_ip(r.host)
        x = unwrap_payload(open(deserialize, r.file, "r+"))
    else
        x = remotecall_fetch(get_worker_at(r.host), r.file) do file
            unwrap_payload(open(deserialize, file, "r+"))
        end
    end
    dref = poolset(x, file=r.file, size=r.size)
    file_to_dref[r.file] = dref
    if enable_who_has_read[]
        remotecall_fetch(1, dref) do dr
            who_has = MemPool.who_has_read
            if !haskey(who_has, r.file)
                who_has[r.file] = DRef[]
            end
            push!(who_has[r.file], dr)
            nothing
        end
    end
    return x
end

function poolget(r::DRef)
    if r.owner == myid()
        _getlocal(r.id, false)
    else
        forwardkeyerror() do
            remotecall_fetch(r.owner, r) do r
                MMWrap(_getlocal(r.id, true))
            end
        end
    end |> unwrap_payload
end

function _getlocal(id, remote)
    state = datastore[id]
    if remote
        if isondisk(state)
            return FileRef(state.file, state.size)
        elseif isinmemory(state)
            #lru_touch(id, state.size)
            return something(state.data)
        end
    else
        # local
        if isinmemory(state)
            #lru_touch(id, state.size)
            return something(state.data)
        elseif isondisk(state)
            # TODO: get memory mapped size and not count it as
            # much as local process size
            #lru_free(state.size)
            data = unwrap_payload(open(deserialize, state.file, "r+"))
            state.data = Some{Any}(data)
            #lru_touch(id, state.size) # load back to memory
            return data
        end
    end
    error("ref id $id not found in memory or on disk!")
end

function datastore_delete(id)
    (id in keys(datastore)) || (return nothing)
    state = datastore[id]
    if isondisk(state)
        f = state.file
        isfile(f) && rm(f; force=true)
    end
    # release
    state.data = nothing
    #delete!(lru_order, id)
    delete!(datastore, id)
    # TODO: maybe cleanup from the who_has_read list?
    nothing
end

function cleanup()
    empty!(file_to_dref)
    empty!(who_has_read)
    map(datastore_delete, collect(keys(datastore)))
    d = default_dir(myid())
    isdir(d) && rm(d; recursive=true)
    nothing
end

function pooldelete(r::DRef)
    if r.owner != myid()
        return remotecall_fetch(pooldelete, r.owner, r)
    end
    datastore_delete(r.id)
end

function pooldelete(r::FileRef)
    isfile(r.file) && rm(r.file)
    (r.file in keys(file_to_dref)) && delete!(file_to_dref, r.file)
    # TODO: maybe cleanup from the who_has_read list?
    nothing
end

function destroyonevict(r::DRef, flag::Bool=true)
    if r.owner == myid()
        datastore[r.id].destroyonevict = flag
    else
        forwardkeyerror() do
            remotecall_fetch(destroyonevict, r.owner, r, flag)
        end
    end
end

default_dir(p) = joinpath(".mempool", "$session-$p")
default_path(r::DRef) = joinpath(default_dir(r.owner), string(r.id))

function movetodisk(r::DRef, path=default_path(r), keepinmemory=false)
    if r.owner != myid()
        return remotecall_fetch(movetodisk, r.owner, r, path)
    end

    mkpath(dirname(path))
    if isfile(path)
        return FileRef(path, r.size)
    end
    open(path, "w") do io
        serialize(io, MMWrap(poolget(r)))
    end
    s = datastore[r.id]
    s.file = path
    fref = FileRef(path, r.size)
    if !keepinmemory
        s.data = nothing
        #delete!(lru_order, r.id)
    end

    return fref
end

copytodisk(r::DRef, path=default_path(r)) = movetodisk(r, path, true)

"""
Allow users to specifically save something to disk.
This does not free the data from memory, nor does it affect
size accounting.
"""
function savetodisk(r::DRef, path)
    if r.owner == myid()
        open(path, "w") do io
            serialize(io, MMWrap(poolget(r)))
        end
        return FileRef(path, r.size)
    else
        return remotecall_fetch(savetodisk, r.owner, r, path)
    end
end

function deletefromdisk(r::DRef, path)
    if r.owner == myid()
        return rm(path)
    else
        return remotecall_fetch(deletefromdisk, r.owner, r, path)
    end
end

#### LRU
#
#using DataStructures
#
#const max_memsize = Ref{UInt}(2^30)
#const lru_order = OrderedDict{Int, UInt}()
#
## marks the ref as most recently used
#function lru_touch(id::Int, sz = datastore[id].size)
#    if haskey(lru_order, id)
#	delete!(lru_order, id)
#	lru_order[id] = sz
#    else
#	lru_order[id] = sz
#    end
#end
#
## returns ref ids that can be evicted to keep
## space below max_memsize
#function lru_evictable(required=0)
#    total_size = UInt(sum(values(lru_order)) + required)
#    if total_size <= max_memsize[]
#	return Int[] # nothing to evict
#    end
#
#    size_to_evict = total_size - max_memsize[]
#    memsum = UInt(0)
#    evictable = Int[]
#
#    for (k, v) in lru_order
#	if total_size - memsum <= max_memsize[]
#	    break
#	end
#	if datastore[k].destroyonevict # remove only recomputable chunks
#	    memsum += v
#	    push!(evictable, k)
#	end
#    end
#    return evictable
#end
#
#const spilltodisk = Ref(false)
#
#function lru_free(sz)
#    list = lru_evictable(sz)
#    for id in list
#	state = datastore[id]
#	ref = DRef(myid(), id, state.size)
#	if state.destroyonevict
#	    pooldelete(ref)
#	else
#	    !spilltodisk[] && continue
#	    #movetodisk(ref)
#	end
#    end
#end


text:
using RayCast.Bresenham: cast_heading
using RoboLib.Util: euclidean, img2grid
using RoboLib.Geom: Pose2D
using LinearAlgebra

export LaserScanModel

@inline amcl_init() = 1
@inline amcl_combine(prob, p) = prob + p^3
@inline amcl_update(weight, prob, inv_squash) = weight *= prob

@noinline std_init() = 1
@noinline std_combine(prob, p) = prob * p
@noinline std_update(weight, prob, inv_squash) = prob ^ inv_squash

@inline log_init() = 0
@inline log_combine(prob, p) = prob + log(p)
@inline log_update(weight, prob, inv_squash) = prob * inv_squash

struct LaserScanModel{D, M, F, F1, F2, F3, S} <: AbstractSensorModel
    beammodel::M
    "rangemethodW(T_WP::Pose2D, theta_W)->(xhit_W, yhit_W)"
    rangemethod::F
    init::F1
    combine::F2
    update::F3
    inv_squash_factor::S
    nposes::Int
    function LaserScanModel{D}(beammodel::M, rangemethod::F, inv_squash_factor::S, nposes::Int, method=:std) where {D,M,F,S}
        if method === :std
            init = std_init
            combine = std_combine
            update = std_update
        elseif method === :amcl
            init = amcl_init
            combine = amcl_combine
            update = amcl_update
        elseif method === :log
            init = log_init
            combine = log_combine
            update = log_update
        else
            error("Type \"$method\" not know")
        end
        new{D,M,F,typeof(init),typeof(combine),typeof(update), S}(beammodel, rangemethod, init, combine, update, inv_squash_factor, nposes)
    end
end

function reweight!(weights_k::Vector{W}, particles_WP_k::Vector{<:Pose2D{T}}, obs_W_k, rng, dt, m::LaserScanModel) where {W,T}
    obsangles_k, obsranges_W_k = obs_W_k

    #totalw = W(0)
    @inbounds Threads.@threads for i in eachindex(particles_WP_k)
        prob = W(m.init())
        T_WP = particles_WP_k[i]
        @inbounds for (obsangle, obsrange_W) in zip(obsangles_k, obsranges_W_k)
            range_exp_W = m.rangemethod(T_WP, -obsangle) #TODO: negative sign?
            p = m.beammodel(obsrange_W, range_exp_W)
            prob = m.combine(prob, p)
        end
        @assert !iszero(prob) && !isnan(prob)
        w = m.update(weights_k[i], prob, m.inv_squash_factor)
        weights_k[i] = w
        #totalw += w
    end
    # Don't actually need to normalize, but keeping it here for now
    # weights_k ./= totalw
    weights_k
end


text:
function mw_rw_cov(t)
  sum_mat = zeros(t,t)
  for i ∈ 1:t
    sum_mat[i, 1:i] = ones(1, i)
  end
  cov = sum_mat*sum_mat'
  return cov
end


text:
"""
Observations from `Normal(μ, 1)` with an unknown `μ`, estimated using a
multivatiate normal auxiliary model (so in this case, the likelihood is exact
and the summary statistics lose no information).
"""
struct NormalMeanModel
    "number of observations"
    N::Int
end

common_random(rng::AbstractRNG, model::NormalMeanModel) = randn(rng, model.N)

common_random!(rng::AbstractRNG, model::NormalMeanModel, ϵ) = randn!(rng, ϵ)

function simulate_data(rng::AbstractRNG, model::NormalMeanModel, θ, ϵ)
    μ = θ[1]
    MvNormalData(reshape(μ .+ ϵ, :, 1))
end

simulate_data(rng::AbstractRNG, ::NormalMeanModel, ::String, ::Any) = nothing

MLE(::NormalMeanModel, data) = MLE(MvNormalModel(), data)

loglikelihood(::NormalMeanModel, data, ϕ) = loglikelihood(MvNormalModel(), data, ϕ)

@testset "indirect likelihood toy problem" begin
    μ₀ = 2.0
    p = simulate_problem(NormalMeanModel(100), zero ∘ first, [μ₀])
    # find the optimum using the same common random numbers
    f(x) = (-indirect_logposterior(p, [x]))[1]
    o = optimize(f, 0.0, 5.0)
    μ₁ = Optim.minimizer(o)
    @test μ₁ ≈ μ₀ atol = 1e-4
    # test common random numbers updating
    @test mean([(p = common_random!(p); mean(p.ϵ))
                for _ in 1:1000]) ≈ 0 atol = 0.01
    # invalid parameters: when data is nothing, indirect log likelihood is -Inf
    @test indirect_logposterior(p, "a fish") == -Inf
    # local jacobian
    parameter_transformation = TransformationTuple((IDENTITY, ))
    @test local_jacobian(p, (2.0, ), parameter_transformation) == [1.0 0.0]'
    # buildup of likelihood
    p2 = IndirectLikelihoodProblem(NormalMeanModel(100), zero ∘ first, p.data)
    θ = [μ₀]
    y = simulate_data(p2.rng, p2.model, θ, p2.ϵ)
    ϕ = MLE(p2.model, y)
    ℓ = loglikelihood(p2.model, p2.data, ϕ)
    @test indirect_logposterior(p2, θ) == ℓ
    @test indirect_logposterior(p2, θ) == indirect_logposterior(p2)(θ)
end

struct ToyModel end

common_random(rng::AbstractRNG, ::ToyModel) = 9.0

simulate_data(rng::AbstractRNG, ::ToyModel, ::Any, ::Any) = 42

@testset "simulation & CRN barebone" begin
    @test common_random(RNG, ToyModel()) == 9.0
    @test common_random!(RNG, ToyModel(), nothing) == 9.0 # fallback
    @test simulate_data(RNG, ToyModel(), :parameter) == 42
    @test simulate_data(ToyModel(), :parameter) == 42
end


text:
module FP

using Distributions
using Printf

# Show floats to 2 digits after the decimal point
Base.show(io::IO, f::Float64) = @printf(io, "%2.2f", f)

Δ = 0.0     # propagation delay upper bound
h = 0.12    # honest mining rate
a = 0.08    # adversarial mining rate
λ = h + a   # total mining rate

K = 30 # confirmation depth

println("Simulation begins.\n")
println(
    "confirmation depth=",
    K,
    "  max delay=",
    Δ,
    "  honest mining rate=",
    h,
    "  adversarial mining rate=",
    a,
)

nmax = round(Int, K * 3 * λ / a)   # In principle the first violation can occur in the arbitrary future.
# However, the probability that it first occurs at a height much larger than K is extremely small.  Here we only simulate up to about 3K adversarial blocks.

ntrials = 1000

violation = 0.0

for nt = 1:ntrials
    xx = rand(Exponential(1.0 / λ), nmax)   # exponential inter-mining times for all blocks
    tt = cumsum(xx)   # the (cumulative) mining times
    bb = rand(Bernoulli(h / λ), nmax)   # each block is honest with probability h/λ
    bb[1] = true    # We assume block 1 is always honest - which is the target of the attack.

    tth = tt[bb]        # honest block mining times
    tta = tt[.!bb]      # adversarial block mining times

    len = min(length(tth), length(tta))

    if any(tth[K:len] .> tta[K:len])
        global violation += 1
    end

    if nt == 1  # print an example of
        println("\nThe following are the mining times in trial #1")
        for n = 1:nmax
            print(bb[n] ? "honest" : "advers")
            println(" block mined at ", tt[n], "    after interval ", xx[n])
            # println( bb[n], " ", xx[n], " ", tt[n] )
        end
        println(
            "\nThe mining times of the first ",
            len,
            " honest blocks:\n",
            tth[1:len],
        )
        println(
            "The mining times of the first ",
            len,
            " adversarial blocks:\n",
            tta[1:len],
        )
        println(
            "When the adversarial is ahead in this trial (1)?  ",
            tth[1:len] .> tta[1:len],
        )
    end
end

println(
    "\nThe probability of block 1's safety violation is about ",
    100.0 * violation / ntrials,
    "% based on ",
    ntrials,
    " trials.",
)
println("Simulation ends.")

end


text:
using ClutteredEnvPathOpt
using Plots

function gen_separator_size(func, generator, range)
    return map(
        i -> begin
            lg = LabeledGraph(generator(i))
            (separator, a, b) = func(lg, 1)
            return length(separator)
        end,
        range
    )
end

function gen_balance(func, generator, range)
    return map(
        i -> begin
            lg = LabeledGraph(generator(i))
            (separator, a, b) = func(lg, 1)
            return length(a) / length(b)
        end,
        range
    )
end

function gen_elapsed(func, generator, range)
    return map(
        i -> begin
            lg = LabeledGraph(generator(i))
            return @elapsed func(lg, 1)
        end,
        range
    )
end

function compare(metric, generator, range, filename)
    data = hcat(
        metric(find_separator_lt, generator, range),
        metric(find_separator_fcs, generator, range),
        metric(find_separator_fcs_best, generator, range),
    )

    plot(range, data, title = filename, label = ["Lipton-Tarjan" "FCS" "FCS Best"])
    savefig(string("stats/", filename, ".png"))
end

# compare(gen_separator_size, ClutteredEnvPathOpt.LightGraphs.ladder_graph, 4:50, "size")
# compare(gen_balance, ClutteredEnvPathOpt.LightGraphs.ladder_graph, 4:50, "balance")
# compare(gen_elapsed, ClutteredEnvPathOpt.LightGraphs.ladder_graph, 4:50, "time")

compare(gen_separator_size, i -> ClutteredEnvPathOpt.LightGraphs.grid([i, i]), 4:50, "size")
compare(gen_balance, i -> ClutteredEnvPathOpt.LightGraphs.grid([i, i]), 4:50, "balance")
compare(gen_elapsed, i -> ClutteredEnvPathOpt.LightGraphs.grid([i, i]), 4:50, "time")


text:
using PixelArt
using Documenter

makedocs(;
    modules=[PixelArt],
    authors="Yuchi Yamaguchi",
    repo="https://github.com/abap34/PixelArt.jl/blob/{commit}{path}#L{line}",
    sitename="PixelArt.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://abap34.github.io/PixelArt.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/abap34/PixelArt.jl",
)


text:
"""
    Reynolds_Number(Tair,pressure,ustar,z0m; constants)

calculates the Roughness Reynolds Number.

# Arguments
- `Tair`      : Air temperature (deg C)
- `pressure`  : Atmospheric pressure (kPa)
- `ustar`     : Friction velocity (m s-1)
- `z0m`       : Roughness length (m)
optional
- `constants=`[`BigleafConstants`](@ref)`()`

# Details
The Roughness Reynolds Number is calculated as in Massman 1999a:
``Re = z0m * ustar / v``, where `v` is the kinematic viscosity (m2 s-1).

# Value
Roughness Reynolds Number (-)

# References
Massman, W_J., 1999a: A model study of kB H- 1 for vegetated surfaces using
'localized near-field' Lagrangian theory. Journal of Hydrology 223, 27-43.

```jldoctest; output = false
Tair,pressure,ustar,z0m = 25.0,100.0,0.5,0.5
R = Reynolds_Number(Tair,pressure,ustar,z0m)
≈(R, 15870, rtol=1e-3)
# output
true
```
"""
function Reynolds_Number(Tair,pressure,ustar,z0m;constants=BigleafConstants())
  v  = kinematic_viscosity(Tair,pressure;constants)
  Re = z0m*ustar/v
end

abstract type RoughnessMethod end
struct Roughness_wind_profile <: RoughnessMethod end
struct RoughnessCanopyHeight <: RoughnessMethod end
struct RoughnessCanopyHeightLAI <: RoughnessMethod end

"""
    roughness_parameters(::RoughnessCanopyHeight()    , zh; frac_d=0.7, frac_z0m=0.1)
    roughness_parameters(::RoughnessCanopyHeightLAI(), zh, LAI; cd=0.2, hs=0.01)
    roughness_parameters(::Roughness_wind_profile()     , ustar, wind, psi_m;
      zh, zr, d = 0.7*zh, constants)

    roughness_parameters(::Roughness_wind_profile(), ustar, wind, Tair, pressure, H;
      zh, zr, d = 0.7*zh, stab_formulation=Dyer1970(), constants)
    roughness_parameters(::Roughness_wind_profile(), df::DFTable; ...)

A approximations of the two roughness parameters displacement height (d)
and roughness length for momentum (z0m).

# Arguments
- `zh`        : Vegetation height (m)
- `constants=`[`BigleafConstants`](@ref)`()`:

By canopy height:
- `frac_d`    : Fraction of displacement height on canopy height (-)
- `frac_z0m`  : Fraction of roughness length on canopy height (-)

By canopy height and LAI
- `LAI`       : Leaf area index (-)
- `cd`        : Mean drag coefficient for individual leaves. Defaults to 0.2.
- `hs`        : roughness length of the soil surface (m).

By wind profile
- `wind`      : Wind speed at height zr (m s-1)
- `ustar`     : Friction velocity (m s-1)
- `zr`        : Instrument (reference) height (m)
- `d`         : Zero-plane displacement height (-)
- `psi_m`     : value of the stability function for heat, see [`stability_correction`](@ref)

Another variant estimates of `psi_m` by [`stability_correction`](@ref), which
requires further input arguments.
For convenience, these arguments can be provided using a DataFrame, however, the result
then is not type stable.

# Details
The two main roughness parameters, the displacement height (d)
and the roughness length for momentum (z0m) can be estimated from simple
empirical relationships with canopy height (zh). If `method = RoughnessCanopyHeight()`,
the following formulas are used:

``d = frac_d * zh``

``z0m = frac_{z0m} * zh``

where ``frac_d`` defaults to 0.7 and ``frac_{z0m}`` to 0.1.

Alternatively, d and z0m can be estimated from both canopy height and LAI
(If `method = RoughnessCanopyHeightLAI()`).
Based on data from Shaw & Pereira 1982, Choudhury & Monteith 1988 proposed
the following semi-empirical relations:

``X = cd * LAI``

``d = 1.1 * zh * ln(1 + X^{1/4})``

``z0m = hs + 0.3 * zh * X^{1/2}``   for ``0 <= X <= 0.2``

``z0m = hs * zh * (1 - d/zh)``   for ``0.2 < X``

If `method = Roughness_wind_profile()`, z0m is estimated by solving
the wind speed profile for z0m:

``z0m = median((zr - d) * exp(-k*wind / ustar - psi_m)``

By default, d in this equation is fixed to 0.7*zh, but can be set to any
other value.

# Value
a NamedTuple with the following components:
- `d`: Zero-plane displacement height (m)
- `z0m`: Roughness length for momentum (m)
- `z0m_se`: Only if `method = wind_profile`: Standard Error of the median for z0m (m)

# References
- Choudhury, B. J., Monteith J_L., 1988: A four-layer model for the heat
  budget of homogeneous land surfaces. Q. J. R. Meteorol. Soc. 114, 373-398.
- Shaw, R. H., Pereira, A., 1982: Aerodynamic roughness of a plant canopy:
  a numerical experiment. Agricultural Meteorology, 26, 51-65.

# See also
[`wind_profile`](@ref)

```jldoctest; output = false
using DataFrames
# estimate d and z0m from canopy height for a dense (LAI=5) and open (LAI=2) canopy
zh = 25.0
roughness_parameters(RoughnessCanopyHeightLAI(),zh,5)
roughness_parameters(RoughnessCanopyHeightLAI(),zh,2)

# fix d to 0.7*zh and estimate z0m from the wind profile
df = DataFrame(
  Tair=[25.0,25.0,25.0],pressure=100.0,wind=[3.0,4.0,5.0],
  ustar=[0.5,0.6,0.65],H=200.0)
roughness_parameters(Roughness_wind_profile(),df;zh,zr=40,d=0.7*zh)

# assume d = 0.8*zh
rp = roughness_parameters(Roughness_wind_profile(),df;zh,zr=40,d=0.8*zh)
≈(rp.z0m, 0.55, rtol=0.1)
# output
true
```
"""
function roughness_parameters(::RoughnessCanopyHeight, zh; frac_d=0.7, frac_z0m=0.1)
  d      = frac_d*zh
  z0m    = frac_z0m*zh
  z0m_se = missing
  (;d, z0m, z0m_se)
end

function roughness_parameters(::RoughnessCanopyHeightLAI, zh, LAI;
  cd=0.2, hs=0.01)
  X = cd * LAI
  d = 1.1 * zh * log(1 + X^(1/4))
  z0m = ifelse(0 <= X <= 0.2, hs + 0.3 * X^(1/2), 0.3 * zh * (1 - d/zh))
  z0m_se = missing
  (;d, z0m, z0m_se)
end

function roughness_parameters(::Roughness_wind_profile, ustar::AbstractVector{FT}, wind, psi_m;
  zh, zr, d = 0.7*zh, constants=BigleafConstants()
  ) where FT
  FT == Missing && return((d = d, z0m = missing, z0m_se = missing))
  z0m_all = allowmissing(@. (nonmissingtype(FT)(zr) - nonmissingtype(FT)(d)) *
    exp(-nonmissingtype(FT)(constants.k)*wind / ustar - psi_m))
  #z0m_all[(z0m_all .> zh)] .= missing # problems with missings
  replace!(x -> !ismissing(x) && x > zh ? missing : x, z0m_all)
  nval = sum(.!ismissing.(z0m_all))
  z0m    = median(skipmissing(z0m_all))
  z0m_se0 = (std(skipmissing(z0m_all)) / sqrt(nval))
  # http://influentialpoints.com/Training/standard_error_of_median.htm
  z0m_se = oftype(z0m_se0,constants.se_median) * z0m_se0
  (;d, z0m, z0m_se)
end

function roughness_parameters(method::Roughness_wind_profile,
  ustar::AbstractVector, wind, Tair, pressure, H; zh, zr, d = 0.7*zh,
  stab_formulation=Dyer1970(), constants=BigleafConstants(), kwargs...
  )
  # psi_m = Tables.Columns(stability_correction.(
  #   zr, d, Tair, pressure, ustar, H; stab_formulation, constants)).psi_m
  psi_m = getindex.(stability_correction.(
    zr, d, Tair, pressure, ustar, H; stab_formulation, constants), :psi_m)
  # https://discourse.julialang.org/t/eltype-changed-during-broadcast-with-missing/71312?u=bgctw
  res_eltype = Union{eltype(ustar),eltype(wind),eltype(Tair),eltype(pressure),eltype(H)}
  psi_m = convert(Vector{res_eltype}, psi_m)::Vector{res_eltype}
  rp = roughness_parameters(method, ustar, wind, psi_m; zh, zr, d, constants, kwargs...)
end

function roughness_parameters(method::Roughness_wind_profile, df::DFTable;
  psi_m = nothing, stab_formulation=Dyer1970(), kwargs...
  )
  !isnothing(psi_m) && return(roughness_parameters(
    method, df.ustar, df.wind, psi_m; kwargs...))
  if stab_formulation isa NoStabilityCorrection
    psi_m = 0.0
    roughness_parameters(method, df.ustar, df.wind, psi_m; kwargs...)
  else
    roughness_parameters(method, df.ustar, df.wind, df.Tair, df.pressure, df.H;
      stab_formulation, kwargs...)
  end
end

"""
    wind_profile(z::Number, ustar, d, z0m, psi_m = zero(z); constants)
    wind_profile(z, df::AbstractDataFrame, d, z0m, psi_m::AbstractVector; constants)
    wind_profile(z, df::DFTable, d, z0m; psi_m = nothing,
      stab_formulation = Dyer1970(), constants)

Wind speed at a given height above the canopy estimated from single-level
measurements of wind speed.

# Arguments
- `z`         : Height above ground for which wind speed is calculated.
- `ustar`     : Friction velocity (m s-1)
- `d`         : Zero-plane displacement height (-)
- `z0m`       : Roughness length (m)
- `constants=`[`BigleafConstants`](@ref)`()`
For DataFrame variant with supplying stability_parameter
- `df`:         : DataFrame with columns
  - `ustar`     : Friction velocity (m s-1)
- `psi_m`     : value of the stability function for heat, see [`stability_correction`](@ref)
  Pass `psi_m = 0.0` to neglect stability correction.
For DataFrame varinat where psi_m is to be estimated
- additional columns in `df`:
  - `Tair`, `pressure`, `H` : see [`stability_correction`](@ref)

# Details
The underlying assumption is the existence of a logarithmic wind profile
above the height d + z0m (the height at which wind speed mathematically reaches zero
according to the Monin-Obhukov similarity theory).
In this case, the wind speed at a given height z is given by:

``u(z) = (ustar/k) * (ln((z - d) / z0m) - \\psi_m``

The roughness parameters zero-plane displacement height (d) and roughness length (z0m)
can be approximated from [`roughness_parameters`](@ref).
``\\psi_m`` is the stability correction. Set it to zero (not recommended) to neglect
statbility correction. By default it is estimated from wind profile using
[`stability_correction`](@ref)

# Note
Note that this equation is only valid for z >= d + z0m, and it is not
meaningful to calculate values closely above d + z0m. All values in `heights`
smaller than d + z0m will return 0.

# Value
wind speed at given height `z`.

# References
- Monteith, J_L., Unsworth, M_H., 2008: Principles of Environmental Physics.
  3rd edition. Academic Press, London.
- Newman, J_F., Klein, P_M., 2014: The impacts of atmospheric stability on
  the accuracy of wind speed extrapolation methods. Resources 3, 81-105.

# See also
[`roughness_parameters`](@ref)

```jldoctest; output = false
using DataFrames
heights = 18:2:40  # heights above ground for which to calculate wind speed
df = DataFrame(
  Tair=25.0,pressure=100.0,wind=[3.0,4,5],ustar=[0.5,0.6,0.65],H=[200.0,230,250])
zr=40;zh=25;d=16
# z0m and MOL are independent of height, compute before
MOL = Monin_Obukhov_length.(df.Tair, df.pressure, df.ustar, df.H)
z0m = roughness_parameters(
  Roughness_wind_profile(), df.ustar, df.wind, df.Tair, df.pressure, df.H; zh, zr).z0m
ws = map(heights) do z
  wind_profile(z,df,d,z0m; MOL)
end
using Plots # plot wind profiles for the three rows in df
plot(first.(ws), heights, ylab = "height (m)", xlab = "wind speed (m/s)", legend=:topleft)
plot!(getindex.(ws, 2), heights)
plot!(getindex.(ws, 3), heights)
nothing
# output
```
"""
function wind_profile(
  z::Number, ustar::Union{Missing,FT}, d, z0m, psi_m; constants=BigleafConstants()
  ) where FT<:Number
  wind_heights = max(FT(0),(ustar / FT(constants.k)) * (log(max(FT(0),(FT(z) -FT(d))) / FT(z0m)) - psi_m))
  wind_heights
end

function wind_profile(z::Number, ustar::Union{Missing,Number}, d, z0m, Tair,pressure,H,
  stab_formulation=Dyer1970(), constants=BigleafConstants())
  psi_m = stability_correction(
    z,d, Tair,pressure,ustar,H; stab_formulation, constants).psi_m
  wind_profile(z, ustar, d, z0m, psi_m)
end

function wind_profile(z::Number, df::DFTable, d, z0m; psi_m = nothing, MOL = nothing,
  stab_formulation = Dyer1970(), constants =BigleafConstants()
  )
  if isnothing(psi_m)
    if isnothing(MOL)
      psi_m = stability_correction(df; z, d, stab_formulation, constants).psi_m
    else
      zeta = stability_parameter.(z,d,MOL)
      psi_m = getindex.(stability_correction.(zeta), :psi_m)
    end
  end
  #wind_profile(df, z, d, z0m, psi_m; constants)
  windz = wind_profile.(z, df.ustar, d, z0m, psi_m; constants)
end


text:
module Compose

using Colors
using IterTools
using DataStructures
using Measures
using Requires
using Dates
using Printf
using Base.Iterators
using Statistics
import JSON

import Base: length, isempty, getindex, setindex!,
    display, show, convert, zero, isless, max, fill, size, copy,
    min, max, abs, +, -, *, /, ==
import Measures: resolve, w, h

export compose, compose!, Context, UnitBox, AbsoluteBoundingBox,
        Rotation, Mirror, Shear,
       ParentDrawContext, context, ctxpromise, table, set_units!, minwidth, minheight,
       text_extents, max_text_extents, polygon, ngon, star, xgon, bezigon,
       line, rectangle, circle, arc, sector, ellipse, text, curve, bitmap,
       stroke, fill, strokedash, strokelinecap, arrow, strokelinejoin,
       linewidth, visible, fillopacity, strokeopacity, clip, points,
       font, fontsize, svgid, svgclass, svgattribute, jsinclude, jscall, Measure,
       inch, mm, cm, pt, px, cx, cy, sx, sy, w, h, hleft, hcenter, hright, vtop, vcenter,
       vbottom, SVG, SVGJS, PGF, PNG, PS, PDF, draw, pad, pad_inner, pad_outer,
       hstack, vstack, gridstack, LineCapButt, LineCapSquare, LineCapRound,
       CAIROSURFACE, introspect, set_default_graphic_size, set_default_jsmode,
       boundingbox, Patchable

abstract type Backend end

"""
Some backends can more efficiently draw forms by batching. If so, they
shuld define a similar method that returns true.
"""
canbatch(::Backend) = false

# Allow users to supply strings without deprecation warnings
parse_colorant(c::Colorant) = c
parse_colorant(str::AbstractString) = parse(Colorant, str)
parse_colorant(c::Union{Tuple,Vector}) = [parse_colorant(x) for x in c]
parse_colorant(c...) = parse_colorant(c)
@deprecate parse_colorant_vec(c...) parse_colorant(c)

include("misc.jl")
include("measure.jl")
include("list.jl")

# Every graphic in Compose consists of a tree.
abstract type ComposeNode end

# Used to mark null child pointers
struct NullNode <: ComposeNode end
nullnode = NullNode()

include("form.jl")
include("property.jl")
include("container.jl")
include("batch.jl")
include("table.jl")
include("stack.jl")

# How large to draw graphics when not explicitly drawing to a backend
default_graphic_width = sqrt(2)*10*cm
default_graphic_height = 10cm

function set_default_graphic_size(width::MeasureOrNumber,
                                  height::MeasureOrNumber)
    global default_graphic_width
    global default_graphic_height
    default_graphic_width = x_measure(width)
    default_graphic_height = y_measure(height)
    nothing
end

default_graphic_format = :html

function set_default_graphic_format(fmt::Symbol)
    fmt in [:html, :png, :svg, :pdf, :ps, :pgf] || error("$(fmt) is not a supported plot format")
    global default_graphic_format
    default_graphic_format = fmt
    nothing
end

# Default means to include javascript dependencies in the SVGJS backend.
default_jsmode = :embed

function set_default_jsmode(mode::Symbol)
    global default_jsmode
    if mode in [:none, :exclude, :embed, :linkabs, :linkrel]
        default_jsmode = mode
    else
        error("$(mode) is not a valid jsmode")
    end
    nothing
end

function default_mime()
    if default_graphic_format == :png
        "image/png"
    elseif default_graphic_format == :svg
        "image/svg+xml"
    elseif default_graphic_format == :html
        "text/html"
    elseif default_graphic_format == :ps
        "application/postscript"
    elseif default_graphic_format == :pdf
        "application/pdf"
    elseif default_graphic_format == :pgf
        "application/x-tex"
    else
        ""
    end
end

# Default property values
default_font_family = "Helvetica Neue,Helvetica,Arial,sans"
default_font_size = 11pt
default_line_width = 0.3mm
default_stroke_color = nothing
default_fill_color = colorant"black"

# If Cairo is not available, throw an error when trying to save with a Cairo backend
missing_cairo_error(backend::String, invocation::String=backend) =
    """
    The Cairo and Fontconfig packages are necessary for saving as $backend.
    Add them with the package manager if necessary, then run `import Cairo,
    Fontconfig` before invoking `$invocation`.
    """
missing_cairo_error(m::MIME) =
    missing_cairo_error(string(m), "show(::IO, ::MIME\"$m\", ::Context)")

for backend in [:PNG, :PS, :PDF]
  docstr = missing_cairo_error(string(backend))
  @eval @doc $docstr $backend(args...; kwargs...) = error(missing_cairo_error(string($backend)))
end

CairoMIME = Union{MIME"image/png", MIME"application/ps", MIME"application/pdf"}
show(io::IO, m::CairoMIME, ctx::Context) = error(missing_cairo_error(m))

include("svg.jl")
include("pgf_backend.jl")

# If available, pango and fontconfig are used to compute text extents and match
# fonts. Otherwise a simplistic pure-julia fallback is used.

include("fontfallback.jl")

const pango_cairo_ctx = Ref{Ptr}(C_NULL)
const pango_cairo_fm = Ref{Ptr}(C_NULL)
const pangolayout = Ref{Any}(nothing)

function link_fontconfig()
    @debug "Loading Fontconfig backend into Compose.jl"
    pango_cairo_ctx[] = C_NULL

    ccall((:g_type_init, libgobject), Cvoid, ())
    pango_cairo_fm[]  = ccall((:pango_cairo_font_map_new, libpangocairo),
                               Ptr{Cvoid}, ())
    pango_cairo_ctx[] = ccall((:pango_font_map_create_context, libpango),
                               Ptr{Cvoid}, (Ptr{Cvoid},), pango_cairo_fm[])
    pangolayout[] = PangoLayout()
end

function link_cairo()
    @debug "Loading Cairo backend into Compose.jl"
    include(joinpath(@__DIR__,"cairo_backends.jl"))
    include(joinpath(@__DIR__,"immerse_backend.jl"))
end

function __init__()
    @require Cairo="159f3aea-2a34-519c-b102-8c37f9878175" link_cairo()
    @require Fontconfig="186bb1d3-e1f7-5a2c-a377-96d770f13627" begin
      include("pango.jl")
      link_fontconfig()
    end
end

show(io::IO, m::MIME"text/html", ctx::Context) =
    draw(SVGJS(io, default_graphic_width, default_graphic_height, false,
               jsmode=default_jsmode), ctx)

show(io::IO, m::MIME"image/svg+xml", ctx::Context) =
    draw(SVG(io, default_graphic_width, default_graphic_height, false), ctx)

function pad_outer(c::Context,
                   left_padding::MeasureOrNumber,
                   right_padding::MeasureOrNumber,
                   top_padding::MeasureOrNumber,
                   bottom_padding::MeasureOrNumber)

    left_padding   = size_measure(left_padding)
    right_padding  = size_measure(right_padding)
    top_padding    = size_measure(top_padding)
    bottom_padding = size_measure(bottom_padding)

    root = context(c.box.x0[1], c.box.x0[1],
                   c.box.a[1] + left_padding + right_padding,
                   c.box.a[2] + top_padding + bottom_padding,
                   minwidth=c.minwidth,
                   minheight=c.minheight)
    c = copy(c)
    c.box = BoundingBox(left_padding, top_padding,
                        1w - left_padding - right_padding,
                        1h - top_padding - bottom_padding)
    return compose!(root, c)
end

pad_outer(c::Context, padding::MeasureOrNumber) =
        pad_outer(c, padding, padding, padding, padding)

pad_outer(cs::Vector{Context},
                   left_padding::MeasureOrNumber,
                   right_padding::MeasureOrNumber,
                   top_padding::MeasureOrNumber,
                   bottom_padding::MeasureOrNumber) =
        map(c -> pad_outer(c, left_padding, right_padding, top_padding, bottom_padding), cs)

pad_outer(cs::Vector{Context}, padding::MeasureOrNumber) =
        pad_outer(cs, padding, padding, padding, padding)

function pad_inner(c::Context,
                   left_padding::MeasureOrNumber,
                   right_padding::MeasureOrNumber,
                   top_padding::MeasureOrNumber,
                   bottom_padding::MeasureOrNumber)

    left_padding   = size_measure(left_padding)
    right_padding  = size_measure(right_padding)
    top_padding    = size_measure(top_padding)
    bottom_padding = size_measure(bottom_padding)

    root = context(c.box.x0[1], c.box.x0[2],
                   c.box.a[1], c.box.a[2],
                   minwidth=c.minwidth,
                   minheight=c.minheight)
    c = copy(c)
    c.box = BoundingBox(left_padding, top_padding,
                        1w - left_padding - right_padding,
                        1h - top_padding - bottom_padding)
    return compose!(root, c)
end

pad_inner(c::Context, padding::MeasureOrNumber) =
        pad_inner(c, padding, padding, padding, padding)

pad_inner(cs::Vector{Context},
                   left_padding::MeasureOrNumber,
                   right_padding::MeasureOrNumber,
                   top_padding::MeasureOrNumber,
                   bottom_padding::MeasureOrNumber) =
        map(c -> pad_inner(c, left_padding, right_padding, top_padding, bottom_padding), cs)

pad_inner(cs::Vector{Context}, padding::MeasureOrNumber) =
        pad_inner(cs, padding, padding, padding, padding)

function gridstack(cs::Matrix{Context})
    m, n = size(cs)
    t = Table(m, n, 1:m, 1:n, x_prop=ones(n), y_prop=ones(m))
    for i in 1:m, j in 1:n
        t[i, j] = [cs[i, j]]
    end
    return compose!(context(), t)
end

const pad = pad_outer

end # module Compose


text:
using Test
using JSON
using MianalyzerBackend
const CorrelationsController = MianalyzerBackend.CorrelationsController

include("../TestSupport.jl")

@testset "CorrelationsController" begin
  @testset ".call" begin
    request = Dict([("records", [[1.0, 100.0], [2.0, 200.0], [3.0, 300.0]])])

    actual = JSON.parse(IOBuffer(CorrelationsController.call(request)))

    @test actual["combs"] == [[1,2]]
    @test actual["correlations"] == [1.0]
    @test actual["spearman_correlations"] == [1.0]
  end
end


text:
import Base:∪, ∩, isless, size

abstract type AbstractGate end

struct Gate{T,V} <: AbstractGate
    low::T
    high::V
    function Gate(;low::T=nothing, high::V=nothing) where {T,V}
        if !isnothing(low) && !isnothing(high)
            @assert low < high "Lower limit must be lower than higher: $low ≮ $high"
        end
        new{T,V}(low, high)
    end
end
#Gate(;low=nothing, high=nothing) = Gate(low, high)
Gate(low, high) = Gate(low=low, high=high)
Gate(mean; width=50) = Gate(mean-width, mean+width)

(gate::Gate{T,V} where {T,V})(x)::Bool        = gate.low < x < gate.high
(gate::Gate{Nothing,V} where {V})(x)::Bool    = x < gate.high
(gate::Gate{T,Nothing} where {T})(x)::Bool    = gate.low < x
function (gate::Gate{Nothing,Nothing})(x)::Bool
    return true
    #throw(ArgumentError("Set the gate limit(s)"))
end

∩(g1::Gate{Nothing, Nothing}, g2::Gate{Nothing, Nothing}) = g1
function ∩(g1::Gate, g2::Gate)
    low = if !isnothing(g1.low) && !isnothing(g2.low)
        g1.low > g2.low ? g1.low : g2.low
        elseif !isnothing(g1.low)
            g1.low
        else
            g2.low
    end
    high = if !isnothing(g1.high) && !isnothing(g2.high)
        g1.high < g2.high ? g1.high : g2.high
        elseif !isnothing(g1.high)
            g1.high
        else
            g2.high
    end
    Gate(low=low, high=high)
end

center(g::Gate) = g.low + 1/2*(g.high - g.low)
isless(g1::Gate, g2::Gate) = center(g1) < center(g2)

struct Point
    x::Float64
    y::Float64
end
Point(p::Union{AbstractVector, Tuple{Number, Number}}) = Point(p[1], p[2])

struct Gate2D <: AbstractGate
    lowerleft::Point
    lowerright::Point
    higherright::Point
    higherleft::Point
    function Gate2D(p1::Point, p2::Point, p3::Point, p4::Point)
        p = [p1, p2, p3, p4]
        if !isrectangle(p)
            throw(ArgumentError("Points do not constitute a rectangle"))
        end
        lessy(p1, p2) = p1.y < p2.y
        lessynx(p1, p2) = p1.y <= p2.y && p1.x < p2.x
        p = sort(sort(p, lt=lessy), lt=lessynx)
        new(p[1], p[2], p[4], p[3])
    end
end

Gate2D(points::AbstractVector{Point}) = Gate2D(points...)

function Gate2D(ll::T, lr::T, hr::T, hl::T) where T <: Union{AbstractVector, Tuple{Number, Number}}
    Gate2D(Point(ll), Point(lr), Point(hr), Point(hl))
end

function Gate2D(center::Point; width=50)
    d = width/2
    Gate2D((center.x - d, center.y-d),
         (center.x + d, center.y - d),
         (center.x + d, center.y + d),
         (center.x - d, center.y + d))
end

function Gate2D(xmin::Number, xmax::Number, ymin::Number, ymax::Number)
    Gate2D((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax))
end

(gate::Gate2D)(p::Point) = gate(p.x, p.y)

function (gate::Gate2D)(x, y)::Bool
    if gate.lowerleft.x < x < gate.lowerright.x && gate.lowerleft.y < y < gate.higherleft.y
        true
    else
        false
    end
end

function isrectangle(points::AbstractVector{Point})::Bool
    if length(points) != 4
        return false
    end

    for i in 1:4
        p = points[i]
        found_x = 0
        found_y = 0
        for j in 1:4
            i == j && continue
            p2 = points[j]
            if p2.x ≈ p.x
                found_x += 1
            end
            if p2.y ≈ p.y
                found_y += 1
            end
        end
        if found_x != 1 || found_y != 1
            return false
        end
    end
    true
end

size(gate::Gate2D) = gate.lowerleft.x, gate.lowerright.x, gate.lowerleft.y, gate.higherleft.y
center(gate::Gate2D) = (gate.lowerleft.x + 1/2*(gate.lowerright.x - gate.lowerleft.x),
                        gate.lowerleft.y + 1/2*(gate.higherright.y - gate.lowerleft.y))

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
#HDF5 Plots: Save/replay plots to/from HDF5
#-------------------------------------------------------------------------------

#==Usage
===============================================================================
Write to .hdf5 file using:
    p = plot(...)
   Plots.hdf5plot_write(p, "plotsave.hdf5")

Read from .hdf5 file using:
    pyplot() #Must first select backend
    pread = Plots.hdf5plot_read("plotsave.hdf5")
    display(pread)
==#

#==TODO
===============================================================================
 1. Support more features.
    - GridLayout known not to be working.
 2. Improve error handling.
    - Will likely crash if file format is off.
 3. Save data in a folder parallel to "plot".
    - Will make it easier for users to locate data.
    - Use HDF5 reference to link data?
 4. Develop an actual versioned file format.
    - Should have some form of backward compatibility.
    - Should be reliable for archival purposes.
 5. Fix construction of plot object with hdf5plot_read.
    - Layout doesn't seem to get transferred well (ex: `Plots._examples[40]`).
    - Not building object correctly when backends do not natively support
      a certain feature (ex: :steppre)
    - No support for CategoricalArrays.* structures. But they appear to be
      brought into `Plots._examples[25,30]` through DataFrames.jl - so we can't
      really reference them in this code.
==#
"""
    _hdf5_implementation

Create module (namespace) for implementing HDF5 "plots".
(Avoid name collisions, while keeping names short)
"""
module _hdf5_implementation #Tools required to implements HDF5 "plots"

import Dates

#Plots.jl imports HDF5 to main:
import ..HDF5
import ..HDF5: Group, Dataset

import ..Colors, ..Colorant
import ..PlotUtils.ColorSchemes.ColorScheme

import ..HDF5Backend
import ..HDF5PLOT_MAP_STR2TELEM, ..HDF5PLOT_MAP_TELEM2STR
import ..HDF5Plot_PlotRef, ..HDF5PLOT_PLOTREF
import ..BoundingBox, ..Extrema, ..Length
import ..RecipesPipeline.datetimeformatter
import ..PlotUtils.ColorPalette, ..PlotUtils.CategoricalColorGradient, ..PlotUtils.ContinuousColorGradient
import ..Surface, ..Shape, ..Arrow
import ..GridLayout, ..RootLayout
import ..Font, ..PlotText, ..SeriesAnnotations
import ..Axis, ..Subplot, ..Plot
import ..AKW, ..KW, ..DefaultsDict
import .._axis_defaults
import ..plot, ..plot!

#Types that already have built-in HDF5 support (just write out natively):
const HDF5_SupportedTypes = Union{Number, String}
#TODO: Types_HDF5Support

#Dispatch types:
struct CplxTuple; end #Identifies a "complex" tuple structure (not merely numbers)
#HDF5 reader will auto-detect type correctly:
struct HDF5_AutoDetect; end #See HDF5_SupportedTypes

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

if length(HDF5PLOT_MAP_TELEM2STR) < 1
    #Possible element types of high-level data types:
    #(Used to add type information as an HDF5 string attribute)
    #(Also used to dispatch appropriate read function through _read_typed())
    _telem2str = Dict{String, Type}(
        "NOTHING" => Nothing,
        "SYMBOL" => Symbol,
        "RGBA" => Colorant, #Write out any Colorant to an #RRGGBBAA string
        "TUPLE" => Tuple,
        "CTUPLE" => CplxTuple, #Tuple of complex structures

        "EXTREMA" => Extrema,
        "LENGTH" => Length,
        "ARRAY" => Array, #Array{Any} (because Array{T<:Union{Number, String}} natively supported by HDF5)

        #Sub-structure types:
        "T_DATETIMEFORMATTER" => typeof(datetimeformatter),

        #Sub-structure types:
        "DEFAULTSDICT" => DefaultsDict,
        "FONT" => Font,
        "BOUNDINGBOX" => BoundingBox,
        "GRIDLAYOUT" => GridLayout,
        "ROOTLAYOUT" => RootLayout,
        "SERIESANNOTATIONS" => SeriesAnnotations,
        "PLOTTEXT" => PlotText,
        "SHAPE" => Shape,
        "ARROW" => Arrow,
        "COLORSCHEME" => ColorScheme,
        "COLORPALETTE" => ColorPalette,
        "CONT_COLORGRADIENT" => ContinuousColorGradient,
        "CAT_COLORGRADIENT" => CategoricalColorGradient,
        "AXIS" => Axis,
        "SURFACE" => Surface,
        "SUBPLOT" => Subplot,
    )
    merge!(HDF5PLOT_MAP_STR2TELEM, _telem2str) #Faster to create than push!()??
    merge!(HDF5PLOT_MAP_TELEM2STR, Dict{Type, String}(v=>k for (k,v) in HDF5PLOT_MAP_STR2TELEM))
end

#==Helper functions
===============================================================================#

h5plotpath(plotname::String) = "plots/$plotname"

#Version info
#NOTE: could cache output, but we seem to not want const declarations in backend files.
function _get_Plots_versionstr()
    #Adds to load up time... Maybe a more efficient way??
    try
        deps = Pkg.dependencies()
        uuid = Base.UUID("91a5bcdd-55d7-5caf-9e0b-520d859cae80") #Plots.jl
        vinfo = deps[uuid].version
        return "Source: Plots.jl v$vinfo"
    catch
        now = string(Dates.now()) #Use time in case it can help recover plot
        return "Source: Plots.jl v? - $now"
    end
end

function _hdf5_merge!(dest::AKW, src::AKW)
    for (k, v) in src
        if isa(v, Axis)
            _hdf5_merge!(dest[k].plotattributes, v.plotattributes)
        else
            dest[k] = v
        end
    end
    return
end

#_type_for_map returns the type to use with HDF5PLOT_MAP_TELEM2STR[], in case it is not concrete:
_type_for_map(::Type{T}) where T = T #Catch-all
_type_for_map(::Type{T}) where T<:BoundingBox = BoundingBox
_type_for_map(::Type{T}) where T<:ColorScheme = ColorScheme
_type_for_map(::Type{T}) where T<:Surface = Surface

#==Read/write things like type name in attributes
===============================================================================#
function _write_datatype_attr(ds::Union{Group, Dataset}, ::Type{T}) where T
    typestr = HDF5PLOT_MAP_TELEM2STR[T]
    HDF5.attributes(ds)["TYPE"] = typestr
end
function _read_datatype_attr(ds::Union{Group, Dataset})
    if !Base.haskey(HDF5.attributes(ds), "TYPE")
        return HDF5_AutoDetect
    end

    typestr = HDF5.read(HDF5.attributes(ds)["TYPE"])
    return HDF5PLOT_MAP_STR2TELEM[typestr]
end

#Type parameter attributes:
function _write_typeparam_attr(ds::Dataset, v::Length{T}) where T
    HDF5.attributes(ds)["TYPEPARAM"] = string(T) #Need to add units for Length
end
_read_typeparam_attr(ds::Dataset) = HDF5.read(HDF5.attributes(ds)["TYPEPARAM"])

function _write_length_attr(grp::Group, v::Vector) #of a vector
    HDF5.attributes(grp)["LENGTH"] = length(v)
end
_read_length_attr(::Type{Vector}, grp::Group) = HDF5.read(HDF5.attributes(grp)["LENGTH"])

function _write_size_attr(grp::Group, v::Array) #of an array
    HDF5.attributes(grp)["SIZE"] = [size(v)...]
end
_read_size_attr(::Type{Array}, grp::Group) = tuple(HDF5.read(HDF5.attributes(grp)["SIZE"])...)

#==_write_typed(): Simple (leaf) datatypes. (Labels with type name.)
===============================================================================#
#= No: write out struct instead!
function _write_typed(grp::Group, name::String, v::T) where T
    tstr = string(T)
    path = HDF5.name(grp) * "/" * name
    @info("Type not supported: $tstr\npath: $path")
    return
end
=#
#Default behaviour: Assumes value is supported by HDF5 format
function _write_typed(grp::Group, name::String, v::HDF5_SupportedTypes)
    grp[name] = v
    return #No need to _write_datatype_attr
end
function _write_typed(grp::Group, name::String, v::Nothing)
    grp[name] = "nothing" #Redundancy check/easier to read HDF5 file
    _write_datatype_attr(grp[name], Nothing)
end
function _write_typed(grp::Group, name::String, v::Symbol)
    grp[name] = String(v)
    _write_datatype_attr(grp[name], Symbol)
end
function _write_typed(grp::Group, name::String, v::Colorant)
    vstr = "#" * Colors.hex(v, :RRGGBBAA)
    grp[name] = vstr
    _write_datatype_attr(grp[name], Colorant)
end
function _write_typed(grp::Group, name::String, v::Extrema)
    grp[name] = [v.emin, v.emax] #More compact than writing struct
    _write_datatype_attr(grp[name], Extrema)
end
function _write_typed(grp::Group, name::String, v::Length)
    grp[name] = v.value
    _write_datatype_attr(grp[name], Length)
    _write_typeparam_attr(grp[name], v)
end
function _write_typed(grp::Group, name::String, v::typeof(datetimeformatter))
    grp[name] = string(v) #Just write something that helps reader
   _write_datatype_attr(grp[name], typeof(datetimeformatter))
end
function _write_typed(grp::Group, name::String, v::Array{T}) where T<:Number #Default for arrays
    grp[name] = v
    return #No need to _write_datatype_attr
end
function _write_typed(grp::Group, name::String, v::AbstractRange)
    _write_typed(grp, name, collect(v)) #For now
end

#== Helper functions for writing complex data structures
===============================================================================#

#Write an array using HDF5 hierarchy (when not using simple numeric eltype):
function _write_harray(grp::Group, name::String, v::Array)
    sgrp = HDF5.create_group(grp, name)
    sz = size(v)
    lidx = LinearIndices(sz)

    for iter in eachindex(v)
        coord = lidx[iter]
        elem = v[iter]
        idxstr = join(coord, "_")
        _write_typed(sgrp, "v$idxstr", elem)
    end

    _write_size_attr(sgrp, v)
end

#Write Dict without tagging with type:
function _write(grp::Group, name::String, d::AbstractDict)
    sgrp = HDF5.create_group(grp, name)
    for (k, v) in d
        kstr = string(k)
        _write_typed(sgrp, kstr, v)
    end
    return
end

#Write out arbitrary `struct`s:
function _writestructgeneric(grp::Group, obj::T) where T
    for fname in fieldnames(T)
       v = getfield(obj, fname)
       _write_typed(grp, String(fname), v)
    end
    return
end

#==_write_typed(): More complex structures. (Labels with type name.)
===============================================================================#

#Catch-all (default behaviour for `struct`s):
function _write_typed(grp::Group, name::String, v::T) where T
    #NOTE: need "name" parameter so that call signature is same with built-ins
    MT = _type_for_map(T)
    try #Check to see if type is supported
        typestr = HDF5PLOT_MAP_TELEM2STR[MT]
    catch
        @warn("HDF5Plots does not yet support structs of type `$MT`\n\n$grp")
        return
    end

    #If attribute is supported and no writer is defined, then this should work:
    objgrp = HDF5.create_group(grp, name)
    _write_datatype_attr(objgrp, MT)
    _writestructgeneric(objgrp, v)
end

function _write_typed(grp::Group, name::String, v::Array{T}) where T
    _write_harray(grp, name, v)
    _write_datatype_attr(grp[name], Array) #{Any}
end

function _write_typed(grp::Group, name::String, v::Tuple, ::Type{ELT}) where ELT<: Number #Basic Tuple
    _write_typed(grp, name, [v...])
    _write_datatype_attr(grp[name], Tuple)
end
function _write_typed(grp::Group, name::String, v::Tuple, ::Type) #CplxTuple
    _write_harray(grp, name, [v...])
    _write_datatype_attr(grp[name], CplxTuple)
end
_write_typed(grp::Group, name::String, v::Tuple) = _write_typed(grp, name, v, eltype(v))

function _write_typed(grp::Group, name::String, v::Dict)
#=
    tstr = string(Dict)
    path = HDF5.name(grp) * "/" * name
    @info("Type not supported: $tstr\npath: $path")
    return
=#
    #No support for structures with Dicts yet
end
function _write_typed(grp::Group, name::String, d::DefaultsDict) #Typically for plot attributes
    _write(grp, name, d)
    _write_datatype_attr(grp[name], DefaultsDict)
end

function _write_typed(grp::Group, name::String, v::Axis)
    sgrp = HDF5.create_group(grp, name)
    #Ignore: sps::Vector{Subplot}
    _write_typed(sgrp, "plotattributes", v.plotattributes)
    _write_datatype_attr(sgrp, Axis)
end

function _write_typed(grp::Group, name::String, v::Subplot)
    #Not for use in main "Plot.subplots[]" hierarchy.  Just establishes reference with subplot_index.
    sgrp = HDF5.create_group(grp, name)
    _write_typed(sgrp, "index", v[:subplot_index])
    _write_datatype_attr(sgrp, Subplot)
    return
end

function _write_typed(grp::Group, name::String, v::Plot)
    #Don't write plot references
end

#==_write(): Write out more complex structures
NOTE: No need to write out type information (inferred from hierarchy)
===============================================================================#

function _write(grp::Group, sp::Subplot{HDF5Backend})
    _write_typed(grp, "attr", sp.attr)

    listgrp = HDF5.create_group(grp, "series_list")
    _write_length_attr(listgrp, sp.series_list)
    for (i, series) in enumerate(sp.series_list)
        #Just write .plotattributes part:
        _write(listgrp, "$i", series.plotattributes)
    end

    return
end

function _write(grp::Group, plt::Plot{HDF5Backend})
    _write_typed(grp, "attr", plt.attr)

    listgrp = HDF5.create_group(grp, "subplots")
    _write_length_attr(listgrp, plt.subplots)
    for (i, sp) in enumerate(plt.subplots)
        sgrp = HDF5.create_group(listgrp, "$i")
        _write(sgrp, sp)
    end

    return
end

function hdf5plot_write(plt::Plot{HDF5Backend}, path::AbstractString; name::String="_unnamed")
    HDF5.h5open(path, "w") do file
        HDF5.write_dataset(file, "VERSION_INFO", _get_Plots_versionstr())
        grp = HDF5.create_group(file, h5plotpath(name))
        _write(grp, plt)
    end
end

#== _read(): Read data, but not type information.
===============================================================================#

#Types with built-in HDF5 support:
_read(::Type{HDF5_AutoDetect}, ds::Dataset) = HDF5.read(ds)

function _read(::Type{Nothing}, ds::Dataset)
    nstr = "nothing"
    v = HDF5.read(ds)
    if nstr != v
        path = HDF5.name(ds)
        throw(Meta.ParseError("_read(::Nothing, ::Group): Read $v != $nstr:\n$path"))
    end
    return nothing
end
_read(::Type{Symbol}, ds::Dataset) = Symbol(HDF5.read(ds))
_read(::Type{Colorant}, ds::Dataset) = parse(Colorant, HDF5.read(ds))
_read(::Type{Tuple}, ds::Dataset) = tuple(HDF5.read(ds)...)
function _read(::Type{Extrema}, ds::Dataset)
    v = HDF5.read(ds)
    return Extrema(v[1], v[2])
end
function _read(::Type{Length}, ds::Dataset)
    TUNIT = Symbol(_read_typeparam_attr(ds))
    v = HDF5.read(ds)
    T = typeof(v)
    return Length{TUNIT,T}(v)
end
_read(::Type{typeof(datetimeformatter)}, ds::Dataset) = datetimeformatter

#== Helper functions for reading in complex data structures
===============================================================================#

#When type is unknown, _read_typed() figures it out:
function _read_typed(grp::Group, name::String)
    ds = grp[name]
    t = _read_datatype_attr(ds)
    return _read(t, ds)
end

#_readstructgeneric: Needs object values to be written out with _write_typed():
function _readstructgeneric(::Type{T}, grp::Group) where T
    vlist = Array{Any}(nothing, fieldcount(T))
    for (i, fname) in enumerate(fieldnames(T))
        vlist[i] = _read_typed(grp, String(fname))
    end
    return T(vlist...)
end

#Read KW from group:
function _read(::Type{KW}, grp::Group)
    d = KW()
    gnames = names(grp)
    for k in gnames
        try
            v = _read_typed(grp, k)
            d[Symbol(k)] = v
        catch e
            @show e
            @show grp
            @warn("Could not read field $k")
        end
    end
    return d
end

#== _read(): More complex structures.
===============================================================================#

#Catch-all (default behaviour for `struct`s):
_read(T::Type, grp::Group) = _readstructgeneric(T, grp)

function _read(::Type{Array}, grp::Group) #Array{Any}
    sz = _read_size_attr(Array, grp)
    if tuple(0) == sz; return []; end
    result = Array{Any}(undef, sz)
    lidx = LinearIndices(sz)

    for iter in eachindex(result)
        coord = lidx[iter]
        idxstr = join(coord, "_")
        result[iter] = _read_typed(grp, "v$idxstr")
    end

    #Hack: Implicitly make Julia detect element type.
    #      (Should probably write it explicitly to file)
    result = [elem for elem in result] #Potentially make more specific
    return reshape(result, sz)
end

_read(::Type{CplxTuple}, grp::Group) = tuple(_read(Array, grp)...)

function _read(::Type{GridLayout}, grp::Group)
    #parent = _read_typed(grp, "parent") #Can't use generic reader
    parent = RootLayout() #TODO: support parent???
    minpad = _read_typed(grp, "minpad")
    bbox = _read_typed(grp, "bbox")
    grid = _read_typed(grp, "grid")
    widths = _read_typed(grp, "widths")
    heights = _read_typed(grp, "heights")
    attr = KW() #TODO support attr: _read_typed(grp, "attr")

    return GridLayout(parent, minpad, bbox, grid, widths, heights, attr)
end
#Defaults depends on context. So: user must constructs with defaults, then read.
function _read(::Type{DefaultsDict}, grp::Group)
    #User should set DefaultsDict.defaults to one of:
    #   _plot_defaults, _subplot_defaults, _axis_defaults, _series_defaults
    path = HDF5.name(ds)
    @warn("Cannot yet read DefaultsDict using _read_typed():\n    $path\nCannot fully reconstruct plot.")
end
function _read(::Type{Axis}, grp::Group)
    #1st arg appears to be ref to subplots. Seems to work without it.
    return Axis([], DefaultsDict(_read(KW, grp["plotattributes"]), _axis_defaults))
end
function _read(::Type{Subplot}, grp::Group)
    #Not for use in main "Plot.subplots[]" hierarchy.  Just establishes reference with subplot_index.
    idx = _read_typed(grp, "index")
    return HDF5PLOT_PLOTREF.ref.subplots[idx]
end

#== _read(): Main plot structures
===============================================================================#

function _read(grp::Group, sp::Subplot)
    listgrp = HDF5.open_group(grp, "series_list")
    nseries = _read_length_attr(Vector, listgrp)

    for i in 1:nseries
        sgrp = HDF5.open_group(listgrp, "$i")
        seriesinfo = _read(KW, sgrp)

        plot!(sp, seriesinfo[:x], seriesinfo[:y]) #Add data & create data structures
        _hdf5_merge!(sp.series_list[end].plotattributes, seriesinfo)
    end

    #Perform after adding series... otherwise values get overwritten:
    agrp = HDF5.open_group(grp, "attr")
    _hdf5_merge!(sp.attr, _read(KW, agrp))

    return sp
end

function _read_plot(grp::Group)
    listgrp = HDF5.open_group(grp, "subplots")
    n = _read_length_attr(Vector, listgrp)

    #Construct new plot, +allocate subplots:
    plt = plot(layout=n)
    HDF5PLOT_PLOTREF.ref = plt #Used when reading "layout"

    agrp = HDF5.open_group(grp, "attr")
    _hdf5_merge!(plt.attr, _read(KW, agrp))

    for (i, sp) in enumerate(plt.subplots)
        sgrp = HDF5.open_group(listgrp, "$i")
        _read(sgrp, sp)
    end

    return plt
end

function hdf5plot_read(path::AbstractString; name::String="_unnamed")
    HDF5.h5open(path, "r") do file
        grp = HDF5.open_group(file, h5plotpath("_unnamed"))
        return _read_plot(grp)
    end
end

end #module _hdf5_implementation

#==Implement Plots.jl backend interface for HDF5Backend
===============================================================================#

is_marker_supported(::HDF5Backend, shape::Shape) = true

# Create the window/figure for this backend.
function _create_backend_figure(plt::Plot{HDF5Backend})
    #Do nothing
end

# ---------------------------------------------------------------------------

# # this is called early in the pipeline, use it to make the plot current or something
# function _prepare_plot_object(plt::Plot{HDF5Backend})
# end

# ---------------------------------------------------------------------------

# Set up the subplot within the backend object.
function _initialize_subplot(plt::Plot{HDF5Backend}, sp::Subplot{HDF5Backend})
    #Do nothing
end

# ---------------------------------------------------------------------------

# Add one series to the underlying backend object.
# Called once per series
# NOTE: Seems to be called when user calls plot()... even if backend
#       plot, sp.o has not yet been constructed...
function _series_added(plt::Plot{HDF5Backend}, series::Series)
    #Do nothing
end

# ---------------------------------------------------------------------------

# When series data is added/changed, this callback can do dynamic updates to the backend object.
# note: if the backend rebuilds the plot from scratch on display, then you might not do anything here.
function _series_updated(plt::Plot{HDF5Backend}, series::Series)
    #Do nothing
end

# ---------------------------------------------------------------------------

# called just before updating layout bounding boxes... in case you need to prep
# for the calcs
function _before_layout_calcs(plt::Plot{HDF5Backend})
    #Do nothing
end

# ----------------------------------------------------------------

# Set the (left, top, right, bottom) minimum padding around the plot area
# to fit ticks, tick labels, guides, colorbars, etc.
function _update_min_padding!(sp::Subplot{HDF5Backend})
    #Do nothing
end

# ----------------------------------------------------------------

# Override this to update plot items (title, xlabel, etc), and add annotations (plotattributes[:annotations])
function _update_plot_object(plt::Plot{HDF5Backend})
    #Do nothing
end

# ----------------------------------------------------------------

# Display/show the plot (open a GUI window, or browser page, for example).
function _display(plt::Plot{HDF5Backend})
    msg = "HDF5 interface does not support `display()` function."
    msg *= "\nUse `Plots.hdf5plot_write(::String)` method to write to .HDF5 \"plot\" file instead."
    @warn(msg)
    return
end

#==Interface actually required to use HDF5Backend
===============================================================================#

hdf5plot_write(plt::Plot{HDF5Backend}, path::AbstractString) = _hdf5_implementation.hdf5plot_write(plt, path)
hdf5plot_write(path::AbstractString) = _hdf5_implementation.hdf5plot_write(current(), path)
hdf5plot_read(path::AbstractString) = _hdf5_implementation.hdf5plot_read(path)

#Last line


text:
#
# To Do:
# - implement assembleLanczosVectors
# - check slice ranges - change end value by 1?
#

function assemble_lanczos_vecs(lanczos_vectors, linear_comb, norm)
  #if length(lanczos_vectors) != length(linear_comb)
  #  @show length(lanczos_vectors)
  #  @show length(linear_comb)
  #end
  xt = norm * linear_comb[1] * lanczos_vectors[1]
  for i in 2:length(lanczos_vectors)
    xt += norm * linear_comb[i] * lanczos_vectors[i]
  end
  return xt
end

struct ApplyExpInfo
  numops::Int
end

function applyexp(H, tau::Number, x0; kwargs...)
  maxiter = get(kwargs, :maxiter, 30)
  tol = get(kwargs, :tol, 1E-12)
  outputlevel = get(kwargs, :outputlevel, 0)
  beta_tol = get(kwargs, :normcutoff, 1E-7)

  # Initialize Lanczos vectors
  v1 = copy(x0)
  nrm = norm(v1)
  v1 /= nrm
  lanczos_vectors = [v1]

  ElT = promote_type(typeof(tau), eltype(x0))

  bigTmat = zeros(ElT, maxiter + 2, maxiter + 2)

  nmatvec = 0

  v0 = nothing
  beta = 0.0
  for iter in 1:maxiter
    tmat_size = iter + 1

    # Matrix-vector multiplication
    w = H(v1)
    nmatvec += 1

    avnorm = norm(w)
    alpha = dot(w, v1)

    bigTmat[iter, iter] = alpha

    w -= alpha * v1
    if iter > 1
      w -= beta * v0
    end
    v0 = copy(v1)

    beta = norm(w)

    # check for Lanczos sequence exhaustion
    if abs(beta) < beta_tol
      # Assemble the time evolved state
      tmat = bigTmat[1:tmat_size, 1:tmat_size]
      tmat_exp = exp(tau * tmat)
      linear_comb = tmat_exp[:, 1]
      xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm)
      return xt, ApplyExpInfo(nmatvec)
    end

    # update next lanczos vector
    v1 = copy(w)
    v1 /= beta
    push!(lanczos_vectors, v1)
    bigTmat[iter + 1, iter] = beta
    bigTmat[iter, iter + 1] = beta

    # Convergence check
    if iter > 0
      # Prepare extended T-matrix for exponentiation
      tmat_ext_size = tmat_size + 2
      tmat_ext = bigTmat[1:tmat_ext_size, 1:tmat_ext_size]

      tmat_ext[tmat_size - 1, tmat_size] = 0.0
      tmat_ext[tmat_size + 1, tmat_size] = 1.0

      # Exponentiate extended T-matrix
      tmat_ext_exp = exp(tau * tmat_ext)

      ϕ1 = abs(nrm * tmat_ext_exp[tmat_size, 1])
      ϕ2 = abs(nrm * tmat_ext_exp[tmat_size + 1, 1] * avnorm)

      if ϕ1 > 10 * ϕ2
        error = ϕ2
      elseif (ϕ1 > ϕ2)
        error = (ϕ1 * ϕ2) / (ϕ1 - ϕ2)
      else
        error = ϕ1
      end

      if outputlevel >= 3
        @printf("  Iteration: %d, Error: %.2E\n", iter, error)
      end

      if ((error < tol) || (iter == maxiter))
        if (iter == maxiter)
          println("warning: applyexp not converged in $maxiter steps")
        end

        # Assemble the time evolved state
        linear_comb = tmat_ext_exp[:, 1]
        xt = assemble_lanczos_vecs(lanczos_vectors, linear_comb, nrm)

        if outputlevel >= 3
          println("  Number of iterations: $iter")
        end

        return xt, ApplyExpInfo(nmatvec)
      end
    end  # end convergence test
  end # iter

  if outputlevel >= 0
    println("In applyexp, number of matrix-vector multiplies: ", nmatvec)
  end

  return x0
end


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

# %% [markdown]
# StaticArrays version: https://github.com/genkuroki/public/blob/main/0018/HMC%20leapfrog%20StaticArrays.ipynb

# %%
module My

using ConcreteStructs: @concrete
using ForwardDiff: derivative
using Parameters: @unpack
using Random: default_rng

@concrete struct LFProblem ϕ; H; F; dt; nsteps end

"""ϕ should be a potential function."""
function LFProblem(ϕ; dt = 0.05, nsteps = 100)
    H(x, v, param) = v^2/2 + ϕ(x, param)
    F(x, param) = -derivative(x -> ϕ(x, param), x)
    LFProblem(ϕ, H, F, dt, nsteps)
end

"""Numerically solve Hamilton's equation of motion with leapfrog method"""
function solve(lf::LFProblem, x, v, param)
    @unpack F, dt, nsteps = lf
    v = v + F(x, param)*dt/2
    x = x + v*dt
    for _ in 2:nsteps
        v = v + F(x, param)*dt
        x = x + v*dt
    end
    v = v + F(x, param)*dt/2
    x, v
end

@inline function _update!(lf::LFProblem{dim}, x, T, param, rng) where dim
    @unpack H = lf
    v = randn(rng, T)
    xnew, vnew = solve(lf, x, v, param)
    dH = H(xnew, vnew, param) - H(x, v, param)
    rand(rng) ≤ exp(-dH) ? xnew : x
end

"""Hamiltonian Monte Carlo"""
function HMC(lf::LFProblem, param = nothing;
        niters = 10^5, nskips = 1, thin = 1, nwarmups = 0, rng = default_rng(),
        x0 = randn(rng))
    T = typeof(x0)
    x = x0
    for _ in 1:nwarmups
        x = _update!(lf, x, T, param, rng)
    end
    sample = Vector{T}(undef, niters)
    for i in 1:niters
        #or _ in 1:thin
            x = _update!(lf, x, T, param, rng)
        #end
        @inbounds sample[i] = x
    end
    sample
end

end

# %%
using Plots
using QuadGK
using BenchmarkTools

# %%
ϕ(x, param=nothing) = x^2/2
lf = My.LFProblem(ϕ)
@time X = My.HMC(lf)
histogram(X; norm=true, alpha=0.3, label="HMC LF sample")
plot!(x -> exp(-x^2/2)/√(2π), -4, 4, label="std Gaussian", lw=2)

# %%
@btime My.HMC($lf);

# %%
ϕ2(x, param=nothing) = x^4 - 4x^2
Z = quadgk(x -> exp(-ϕ2(x)), -Inf, Inf)[1]
lf2 = My.LFProblem(ϕ2)
@time X2 = My.HMC(lf2)
histogram(X2; norm=true, alpha=0.3, label="HMC LF sample", bin=100)
plot!(x -> exp(-ϕ2(x))/Z, -2.5, 2.5; label="exp(-ϕ2(x))/Z", lw=2)
plot!(; legend=:outertop)

# %%
@btime My.HMC($lf2);

# %% [markdown]
# https://github.com/genkuroki/public/blob/main/0019/HMC%20leapfrog.ipynb

# %%
ϕ3(x, param=nothing) = x^2*(x - 3)^2/2
Z = quadgk(x -> exp(-ϕ3(x)), -Inf, Inf)[1]
lf3 = My.LFProblem(ϕ3)
@time X3 = My.HMC(lf3; niters=10^6)
histogram(X3; norm=true, alpha=0.3, label="HMC LF sample", bin=100)
plot!(x -> exp(-ϕ3(x))/Z, -1.5, 4.5; label="exp(-ϕ3(x))/Z", lw=2)
plot!(; legend=:outertop)

# %%
@btime My.HMC($lf3; niters=10^6);

# %%
plot(X[1:10000]; label="", size=(800, 200), lw=0.5)

# %%
plot(X2[1:10000]; label="", size=(800, 200), lw=0.5)

# %%
plot(X3[1:10000]; label="", size=(800, 200), lw=0.5)

# %%
ϕ4(x, a) = a * (x^2 - 1)^2
a = [3, 4, 5, 6, 7, 8]
XX = Vector{Float64}[]
ZZ = Float64[]
PP = []
for i in eachindex(a)
    Z = quadgk(x -> exp(-ϕ4(x, a[i])), -Inf, Inf)[1]
    push!(ZZ, Z)
    lf = My.LFProblem(ϕ4)
    @time X = My.HMC(lf, a[i])
    flush(stdout)
    push!(XX, X)
    P = plot()
    histogram!(X; norm=true, alpha=0.3, label="HMC LF sample", bin=100, c=i)
    plot!(x -> exp(-ϕ4(x, a[i]))/Z, -2, 2; label="exp(-ϕ2(x))/Z", lw=2, c=i)
    plot!(; legend=false, xtick=-2:0.5:2)
    title!("ϕ(x) = a(x² - 1)²  for  a = $(a[i])", titlefontsize=9)
    push!(PP, P)
end
plot(PP...; size=(800, 450), layout=(3, 2))

# %%
QQ = []
for i in eachindex(a)
    Q = plot(XX[i][1:10000]; ylim=(-1.5, 1.5), label="", c=i, lw=0.5)
    title!("ϕ(x) = a(x² - 1)²  for  a = $(a[i])", titlefontsize=9)
    push!(QQ, Q)
end
plot(QQ...; size=(800, 900), layout=(length(a), 1))

# %%
# https://github.com/moruten/julia-code/blob/c0dfc2443d6b74256364e698b3edb37f98214ce7/Test-2021-9-3-no4.ipynb

#HMC method
using Plots
using StatsBase
using Random
using Printf
using Distributions
using BenchmarkTools
using QuadGK
using SpecialFunctions
using Zygote

# Generating Gaussian distribution (1-dim) by using HMC method

# # of configurations
n_config = 1_00_000

# size of a configuration
n_size = 1

# # of step of Leapfrog method
n_tau = 100

# time period of Hamilton equation
τ = 5.0

function S(x)
    sum(0.5*x[i]^2 for i=1:n_size)
end

function dSdx(x,k)
    return x
end

function H(x,p)
    K = sum(p[i]^2 for i=1:n_size)
    return S(x) + 0.5*K
end

#１変数更新のHMC法
function HMC_1dim(x)
        #初期化
        p_ini = rand(Normal(0,1),n_size)
        p_fin = zeros(n_size)
        x_ini = zeros(n_size)
        x_fin = zeros(n_size)

        #初期配位を保存
        for i=1:n_size x_ini[i] = x[i] end

        #初期配位のハミルトニアンを計算
        H_ini = H(x_ini,p_ini)

        #リープフロッグで時間発展
        leapfrog!(x,p_ini)

        #時間発展後の配位を保存
        for i=1:n_size x_fin[i] = x[i]; p_fin[i] = p_ini[i] end

        #時間発展後のハミルトニアンを計算
        H_fin = H(x_fin,p_fin)

        #Metropolis check
        r=rand()
        if r<exp(H_ini - H_fin)
            for i=1:n_size x[i] = x_fin[i] end
        else
            for i=1:n_size x[i] = x_ini[i] end
        end
end

#x,pをリープフロッグ法で更新する．
function leapfrog!(x,p)
    Δτ  = τ/n_tau
    p0  = 0.0
    p1  = 0.0
    x05 = 0.0
    x15 = 0.0
    for k=1:n_size
        #initial setting : τ0
        p0  = p[k]
        x05 = x[k] + p0 * 0.5*Δτ

        #(n-1)step : τ0 -> τn-1
        for i=1:n_tau-1
            p1  = p0 - dSdx(x05,k)*Δτ
            x15 = x05 + p1*Δτ

            p0  = p1
            x05 = x15
        end

        #final step : τn-1 -> τn
        p_end = p1 - dSdx(x15,k)*Δτ
        x_end = x15 + p_end*0.5*Δτ

        x[k]  = x_end
        p[k]  = p_end
    end
end

""" main """
#set of configurations
Ω = zeros(n_config)

#configuration
x = zeros(n_size)

@time for j=1:n_config
    HMC_1dim(x)
    Ω[j] = x[1]
end

Ω

gauss(x) = exp(-0.5*x^2)/sqrt(2*pi)
xr = range(-4,4,length=10001)

plot(Ω,st=:histogram,nbins=50,norm=:pdf,label="HMC's samples")
plot!(xr,gauss.(xr),label="Gaussian") |> display

# %%
function HMC_orig!()
    for j=1:n_config
        HMC_1dim(x)
        Ω[j] = x[1]
    end
end

""" main """
#set of configurations
Ω = zeros(n_config)

#configuration
x = zeros(n_size)

@btime HMC_orig!()

# %%
length(Ω), length(X)

# %%
1.3/0.057

# %%

max_stars_repo_path: 0018/HMC leapfrog.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0018/HMC leapfrog.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0018/HMC leapfrog.jl

max_forks_repo_name: genkuroki/public