text:
const sounds = (
    fire        = Sound(asset"fire.wav",      volume = 20),
    explosion   = Sound(asset"explosion.wav", volume = 20),
    gameover    = Sound(asset"game-over.wav", volume = 50)
)

@setup function()
    ship = Ship()
    asteroids = Asteroid[]
    lasers = Laser[]
    points = 0
    life = 100
    time = 0.0
    waiting = true
    wave = 0
    paused = false
end

onkeypress(@window) do event
    key = keyname(event) |> uppercase
    if key == "D"
        setrotation!(ship, 0.1)
    elseif key == "A"
        setrotation!(ship, -0.1)
    elseif key == "W"
        boosting!(ship, true)
    elseif key == "ESCAPE"
        global paused = !paused
    end
end

onmousedown(@canvas) do event
    if !paused
        push!(lasers, Laser(ship.pos..., ship.heading))
        stop(sounds.fire)
        play(sounds.fire)
    end
end

onmouseup(@canvas) do event
    boosting!(ship, false)
end

onevent("key-release-event", @window) do args...
    setrotation!(ship, 0)
end

function update(dt)
    w, h = size(@canvas)
    @layer begin
        background("#111")
        renderhud()
        if paused
            sethue("yellow")
            fontsize(28)
            text("Game paused", w/2, h/2, halign = :center, valign = :center)
            return;
        end

        if life <= 0
            return endgame()
        end

        for asteroid in asteroids
            render(asteroid)
            update!(asteroid)
            edges(asteroid)

            if hits(ship, asteroid)
                global life -= 1
            end
        end

        for i in length(lasers):-1:1
            render(lasers[i])
            update!(lasers[i])
            for j in length(asteroids):-1:1
                if lasers[i] --> asteroids[j]
                    new = breakup!(asteroids[j])
                    !isempty(new) && append!(asteroids, new)

                    stop(sounds.explosion)
                    play(sounds.explosion)

                    splice!(asteroids, j)
                    splice!(lasers, i)

                    global points += 1
                    break
                end
            end
            cleanlasers(i)
        end

        render(ship)
        turn!(ship)
        update!(ship)
        edges(ship)

        if isempty(asteroids) && !waiting
            global waiting = true
            global time = 0
        end

        newwave()

        global time += dt
    end
end

function cleanlasers(i::Int)
    if !isempty(lasers) && i < length(lasers)
        x, y = lasers[i].pos
        if x < 0 || x > @width() || y < 0 || y > @height()
            splice!(lasers, i)
        end
    end
end

function renderhud()
    sethue("white")
    fontsize(16)
    text("Points: $points", 30, 30)
    text("Wave: $wave", 120, 30)

    enemies = length(asteroids)
    text("Enemies remaining: $enemies", 200, 30)

    lp = life / 100
    hue = mapr(lp, 0, 1, 0, 90)

    sethue(HSL(hue, .65, .6))
    barh = 10
    rect(0, @height() - barh, @width() * lp, barh, :fill)
    text("Life: $(floor(Int, 100 * lp))%", 10, @height() - 20)
end

function endgame()
    @layer begin
        background("red")
        origin()
        fontsize(28)
        sethue("#111")
        text("Game Over", O, halign = :center, valign = :middle)
    end
    renderhud()
    play(sounds.gameover)
    loop!(false)
end

function showtimer()
    t = floor(Int, time)
    @layer begin
        origin()
        fontsize(28)
        sethue("white")
        text("Wave starting at $(5 - t)s", 0, -120, halign = :center, valign = :middle)
    end
end

function newwave()
    global waiting, time, wave, life;
    if waiting
        if time >= 5
            waiting = false
            wave += 1
            time = 0
            if life < 100
                life += 20
            end
            append!(asteroids, [Asteroid() for i in 1:3wave])
            @info "Wave started!"
        else
            showtimer()
        end
    end
end

onevent(:destroy, @window) do args...
    for sound in sounds
        Sounds.destroy!(sound)
    end
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
ENV["PYTHON_JL_RUNTIME_PYTHON"] = Sys.which("python")
dir(x) = fieldnames(typeof(x))
singlify(x::Float64) = Float32(x)

using Pkg
using PyCall
using OrderedCollections
using LinearAlgebra
using SpikeSynchrony
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
SNN.@load_units
using Distributed
import DataStructures
using JLD
using Metaheuristics
using Plots

const ngt_spikes
const opt_vec
const extremum
const extremum_param
const ngt_spikes
const fitness

PARALLEL = false
THREADS = false
if PARALLEL
    if nprocs() == 1
        addprocs(8)
    else
        @everywhere include("../utils.jl")
        @everywhere include("../current_search.jl")
        eval(macroexpand(quote
            @everywhere using SpikingNeuralNetworks
        end))
        eval(macroexpand(quote
            @everywhere using PyCall
        end))

    end
else
    include("../utils.jl")
end

###
const cell_type = "ADEXP"
const vecp = false
###

(vmgtv, vmgtt, ngt_spikes, ground_spikes, julia_version_vm) = get_data()
const simdur = last(vmgtt) * 1000

println("Ground Truth")
plot(vmgtt[:], vmgtv[:]) |> display

lower, upper = get_adexp_ranges()
ALLEN_DURATION = 2000 * ms
ALLEN_DELAY = 1000 * ms

function loss(E, ngt_spikes, ground_spikes)
    spikes = get_spikes(E)
    spikes = [s / 1000.0 for s in spikes]

    maxt = findmax(sort!(unique(vcat(spikes, ground_spikes))))[1]
    if size(spikes)[1] > 1
        t, S = SPIKE_distance_profile(spikes, ground_spikes; t0 = 0, tf = maxt)
        spkdistance = SpikeSynchrony.trapezoid_integral(t, S) / (t[end] - t[1]) # == SPIKE_distance(y1, y2)
    else
        spkdistance = 10.0
    end
    if length(spikes) > 1
        custom_raster2(spikes, ground_spikes)
        custom_raster(spikes, ground_spikes)
    end
    spkdistance *= spkdistance

    delta = abs(size(spikes)[1] - ngt_spikes)
    return spkdistance + delta

end

function loss(param)
    current = current_search(cell_type, param, ngt_spikes)
    if cell_type == "IZHI"
        pp = SNN.IZParameter(; a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(; N = 1, param = pp)
    end
    if cell_type == "ADEXP"
        adparam = SNN.ADEXParameter(;
            a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10],
        )

        E = SNN.AD(; N = 1, param = adparam)
    end

    SNN.monitor(E, [:v])
    SNN.monitor(E, [:fire])
    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms

    E.I = [current * nA]
    #@show(simdur)
    SNN.sim!(
        [E];
        dt = 1 * ms,
        delay = ALLEN_DELAY,
        stimulus_duration = ALLEN_DURATION,
        simulation_duration = simdur,
    )
    vecplot(E, :v) |> display
    error = loss(E, ngt_spikes, ground_spikes)
    @show(error)

    error
end

ɛ = 0.125
options = GA(
    populationSize = 10,
    ɛ = ɛ,
    selection = ranklinear(1.5),
    crossover = intermediate(1.0),
    mutation = uniform(1.0),
)
@time result = Evolutionary.optimize(
    loss,
    lower,
    upper,
    initd,
    options,
    Evolutionary.Options(
        iterations = 140,
        successive_f_tol = 25,
        show_trace = true,
        store_trace = true,
    ),#,parallelisation=:thread)
)

fitness = minimum(result)
println("GA: ɛ=$ɛ) => F: $(minimum(result))")
extremum_param = Evolutionary.minimizer(result)

using SignalAnalysis
opt_vec, opt_spikes = checkmodel(extremum_param, cell_type, ngt_spikes)
trace = result.trace
filename = string("../JLD/PopulationScatter_adexp.jld")
save(
    filename,
    "trace",
    trace,
    "opt_vec",
    opt_vec,
    "extremum_param",
    extremum_param,
    "vmgtt",
    vmgtt,
    "vmgtv",
    vmgtv,
    "ground_spikes",
    ground_spikes,
    "opt_spikes",
    opt_spikes,
)
evo_loss = [t.value for t in trace]
loss_evolution = []
for (i, l) in enumerate(evo_loss)
    if i > 1
        append!(loss_evolution, l)
    end
end
display(plot(loss_evolution))
println("Fittest model parameters")
println(extremum_param)

opt_vec = [i[1] for i in opt_vec]

s_a = signal(opt_vec, length(opt_vec) / last(vmgtt))
s_b = signal(vmgtv, length(vmgtt) / last(vmgtt))
p = plot(s_a)
p2 = plot!(p, s_b)
display(p2)
savefig(p2, "../imgs/aligned_VMs_adexp.png")

crp = custom_raster(opt_spikes, ground_spikes)
savefig(crp, "../imgs/aligned_VMs_adexp.png")
display(crp)

@show(result.minimizer)
@show(fitness)

ext: jl


text:
using CharFuncPricing, Test, Printf
using StatsBase

# bench test from AndersenLake (goal: verify accuracy of Cos & Flinn)
@testset "HestonALBench" begin
    forward = 100.0
    strikes = [100.0001, 101.0, 110.0, 200.0, 1000.0, 10000.0]
    τs = [0.0025, 0.1, 0.5, 2.0, 10.0, 30.0]
    v0s = [0.0001, 0.0025, 0.04, 0.25, 1.0]
    θs = [0.0001, 0.0025, 0.04, 0.25, 1.0]
    κs = [0.01, 0.1, 0.5, 2.0]
    σs = [0.0001, 0.1, 0.5, 1.0, 3.0]
    ρs = [-0.95, -0.5, -0.1, 0.0, 0.1, 0.5, 0.95]
    paramset = Vector(undef, 0)
    for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
    params = HestonParams(v0, κ, θ, ρ, σ)
    for strike in strikes
        push!(paramset, (τ, strike, params))
    end
end

    pricesA1 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                pricesA1[counter] = price
            end
        end
    end

    pricesA2 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, n=1000)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                pricesA2[counter] = price
            end
        end
    end

    pricesA3 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf,  n=2000)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                pricesA3[counter] = price
            end
        end
    end

    pricesA0GL = zeros(length(paramset));
    qgl = ModlobQuadrature(1e-10)
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, qgl)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward,τ, 1.0)
                pricesA0GL[counter] = price
            end
        end
    end

    pricesA0L = zeros(length(paramset));
    qde = DEQuadrature(1e-10)
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, qde)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward,τ, 1.0)
                pricesA0L[counter] = price
            end
        end
    end
    pricesC1 = Vector{Float64}(undef, 0);
    elapsed = @elapsed begin
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
        params = HestonParams(v0, κ, θ, ρ, σ)
        cf = DefaultCharFunc(params)
        refPricer = makeCosCharFuncPricer(cf, τ, 200 , 8)
        for strike in strikes
            price = priceEuropean(refPricer, false, strike, forward, 1.0)
            push!(pricesC1, price)
        end
    end
    end

    pricesC2 = Vector{Float64}(undef, 0);
    elapsed = @elapsed begin
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
        params = HestonParams(v0, κ, θ, ρ, σ)
        cf = DefaultCharFunc(params)
        refPricer = makeCosCharFuncPricer(cf, τ, 8)
        for strike in strikes
            price = priceEuropean(refPricer, false, strike, forward, 1.0)
            push!(pricesC2, price)
        end
    end
end

    pricesF1 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = CharFuncPricing.AdaptiveFlinnCharFuncPricer(cf, τ)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                pricesF1[counter] = price
            end
        end
    end

    pricesF2 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = CharFuncPricing.AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-9)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                pricesF2[counter] = price
            end
        end
    end

    pricesF3 = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = CharFuncPricing.AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-10)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                pricesF3[counter] = price
            end
        end
    end

    pricesF2B = zeros(length(paramset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = FlinnCharFuncPricer(cf, τ, tTol=1e-6, qTol=1e-10)
            for strike in strikes
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                pricesF2B[counter] = price
            end
        end
    end

    allp = [pricesC1, pricesC2, pricesF1, pricesF2, pricesF3, pricesF2B, pricesA1, pricesA2, pricesA0L]
    names = ["Cos","Cos-Adaptive","Flinn-Transf.","Flinn-Transf. (1e-9)", "Flinn-Transf. (1e-10)","Flinn (1e-10)","Andersen-Lake", "Andersen-Lake 1000", "Andersen-Lake GL"]
    for (p,name) in zip(allp,names)
        @printf("%s & %2e & %2e & %2e & %2e\n",name, maxad(p, pricesA3), rmsd(p, pricesA3),  maxad(p ./ pricesA3, ones(length(pricesA3))), rmsd(p ./ pricesA3, ones(length(pricesA3))))
    end

    allp = [pricesA1, pricesA2, pricesA0L]
    names = ["Andersen-Lake", "Andersen-Lake 1000", "Andersen-Lake GL"]
    for (p,name) in zip(allp,names)
        p3 = pricesA3
        @printf("%s & %.1e & %.1e & %.1e & %.1e\n",name, maxad(p, p3), rmsd(p, p3),  maxad(p ./ pricesA3, ones(length(p3))), rmsd(p ./ p3, ones(length(p3))))
    end
end
#
@testset "HestonALBenchForward" begin
    strike = 100.0
    forwards = [100.0,100.0001, 101.0, 110.0, 200.0, 1000.0, 10000.0]
    τs = [0.0025, 0.1, 0.5, 2.0, 10.0, 30.0]
    v0s = [0.0001, 0.0025, 0.04, 0.25, 1.0]
    θs = [0.0001, 0.0025, 0.04, 0.25, 1.0]
    κs = [0.01, 0.1, 0.5, 2.0]
    σs = [0.0001, 0.1, 0.5, 1.0, 3.0]
    ρs = [-0.95, -0.5, -0.1, 0.0, 0.1, 0.5, 0.95]
    fparamset = Vector(undef, 0)
    for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
    params = HestonParams(v0, κ, θ, ρ, σ)
    for forward in forwards
        push!(fparamset, (τ, forward, params))
    end
end

fpricesC1 = Vector{Float64}(undef, 0);
elapsed = @elapsed begin
    for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
    params = HestonParams(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    refPricer = makeCosCharFuncPricer(cf, τ, 200 , 8)
    for forward in forwards
        price = priceEuropean(refPricer, false, strike, forward, 1.0)
        push!(fpricesC1, price)
    end
end
end

fpricesC2 = Vector{Float64}(undef, 0);
elapsed = @elapsed begin
    for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
    params = HestonParams(v0, κ, θ, ρ, σ)
    cf = DefaultCharFunc(params)
    refPricer = makeCosCharFuncPricer(cf, τ, 8)
    for forward in forwards
        price = priceEuropean(refPricer, false, strike, forward, 1.0)
        push!(fpricesC2, price)
    end
end
end
    fpricesF1 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = AdaptiveFlinnCharFuncPricer(cf, τ)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                fpricesF1[counter] = price
            end
        end
    end

    fpricesF2 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = CharFuncPricing.AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-9)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                fpricesF2[counter] = price
            end
        end
    end

    fpricesF3 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = CharFuncPricing.AdaptiveFlinnCharFuncPricer(cf, τ, qTol=1e-10)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                fpricesF3[counter] = price
            end
        end
    end
    fpricesF2B = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = FlinnCharFuncPricer(cf, τ, tTol=1e-10, qTol=1e-10)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, 1.0)
                fpricesF2B[counter] = price
            end
        end
    end

    quadrature = TanhSinhQuadrature(200,eps())
    fpricesA1 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf,quadrature)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                fpricesA1[counter] = price
            end
        end
    end

    quadrature = TanhSinhQuadrature(1000,eps())
fpricesA2 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, quadrature)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                fpricesA2[counter] = price
            end
        end
    end

    quadrature = TanhSinhQuadrature(2000,eps())
    fpricesA3 = zeros(length(fparamset));
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, quadrature)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                fpricesA3[counter] = price
            end
        end
    end

    quadrature = ModlobQuadrature(1e-10)
    fpricesA0L = zeros(length(fparamset))
    elapsed = @elapsed begin
        local counter = 0
        for τ in τs, v0 in v0s, θ in θs, κ in κs, σ in σs, ρ in ρs
            params = HestonParams{Float64}(v0, κ, θ, ρ, σ)
            cf = DefaultCharFunc(params)
            pricer = ALCharFuncPricer(cf, quadrature)
            for forward in forwards
                counter += 1
                price = priceEuropean(pricer, false, strike, forward, τ, 1.0)
                fpricesA0L[counter] = price
            end
        end
    end

    allp = [fpricesC1, fpricesC2, fpricesF1, fpricesF2, fpricesF2B, fpricesA1, fpricesA2, fpricesA0L]
    names = ["Cos","Cos-Adaptive","Flinn-Transf.","Flinn-Transf. (1e-10)","Flinn (1e-10)","Andersen-Lake", "Andersen-Lake 1000", "Andersen-Lake GL"]
    for (p,name) in zip(allpf,names)
        @printf("%s & %2e & %2e & %2e & %2e\n",name, maxad(p, fpricesA3), rmsd(p, fpricesA3),  maxad(p ./ fpricesA3, ones(length(fpricesA3))), rmsd(p ./ fpricesA3, ones(length(fpricesA3))))
    end

    allp = [pricesC1, pricesC2, pricesF1, pricesF3, pricesF2B, pricesA1, pricesA2, pricesA0L]
    allpf = [fpricesC1, fpricesC2, fpricesF1, fpricesF3, fpricesF2B, fpricesA1, fpricesA2, fpricesA0L]
    names = ["Cos","Cos-Adaptive","Flinn-Transf.","Flinn-Transf. (1e-10)","Flinn (1e-10)","Andersen-Lake", "Andersen-Lake 1000", "Andersen-Lake GL"]
    for (pk, pf, name) in zip(allp, allpf, names)
        p3 = vcat(pricesA3,fpricesA3)
        ind = findall(z -> z > 1e-90, p3)
        p = vcat(pk, pf)
        @printf("%s & %.1e & %.1e & %.1e & %.1e\n",name, maxad(p[ind], p3[ind]), rmsd(p[ind], p3[ind]),  maxad(p[ind] ./ p3[ind], ones(length(p3[ind]))), rmsd(p[ind] ./ p3[ind], ones(length(p3[ind]))))
    end

    allpf = [fpricesA1, fpricesA2]
    names = ["Andersen-Lake", "Andersen-Lake 1000"]
    for (pf, name) in zip(allpf, names)
        p3 = fpricesA3
        ind = findall(z -> z > 1e-90, p3)
        p =  pf
        @printf("%s & %.1e & %.1e & %.1e & %.1e\n",name, maxad(p[ind], p3[ind]), rmsd(p[ind], p3[ind]),  maxad(p[ind] ./ p3[ind], ones(length(p3[ind]))), rmsd(p[ind] ./ p3[ind], ones(length(p3[ind]))))
    end
end


text:
# Use Chemfiles_jll & the artifcat system on julia >= 1.3
if VERSION >= v"1.3.0"
    Sys.exit()
end

include("build_Chemfiles.v0.9.3.jl")


text:
importall Bukdu
importall Bukdu.Octo
import Bukdu.Octo: |>

type User
    name::String
    username::String
end

using Base.Test

model = default(User)
params = Dict("user[name]"=>"José", "user[username]"=>"josevalim")

@test_throws MethodError validates(model, params)

function validates(model::User, params)
    model |>
    cast(params, [:name, :username]) |>
    validate_length(:username, min= 1, max= 20)
end

@test isa(validates(model, params), Changeset)


text:
using Test
using ClimateMachine
import ClimateMachine.BalanceLaws
import ClimateMachine.BalanceLaws: boundary_state!
using ClimateMachine.DGMethods: ESDGModel, init_ode_state
using ClimateMachine.DGMethods.NumericalFluxes:
    numerical_volume_flux_first_order!
using ClimateMachine.Mesh.Topologies: BrickTopology
using ClimateMachine.Mesh.Grids: DiscontinuousSpectralElementGrid
using ClimateMachine.VariableTemplates: varsindex
using StaticArrays: MArray, @SVector
using KernelAbstractions: wait
using Random
using DoubleFloats
using MPI
using ClimateMachine.MPIStateArrays

Random.seed!(7)

include("DryAtmos.jl")

boundary_state!(::Nothing, _...) = nothing

struct TestProblem <: AbstractDryAtmosProblem end

# Random initialization function
function init_state_prognostic!(
    ::DryAtmosModel,
    ::TestProblem,
    state_prognostic,
    state_auxiliary,
    _...,
) where {dim}
    FT = eltype(state_prognostic)
    ρ = state_prognostic.ρ = rand(FT) + 1
    ρu = state_prognostic.ρu = 2 * (@SVector rand(FT, 3)) .- 1
    p = rand(FT) + 1
    Φ = state_auxiliary.Φ
    state_prognostic.ρe = totalenergy(ρ, ρu, p, Φ)
end

function check_operators(FT, dim, mpicomm, N, ArrayType)
    # Create a warped mesh so the metrics are not constant
    Ne = (8, 9, 10)
    brickrange = (
        range(FT(-1); length = Ne[1] + 1, stop = 1),
        range(FT(-1); length = Ne[2] + 1, stop = 1),
        range(FT(-1); length = Ne[3] + 1, stop = 1),
    )
    topl = BrickTopology(
        mpicomm,
        ntuple(k -> brickrange[k], dim);
        periodicity = ntuple(k -> true, dim),
    )
    warpfun =
        (x1, x2, x3) -> begin
            α = (4 / π) * (1 - x1^2) * (1 - x2^2) * (1 - x3^2)
            # Rotate by α with x1 and x2
            x1, x2 = cos(α) * x1 - sin(α) * x2, sin(α) * x1 + cos(α) * x2
            # Rotate by α with x1 and x3
            if dim == 3
                x1, x3 = cos(α) * x1 - sin(α) * x3, sin(α) * x1 + cos(α) * x3
            end
            return (x1, x2, x3)
        end
    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
        meshwarp = warpfun,
    )

    # Orientation does not matter since we will be setting the geopotential to a
    # random field
    model = DryAtmosModel{dim}(FlatOrientation(), TestProblem())

    ##################################################################
    # check that the volume terms lead to only surface contributions #
    ##################################################################
    # Create the ES model
    esdg = ESDGModel(
        model,
        grid;
        volume_numerical_flux_first_order = EntropyConservative(),
        surface_numerical_flux_first_order = nothing,
    )

    # Make the Geopotential random
    esdg.state_auxiliary .= ArrayType(2rand(FT, size(esdg.state_auxiliary)))
    start_exchange = MPIStateArrays.begin_ghost_exchange!(esdg.state_auxiliary)
    end_exchange = MPIStateArrays.end_ghost_exchange!(
        esdg.state_auxiliary,
        dependencies = start_exchange,
    )
    wait(end_exchange)

    # Create a random state
    state_prognostic = init_ode_state(esdg; init_on_cpu = true)

    # Storage for the tendency
    volume_tendency = similar(state_prognostic)

    # Compute the tendency function
    esdg(volume_tendency, state_prognostic, nothing, 0)

    # Check that the volume terms only lead to surface integrals of
    #    ∑_{j} n_j ψ_j
    # where Ψ_j = β^T f_j - ζ_j = ρu_j
    Np, K = prod(N .+ 1), length(esdg.grid.topology.realelems)

    # Get the mass matrix on the host
    _M = ClimateMachine.Grids._M
    M = Array(grid.vgeo[:, _M:_M, 1:K])

    # Get the state, tendency, and aux on the host
    Q = Array(state_prognostic.data[:, :, 1:K])
    dQ = Array(volume_tendency.data[:, :, 1:K])
    A = Array(esdg.state_auxiliary.data[:, :, 1:K])

    # Compute the entropy variables
    β = similar(Q, Np, number_states(model, Entropy()), K)
    @views for e in 1:K
        for i in 1:Np
            state_to_entropy_variables!(
                model,
                β[i, :, e],
                Q[i, :, e],
                A[i, :, e],
            )
        end
    end

    # Get the unit normals and surface mass matrix
    sgeo = Array(grid.sgeo)
    n1 = sgeo[ClimateMachine.Grids._n1, :, :, 1:K]
    n2 = sgeo[ClimateMachine.Grids._n2, :, :, 1:K]
    n3 = sgeo[ClimateMachine.Grids._n3, :, :, 1:K]
    sM = sgeo[ClimateMachine.Grids._sM, :, :, 1:K]

    num_state = number_states(model, Prognostic())
    volume = sum(β[:, 1:num_state, :] .* M .* dQ, dims = (1, 2))[:]

    ###########################################
    # check that the volume and surface match #
    ###########################################
    esdg = ESDGModel(
        model,
        grid;
        state_auxiliary = esdg.state_auxiliary,
        volume_numerical_flux_first_order = nothing,
        surface_numerical_flux_first_order = EntropyConservative(),
    )

    surface_tendency = similar(state_prognostic)

    # Compute the tendency function
    esdg(surface_tendency, state_prognostic, nothing, 0)

    # Surface integral should be equal and opposite to the volume integral
    dQ = Array(surface_tendency.data[:, :, 1:K])
    volume_integral = MPI.Allreduce(sum(volume), +, mpicomm)
    surface_integral =
        MPI.Allreduce(sum(β[:, 1:num_state, :] .* M .* dQ), +, mpicomm)
    @test volume_integral ≈ -surface_integral

    ########################################################
    # check that the full tendency is entropy conservative #
    ########################################################
    esdg = ESDGModel(
        model,
        grid;
        state_auxiliary = esdg.state_auxiliary,
        volume_numerical_flux_first_order = EntropyConservative(),
        surface_numerical_flux_first_order = EntropyConservative(),
    )

    tendency = similar(state_prognostic)

    # Compute the tendency function
    esdg(tendency, state_prognostic, nothing, 0)

    # Check for entropy conservation
    dQ = Array(tendency.data[:, :, 1:K])
    integral = MPI.Allreduce(sum(β[:, 1:num_state, :] .* M .* dQ), +, mpicomm)
    @test isapprox(integral, 0, atol = sqrt(eps(sum(volume))))
end

let
    model = DryAtmosModel{3}(FlatOrientation(), TestProblem())
    num_state = number_states(model, Prognostic())
    num_aux = number_states(model, Auxiliary())
    num_entropy = number_states(model, Entropy())

    @testset "state to entropy variable transforms" begin
        for FT in (Float32, Float64)
            state_in =
                [1, 2, 2, 2, 1] .* rand(FT, num_state) + [3, -1, -1, -1, 100]
            aux_in = rand(FT, num_aux)
            state_out = similar(state_in)
            aux_out = similar(aux_in)
            entropy = similar(state_in, num_entropy)

            state_to_entropy_variables!(model, entropy, state_in, aux_in)
            entropy_variables_to_state!(model, state_out, aux_out, entropy)

            @test all(state_in .≈ state_out)
            #@test all(aux_in .≈ aux_out)
        end
    end

    @testset "test numerical flux for Tadmor shuffle" begin
        for FT in (Float32, Float64)
            # Create some random states
            state_1 =
                [1, 2, 2, 2, 1] .* rand(FT, num_state) + [3, -1, -1, -1, 100]
            aux_1 = 0 * rand(FT, num_aux)

            state_2 =
                [1, 2, 2, 2, 1] .* rand(FT, num_state) + [3, -1, -1, -1, 100]
            aux_2 = 0 * rand(FT, num_aux)

            # Get the entropy variables for the two states
            entropy_1 = similar(state_1, num_entropy)
            state_to_entropy_variables!(model, entropy_1, state_1, aux_1)

            entropy_2 = similar(state_1, num_entropy)
            state_to_entropy_variables!(model, entropy_2, state_2, aux_2)

            # Get the values of Ψ_j = β^T f_j - ζ_j = ρu_j where β is the
            # entropy variables, f_j is the conservative flux, and ζ_j is the
            # entropy flux. For conservation laws this is the entropy potential.
            Ψ_1 = Vars{vars_state(model, Prognostic(), FT)}(state_1).ρu
            Ψ_2 = Vars{vars_state(model, Prognostic(), FT)}(state_2).ρu

            # Evaluate the flux with both orders of the two states
            H_12 = fill!(MArray{Tuple{3, num_state}, FT}(undef), -zero(FT))
            numerical_volume_flux_first_order!(
                EntropyConservative(),
                model,
                H_12,
                state_1,
                aux_1,
                state_2,
                aux_2,
            )

            H_21 = fill!(MArray{Tuple{3, num_state}, FT}(undef), -zero(FT))
            numerical_volume_flux_first_order!(
                EntropyConservative(),
                model,
                H_21,
                state_2,
                aux_2,
                state_1,
                aux_1,
            )

            # Check that we satisfy the Tadmor shuffle
            @test all(
                H_12 * entropy_1[1:num_state] - H_21 * entropy_2[1:num_state] .≈
                Ψ_1 - Ψ_2,
            )
        end
    end

    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()
    mpicomm = MPI.COMM_WORLD
    polynomialorders = ((2, 2, 1), (1, 1, 2))
    test_types = (Float32, Float64)
    for FT in test_types
        for dim in (3,)
            for polynomialorder in polynomialorders
                dim = 3
                @testset "check ESDGMethods relations for dim = $dim, FT = $FT, and polynomial order $polynomialorder" begin
                    check_operators(
                        FT,
                        dim,
                        mpicomm,
                        polynomialorder,
                        ArrayType,
                    )
                end
            end
        end
    end
end


text:
# OedipusRex.jl
# - Polaron optical absorption
#     A WORK IN PROGRESS!

#using QuadGK
"""
    ReX(v=20,w=20)

Codes to implement 'Optical Absorption of Polarons in the
Feynman-Hellwarth-Iddings-Platzman Approximation',
https://doi.org/10.1103/PhysRevB.5.2367 ;

In particular we are trying to calculate Re[χ] , the real part of the polaron optical absorption.
# (13) - defn of Re[χ]

"""
function ReX(v=20,w=20)
    R=(v^2-w^2)/(w^2*v)

    Reχintegrand(u,Ω) =  (1-cos(Ω*u)*exp(im*u))/(R*(1-exp(im*v*u))-im*u)^(3/2)

    Ω=1
    #Reχ = quadgk(u->Reχintegrand(u,Ω), 0.0, Inf)
# OK, problematic as this is a complex (multi valued integration!)

    [  Reχintegrand(u,Ω) for u=0:20 ]
    # Oh yikes, and it explodes as u->0 takes it to ->0 in the denominator
end


text:
"""
TODO
"""
struct HeptaOrder{T} <: FieldVector{6, T}
  front::Union{T,Void}
  back::Union{T,Void}
  left::Union{T,Void}
  right::Union{T,Void}
  top::Union{T,Void}
  bottom::Union{T,Void}
end

"""
TODO
"""
type Block
  typ::Int32
  pos::Vec3f
  flags::UInt32
  active::Bool
  visible::Bool
  surrounded::Bool
  sides::Array{UInt32,1}
  friends::HeptaOrder{Block}
end

const BlockOrder = HeptaOrder{Block}

"""
TODO
"""
BlockOrder() = HeptaOrder{Block}(nothing,nothing,nothing,nothing,nothing,nothing)

"""
TODO
"""
Block(pos=Vec3f,typ=0) = Block(typ,pos,0,true,true,false,resetSides(),BlockOrder()) #,zeros(Mat4x4f),zeros(Mat4x4f))

"""
TODO
"""
type Chunk
  active::Bool
  len::UInt32
  childs::Array{Block,3}
  filtered::Array{Block,1}
  count::UInt32
  fileredCount::UInt32
end

"""
TODO
"""
function Chunk(len::Integer)
  this = Chunk(true,len,Array{Block,3}(len,len,len),Block[],0,0)
  for x=1:len; for y=1:len; for z=1:len; this.childs[x,y,z] = Block(Vec3f(x,y,z)); end; end; end
  this
end

"""
TODO
"""
function clean(this::Union{Void,Chunk})
  if this == nothing return end
  this.childs = Array{Block,3}(0,0,0)
  this.filtered = Block[]
  gc() # force garbage collection, free memory
end

"""
TODO
"""
isType(this::Block, typ) = this.typ == typ

"""
TODO
"""
isValid(this::Chunk) = this.count > 0 && this.fileredCount > 0

"""
TODO
"""
isSeen(this::Block) = isActive(this) && isVisible(this) && this.typ > 0

MAX_SIDES = 6

LEFT_SIDE = 1 #0x1
RIGHT_SIDE = 2 #0x2
TOP_SIDE = 3 #0x4
BOTTOM_SIDE = 4 #0x8
FRONT_SIDE = 5 #0x10
BACK_SIDE = 6 #0x20

"""
TODO
"""
resetSides() = fill(UInt32(1),MAX_SIDES)

"""
TODO
"""
resetSides(this::Block) = this.sides=resetSides()

"""
TODO
"""
function hideUnseen(this::Chunk)
  a = this.childs
  len = this.len

  for x=1:len; for y=1:len; for z=1:len;

    b=a[x,y,z]

    if !isSeen(b) continue end # skip invisible
    #if  x==1 || y==1 || z==1 || x==len || y==len || z==len continue end  # skip border

    blocked=0

    if x>1
      r=b.sides[LEFT_SIDE]=!isSeen(a[x-1,y,z])
      if !r blocked+=1 end
    end

    if x<len
      r=b.sides[RIGHT_SIDE]=!isSeen(a[x+1,y,z])
      if !r blocked+=1 end
    end

    if y>1
      r=b.sides[BOTTOM_SIDE]=!isSeen(a[x,y-1,z])
      if !r blocked+=1 end
    end

    if y<len
      r=b.sides[TOP_SIDE]=!isSeen(a[x,y+1,z])
      if !r blocked+=1 end
    end

    if z>1
      r=b.sides[BACK_SIDE]=!isSeen(a[x,y,z-1])
      if !r blocked+=1 end
    end

    if z<len
      r=b.sides[FRONT_SIDE]=!isSeen(a[x,y,z+1])
      if !r blocked+=1 end
    end

    setSurrounded(b,blocked >= 6)

  end; end; end
end

"""
TODO
"""
function setFlag(this::Block, flag::Unsigned, add::Bool)
    if add this.flags |= flag
    else this.flags ⊻= flag #⊻ = xor and is not shown on windows...
    end
end

"""
TODO
"""
isActive(this::Block) = this.active

"""
TODO
"""
isVisible(this::Block) = this.visible

"""
TODO
"""
isSurrounded(this::Block) = this.surrounded

"""
TODO
"""
isValid(this::Block) = isActive(this) && isVisible(this) && !isSurrounded(this) && this.typ > 0

"""
TODO
"""
setActive(this::Block, active::Bool) = this.active = active

"""
TODO
"""
setVisible(this::Block, visible::Bool) = this.visible = visible

"""
TODO
"""
setSurrounded(this::Block, surrounded::Bool) = this.surrounded = surrounded

"""
TODO
"""
hideType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setVisible(b,false); end; end

"""
TODO
"""
removeType(this::Chunk, typ::Integer) = for b in this.childs; if b.typ == typ; setActive(b,false); end; end

"""
TODO
"""
showAll(this::Chunk) = for b in this.childs
  setVisible(b,true)
  setSurrounded(b,false)
  resetSides(b)
end

"""
TODO
"""
checkInFrustum(this::Chunk, fstm::Frustum) = for b in this.childs setVisible(b,checkSphere(fstm, b.pos, 1.5) != 0) end

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

"""
TODO
"""
setFilteredChilds(this::Chunk, r::Array{Block,1}) = begin this.filtered = r; this.fileredCount = length(r); r end

"""
TODO
"""
getFilteredChilds(this::Chunk) = this.filtered

"""
TODO
"""
getActiveChilds(this::Chunk) = filter(b->isActive(b),this.childs)

"""
TODO
"""
getVisibleChilds(this::Chunk) = filter(b->isVisible(b),this.childs)

"""
TODO
"""
getValidChilds(this::Chunk) = filter(b->isValid(b),this.childs)

"""
TODO
"""
function getData(this::Block)

  sides=0
  if this.sides[LEFT_SIDE] > 0 sides |= 0x1 end
  if this.sides[RIGHT_SIDE] > 0 sides |= 0x2 end
  if this.sides[TOP_SIDE] > 0 sides |= 0x4 end
  if this.sides[BOTTOM_SIDE] > 0 sides |= 0x8 end
  if this.sides[FRONT_SIDE] > 0 sides |= 0x10 end
  if this.sides[BACK_SIDE] > 0 sides |= 0x20 end

  # get visible sides
  #i=0; sides=0
  #for side in this.sides
  #  if side > 0 sides |= (0x1 << i) end
  #  i+=1
  #end

  #count = MAX_SIDES*5
  #a=fill(0f0,count)
  #i=1

  #for side in this.sides
  #  for p in this.pos a[i] = p; i+=1 end
  #  a[i] = this.typ; i +=1
  #  a[i] = side; i +=1
  #end

  #SVector(Float32[
  #  this.pos...,this.typ,0,sides
  #  this.pos...,this.typ,1,sides
  #  this.pos...,this.typ,2,sides
  #  this.pos...,this.typ,3,sides
  #  this.pos...,this.typ,4,sides
  #  this.pos...,this.typ,5,sides
  #]...)

  SVector(Float32[this.pos...,this.typ,sides]...)
end

"""
TODO
"""
getData(this::Chunk) = isValid(this)?vec((b->getData(b)).(getFilteredChilds(this))):Float32[]

"""
TODO
"""
function update(this::Chunk)
  this.count = length(this.childs)
  setFilteredChilds(this, getValidChilds(this))

  #fill copies references
  #refblocks = Ptr{Void}[]
  #for b in this.childs push!(refblocks, pointer_from_objref(b)) end
end

"""
TODO
"""
function createSingle(this::Chunk)
  this.childs[1,1,1] = Block(Vec3f(0, 0, -10),1)
end

"""
TODO
"""
function createExample(this::Chunk)
  const DIST = Vec3f(2,2,2); #Vec3f(3,5,3) #r = 1f0/30 #0.005f0
  const START = Vec3f(-(this.len*DIST.x) / 2.0f0, -(this.len*DIST.y) / 2.0f0, (this.len*DIST.z) / 2.0f0)

  translate(x,y,z) = Vec3f(START.x+x*DIST.x, START.y+y*DIST.y, START.z-z*DIST.z)

  x=1; y=1; z=1;

  for i=1:(this.len^3)
    this.childs[i] = b = Block(translate(x-1,y-1,z-1),rand(1:16))
    #model = Mat4x4f(translation(Array(b.pos)))
    #push!(refblocks, pointer_from_objref(b))

    x += 1; if x > this.len
      y += 1; x=1;
      if y > this.len z += 1; y=1;
        if z > this.len break; end
      end
    end
  end
end

"""
TODO
"""
function createLandscape(this::Chunk)
  #Texture* heightTexture = m_pRenderer->GetTexture(m_pChunkManager->GetHeightMapTexture());

  rgbImage=Images.load("heightmap.jpg")
  (imgwidth, imgheight) = size(rgbImage)

  grayImage = Gray.(rgbImage)
  imgMatrix = reinterpret.(channelview(grayImage)) #convert(Array{Int32},raw(grayImage))

  w = imgwidth/this.len
  h = imgheight/this.len

  const DIST = Vec3f(2,2,2)
  const START = Vec3f(-(this.len*DIST.x) / 2.0f0, -(this.len*DIST.y) / 2.0f0, (this.len*DIST.z) / 2.0f0)

  translate(x,y,z) = Vec3f(START.x+x*DIST.x, START.y+y*DIST.y, START.z-z*DIST.z)

  for x=1:this.len
    for z=1:this.len
      # Use the height map texture to get the height value of x, z
      height = (imgMatrix[UInt32(trunc(z*h)), UInt32(trunc(x*w))] / 0xFF)* 1.25 * this.len
      if height >= this.len height=this.len end

      level_air = height * 0.95
      level_grass = height * 0.9
      level_dirt = height * 0.8
      level_stonebricks = height * 0.7
      level_stone = height * 0.4
      level_lava = height * 0

      for y=1:UInt32(round(height))
        #id = y*this.len^2+z*this.len+x
        #c = this.childs[id]
        #c.active = true
        b = Block(translate(x-1,y-1,z-1))
        #if y >= level_air b.typ = 0 # air or nothing
        if y >= level_grass  b.typ = 2 #grass
        elseif y >= level_dirt b.typ = 1 #dirt
        elseif y >= level_stonebricks b.typ = 5 #stonebricks
        elseif y >= level_stone b.typ = 4 #stone
        elseif y >= level_lava b.typ = 15 #lava
        end
        if b.typ > 0 this.childs[x,y,z]=b end #push!(this.childs,b) end
      end
    end
  end

  h=1 #UInt32(trunc(this.len/2))
  this.childs[h,h,h] = Block(Vec3f(0,0,0),7)
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
using Colors
using ShaderAbstractions: InstancedProgram, Program
using AbstractPlotting: Key, plotkey
using Colors: N0f8

AbstractPlotting.plotkey(::Nothing) = :scatter

function lift_convert(key, value, plot)
    val = lift(value) do value
        return wgl_convert(value, Key{key}(), Key{plotkey(plot)}())
    end
    if key == :colormap && val[] isa AbstractArray
        return ShaderAbstractions.Sampler(val)
    else
        return val
    end
end

function Base.pairs(mesh::GeometryBasics.Mesh)
    return (kv for kv in GeometryBasics.attributes(mesh))
end

function GeometryBasics.faces(x::VertexArray)
    return GeometryBasics.faces(getfield(x, :data))
end

tlength(T) = length(T)
tlength(::Type{<:Real}) = 1

serialize_three(val::Number) = val
serialize_three(val::Vec2f0) = convert(Vector{Float32}, val)
serialize_three(val::Vec3f0) = convert(Vector{Float32}, val)
serialize_three(val::Vec4f0) = convert(Vector{Float32}, val)
serialize_three(val::Quaternion) = convert(Vector{Float32}, collect(val.data))
serialize_three(val::RGB) = Float32[red(val), green(val), blue(val)]
serialize_three(val::RGBA) = Float32[red(val), green(val), blue(val), alpha(val)]
serialize_three(val::Mat4f0) = vec(val)
serialize_three(val::Mat3) = vec(val)

function serialize_three(observable::Observable)
    return Dict(:type => "Observable", :id => observable.id,
                :value => serialize_three(observable[]))
end

function serialize_three(array::AbstractArray)
    return serialize_three(flatten_buffer(array))
end

function serialize_three(array::Buffer)
    return serialize_three(flatten_buffer(array))
end

function serialize_three(array::AbstractArray{UInt8})
    return Dict(:type => "Uint8Array", :data => array)
end

function serialize_three(array::AbstractArray{Int32})
    return Dict(:type => "Int32Array", :data => array)
end

function serialize_three(array::AbstractArray{UInt32})
    return Dict(:type => "Uint32Array", :data => array)
end

function serialize_three(array::AbstractArray{Float32})
    return Dict(:type => "Float32Array", :data => array)
end

function serialize_three(array::AbstractArray{Float16})
    return Dict(:type => "Float32Array", :data => array)
end

function serialize_three(array::AbstractArray{Float64})
    return Dict(:type => "Float64Array", :data => array)
end

function serialize_three(color::Sampler{T,N}) where {T,N}
    tex = Dict(:type => "Sampler", :data => serialize_three(color.data),
               :size => [size(color.data)...], :three_format => three_format(T),
               :three_type => three_type(eltype(T)),
               :minFilter => three_filter(color.minfilter),
               :magFilter => three_filter(color.magfilter),
               :wrapS => three_repeat(color.repeat[1]), :anisotropy => color.anisotropic)
    if N > 1
        tex[:wrapT] = three_repeat(color.repeat[2])
    end
    if N > 2
        tex[:wrapR] = three_repeat(color.repeat[3])
    end
    return tex
end

function serialize_uniforms(dict::Dict)
    result = Dict{Symbol,Any}()
    for (k, v) in dict
        result[k] = serialize_three(to_value(v))
    end
    return result
end

three_format(::Type{<:Real}) = "RedFormat"
three_format(::Type{<:RGB}) = "RGBFormat"
three_format(::Type{<:RGBA}) = "RGBAFormat"

three_type(::Type{Float16}) = "FloatType"
three_type(::Type{Float32}) = "FloatType"
three_type(::Type{N0f8}) = "UnsignedByteType"

function three_filter(sym)
    sym == :linear && return "LinearFilter"
    return sym == :nearest && return "NearestFilter"
end

function three_repeat(s::Symbol)
    s == :clamp_to_edge && return "ClampToEdgeWrapping"
    s == :mirrored_repeat && return "MirroredRepeatWrapping"
    return s == :repeat && return "RepeatWrapping"
end

"""
    flatten_buffer(array::AbstractArray)
Flattens `array` array to be a 1D Vector of Float32 / UInt8.
If presented with AbstractArray{<: Colorant/Tuple/SVector}, it will flatten those
to their element type.
"""
function flatten_buffer(array::AbstractArray{<: Number})
    return array
end

function flatten_buffer(array::Buffer)
    return flatten_buffer(getfield(array, :data))
end

function flatten_buffer(array::AbstractArray{T}) where {T<:N0f8}
    return reinterpret(UInt8, array)
end

function flatten_buffer(array::AbstractArray{T}) where {T}
    return flatten_buffer(reinterpret(eltype(T), array))
end

lasset(paths...) = read(joinpath(@__DIR__, "..", "assets", paths...), String)

isscalar(x::StaticArrays.StaticArray) = true
isscalar(x::AbstractArray) = false
isscalar(x::Observable) = isscalar(x[])
isscalar(x) = true

function ShaderAbstractions.type_string(::ShaderAbstractions.AbstractContext,
                                        ::Type{<:AbstractPlotting.Quaternion})
    return "vec4"
end

function ShaderAbstractions.convert_uniform(::ShaderAbstractions.AbstractContext,
                                            t::Quaternion)
    return convert(Quaternion, t)
end

function wgl_convert(value, key1, key2)
    val = AbstractPlotting.convert_attribute(value, key1, key2)
    return if val isa AbstractArray{<:Float64}
        return AbstractPlotting.el32convert(val)
    else
        return val
    end
end

function wgl_convert(value::AbstractMatrix, ::key"colormap", key2)
    return ShaderAbstractions.Sampler(value)
end

function serialize_buffer_attribute(buffer::AbstractVector{T}) where {T}
    return Dict(:flat => serialize_three(buffer), :type_length => tlength(T))
end

function serialize_named_buffer(buffer)
    return Dict(map(pairs(buffer)) do (name, buff)
                    return name => serialize_buffer_attribute(buff)
                end)
end

function register_geometry_updates(update_buffer::Observable, named_buffers)
    for (name, buffer) in pairs(named_buffers)
        if buffer isa Buffer
            on(ShaderAbstractions.updater(buffer).update) do (f, args)
                # update to replace the whole buffer!
                if f === (setindex!) && args[1] isa AbstractArray && args[2] isa Colon
                    new_array = args[1]
                    flat = flatten_buffer(new_array)
                    update_buffer[] = [name, serialize_three(flat), length(new_array)]
                end
                return
            end
        end
    end
    return update_buffer
end

function register_geometry_updates(update_buffer::Observable, program::Program)
    return register_geometry_updates(update_buffer, program.vertexarray)
end

function register_geometry_updates(update_buffer::Observable, program::InstancedProgram)
    return register_geometry_updates(update_buffer, program.per_instance)
end

function uniform_updater(uniforms::Dict)
    updater = Observable(Any[:none, []])
    for (name, value) in uniforms
        if value isa Sampler
            on(ShaderAbstractions.updater(value).update) do (f, args)
                if args[2] isa Colon && f == setindex!
                    updater[] = [name, serialize_three(args[1])]
                end
                return
            end
        else
            value isa Observable || continue
            on(value) do value
                updater[] = [name, serialize_three(value)]
                return
            end
        end
    end
    return updater
end

function serialize_three(ip::InstancedProgram)
    program = serialize_three(ip.program)
    program[:instance_attributes] = serialize_named_buffer(ip.per_instance)
    register_geometry_updates(program[:attribute_updater], ip)
    return program
end

function serialize_three(program::Program)
    indices = GeometryBasics.faces(program.vertexarray)
    indices = reinterpret(UInt32, indices)
    uniforms = serialize_uniforms(program.uniforms)
    attribute_updater = Observable(["", [], 0])
    register_geometry_updates(attribute_updater, program)
    return Dict(:vertexarrays => serialize_named_buffer(program.vertexarray),
                :faces => indices, :uniforms => uniforms,
                :vertex_source => program.vertex_source,
                :fragment_source => program.fragment_source,
                :uniform_updater => uniform_updater(program.uniforms),
                :attribute_updater => attribute_updater)
end

function serialize_scene(scene::Scene, serialized_scenes=[])
    hexcolor(c) = "#" * hex(Colors.color(to_color(c)))
    pixel_area = lift(area -> [minimum(area)..., widths(area)...], pixelarea(scene))
    cam_controls = cameracontrols(scene)
    cam3d_state = if cam_controls isa Camera3D
        fields = (:lookat, :upvector, :eyeposition, :fov, :near, :far)
        Dict((f => serialize_three(getfield(cam_controls, f)[]) for f in fields))
    else
        nothing
    end
    serialized = Dict(:pixelarea => pixel_area,
                      :backgroundcolor => lift(hexcolor, scene.backgroundcolor),
                      :clearscene => scene.clear,
                      :camera => serialize_camera(scene),
                      :plots => serialize_plots(scene, scene.plots),
                      :cam3d_state => cam3d_state,
                      :visible => scene.visible,
                      :uuid => js_uuid(scene))

    push!(serialized_scenes, serialized)
    foreach(child -> serialize_scene(child, serialized_scenes), scene.children)

    return serialized_scenes
end

function serialize_plots(scene::Scene, plots::Vector{T}, result=[]) where {T<:AbstractPlot}
    for plot in plots
        # if no plots inserted, this truely is an atomic
        if isempty(plot.plots)
            plot_data = serialize_three(scene, plot)
            plot_data[:uuid] = js_uuid(plot)
            push!(result, plot_data)
        else
            serialize_plots(scene, plot.plots, result)
        end
    end
    return result
end

function serialize_three(scene::Scene, plot::AbstractPlot)
    program = create_shader(scene, plot)
    mesh = serialize_three(program)
    mesh[:name] = string(AbstractPlotting.plotkey(plot)) * "-" * string(objectid(plot))
    mesh[:visible] = plot.visible
    mesh[:uuid] = js_uuid(plot)
    uniforms = mesh[:uniforms]
    updater = mesh[:uniform_updater]

    delete!(uniforms, :lightposition)

    if haskey(plot, :lightposition)
        eyepos = scene.camera.eyeposition
        lightpos = lift(plot.lightposition, eyepos,
                        typ=Vec3f0) do pos, eyepos
            return ifelse(pos == :eyeposition, eyepos, pos)::Vec3f0
        end
        uniforms[:lightposition] = serialize_three(lightpos[])
        on(lightpos) do value
            updater[] = [:lightposition, serialize_three(value)]
            return
        end
    end
    if haskey(plot, :space)
        mesh[:space] = plot.space[]
    end

    return mesh
end

function serialize_camera(scene::Scene)
    cam = scene.camera
    return lift(cam.view, cam.projection, cam.resolution) do v, p, res
        # projectionview updates with projection & view
        pv = cam.projectionview[]
        # same goes for eyeposition, since an eyepos change will trigger
        # a view matrix change!
        ep = cam.eyeposition[]
        pixel_space = cam.pixel_space[]
        return [serialize_three.((v, p, pv, res, ep, pixel_space))...]
    end
end


text:
#https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/09.%20Julia%20is%20fast.ipynb

using PyCall
using Libdl

py"""
def pysum(arr):
    s=0
    for a in arr:
        s+=a
    return s
"""

pysum=py"pysum"

@show pysum([1,2,3,4,5])

py"""
class Goma():
    def __init__(self,msg):
        self.msg=msg
    def say_msg(self):
        return " ".join([self.msg,'Goma goma kyukkyu'])
"""

Goma = py"Goma"
goma = Goma("Hi")
@show goma[:say_msg]()

C_code = """
#include<stddef.h>
double c_sum(size_t n, double *X){
    double s = 0.0;
    for(size_t i=0; i<n; ++i){
        s += X[i];
    }
    return s;
}
"""

const C_lib = tempname()
open(`gcc -fPIC -O3 -msse3 -xc -shared -o $(C_lib * "." * dlext) -`, "w") do f
    print(f, C_code)
end

c_sum(X::Array{Float64}) = ccall(("c_sum", C_lib), Float64, (Csize_t, Ptr{Float64}), length(X), X)
@show c_sum([1.0,2.0,3.0,4.0,5.0])

using InteractiveUtils
versioninfo()

#= output
pysum([1, 2, 3, 4, 5]) = 15
(goma[:say_msg])() = "Hi Goma goma kyukkyu"
c_sum([1.0, 2.0, 3.0, 4.0, 5.0]) = 15.0
Julia Version 0.7.0
Commit a4cb80f3ed (2018-08-08 06:46 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) M-5Y51 CPU @ 1.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)
[Finished in 9.8s]
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
__precompile__()

module Contour

using StaticArrays

export
    ContourLevel,
    Curve2,
    contour,
    contours,
    level,
    levels,
    lines,
    coordinates

import Base: push!, length, eltype, show

struct Curve2{T}
    vertices::Vector{SVector{2,T}}
end
Curve2(::Type{T}) where {T} = Curve2(SVector{2,T}[])
show(io::IO, ::MIME"text/plain", c2::Curve2) = write(io, "$(typeof(c2))\n  with $(length(c2.vertices)-1) vertices")
show(io::IO, ::MIME"text/plain", c2s::Vector{TC}) where {TC <: Curve2} = write(io, "$(typeof(c2s))\n  $(length(c2s)) contour line(s)")

struct ContourLevel{T}
    level::T
    lines::Vector{Curve2{T}}
end
ContourLevel(h::T) where {T <: AbstractFloat} = ContourLevel(h, Curve2{T}[])
ContourLevel(h::T) where {T} = ContourLevel(Float64(h))
show(io::IO, ::MIME"text/plain", cl::ContourLevel) = write(io, "$(typeof(cl))\n  at $(level(cl)) with $(length(lines(cl))) line(s)")
show(io::IO, ::MIME"text/plain", cls::Vector{CL}) where {CL <: ContourLevel} = write(io, "$(typeof(cls))\n  $(length(cls)) contour level(s)")
"""
`lines(c)` Extracts an iterable collection of isolines from a contour level.
Use [`coordinates`](@ref) to get the coordinates of a line.
"""
lines(cl::ContourLevel) = cl.lines
"""
`level(c)` Indicates the `z`-value at which the contour level `c` was traced.
"""
level(cl::ContourLevel) = cl.level

struct ContourCollection{Tlevel<:ContourLevel}
    contours::Vector{Tlevel}
end
ContourCollection() = ContourCollection(Float64)
ContourCollection(::Type{Tlevel}) where {Tlevel} = ContourCollection(ContourLevel{Tlevel}[])
show(io::IO, ::MIME"text/plain", cc::ContourCollection) = write(io, "$(typeof(cc))\n with $(length(levels(cc))) level(s).")

"""
Turns the output of [`contours`](@ref) into an iterable with each of the traced
contour levels. Each of the objects support [`level`](@ref) and
[`coordinates`](@ref).
"""
levels(cc::ContourCollection) = cc.contours

"""
`contour(x, y, z, level::Number)` Trace a single contour level, indicated by the
argument `level`.

You'll usually call [`lines`](@ref) on the output of `contour`, and then iterate
over the result.
"""
function contour(x, y, z, level::Number)
    # Todo: size checking on x,y,z
    trace_contour(x, y, z, level, get_level_cells(z, level))
end

"""
`contours` returns a set of isolines.

You'll usually call [`levels`](@ref) on the output of `contours`.
"""
contours(::Any...) = error("This method exists only for documentation purposes")

"""
`contours(x,y,z,levels)` Trace the contour levels indicated by the `levels`
argument.
"""
contours(x, y, z, levels) = ContourCollection([contour(x, y, z, l) for l in levels])

"""
`contours(x,y,z,Nlevels::Integer)` Trace `Nlevels` contour levels at heights
chosen by the library (using the  [`contourlevels`](@ref) function).
"""
function contours(x, y, z, Nlevels::Integer)
    contours(x, y, z, contourlevels(z, Nlevels))
end

"""
`contours(x,y,z)` Trace 10 automatically chosen contour levels.
"""
contours(x, y, z) = contours(x, y, z, 10)

"""
`contourlevels(z,n)` Examines the values of `z` and chooses `n` evenly spaced
levels to trace.
"""
function contourlevels(z, n)
    zmin, zmax = extrema(z)
    dz = (zmax - zmin) / (n + 1)
    range(zmin + dz; step = dz, length = n)
end

"""
`coordinates(c)` Returns the coordinates of the vertices of the contour line as
a tuple of lists.
"""
function coordinates(c::Curve2{T}) where {T}
    N = length(c.vertices)
    xlist = Vector{T}(undef, N)
    ylist = Vector{T}(undef, N)

    for (i, v) in enumerate(c.vertices)
        xlist[i] = v[1]
        ylist[i] = v[2]
    end
    xlist, ylist
end

# The marching squares algorithm defines 16 cell types
# based on the edges that a contour line enters and exits
# through. The edges of the cells are identified using
# compass directions, while the vertices are ordered as
# follows:
#
#      N
#  4 +---+ 3
# W  |   |  E
#  1 +---+ 2
#      S
#
# Each cell type is identified with 4 bits, with each
# bit corresponding to a vertex (MSB -> 4, LSB -> 1).
# A bit is set for vertex v_i is set if z(v_i) > h. So a cell
# where a contour line only enters from the W edge and exits
# through the N edge will have the cell type: 0b0111
# Note that there are two cases where there are two
# lines crossing through the same cell: 0b0101, 0b1010.
const N, S, E, W = (UInt8(1)), (UInt8(2)), (UInt8(4)), (UInt8(8))
const NS, NE, NW = N|S, N|E, N|W
const SN, SE, SW = S|N, S|E, S|W
const EN, ES, EW = E|N, E|S, E|W
const WN, WS, WE = W|N, W|S, W|E

const dirStr = ["N", "S", "NS", "E", "NE", "NS", "Invalid crossing",
                "W", "NW", "SW", "Invalid crossing", "WE"]

# The way a contour crossing goes through a cell is labeled
# by combining compass directions (e.g. a NW crossing connects
# the N edge and W edges of the cell).  The Cell type records
# the type of crossing that a cell contains.  While most
# cells will have only one crossing, cell type 5 and 10 will
# have two crossings.
struct Cell
    crossings::Vector{UInt8}
end

function get_next_edge!(cell::Cell, entry_edge::UInt8)
    for (i, edge) in enumerate(cell.crossings)
        if edge & entry_edge != 0
            next_edge = edge ⊻ entry_edge
            deleteat!(cell.crossings, i)

            return next_edge
        end
    end
    error("There is no edge containing ", entry_edge)
end

# Maps cell type to crossing types for non-ambiguous cells
const edge_LUT = [SW, SE, EW, NE, 0, NS, NW, NW, NS, 0, NE, EW, SE, SW]

function get_level_cells(z, h::Number)
    cells = Dict{(Tuple{Int,Int}),Cell}()
    xi_max, yi_max = size(z)

    local case::Int8

    for xi in 1:xi_max - 1
        for yi in 1:yi_max - 1
            case = 1(z[xi, yi] > h)     |
                   2(z[xi + 1, yi] > h)   |
                   4(z[xi + 1, yi + 1] > h) |
                   8(z[xi, yi + 1] > h)

            # Contour does not go through these cells
            if case == 0 || case == 15
                continue
            end

            # Process ambiguous cells (case 5 and 10) using
            # a bilinear interpolation of the cell-center value.
            if case == 5
                if 0.25(z[xi, yi] + z[xi, yi + 1] + z[xi + 1, yi] + z[xi + 1, yi + 1]) >= h
                    cells[(xi, yi)] = Cell([NW, SE])
                else
                    cells[(xi, yi)] = Cell([NE, SW])
                end
            elseif case == 10
                if 0.25(z[xi, yi] + z[xi, yi + 1] + z[xi + 1, yi] + z[xi + 1, yi + 1]) >= h
                    cells[(xi, yi)] = Cell([NE, SW])
                else
                    cells[(xi, yi)] = Cell([NW, SE])
                end
            else
                cells[(xi, yi)] = Cell([edge_LUT[case]])
            end
        end
    end

    return cells
end

# Some constants used by trace_contour

const fwd, rev = (UInt8(0)), (UInt8(1))

function add_vertex!(curve::Curve2{T}, pos::(Tuple{T,T}), dir::UInt8) where {T}
    if dir == fwd
        push!(curve.vertices, SVector{2,T}(pos...))
    else
        pushfirst!(curve.vertices, SVector{2,T}(pos...))
    end
end

# Given the row and column indices of the lower left
# vertex, add the location where the contour level
# crosses the specified edge.
function interpolate(x, y, z::AbstractMatrix{T}, h::Number, xi::Int, yi::Int, edge::UInt8) where {T <: AbstractFloat}
    if edge == W
        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi, yi]) / (z[xi, yi + 1] - z[xi, yi])
        x_interp = x[xi]
    elseif edge == E
        y_interp = y[yi] + (y[yi + 1] - y[yi]) * (h - z[xi + 1, yi]) / (z[xi + 1, yi + 1] - z[xi + 1, yi])
        x_interp = x[xi + 1]
    elseif edge == N
        y_interp = y[yi + 1]
        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi + 1]) / (z[xi + 1, yi + 1] - z[xi, yi + 1])
    elseif edge == S
        y_interp = y[yi]
        x_interp = x[xi] + (x[xi + 1] - x[xi]) * (h - z[xi, yi]) / (z[xi + 1, yi] - z[xi, yi])
    end

    return x_interp, y_interp
end

# Given a cell and a starting edge, we follow the contour line until we either
# hit the boundary of the input data, or we form a closed contour.
function chase!(cells, curve, x, y, z, h, xi_start, yi_start, entry_edge, xi_max, yi_max, dir)

    xi, yi = xi_start, yi_start

    # When the contour loops back to the starting cell, it is possible
    # for it to not intersect with itself.  This happens if the starting
    # cell contains a saddle-point. So a loop is only closed if the
    # contour returns to the starting edge of the starting cell
    loopback_edge = entry_edge

    while true
        cell = cells[(xi, yi)]
        exit_edge = get_next_edge!(cell, entry_edge)
        if length(cell.crossings) == 0
            delete!(cells, (xi, yi))
        end

        add_vertex!(curve, interpolate(x, y, z, h, xi, yi, exit_edge), dir)

        if exit_edge == N
            yi += 1
            entry_edge = S
        elseif exit_edge == S
            yi -= 1
            entry_edge = N
        elseif exit_edge == E
            xi += 1
            entry_edge = W
        elseif exit_edge == W
            xi -= 1
            entry_edge = E
        end
        !((xi, yi, entry_edge) != (xi_start, yi_start, loopback_edge) &&
           0 < yi < yi_max && 0 < xi < xi_max) && break
    end

    return xi, yi
end

function trace_contour(x, y, z, h::Number, cells::Dict{(Tuple{Int,Int}),Cell})

    contours = ContourLevel(h)

    local yi::Int
    local xi::Int
    local xi_0::Int
    local yi_0::Int

    local xi_max::Int
    local yi_max::Int

    (xi_max, yi_max) = size(z)

    # When tracing out contours, this algorithm picks an arbitrary
    # starting cell, then first follows the contour in one direction
    # until it either ends up where it started # or at one of the boundaries.
    # It then tries to trace the contour in the opposite direction.

    while length(cells) > 0
        contour = Curve2(promote_type(map(eltype, (x, y, z))...))

        # Pick initial box
        (xi_0, yi_0), cell = first(cells)
        (xi, yi) = (xi_0, yi_0)

        # Pick a starting edge
        crossing = first(cell.crossings)
        starting_edge = UInt8(0)
        for edge in [N, S, E, W]
            if edge & crossing != 0
                starting_edge = edge
                break
            end
        end

        # Add the contour entry location for cell (xi_0,yi_0)
        add_vertex!(contour, interpolate(x, y, z, h, xi_0, yi_0, starting_edge), fwd)

        # Start trace in forward direction
        (xi_end, yi_end) = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, fwd)

        if (xi_end, yi_end) == (xi_0, yi_0)
            push!(contours.lines, contour)
            continue
        end

        if starting_edge == N
            yi = yi_0 + 1
            starting_edge = S
        elseif starting_edge == S
            yi = yi_0 - 1
            starting_edge = N
        elseif starting_edge == E
            xi = xi_0 + 1
            starting_edge = W
        elseif starting_edge == W
            xi = xi_0 - 1
            starting_edge = E
        end

        if !(0 < yi < yi_max && 0 < xi < xi_max)
            push!(contours.lines, contour)
            continue
        end

        # Start trace in reverse direction
        (xi, yi) = chase!(cells, contour, x, y, z, h, xi, yi, starting_edge, xi_max, yi_max, rev)
        push!(contours.lines, contour)
    end

    return contours

end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3f57a160-9042-11eb-2981-a58c990f6dfe
begin
    import Pkg
    Pkg.activate(mktempdir())

    Pkg.add([
        Pkg.PackageSpec(name="Pluto"),
        Pkg.PackageSpec(name="PlutoUI"),
        Pkg.PackageSpec(name="CSV"),
        Pkg.PackageSpec(name="DataFrames"),
        Pkg.PackageSpec(name="FreqTables"),
        Pkg.PackageSpec(name="Plots"),
        Pkg.PackageSpec(name="StatsPlots"),
        Pkg.PackageSpec(name="Pipe"),
        Pkg.PackageSpec(name="Statistics"),
        Pkg.PackageSpec(name="DisplayAs"),
		Pkg.PackageSpec(name="MixedModels", rev="main"),
    ])

    using Pluto
	using PlutoUI
	using CSV
	using DataFrames
	using FreqTables
	using Plots
	using StatsPlots
	using Pipe
	using Statistics
	using MixedModels
	using DisplayAs
end

# ╔═╡ 92256604-8f75-11eb-0e65-4df8abfdc276
TableOfContents()

# ╔═╡ e385fbc8-8f75-11eb-08bf-232a909b8db2
md"""
# 0. Load Packages and Set WD
"""

# ╔═╡ ecb582ce-8fa8-11eb-0a7a-6398e820bfa4
theme(:wong)

# ╔═╡ 35723342-8d59-11eb-28ea-3f5af50234e9
cd("/Users/lzhan/Desktop/ITOM/Data"); pwd()

# ╔═╡ 2e42988c-903d-11eb-1f06-9f74eaa45252
md"""
**Levels of `Agent` variable**

- Ex1: LN - Left Normal; RN - Right Normal
- Ex2: LN/M - Left Normal/Mirror; RN/M - Right Normal/Mirror
- Ex3: LA/H - Left Arrow / Human; RA/H - Right Arrow / Human
"""

# ╔═╡ e1eb18d4-8de3-11eb-1d26-4b637d2ab921
md"""
# 1. Experiment One
"""

# ╔═╡ 6299cd52-8ecc-11eb-1558-27598f1704f1
fm1 = @formula(rt ~ 1 + agent + (1 + agent | participant) + (1 + agent | item));

# ╔═╡ e4516da4-8ecd-11eb-274c-2d51c0dd2ad3
contr1 = Dict(:agent => DummyCoding(base = "__"));

# ╔═╡ 79c1fa5c-8ed3-11eb-02ec-ab660761b467
md"""
# 2. Experiment Two
"""

# ╔═╡ 4c56aec6-902e-11eb-3fa7-e5b250746180
contr21 = Dict(:agent => DummyCoding(base = "LM"));

# ╔═╡ aa510800-902e-11eb-2db0-3bf83a5d517f
contr22 = Dict(:agent => DummyCoding(base = "RM"));

# ╔═╡ be05135c-8eec-11eb-3e6a-c3bd9c8a1950
md"""
# 3. Experiment Three
"""

# ╔═╡ 40e48996-902b-11eb-1699-a38768e1b11a
contr3 = Dict(:agent => DummyCoding(base = "LH"));

# ╔═╡ e8045e32-8f74-11eb-08d7-0781a63175a3
md"""
# 4. Appendixes
"""

# ╔═╡ 6110d00c-903e-11eb-2b38-99a816047bc3
md"""
**Column names and their meanings in the data**
- participant: Participant Number;
- type: Letter Type, N - Normal; M - Mirror;
- Orientation: Degree of rotation;
- human (Exp3 - face_orientation): L - Left; R -  Right; _ - No human;
- face_type: N - faced table; M - back(Ex2); _ - No human;
- Condition(Ex3): H - Human(faced_table); A - Arrow; __ - No;
- key\_JUDGE.keys (Ex3: key\_judge.keys): left - normal; right - mirror
- key\_JUDGE.rt (Ex3: key\_judge.rt): Response latency
- Number: Image Number
"""

# ╔═╡ 692a8fc4-8f96-11eb-0614-b1fd9fe4fe56
function key_2_TF(letter_type, key)
	if letter_type == "M"
		return key ∈ ["down", "right"] ? 1 : 0
	elseif letter_type == "N"
		return key ∈ ["up",    "left"] ? 1 : 0
	end
end

# ╔═╡ a110c5fe-8f94-11eb-0f0f-b71e5b9353a6
function read_rename_select(directory, rename_columns, select_columns)
	files = filter(x -> endswith(x, ".csv"), readdir(directory, join = true))
	list = [DataFrame(CSV.File(f))[1:(end-1), :] for f in files]
	for (ind, df) in enumerate(list)
		rename!(df, rename_columns)
		df.participant = string.(df.participant)
		df.item        = string.(df.participant)
		ncol(df) > 40 && (df.rt = df.rt - df.stimulus_duration)
		insertcols!(df, 1, :tf      => key_2_TF.(df.letter_type, df.key));
		insertcols!(df, 1, :agent   => df.agent_position .* df.agent_type);
		list[ind] = df[:, select_columns]
	end
	df = reduce(vcat, list)
	return df
end

# ╔═╡ 89265a5a-8f95-11eb-1d1d-0b9c84742e9b
function filter_data(data, correct_rate, rt_lower, rt_upper)
	participants = @pipe data |>
		groupby(_, :participant) |>
		combine(_, :tf => mean)  |>
		filter(:tf_mean => >=(correct_rate), _) |>
		select(_, :participant)
	dt = @pipe data |>
		dropmissing(_, :rt) |>
		filter(:participant => ∈(participants.participant), _) |>
		# filter(:tf => ==(1), _) |>
		filter(:rt => x -> x >= rt_lower && x <= rt_upper, _)
	return dt
end

# ╔═╡ 00147106-8f98-11eb-34e2-a98d0f6c046a
function summarize_4_plot(dt)
	dtm = @pipe dt |>
		groupby(_, [:letter_rotation, :agent]) |>
		combine(_, :rt => mean)
	new_rows = insertcols!(dtm[dtm.letter_rotation .== 0, 2:end], 1,
	:letter_rotation => 360)
	dtmn = vcat(dtm, new_rows)
	return dtmn
end

# ╔═╡ 10acbb8a-8fad-11eb-2d8a-913e1c002326
function polar_plot(dt; lower = 1, upper = 1.4)
	plot(
		(dt.letter_rotation) / 360 * 2π,
		dt.rt_mean,
		group = dt.agent,
		ylim = (lower, upper),
		seriestype = :line,
		projection = :polar,
		theta_direction = -2,
		set_theta_zero_location = "S",
		linewidth = 1.5,
		marker = 4)
end

# ╔═╡ 083f51ca-8de0-11eb-05d7-513b912438ed
## Experiments 1 and 2
rename1 = Dict(
	:Letter      => "letter",
	:type        => "letter_type",
	:Orientation => "letter_rotation",
	:human       => "agent_position",
	:face_type   => "agent_type",
	Symbol("trials.thisIndex") => "item",
	Symbol("key_JUDGE.keys")   => "key",
	Symbol("key_JUDGE.rt")     => "rt"
);

# ╔═╡ e8ca0eb0-8dd9-11eb-20a5-e1d5e04fa2e8
columns1 = vcat(collect(values(rename1)),
	"participant", "tf", "agent"); # [[1:4..., 6:7...]]

# ╔═╡ 35b40010-8f93-11eb-3d7a-bbd6941f95dc
ex1_original = read_rename_select("Experiment_1", rename1, columns1);

# ╔═╡ e876a08c-8ef7-11eb-3650-31ae7ed843b7
freqtable(ex1_original, :agent, :letter_rotation)

# ╔═╡ 07386f3e-8f96-11eb-20f4-09d27ddfa926
ex1 = filter_data(ex1_original, 0.8, 0.15, 2);

# ╔═╡ 5ef3d9c0-9031-11eb-3839-e7db22fbf21b
length(unique(ex1.participant))

# ╔═╡ 6ea41f1c-8fad-11eb-193d-dbc8b73fadc2
polar_plot(summarize_4_plot(ex1))

# ╔═╡ 2b3d9e00-8ecd-11eb-3d55-b92db51bd635
fit(MixedModel, fm1, ex1[ex1.letter_rotation .== 135, :], contrasts = contr1)

# ╔═╡ 5503914a-8eeb-11eb-3dca-b96cf673567d
fit(MixedModel, fm1, ex1[ex1.letter_rotation .== 225, :], contrasts = contr1)

# ╔═╡ 56a57220-8eeb-11eb-2cce-ed6739c33c9d
fit(MixedModel, fm1, ex1[ex1.letter_rotation .== 315, :], contrasts = contr1)

# ╔═╡ 308ef3b4-8f94-11eb-293f-2d2ddb3ace9f
ex2_original = read_rename_select("Experiment_2", rename1, columns1);

# ╔═╡ 83cb4264-8ef7-11eb-33be-69bced1d6db1
freqtable(ex2_original, :agent, :letter_rotation)

# ╔═╡ 2b93a9e8-8f96-11eb-3c59-5155ec37e662
ex2 = filter_data(ex2_original, 0.80, 0.15, 2);

# ╔═╡ 7de82660-9036-11eb-1c3f-49b8e64e495b
length(unique(ex2.participant))

# ╔═╡ ff604e88-902d-11eb-2481-751fb18eed11
ex2L = ex2[findall(∈(["LM", "LN"]), ex2.agent), :];

# ╔═╡ 2a97c1e8-8f56-11eb-3c7d-bdee7918a422
fit(MixedModel, fm1, ex2L[ex2L.letter_rotation .== 90, :], contrasts = contr21)

# ╔═╡ 67951164-902e-11eb-024f-eb45dfaeac7c
fit(MixedModel, fm1, ex2L[ex2L.letter_rotation .== 270, :], contrasts = contr21)

# ╔═╡ 8f26ab3e-902e-11eb-1b83-ed1c24af3325
fit(MixedModel, fm1, ex2L[ex2L.letter_rotation .== 315, :], contrasts = contr21)

# ╔═╡ 0f1791ec-902e-11eb-1927-2bdc219c0976
ex2R = ex2[findall(∈(["RM", "RN"]), ex2.agent), :];

# ╔═╡ c798e31c-8faa-11eb-22c3-69269a304f67
plot(polar_plot(summarize_4_plot(ex2L)),
	 polar_plot(summarize_4_plot(ex2R)), layout = 2)

# ╔═╡ 3a63a3ec-8f55-11eb-390f-b5262e5f8646
fit(MixedModel, fm1, ex2R[ex2R.letter_rotation .== 270, :], contrasts = contr22)

# ╔═╡ 9c771fc6-8f57-11eb-2408-33b89126fa4d
## Experiments 3
rename3 = Dict(
	:Letter      => "letter",
	:type        => "letter_type",
	:Orientation => "letter_rotation",
	Symbol("trials.thisIndex") => "item",
	Symbol("key_judge.keys")   => "key",
	Symbol("key_judge.rt")     => "rt",
	Symbol("face_orientation") => "agent_position",
	:Condition   => "agent_type"
);

# ╔═╡ fb0dfa1e-8f57-11eb-306d-bf0057bc7e12
columns3 = vcat(collect(values(rename3)), "participant", "tf", "agent"); # [[1,2,4,5,7,8]]

# ╔═╡ 53c5a418-8f94-11eb-21ed-6905d9a227ae
ex3_original = read_rename_select("Experiment_3", rename3, columns3);

# ╔═╡ 9b1b7fcc-8f58-11eb-1b2a-e1104e001cb7
freqtable(ex3_original, :agent, :letter_rotation)

# ╔═╡ 3a73aeae-8f96-11eb-0a03-f147003a4128
ex3 = filter_data(ex3_original, 0.8, 0.15, 2);

# ╔═╡ 66990d2e-9031-11eb-2cc7-09cdde9af93b
length(unique(ex3.participant))

# ╔═╡ 9771d9fa-902c-11eb-3ff5-0fec27d703f1
ex3L = ex3[findall(∈(["LA", "LH"]), ex3.agent), :];

# ╔═╡ 634aeff8-902d-11eb-31aa-4121d2a7c318
fit(MixedModel, fm1, ex3L[ex3L.letter_rotation .== 0, :], contrasts = contr3)

# ╔═╡ c4537ec6-8f62-11eb-25b5-b5a7ea79f741
fit(MixedModel, fm1, ex3L[ex3L.letter_rotation .== 315, :], contrasts = contr3)

# ╔═╡ 85ae5cde-902c-11eb-282c-79a93395e5f0
ex3R = ex3[findall(∈(["RA", "RH"]), ex3.agent), :];

# ╔═╡ d74bde4c-8fad-11eb-20d7-41b985c55506
plot(polar_plot(summarize_4_plot(ex3L), lower = 0.5, upper = 0.9),
	 polar_plot(summarize_4_plot(ex3R), lower = 0.5, upper = 0.9),
	 layout = 2)

# ╔═╡ Cell order:
# ╟─92256604-8f75-11eb-0e65-4df8abfdc276
# ╟─e385fbc8-8f75-11eb-08bf-232a909b8db2
# ╠═3f57a160-9042-11eb-2981-a58c990f6dfe
# ╠═ecb582ce-8fa8-11eb-0a7a-6398e820bfa4
# ╠═35723342-8d59-11eb-28ea-3f5af50234e9
# ╟─2e42988c-903d-11eb-1f06-9f74eaa45252
# ╟─e1eb18d4-8de3-11eb-1d26-4b637d2ab921
# ╠═35b40010-8f93-11eb-3d7a-bbd6941f95dc
# ╠═e876a08c-8ef7-11eb-3650-31ae7ed843b7
# ╠═07386f3e-8f96-11eb-20f4-09d27ddfa926
# ╠═5ef3d9c0-9031-11eb-3839-e7db22fbf21b
# ╠═6ea41f1c-8fad-11eb-193d-dbc8b73fadc2
# ╠═6299cd52-8ecc-11eb-1558-27598f1704f1
# ╠═e4516da4-8ecd-11eb-274c-2d51c0dd2ad3
# ╠═2b3d9e00-8ecd-11eb-3d55-b92db51bd635
# ╠═5503914a-8eeb-11eb-3dca-b96cf673567d
# ╠═56a57220-8eeb-11eb-2cce-ed6739c33c9d
# ╟─79c1fa5c-8ed3-11eb-02ec-ab660761b467
# ╠═308ef3b4-8f94-11eb-293f-2d2ddb3ace9f
# ╠═83cb4264-8ef7-11eb-33be-69bced1d6db1
# ╠═2b93a9e8-8f96-11eb-3c59-5155ec37e662
# ╠═7de82660-9036-11eb-1c3f-49b8e64e495b
# ╠═ff604e88-902d-11eb-2481-751fb18eed11
# ╠═0f1791ec-902e-11eb-1927-2bdc219c0976
# ╠═c798e31c-8faa-11eb-22c3-69269a304f67
# ╠═4c56aec6-902e-11eb-3fa7-e5b250746180
# ╠═2a97c1e8-8f56-11eb-3c7d-bdee7918a422
# ╠═67951164-902e-11eb-024f-eb45dfaeac7c
# ╠═8f26ab3e-902e-11eb-1b83-ed1c24af3325
# ╠═aa510800-902e-11eb-2db0-3bf83a5d517f
# ╠═3a63a3ec-8f55-11eb-390f-b5262e5f8646
# ╟─be05135c-8eec-11eb-3e6a-c3bd9c8a1950
# ╠═53c5a418-8f94-11eb-21ed-6905d9a227ae
# ╠═9b1b7fcc-8f58-11eb-1b2a-e1104e001cb7
# ╠═3a73aeae-8f96-11eb-0a03-f147003a4128
# ╠═66990d2e-9031-11eb-2cc7-09cdde9af93b
# ╠═9771d9fa-902c-11eb-3ff5-0fec27d703f1
# ╠═85ae5cde-902c-11eb-282c-79a93395e5f0
# ╠═d74bde4c-8fad-11eb-20d7-41b985c55506
# ╠═40e48996-902b-11eb-1699-a38768e1b11a
# ╠═634aeff8-902d-11eb-31aa-4121d2a7c318
# ╠═c4537ec6-8f62-11eb-25b5-b5a7ea79f741
# ╟─e8045e32-8f74-11eb-08d7-0781a63175a3
# ╟─6110d00c-903e-11eb-2b38-99a816047bc3
# ╠═a110c5fe-8f94-11eb-0f0f-b71e5b9353a6
# ╠═692a8fc4-8f96-11eb-0614-b1fd9fe4fe56
# ╠═89265a5a-8f95-11eb-1d1d-0b9c84742e9b
# ╠═00147106-8f98-11eb-34e2-a98d0f6c046a
# ╠═10acbb8a-8fad-11eb-2d8a-913e1c002326
# ╠═083f51ca-8de0-11eb-05d7-513b912438ed
# ╠═e8ca0eb0-8dd9-11eb-20a5-e1d5e04fa2e8
# ╠═9c771fc6-8f57-11eb-2408-33b89126fa4d
# ╠═fb0dfa1e-8f57-11eb-306d-bf0057bc7e12

max_stars_repo_name: likanzhan/statistics

max_issues_repo_name: likanzhan/statistics

max_forks_repo_name: likanzhan/statistics


text:
"""
These are two ways to calculate the total burden of a population, attributed to
each cause, by short-circuiting the complete sum. The complete sum has 2^n terms,
where n is the number of causes. These examples try to cut off after burden
terms get small or the number of comorbidities gets high. Neither works super well.
"""

function nextburden(state, leaf = false)
    which = copy(state)
    cstate = cumsum(which)
    depth = cstate[end]
    if depth == 0
        which[1] = true
        depth = 1
    else
        idx = findfirst(cstate .== depth)
        if idx == length(which)
            if depth == 1
                fill!(which, false)
            else
                depth -= 1
                which[idx] = false
                up_one = findfirst(cstate .== depth)
                which[up_one] = false
                which[up_one + 1] = true
            end
        elseif leaf
            which[idx] = false
            which[idx + 1] = true
        else
            which[idx + 1] = true
        end
    end
    which
end

function show_iterate_burdens(n)
    state = zeros(Bool, n)
    for i in 1:2^n
        state = nextburden(state)
        println(state)
    end
end

"""
  treeburden(weights, prevalences)

Calculate the burden by searching term-by-term until terms get too small.
"""
function treeburden(weights, prevalences, epsilon)
    n = length(prevalences)
    b = zeros(Float64, n)
    which = nextburden(zeros(Bool, n))
    iter_count = 0
    while any(which)
        term = exact_burden_term(weights, prevalences, which)
        b += term
        if sum(term) >= epsilon
            which = nextburden(which)
        else
            which = nextburden(which, true)
        end
        iter_count += 1
    end
    iterfrac = iter_count / (2^n - 1)
    println("treeburden $(iterfrac)")
    b
end

function depthburden(weights, prevalences, maxdepth)
    n = length(prevalences)
    b = zeros(Float64, n)
    which = nextburden(zeros(Bool, n))
    iter_count = 0
    while any(which)
        term = exact_burden_term(weights, prevalences, which)
        b += term
        if sum(which) < maxdepth
            which = nextburden(which)
        else
            which = nextburden(which, true)
        end
        iter_count += 1
    end
    iterfrac = iter_count / (2^n - 1)
    println("treeburden $(iterfrac)")
    b
end

function tree_burden_random(cnt, epsilon)
    w = rand(cnt)
    p = rand(cnt)
    exact = exact_burden(w, p)
    tree = treeburden(w, p, epsilon)
    [sum(tree), total_burden(w, p), sum(exact), maximum(abs.(tree .- exact))]
end

"""
Using an arbitrary depth cutoff performs pretty poorly too. Even evaluating
50% of the total possible terms results in 30% error.
"""
function depth_burden_random(cnt, depth)
    w = rand(cnt)
    p = rand(cnt)
    exact = exact_burden(w, p)
    tree = depthburden(w, p, depth)
    total = sum(exact)
    trial = tree * total / sum(tree)
    [sum(tree), total_burden(w, p), sum(exact), maximum(abs.((trial .- exact) ./ exact))]
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "util_linux"
version = v"2.35.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://mirrors.edge.kernel.org/pub/linux/utils/util-linux/v$(version.major).$(version.minor)/util-linux-$(version.major).$(version.minor).tar.gz",
                  "98acab129a8490265052e6c1e033ca96d68758a13bb7fcd232c06bf16cc96238")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/util-linux-*
export CPPFLAGS="-I${includedir}"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-makeinstall-chown
make -j${nproc}
make install
"""

# Build only for Linux
platforms = filter(Sys.islinux, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("chmem", :chmem),
    ExecutableProduct("partx", :partx, "sbin"),
    ExecutableProduct("fallocate", :fallocate),
    ExecutableProduct("lsipc", :lsipc),
    ExecutableProduct("findfs", :findfs, "sbin"),
    ExecutableProduct("whereis", :whereis),
    ExecutableProduct("lslogins", :lslogins),
    ExecutableProduct("fsfreeze", :fsfreeze, "sbin"),
    ExecutableProduct("ctrlaltdel", :ctrlaltdel, "sbin"),
    ExecutableProduct("lsns", :lsns),
    ExecutableProduct("sfdisk", :sfdisk, "sbin"),
    ExecutableProduct("look", :look),
    ExecutableProduct("taskset", :taskset),
    ExecutableProduct("mkfs.bfs", :mkfs_bfs, "sbin"),
    ExecutableProduct("chrt", :chrt),
    ExecutableProduct("resizepart", :resizepart, "sbin"),
    ExecutableProduct("rfkill", :rfkill, "sbin"),
    ExecutableProduct("uuidparse", :uuidparse),
    ExecutableProduct("umount", :umount),
    ExecutableProduct("fsck.minix", :fsck_minix, "sbin"),
    ExecutableProduct("hwclock", :hwclock, "sbin"),
    ExecutableProduct("prlimit", :prlimit),
    ExecutableProduct("switch_root", :switch_root, "sbin"),
    LibraryProduct("libmount", :libmount),
    ExecutableProduct("mcookie", :mcookie),
    ExecutableProduct("fstrim", :fstrim, "sbin"),
    ExecutableProduct("delpart", :delpart, "sbin"),
    ExecutableProduct("fsck.cramfs", :fsck_cramfs, "sbin"),
    ExecutableProduct("fincore", :fincore),
    ExecutableProduct("hardlink", :hardlink),
    LibraryProduct("libuuid", :libuuid),
    ExecutableProduct("ipcs", :ipcs),
    ExecutableProduct("eject", :eject),
    ExecutableProduct("uuidd", :uuidd, "sbin"),
    ExecutableProduct("nologin", :nologin, "sbin"),
    ExecutableProduct("readprofile", :readprofile, "sbin"),
    ExecutableProduct("fdformat", :fdformat, "sbin"),
    ExecutableProduct("lscpu", :lscpu),
    ExecutableProduct("col", :col),
    ExecutableProduct("addpart", :addpart, "sbin"),
    ExecutableProduct("sulogin", :sulogin, "sbin"),
    ExecutableProduct("getopt", :getopt),
    ExecutableProduct("rename", :rename_bin),
    ExecutableProduct("blockdev", :blockdev, "sbin"),
    ExecutableProduct("nsenter", :nsenter),
    ExecutableProduct("scriptreplay", :scriptreplay),
    ExecutableProduct("mkfs.cramfs", :mkfs_cramfs, "sbin"),
    ExecutableProduct("logger", :logger),
    ExecutableProduct("wall", :wall),
    ExecutableProduct("renice", :renice),
    ExecutableProduct("blkid", :blkid, "sbin"),
    LibraryProduct("libblkid", :libblkid),
    LibraryProduct("libfdisk", :libfdisk),
    ExecutableProduct("column", :column),
    ExecutableProduct("rtcwake", :rtcwake, "sbin"),
    ExecutableProduct("pivot_root", :pivot_root, "sbin"),
    ExecutableProduct("ipcrm", :ipcrm),
    ExecutableProduct("flock", :flock),
    ExecutableProduct("utmpdump", :utmpdump),
    ExecutableProduct("ionice", :ionice),
    ExecutableProduct("chcpu", :chcpu, "sbin"),
    ExecutableProduct("mesg", :mesg),
    ExecutableProduct("mkfs", :mkfs, "sbin"),
    ExecutableProduct("hexdump", :hexdump),
    ExecutableProduct("agetty", :agetty, "sbin"),
    ExecutableProduct("blkdiscard", :blkdiscard, "sbin"),
    ExecutableProduct("script", :script),
    ExecutableProduct("mountpoint", :mountpoint),
    ExecutableProduct("scriptlive", :scriptlive),
    ExecutableProduct("findmnt", :findmnt),
    ExecutableProduct("mkfs.minix", :mkfs_minix, "sbin"),
    ExecutableProduct("blkzone", :blkzone, "sbin"),
    ExecutableProduct("swapoff", :swapoff, "sbin"),
    ExecutableProduct("lslocks", :lslocks),
    ExecutableProduct("cal", :cal),
    ExecutableProduct("setarch", :setarch),
    ExecutableProduct("lsmem", :lsmem),
    ExecutableProduct("losetup", :losetup, "sbin"),
    ExecutableProduct("raw", :raw, "sbin"),
    ExecutableProduct("colrm", :colrm),
    ExecutableProduct("dmesg", :dmesg),
    ExecutableProduct("rev", :rev),
    ExecutableProduct("fsck", :fsck, "sbin"),
    ExecutableProduct("zramctl", :zramctl, "sbin"),
    ExecutableProduct("mkswap", :mkswap, "sbin"),
    ExecutableProduct("wipefs", :wipefs, "sbin"),
    ExecutableProduct("choom", :choom),
    LibraryProduct("libsmartcols", :libsmartcols),
    ExecutableProduct("namei", :namei),
    ExecutableProduct("mount", :mount),
    ExecutableProduct("setsid", :setsid),
    ExecutableProduct("swapon", :swapon, "sbin"),
    ExecutableProduct("ldattach", :ldattach, "sbin"),
    ExecutableProduct("swaplabel", :swaplabel, "sbin"),
    ExecutableProduct("ipcmk", :ipcmk),
    ExecutableProduct("unshare", :unshare),
    ExecutableProduct("isosize", :isosize),
    ExecutableProduct("colcrt", :colcrt),
    ExecutableProduct("last", :last_bin),
    ExecutableProduct("fdisk", :fdisk, "sbin"),
    ExecutableProduct("uuidgen", :uuidgen),
    ExecutableProduct("kill", :kill_bin),
    ExecutableProduct("lsblk", :lsblk),
    ExecutableProduct("wdctl", :wdctl)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Gettext_jll", uuid="78b55507-aeef-58d4-861c-77aaff3498b1")),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
module PowerAnalysis
    using Distributions
    using HypothesisTests
    using Roots
    using QuadGK

    export
        power,
        effect_size,
        sample_size,
        type_1,
        type_2,
        type_s,
        type_m,
        exaggeration_factor

    include("utils.jl")
    include("generics.jl")
    include("t_test.jl")
    include("z_test.jl")
end


text:
@testset "arbcall" begin
    Mag = Arblib.Mag
    mag_struct = Arblib.mag_struct
    arf_struct = Arblib.arf_struct
    arb_struct = Arblib.arb_struct
    acb_struct = Arblib.acb_struct

    @testset "Carg" begin
        # Supported types
        for (str, name, isconst, jltype, ctype) in (
            ("mag_t res", :res, false, Arblib.MagLike, Ref{mag_struct}),
            ("arf_t res", :res, false, Arblib.ArfLike, Ref{arf_struct}),
            ("arb_t res", :res, false, Arblib.ArbLike, Ref{arb_struct}),
            ("acb_t res", :res, false, Arblib.AcbLike, Ref{acb_struct}),
            ("const mag_t x", :x, true, Arblib.MagLike, Ref{mag_struct}),
            ("const arf_t x", :x, true, Arblib.ArfLike, Ref{arf_struct}),
            ("const arb_t x", :x, true, Arblib.ArbLike, Ref{arb_struct}),
            ("const acb_t x", :x, true, Arblib.AcbLike, Ref{acb_struct}),
            (
                "arf_rnd_t rnd",
                :rnd,
                false,
                Union{Arblib.arb_rnd,RoundingMode},
                Arblib.arb_rnd,
            ),
            ("mpfr_t x", :x, false, BigFloat, Ref{BigFloat}),
            (
                "mpfr_rnd_t rnd",
                :rnd,
                false,
                Union{Base.MPFR.MPFRRoundingMode,RoundingMode},
                Base.MPFR.MPFRRoundingMode,
            ),
            ("mpz_t x", :x, false, BigInt, Ref{BigInt}),
            ("int flags", :flags, false, Integer, Cint),
            ("slong x", :x, false, Integer, Int),
            ("ulong x", :x, false, Unsigned, UInt),
            ("double x", :x, false, Base.GMP.CdoubleMax, Cdouble),
            ("slong * x", :x, false, Vector{<:Integer}, Ref{Int}),
            ("ulong * x", :x, false, Vector{<:Unsigned}, Ref{UInt}),
            ("const char * inp", :inp, true, AbstractString, Cstring),
            ("arb_ptr v", :v, false, Arblib.ArbVectorLike, Ptr{arb_struct}),
            ("arb_srcptr res", :res, true, Arblib.ArbVectorLike, Ptr{arb_struct}),
            ("acb_ptr v", :v, false, Arblib.AcbVectorLike, Ptr{acb_struct}),
            ("acb_srcptr res", :res, true, Arblib.AcbVectorLike, Ptr{acb_struct}),
        )
            arg = Arblib.Carg(str)
            @test Arblib.name(arg) == name
            @test Arblib.isconst(arg) == isconst
            @test Arblib.jltype(arg) == jltype
            @test Arblib.ctype(arg) == ctype
        end

        # Unsupported types
        for str in (
            "FILE * file",
            "fmpr_t x",
            "fmpr_rnd_t rnd",
            "flint_rand_t state",
            "bool_mat_t mat",
        )
            @test_throws Arblib.UnsupportedArgumentType arg = Arblib.Carg(str)
        end

        # Unimplemented types
        for str in (
            "fmpz_t x",
            "fmpq_t x",
            "mag_ptr res",
            "const fmpz_t x",
            "const fmpq_t x",
            "mag_srcptr res",

            # Internal types
            "mp_limb_t lo",
            "mp_bitcnt_t r",
            "mp_ptr ycos",
            "mp_srcptr x",
            "mp_limb_t * error",
            "mp_bitcnt_t * Qexp",
        )
            @test_throws KeyError arg = Arblib.Carg(str)
        end

        # Parse errors
        for str in ()
            @test_throws ArgumentError arg = Arblib.Carg(str)
        end
    end

    @testset "jlfname" begin
        for (arbfname, name) in (# Supported types
            ("mag_set", :set),
            ("mag_set_d", :set),
            ("mag_set_ui", :set),
            ("arf_set", :set),
            ("arf_set_ui", :set),
            ("arf_set_si", :set),
            ("arf_set_d", :set),
            ("arb_set", :set),
            ("arf_set_mpz", :set),
            ("arf_set_mpfr", :set),
            ("arb_set_arf", :set),
            ("arb_set_si", :set),
            ("arb_set_ui", :set),
            ("arb_set_d", :set),
            ("acb_set", :set),
            ("acb_set_ui", :set),
            ("acb_set_si", :set),
            ("acb_set_d", :set),
            ("acb_set_arb", :set),
            ("acb_set_si_si", :set),
            ("acb_set_d_d", :set),
            ("acb_set_arb_arb", :set),

            # Unsupported types
            ("arf_set_fmpz", :set_fmpz),
            ("acb_set_fmpq", :set_fmpq),
            ("arb_bin_uiui", :bin_uiui),

            # Deprecated types
            ("arf_set_fmpr", :set_fmpr),

            # Underscore methods
            ("_acb_vec_set", :set),

            # Removing scalar
            ("_arb_vec_scalar_mul", :mul),

            # Some special cases to be aware of and maybe change
            ("mag_set_d_lower", :set_d_lower),
            ("arb_ui_div", :ui_div),
            ("arb_rising_ui_rec", :rising_ui_rec),
            ("arb_zeta_ui_vec", :zeta_ui_vec),
        )
            @test Arblib.jlfname(arbfname) == name
            @test Arblib.jlfname(arbfname, inplace = true) == Symbol(name, :!)
        end
    end

    @testset "returntype" begin
        # Supported return types
        for (str, T) in (
            ("void arb_init(arb_t x)", Nothing),
            ("slong arb_rel_error_bits(const arb_t x)", Int),
            ("int arb_is_zero(const arb_t x)", Int32),
            ("double arf_get_d(const arf_t x, arf_rnd_t rnd)", Float64),
            ("acb_ptr _acb_vec_init(slong n)", AcbVector),
        )
            @test Arblib.returntype(Arblib.Arbfunction(str)) == T
        end

        # Unsupported return types
        for str in ("mag_ptr _mag_vec_init(slong n)",)
            @test_throws KeyError Arblib.Arbfunction(str)
        end

        # Return types with parse errors
        for str in ()
            @test_throws ArgumentError Arblib.Arbfunction(str)
        end
    end

    @testset "Length argument detection" begin
        len_keywords = Set{Symbol}()
        prev_carg = Arblib.Carg("acb_srcptr A")
        prev_carg_bad = Arblib.Carg("acb_t x")
        carg1 = Arblib.Carg("slong lenA")
        carg2 = Arblib.Carg("slong len")
        carg3 = Arblib.Carg("slong n")
        carg4 = Arblib.Carg("slong m")

        @test Arblib.is_length_argument(carg1, prev_carg, len_keywords)
        @test Arblib.is_length_argument(carg2, prev_carg, len_keywords)
        @test Arblib.is_length_argument(carg3, prev_carg, len_keywords)
        @test !Arblib.is_length_argument(carg4, prev_carg, len_keywords)
        @test !Arblib.is_length_argument(carg1, prev_carg_bad, len_keywords)
        @test !Arblib.is_length_argument(carg2, prev_carg_bad, len_keywords)
        @test !Arblib.is_length_argument(carg3, prev_carg_bad, len_keywords)
        @test !Arblib.is_length_argument(carg4, prev_carg_bad, len_keywords)

        kwargs = []
        Arblib.extract_length_argument!(kwargs, len_keywords, carg1, prev_carg)
        @test kwargs[1] == :($(Expr(:kw, :(lenA::Integer), :(length(A)))))
        @test :lenA ∈ len_keywords
        @test !Arblib.is_length_argument(carg1, prev_carg, len_keywords)
    end

    @testset "Predicate detection" begin
        for sig in (
            "int acb_is_zero(const acb_t z)",
            "int acb_is_int_2exp_si(const acb_t x, slong e)",
            "int acb_equal_si(const acb_t x, slong y)",
            "int arb_ne(const arb_t x, const arb_t y)",
            "int acb_overlaps(const acb_t x, const acb_t y)",
            "int acb_contains(const acb_t x, const acb_t y)",
            "int acb_is_real(const acb_t x)",
            "int _acb_vec_is_zero(acb_srcptr vec, slong len)",
            "int _arb_vec_is_finite(arb_srcptr x, slong len)",
        )
            @test Arblib.ispredicate(Arblib.Arbfunction(sig))
        end

        for sig in (
            "int acb_mat_lu_classical(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)",
            "int acb_mat_lu_recursive(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)",
            "int acb_mat_lu(slong * perm, acb_mat_t LU, const acb_mat_t A, slong prec)",
            "int acb_mat_inv(acb_mat_t X, const acb_mat_t A, slong prec)",
            "int arb_sgn_nonzero(const arb_t x)",
            "int arf_set_round(arf_t res, const arf_t x, slong prec, arf_rnd_t rnd)",
            "int arf_cmp(const arf_t x, const arf_t y)",
        )
            @test !Arblib.ispredicate(Arblib.Arbfunction(sig))
        end
    end

    @testset "jlargs" begin
        for (str, args, kwargs) in (
            ("void arb_init(arb_t x)", [:(x::$(Arblib.ArbLike))], Expr[]),
            (
                "void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)",
                [:(z::$(Arblib.ArbLike)), :(x::$(Arblib.ArbLike)), :(y::$(Arblib.ArbLike))],
                [Expr(:kw, :(prec::Integer), :(_precision(z)))],
            ),
            (
                "int arf_add(arf_t res, const arf_t x, const arf_t y, slong prec, arf_rnd_t rnd)",
                [
                    :(res::$(Arblib.ArfLike)),
                    :(x::$(Arblib.ArfLike)),
                    :(y::$(Arblib.ArfLike)),
                ],
                [
                    Expr(:kw, :(prec::Integer), :(_precision(res))),
                    Expr(
                        :kw,
                        :(rnd::Union{$(Arblib.arb_rnd),RoundingMode}),
                        :(RoundNearest),
                    ),
                ],
            ),
            (
                "int _acb_vec_is_zero(acb_srcptr vec, slong len)",
                [:(vec::$(Arblib.AcbVectorLike))],
                [:($(Expr(:kw, :(len::Integer), :(length(vec)))))],
            ),
        )
            (a, k) = Arblib.jlargs(Arblib.Arbfunction(str))
            @test a == args
            @test k == kwargs
        end
    end

    @testset "arbsignature" begin
        for str in (
            "void arb_init(arb_t x)",
            "slong arb_rel_error_bits(const arb_t x)",
            "int arb_is_zero(const arb_t x)",
            "void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)",
            "void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec)",
            "void arb_add_ui(arb_t z, const arb_t x, ulong y, slong prec)",
            "void arb_add_si(arb_t z, const arb_t x, slong y, slong prec)",
            "void arb_sin(arb_t s, const arb_t x, slong prec)",
            "void arb_cos(arb_t c, const arb_t x, slong prec)",
            "void arb_sin_cos(arb_t s, arb_t c, const arb_t x, slong prec)",
            # Pointer
            "char * arb_get_str(const arb_t x, slong n, ulong flags)",
        )
            @test Arblib.arbsignature(Arblib.Arbfunction(str)) == str
        end
    end
end