text:
############################
# Define BCs and Operators #
############################

# Get boundaries of cylinders

locs = [(10.5, 0.0), (5.5, 0.0), (0.5, 0.0)]
cyl = vcat(map(locs) do loc
    findall(p -> ( ((p[1] - loc[1])^2 + (p[2] - loc[2])^2) <= 0.5^2 + 1e-4),s[:point])
end...)
fuzzy_bound = unique(vcat(incident(s, cyl, :∂v1)..., incident(s, cyl, :∂v0)...))
cyl_edge = filter(e -> (s[e, :∂v1] ∈ cyl)&&(s[e, :∂v0] ∈ cyl), fuzzy_bound)

cyl_bound = vcat(map(locs) do loc
    findall(p -> (0.5^2 - 1e-4 <= ((p[1] - loc[1])^2 + (p[2] - loc[2])^2) <= 0.5^2 + 1e-4),s[:point])
end...)
fuzzy_boundary = unique(vcat(incident(s, cyl_bound, :∂v1)..., incident(s, cyl_bound, :∂v0)...))
cyl_bound_edge = filter(e -> (s[e, :∂v1] ∈ cyl_bound)&&(s[e, :∂v0] ∈ cyl_bound), fuzzy_boundary)
cyl_inner = filter(p -> !(p ∈ cyl_bound), cyl)
slip_edge = filter!(p -> !(p ∈ cyl_bound_edge), cyl_edge)

k_col = fill(k₁, ne(s))
k_col[cyl_edge] .= k₂
k = diagm(k_col)

# Get other boundaries
funcs = sym2func(sd)
∂₀ = Examples.boundary_inds(Val{0}, s)
∂₁ = Examples.boundary_inds(Val{1}, s)

∂ₒ₀ = ∂₀[findall(p-> -9 <= p[1] <= 20 && -9 <= p[2] <= 9, s[∂₀, :point])]

lx = -10.0
rx = 21.0
ty = 15.5
by = -15.5

∂ₗ₀ = ∂₀[findall(p-> p[1] <= lx + 1e-4, s[∂₀, :point])]
∂ᵣ₀ = ∂₀[findall(p-> p[1] >= rx - 1e-4, s[∂₀, :point])]
∂ₜ₀ = ∂₀[findall(p-> p[2] >= ty - 1e-4, s[∂₀, :point])]
∂ᵦ₀ = ∂₀[findall(p-> p[2] <= by + 1e-4, s[∂₀, :point])]
∂ₑ₀ = vcat(∂ₗ₀, ∂ᵣ₀, ∂ₜ₀, ∂ᵦ₀)

∂ₗ₁ = Examples.bound_edges(s, ∂ₗ₀)
∂ᵣ₁ = Examples.bound_edges(s, ∂ᵣ₀)
∂ₑ₁ = Examples.bound_edges(s, ∂ₑ₀)

# Just induced edges
#∂₁₊ = Examples.bound_edges(s, ∂₀)

# Full adjacent boundary
#∂₁₊ = Examples.adj_edges(s, ∂₀)
# Just edges and corners
#=
∂_corners = Examples.adj_edges(s,
                vcat(intersect(∂ₗ₀, ∂ₜ₀), intersect(∂ₗ₀, ∂ᵦ₀),
                     intersect(∂ᵣ₀, ∂ₜ₀), intersect(∂ᵣ₀, ∂ᵦ₀)))
@show ∂_corners
∂₁₊ = collect(unique(vcat(∂ₑ₁, ∂_corners)))
=#
# Full triangle on boundary

∂₁₊ = Examples.adj_edges(s, ∂₀)
#∂_points = unique(vcat(s[∂₁₊, :∂v0], s[∂₁₊, :∂v1]))
#∂₁₊ = Examples.bound_edges(s, ∂_points)

∂ₗ₁₊ = Examples.adj_edges(s, ∂ₗ₀)
∂ᵣ₁₊ = Examples.adj_edges(s, ∂ᵣ₀)
∂ₑ₁₊ = Examples.adj_edges(s, ∂ₑ₀)
∂_points = unique(vcat(s[∂ₑ₁₊, :∂v0], s[∂ₑ₁₊, :∂v1]))
∂ₑ₁₊ = Examples.bound_edges(s, ∂_points)

∂ₗ₀₊ = unique(vcat(s[∂ₗ₁₊, :∂v1], s[∂ₗ₁₊, :∂v0]))
∂ᵣ₀₊ = unique(vcat(s[∂ᵣ₁₊, :∂v1], s[∂ᵣ₁₊, :∂v0]))
∂ₑ₀₊ = unique(vcat(s[∂ₑ₁₊, :∂v1], s[∂ₑ₁₊, :∂v0]))

max_stars_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_issues_repo_name: AlgebraicJulia/DECAPODES-Benchmarks

max_forks_repo_name: AlgebraicJulia/DECAPODES-Benchmarks


text:
using ExtensibleUnions
using Test

abstract type A1 end
@test_throws ArgumentError extensibleunion!(A1)

@test_throws ArgumentError extensibleunion!(Int)

struct S1
    x
end
@test_throws ArgumentError extensibleunion!(S1)

mutable struct S2 end
@test_throws ArgumentError extensibleunion!(S2)

abstract type A2 end
struct S3 <: A2
end
@test_throws ArgumentError extensibleunion!(S3)

struct S4 end
@test !isextensibleunion(S4)
extensibleunion!(S4)
@test isextensibleunion(S4)
extensibleunion!(S4)
@test isextensibleunion(S4)
extensibleunion!(S4)
@test isextensibleunion(S4)

struct S6 end
extensibleunion!(S6)
@test !ExtensibleUnions.unioncurrentlycontains(S6, String)
addtounion!(S6, String)
@test ExtensibleUnions.unioncurrentlycontains(S6, String)

struct S7 end
extensibleunion!(S7)
@test !ExtensibleUnions.unioncurrentlycontains(S7, String)
addtounion!(S7, String)
@test ExtensibleUnions.unioncurrentlycontains(S7, String)

struct S8 end
@test_throws ArgumentError addtounion!(S8)
@test_throws ArgumentError ExtensibleUnions.unioncurrentlycontains(S8, String)

let
    struct MyUnion end
    extensibleunion!(MyUnion)

    struct Foo end
    addtounion!(MyUnion, Foo)

    foo(x::Tuple{Int, T}) where {T} = x[1]

    foo(::MyUnion) = "boo!"
    extensiblefunction!(foo, MyUnion)

    @test foo(Foo()) == "boo!"
end

let
    a = Vector{T} where T
    ExtensibleUnions._replace_types(a)
    ExtensibleUnions._replace_types(a.var)
    ExtensibleUnions._replace_types(a, nothing => nothing)
    ExtensibleUnions._replace_types(a.var, nothing => nothing)
end

ext: jl


text:
using Pkg
using Test
using StructDatabaseMapping
@testset "Connections" begin
     @testset "Queue" begin
        include("TestQueue.jl")
    end
    @testset "Pool" begin
        include("TestPool.jl")
    end
end
@testset "Mapper" begin
    include("TestMapper.jl")
end

@testset "Redis" begin
    Pkg.add(PackageSpec(url="https://github.com/JuliaDatabases/Redis.jl.git"))
    include(joinpath(@__DIR__, "redis",  "TestRedis.jl"))

    try
        TestRedis.test()
    finally
        TestRedis.cleanup()
    end
    Pkg.rm("Redis")
end
# I have to do  this because LibPQ and SQlite cannot be installed
# Because there is some incompatibility with Tables.jl
@testset "SQLite" begin
    Pkg.add("SQLite")
    include(joinpath(@__DIR__, "sqlite",  "TestSQLite.jl"))
    try
        TestSQLite.test()
        TestSQLite.test_cleanup()
    finally
        TestSQLite.cleanup()
    end
    Pkg.rm("SQLite")
end
@testset "PostgreSQL" begin
    Pkg.add("LibPQ")
    include(joinpath(@__DIR__, "postgresql",  "TestLibPQ.jl"))
    try
        TestLibPQ.test()
    finally
        TestLibPQ.cleanup()
    end
    Pkg.rm("LibPQ")
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
struct Summary <: Output
  filename::String
  last_infections::Vector{Float32}
  num_infections::Vector{UInt32}
  Summary(filename::AbstractString, num_trajectories::Integer) = new(filename, sizehint!(Vector{Float32}(), num_trajectories), sizehint!(Vector{UInt32}(), num_trajectories))
end

function beforetrajectories(s::Summary, params::Simulation.SimParams)
  sz = Simulation.numindividuals(params)
  sizehint!(s.last_infections, sz)
  sizehint!(s.num_infections, sz)
end

function pushtrajectory!(s::Summary, state::Simulation.SimState, ::Simulation.SimParams, ::DetectionCallback)
  push!(s.last_infections, maximum(Simulation.time, state.forest.inedges))
  push!(s.num_infections, count(Simulation.istransmission, state.forest.inedges))
  nothing
end

function aftertrajectories(s::Summary, ::Simulation.SimParams)
  file = jldopen(s.filename, "w")
  try
    file["num_infections"] = s.num_infections
    file["last_infections"] = s.last_infections
  finally
    close(file)
  end
  nothing
end


text:
using Cassette
macro trace(f,ex)
    Cx = gensym("C")
    Fx = gensym("F")
    quote
        Cassette.@context $Cx
        Cassette.prehook(::$Cx, ::typeof($f), args...) = println(($f,args...))
        Cassette.posthook(::$Cx, output, ::typeof($f), args...) = println(($f,args...,"=>",output))
        $Fx() = $ex
        Cassette.overdub($Cx(), $Fx)
    end
end


text:
function jumpcount(eventcount)
    c = 0
    [x ? c += 1 : 0 for x in eventcount]
end

function ispokeout(pokeID,event)
    event == pokeID*"_out"
end

function findpokesout!(df)
    idx = findall(df.PokeIn)
    df[!,:PokeOut] .= false
    df[!,:PokeOut_count] .= 0
    for (count,i) in enumerate(idx)
        c = df[i,:PokeIn_count]
        pokeID = df[i,:Event]
        find_po = findfirst(e -> ispokeout(pokeID,e),df[i+1:end,:Event])
        if !isnothing(find_po)
            po_idx = find_po + i
        else
            println("Poke out not found for poke $c out of $(length(idx)) pokes at index $i out of $(nrow(df))")
            if count == length(idx)
                println("setting poke out index equal to last event")
                po_idx = nrow(df)
            else
                println("no solution available, deleting poke $c with PokeID = $pokeID")
                df[i,:PokeIn] = false
                df[i,:PokeIn_count] = 0
                continue
            end
        end
        df[po_idx,:PokeOut] = true
        df[po_idx,:PokeOut_count] = c
    end
end

# function correct_pokes!(df)
#     ins = findall(df.PokeIn)
#     outs = findall(df.PokeOut)
#     length(ins) != length(outs) && error("poke_in and poke_out count not matching")
#     wrongs_idx = findall(outs[1:end-1] .> ins[2:end])
#     for i in wrongs_idx
#         i_out = outs[i]
#         df[i_out,:PokeOut] = false
#         i_in = ins[2:end][i]
#         df[i_in,:PokeIn] = false
#     end
# end
#
# """
#     count_pokes(event_vec::AbstractVector)
#
# Using Event column resulting from `table_raw_data(session_path::String)`
# reads the event string to identify the beginning a new poke
# and counts them to use this as a grouping factor for events related to the same
# poke. Crucially it only begin counting from the first poke_in event excluding
# beginning of sessions where not the entire poke was registered: meaning Poke=0
# Poke beginning is identified with regex: a string containing exactly "poke_"
# followed by any single digit (/d) and exactly concluding after it (\$)
# """
# function count_pokes!(prov::AbstractDataFrame)
#     prov[!,:PokeIn] = ismatch.(r"poke_\d$", prov.Event)
#     prov[!,:PokeIn_count] = jumpcount(prov.PokeIn)
#     findpokesout!(prov)
#     correct_pokes!(prov)
#     # pokes_in = [ismatch(r"^poke_\d$",ev) for ev in event_vec]
#     prov[!,:PreCount] = accumulate(+,prov.PokeIn,init = 0)
#     # pre = DataFrame(ev = event_vec, p = accumulate(+,pokes_in,init = 0))
#     transform!(groupby(prov,:PreCount), :Event => adjust_poke_count => :CleanPoke)
#     prov[!,:Poke] = [c ? p : missing for (c,p) in zip(prov.CleanPoke,prov.PreCount)]
#
# end
#
# function adjust_poke_count(ev_vec)
#     rm = .!isnothing.(match.(r"(^poke_\d_out$)|(^out_((left)|(right))$)|(^travel_out_of_poke$)", ev_vec))
#     if any(rm)
#         if findlast(rm) == length(ev_vec)
#             return trues(length(ev_vec))
#         else
#             return vcat(trues(findlast(rm)), falses(length(ev_vec)-findlast(rm)))
#         end
#     else
#         return falses(length(ev_vec))
#     end
# end


text:
## helpers.jl - exported functions used in generated code

function ungetindex!(dx::AbstractArray, x::AbstractArray, ds, i...)
    dx[[i...]] .= ds
    return dx
end

function ungetindex!(dx::AbstractArray, x::AbstractArray, ds, i::AbstractArray{Int})
    dx[i] .= ds
    return dx
end

function ungetindex(x::AbstractArray, ds, i...)
    dx = zero(x)
    return ungetindex!(dx, x, ds, i...)
end

function sum_grad(x::AbstractArray, ds; opts...)
    return ones(size(x)) .* ds
end


text:
include("../src/tenedyn_module.jl")
include("../../__lib__/math/ode/src/ode_module.jl")

module modeldef
using ..tenedyn

function create_model_pi_controller_pt1(
    setpoint_amplitude,
    setpoint_time_delay,
    ctrl_kp,
    ctrl_ki,
    ctrl_integral_init_value,
    plant_T,
    plant_integral_init_value)

    # model definition
    setpoint = tenedyn.Step(setpoint_amplitude, setpoint_time_delay)

    ctrl_minus      = tenedyn.Minus()
    ctrl_prop_gain  = tenedyn.Gain(ctrl_kp)
    ctrl_int_gain   = tenedyn.Gain(ctrl_ki)
    ctrl_plus       = tenedyn.Plus()
    ctrl_integrator = tenedyn.Integrator(ctrl_integral_init_value)

    plant_pt1_minus      = tenedyn.Minus()
    plant_pt1_mult1T     = tenedyn.Gain(1/plant_T)
    plant_pt1_integrator = tenedyn.Integrator(plant_integral_init_value)

    sensor = tenedyn.Identity()

    model = tenedyn.Model()

    # add elements to the model
    push!(model, setpoint,              "ref")                  # 1

    push!(model, ctrl_minus,            "ctrl_minus")           # 2
    push!(model, ctrl_prop_gain,        "ctrl_prop_gain")       # 3
    push!(model, ctrl_int_gain,         "ctrl_int_gain")        # 4
    push!(model, ctrl_plus,             "ctrl_plus")            # 5
    push!(model, ctrl_integrator,       "ctrl_integrator")      # 6

    push!(model, plant_pt1_minus,       "plant_pt1_minus")      # 7
    push!(model, plant_pt1_mult1T,      "plant_pt1_mult1T")     # 8
    push!(model, plant_pt1_integrator,  "plant_pt1_integrator") # 9

    push!(model, sensor,                "sensor")               # 10

    # topology
    tenedyn.connect!(model, setpoint.y,             ctrl_minus.x1)

    tenedyn.connect!(model, ctrl_minus.y,           ctrl_prop_gain.x)
    tenedyn.connect!(model, ctrl_prop_gain.y,       ctrl_plus.x1)

    tenedyn.connect!(model, ctrl_minus.y,           ctrl_int_gain.x)
    tenedyn.connect!(model, ctrl_int_gain.y,        ctrl_integrator.x)
    tenedyn.connect!(model, ctrl_integrator.y,      ctrl_plus.x2)

    tenedyn.connect!(model, ctrl_plus.y,            plant_pt1_minus.x1)
    tenedyn.connect!(model, plant_pt1_integrator.y, plant_pt1_minus.x2)

    tenedyn.connect!(model, plant_pt1_minus.y,      plant_pt1_mult1T.x)
    tenedyn.connect!(model, plant_pt1_mult1T.y,     plant_pt1_integrator.x)

    tenedyn.connect!(model, plant_pt1_integrator.y, sensor.x)
    tenedyn.connect!(model, sensor.y,               ctrl_minus.x2)

    return model
end

function create_model_pi_controller_pt2(
    setpoint_time,
    setpoint_value,
    ctrl_kp,
    ctrl_ki,
    ctrl_integral_init_value,
    plant_T,
    plat_D,
    plant_x_integral_init_value,
    plant_v_integral_init_value)

    # model definition
    # setpoint = tenedyn.Step(setpoint_amplitude, setpoint_time_delay)
    setpoint = tenedyn.Curve(setpoint_time, setpoint_value)

    ctrl_minus      = tenedyn.Minus()
    ctrl_prop_gain  = tenedyn.Gain(ctrl_kp)
    ctrl_int_gain   = tenedyn.Gain(ctrl_ki)
    ctrl_plus       = tenedyn.Plus()
    ctrl_integrator = tenedyn.Integrator(ctrl_integral_init_value)

    plant_pt2_minus      = tenedyn.Minus()
    plant_pt2_plus       = tenedyn.Plus()

    plant_pt2_mult1T2    = tenedyn.Gain(1/plant_T^2)
    plant_pt2_mult2D1T   = tenedyn.Gain(-2 * plat_D / plant_T)

    plant_pt2_v_integrator = tenedyn.Integrator(plant_v_integral_init_value)
    plant_pt2_x_integrator = tenedyn.Integrator(plant_x_integral_init_value)

    sensor = tenedyn.Identity()

    model = tenedyn.Model()

    # add elements to the model
    push!(model, setpoint,                "ref")                      # 1

    push!(model, ctrl_minus,              "ctrl_minus")               # 2
    push!(model, ctrl_prop_gain,          "ctrl_prop_gain")           # 3
    push!(model, ctrl_int_gain,           "ctrl_int_gain")            # 4
    push!(model, ctrl_plus,               "ctrl_plus")                # 5
    push!(model, ctrl_integrator,         "ctrl_integrator")          # 6

    push!(model, plant_pt2_minus,         "plant_pt2_minus")          # 7
    push!(model, plant_pt2_plus,          "plant_pt2_plus")           # 8
    push!(model, plant_pt2_mult1T2,       "plant_pt2_mult1T2")        # 9
    push!(model, plant_pt2_mult2D1T,      "plant_pt2_mult2D1T")       # 10
    push!(model, plant_pt2_v_integrator,  "plant_pt2_v_integrator")   # 11
    push!(model, plant_pt2_x_integrator,  "plant_pt2_x_integrator")   # 12

    push!(model, sensor,                  "sensor")                   # 13

    # topology
    tenedyn.connect!(model, setpoint.y,                 ctrl_minus.x1)

    tenedyn.connect!(model, ctrl_minus.y,               ctrl_prop_gain.x)
    tenedyn.connect!(model, ctrl_prop_gain.y,           ctrl_plus.x1)
    tenedyn.connect!(model, ctrl_minus.y,               ctrl_int_gain.x)
    tenedyn.connect!(model, ctrl_int_gain.y,            ctrl_integrator.x)
    tenedyn.connect!(model, ctrl_integrator.y,          ctrl_plus.x2)

    tenedyn.connect!(model, ctrl_plus.y,                plant_pt2_minus.x1)

    tenedyn.connect!(model, plant_pt2_x_integrator.y,   plant_pt2_minus.x2)
    tenedyn.connect!(model, plant_pt2_minus.y,          plant_pt2_mult1T2.x)
    tenedyn.connect!(model, plant_pt2_mult1T2.y,        plant_pt2_plus.x1)
    tenedyn.connect!(model, plant_pt2_plus.y,           plant_pt2_v_integrator.x)
    tenedyn.connect!(model, plant_pt2_v_integrator.y,   plant_pt2_mult2D1T.x)
    tenedyn.connect!(model, plant_pt2_mult2D1T.y,       plant_pt2_plus.x2)
    tenedyn.connect!(model, plant_pt2_v_integrator.y,   plant_pt2_x_integrator.x)

    tenedyn.connect!(model, plant_pt2_x_integrator.y,   sensor.x)

    tenedyn.connect!(model, sensor.y,                   ctrl_minus.x2)

    return model
end

end

module modelsim
using ..tenedyn
using ..ode
using ..modeldef

function simulate(
    ode_fcn,
    init_cond;
    time_interval = ode.TimeInterval(1.0),
    solver = ode.Tsitouras45(),
    abs_tol = 1e-5,
    rel_tol = 1e-5,
    max_iter = 1000,
    step_size_controller = ode.ClassicalStepSizeControl(),
    max_step_size = 0.1)

    ode_problem = ode.OdeProblem(
        ode_fcn,
        time_interval,
        init_cond)

    basic_options = ode.BasicOdeOptions(
        ode_problem;
        abs_tol = abs_tol,
        rel_tol = rel_tol,
        max_iter = max_iter)

    stepper_options = ode.StepperOptions(;
        max_step_size = max_step_size)

    @time (ode_res, step_stat) = ode.solve(
        ode_problem,
        solver;
        basic_options        = basic_options,
        step_size_controller = step_size_controller,
        stepper_options      = stepper_options) #

    return (ode_res, step_stat)
end

end

module odeplot
using ..ode
using ..modelsim

using PyPlot
PyPlot.pygui(true)

function plot_ode_res(ode_res, indices = (1,))

    t = ode.time(ode_res)
    q = ode.variables(ode_res)

    L = length(indices)

    figh = PyPlot.figure()

    a1 = PyPlot.subplot(L, 1, 1)
    PyPlot.grid()
    PyPlot.plot( t, q[indices[1], :] )

    for ii = 2:L
        PyPlot.subplot(L, 1, ii, sharex = a1)
        PyPlot.grid()
        PyPlot.plot( t, q[indices[ii], :] )
    end

    return figh
end

end

module simpipeline
using ..ode
using ..tenedyn

using ..modeldef
using ..modelsim
using ..odeplot

function create_parametrized_model_pi_controller_pt1()
    # model parameters
    setpoint_amplitude  = 2.0
    setpoint_time_delay = 0.5

    ctrl_kp = 10.0
    ctrl_ki = 18.0
    ctrl_integral_init_value = 0.0

    plant_T = 0.2
    plant_integral_init_value = 0.5

    ode_fcn = tenedyn.create_ode_function_pipeline(
        modeldef.create_model_pi_controller_pt1,
        setpoint_amplitude,
        setpoint_time_delay,
        ctrl_kp,
        ctrl_ki,
        ctrl_integral_init_value,
        plant_T,
        plant_integral_init_value)

    return ode_fcn
end

function run_simulation_pt1(ode_fcn)

    init_cond = [0.0, 0.0]

    (ode_res, step_stat) = modelsim.simulate(
            ode_fcn,
            init_cond;
            time_interval = ode.TimeInterval(5.0),
            abs_tol = 1e-7,
            rel_tol = 1e-7,
            max_iter = 1000)

    @show step_stat

    odeplot.plot_ode_res(ode_res, (1,2))

    return ode_res
end

function create_parametrized_model_pi_controller_pt2()
    # model parameters
    setpoint_amplitude  = 2.0
    setpoint_time_delay = 0.5

    setpoint_time  = [0.0, 0.5, 0.501, 1.0, 3.0, 3.01, 4.0]
    setpoint_value = [0.0, 0.0, 2.0,   2.0, 5.0, 0.0,  0.0]

    ctrl_kp = 3.0
    ctrl_ki = 15.0
    ctrl_integral_init_value = 0.0

    plant_T = 0.05
    plant_D = 0.35

    plant_integral_x_init_value = 0.0
    plant_integral_y_init_value = 0.0

    ode_fcn = tenedyn.create_ode_function_pipeline(
        modeldef.create_model_pi_controller_pt2,
        setpoint_time,
        setpoint_value,
        ctrl_kp,
        ctrl_ki,
        ctrl_integral_init_value,
        plant_T,
        plant_D,
        plant_integral_x_init_value,
        plant_integral_y_init_value)

    return ode_fcn
end

function run_simulation_pt2(ode_fcn)

    init_cond = [0.0, 0.0, 1.0]

    (ode_res, step_stat) = modelsim.simulate(
            ode_fcn,
            init_cond;
            time_interval = ode.TimeInterval(15.0),
            abs_tol = 1e-6,
            rel_tol = 1e-6,
            max_iter = 10000)

    @show step_stat

    odeplot.plot_ode_res(ode_res, (1, 2, 3))

    return ode_res
end

# ode_fcn = create_parametrized_model_pi_controller_pt1()
# ode_res = run_simulation_pt1(ode_fcn)

ode_fcn = create_parametrized_model_pi_controller_pt2()
ode_res = run_simulation_pt2(ode_fcn)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

################## Load Operator ##################
if Sys.islinux()
py"""
import tensorflow as tf
libSparseLeastSquare = tf.load_op_library('build/libSparseLeastSquare.so')
@tf.custom_gradient
def sparse_least_square(ii,jj,vv,ff,n):
    u = libSparseLeastSquare.sparse_least_square(ii,jj,vv,ff,n)
    def grad(dy):
        return libSparseLeastSquare.sparse_least_square_grad(dy, u, ii,jj,vv,ff,n)
    return u, grad
"""
elseif Sys.isapple()
py"""
import tensorflow as tf
libSparseLeastSquare = tf.load_op_library('build/libSparseLeastSquare.dylib')
@tf.custom_gradient
def sparse_least_square(ii,jj,vv,ff,n):
    u = libSparseLeastSquare.sparse_least_square(ii,jj,vv,ff,n)
    def grad(dy):
        return libSparseLeastSquare.sparse_least_square_grad(dy, u, ii,jj,vv,ff,n)
    return u, grad
"""
elseif Sys.iswindows()
py"""
import tensorflow as tf
libSparseLeastSquare = tf.load_op_library('build/libSparseLeastSquare.dll')
@tf.custom_gradient
def sparse_least_square(ii,jj,vv,ff,n):
    u = libSparseLeastSquare.sparse_least_square(ii,jj,vv,ff,n)
    def grad(dy):
        return libSparseLeastSquare.sparse_least_square_grad(dy, u, ii,jj,vv,ff,n)
    return u, grad
"""
end

sparse_least_square = py"sparse_least_square"
################## End Load Operator ##################

# TODO: specify your input parameters
ii = Int32[1;1;2;2;3;3]
jj = Int32[1;2;1;2;1;2]
vv = Float64[1;2;3;4;5;6]
ff = Float64[1;1;1]
A = Float64[1 2;3 4;5 6]; f = Float64[1;1;1]
sol = A\f
n = 2
# error()
u = sparse_least_square(ii,jj,vv,ff,constant(n, dtype=Int32))
sess = Session()
init(sess)
@show run(sess, u)-sol
# error()

# TODO: change your test parameter to `m`
# gradient check -- v
function scalar_function(m)
    return sum(tanh(sparse_least_square(ii,jj,m,ff,n)))
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(vv)
v_ = rand(6)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session()
init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
printstyled("  rescale\n", color=:light_green)
g = 32.0
g1 = 7.0
g2 = 11.0
segy_nodal = string(path, "/SampleFiles/SEGY/FORGE_78-32_iDASv3-P11_UTC190428135038.sgy")
U = read_data("segy", segy_nodal, ll=0x01)
@test U.n == 33
S = U[1:7]
fill!(S.gain, 1.0)
C = deepcopy(S[1])
U = deepcopy(S)
D = deepcopy(C)

printstyled("    rescale!(S)\n", color=:light_green)
S = deepcopy(U)
S.gain[1] = g1
S2 = rescale(S)
rescale!(S)
@test S == S2
G = S.gain
@test all([getindex(G, i) == g1 for i in 1:S.n])
for i in 2:S.n
  @test isapprox(S.x[i]/g1, U.x[i])
end

printstyled("    rescale!(S, g)\n", color=:light_green)
S = deepcopy(U)
S2 = rescale(S, g)
rescale!(S, g)
@test S == S2
G = S.gain
X = S.x
@test all([getindex(G, i) == g for i in 1:S.n])
@test all([isapprox(X[i], U.x[i].*g) for i in 1:S.n])

printstyled("      channel range\n", color=:light_green)
S = deepcopy(U)
chans = collect(2:2:6)
rest  = setdiff(1:S.n, chans)
S2 = rescale(S, g, chans=chans)
rescale!(S, g, chans=chans)
@test S == S2
G = S.gain
X = S.x
@test all([getindex(G, i) == g for i in chans])
@test any([getindex(G, i) == g for i in rest]) == false
@test all([isapprox(X[i], U.x[i].*g) for i in chans])
@test any([isapprox(X[i], U.x[i].*g) for i in rest]) == false

printstyled("    rescale!(S_targ, S_src)\n", color=:light_green)
S = deepcopy(U)
fill!(U.gain, g2)
S2 = rescale(S, U)
rescale!(S, U)
@test S == S2
G = S.gain
X = S.x
@test all([getindex(G, i) == g2 for i in 1:S.n])
@test all([isapprox(X[i], U.x[i].*g2) for i in 1:S.n])

printstyled("    rescale!(C, g)\n", color=:light_green)
B = rescale(C, g)
rescale!(C, g)
@test B == C
@test C.gain == g
@test isapprox((C.x).*g, D.x)

printstyled("    rescale!(C, D)\n", color=:light_green)
C = deepcopy(D)
D.gain = g2
B = rescale(C, D)
rescale!(C, D)
@test B == C
@test isapprox(C.x, (D.x).*g2)
@test eltype(C.x) == Float32


text:
function from_default_theme(attr)
    theme = default_styles()
    return get(theme, attr) do
        AbstractPlotting.current_default_theme()[attr]
    end
end

line_element(;
             color=from_default_theme(:color),
             linestyle=from_default_theme(:linestyle),
             linewidth=from_default_theme(:linewidth),
             kwargs...) =
    LineElement(; color, linestyle, linewidth, kwargs...)

marker_element(;
               color=from_default_theme(:color),
               marker=from_default_theme(:marker),
               strokecolor=from_default_theme(:strokecolor),
               markerpoints=[Point2f0(0.5, 0.5)],
               kwargs...) =
    MarkerElement( ; color, marker, strokecolor, markerpoints, kwargs...)

poly_element(;
             color=from_default_theme(:color),
             strokecolor=:transparent,
             kwargs...) =
    PolyElement(; color, strokecolor, kwargs...)

legend_elements(::Type{Scatter}; kwargs...) = [marker_element(; kwargs...)]
legend_elements(::Type{Lines}; kwargs...) = [line_element(; kwargs...)]
legend_elements(::Type{Contour}; kwargs...) = [line_element(; kwargs...)]

function legend_elements(::Type{LinesFill}; color=from_default_theme(:color), fillalpha=0.15, kwargs...)
	meshcolor = to_color((color, fillalpha))
	return [poly_element(; color=meshcolor, kwargs...), line_element(; color, kwargs...)]
end

legend_elements(::Any; linewidth=0, strokecolor=:transparent, kwargs...) =
	[poly_element(; linewidth, kwargs...)]


text:
"""
    OneSidedBook{Sz,Px,Oid,Aid}

One-Sided book with order-id type Oid, account-id type Aid,
size type Sz and price type Px.

OneSidedBook is a one-sided book (i.e. :BID or :ASK) of order queues at
varying prices.

OrderQueues are stored in an AVLTree (.book) indexed
either by price (ASK side) or -price (BID side)

The book keeps track of various statistics such as the current best price,
total share and price volume, as well as total contained number of orders.

"""
@kwdef mutable struct OneSidedBook{Sz<:Real,Px<:Real,Oid<:Integer,Aid<:Integer}
    is_bid_side::Bool
    book::AVLTree{Px,OrderQueue{Sz,Px,Oid,Aid}} = AVLTree{Px,OrderQueue{Sz,Px,Oid,Aid}}()
    total_volume::Sz = 0 # Total volume available in shares
    total_volume_funds::Float64 = 0.0 # Total volume available in underlying currency
    num_orders::Int32 = Int32(0) # Number of orders in the book
    best_price::Union{Px,Nothing} = nothing # best bid or ask
end

isbidbook(sb::OneSidedBook) = sb.is_bid_side
isaskbook(sb::OneSidedBook) = !sb.is_bid_side
Base.isempty(sb::OneSidedBook) = isempty(sb.book)

"Compute total volume available below limit price"
function _size_available(sb::OneSidedBook{Sz,Px,Oid,Aid},limit_price::Px) where {Sz,Px,Oid,Aid}
    t = Sz(0)
    if isbidbook(sb)
        for q in sb.queue
            (q.price >= limit_price) ? (t += q.total_volume) : break
        end
    else
        for q in sb.queue
            (q.price <= limit_price) ? (t += q.total_volume) : break
        end
    end
    return t
end

"Compute currency volume available below limit price"
function _funds_available(sb::OneSidedBook{Sz,Px,Oid,Aid},limit_price::Px) where {Sz,Px,Oid,Aid}
    t = Int64(0.0)
    if isbidbook(sb)
        for q in sb.queue
            (q.price >= limit_price) ? (t += Int64(q.total_volume*q.price)) : break
        end
    else
        for q in sb.queue
            (q.price <= limit_price) ? (t += Int64(q.total_volume*q.price)) : break
        end
    end
    return t

end

_size_available(sb::OneSidedBook,::Nothing) = sb.total_volume
_funds_available(sb::OneSidedBook,::Nothing) = sb.total_volume_funds

"Updates the latest best price in a Sidebook (either :BID or :ASK book)."
function _update_next_best_price!(sb::OneSidedBook)
    sb.best_price = isempty(sb.book) ? nothing : abs(first(sb.book)[1])
    return nothing
end

"Retrieve order queue from OneSidedBook at given price"
function _get_price_queue(
    sb::OneSidedBook{<:Real,Px,<:Integer,<:Integer}, price::Px
) where {Px}
    pricekey = isaskbook(sb) ? price : -price
    return AVLTrees.findkey(sb.book, pricekey) # Return the price queue
end

"Delete entire queue associated with given price from OneSidedBook and track stats"
function _popat_queue!(
    sb::OneSidedBook{Sz,Px,Oid,Aid}, price::Px
) where {Sz,Px,Oid,Aid}
    pricekey = isaskbook(sb) ? price : -price
    price_queue = pop!(sb.book, pricekey) # delete price queue
    # update book stats
    (price_queue.price == sb.best_price) && _update_next_best_price!(sb) # Update price only if best price was changed
    sb.num_orders -= price_queue.num_orders[]
    sb.total_volume -= price_queue.total_volume[]
    sb.total_volume_funds -= Float64(price_queue.total_volume[]*price_queue.price)
    return price_queue
end

"Pop first queue and track stats"
function _popfirst_queue!(
    sb::OneSidedBook{Sz,Px,Oid,Aid}
) where {Sz,Px,Oid,Aid}
    price_queue = popfirst!(sb.book) # pop price queue at best price
    # update book stats
    _update_next_best_price!(sb) # Update best price
    sb.num_orders -= price_queue.num_orders[]
    sb.total_volume -= price_queue.total_volume[]
    sb.total_volume_funds -= Float64(price_queue.total_volume[]*price_queue.price)
    return price_queue
end

"Insert entire price queue into OneSidedBook and track stats"
function _insert_queue!(
    sb::OneSidedBook{Sz,Px,Oid,Aid},
    price_queue::OrderQueue
) where {Sz,Px,Oid,Aid}
    pricekey = isaskbook(sb) ? price_queue.price : -price_queue.price
    insert!(sb.book,pricekey,price_queue)
    # update book stats
    (isnothing(sb.best_price) || (price_queue.price === sb.best_price)) && _update_next_best_price!(sb) # Update price only if best price was changed
    sb.num_orders += price_queue.num_orders[]
    sb.total_volume += price_queue.total_volume[]
    sb.total_volume_funds += Float64(price_queue.total_volume[]*price_queue.price)
    return nothing
end

"Insert new_order into OneSidedBook at given price, create new price queue if needed"
function insert_order!(
    sb::OneSidedBook{Sz,Px,Oid,Aid}, new_order::Order{Sz,Px,Oid,Aid}
) where {Sz,Px,Oid,Aid}
    pricekey = isaskbook(sb) ? new_order.price : -new_order.price
    # search for order queue at price
    order_queue = AVLTrees.findkey(sb.book, pricekey)
    if isnothing(order_queue) # If key not present (price doesnt exist in book)
        new_queue = OrderQueue{Sz,Px,Oid,Aid}(new_order.price) # Create new price queue
        push!(new_queue, new_order) # Add order to new price queue
        insert!(sb.book, pricekey, new_queue) # add new price queue to OneSidedBook
        # Update new best price depending on bid/ask
        if isnothing(sb.best_price)
            sb.best_price = new_order.price
        elseif isbidbook(sb)
            sb.best_price = max(new_order.price, sb.best_price)
        else
            sb.best_price = min(new_order.price, sb.best_price)
        end
    else # If order queue present, retrieve queue and insert new order
        push!(order_queue, new_order)
    end

    # Update Onesidedbook info
    sb.num_orders += 1
    sb.total_volume += new_order.size
    sb.total_volume_funds += Float64(new_order.price*new_order.size)
    return nothing
end

# "Bulk insert orders into `OneSidedBook`"
# function insert_order!(
#     sb::OneSidedBook{Sz,Px,Oid,Aid}, new_orders::AbstractVector{Order{Sz,Px,Oid,Aid}}
# ) where {Sz,Px,Oid,Aid}
#     for ord in new_orders
#         insert_order!(sb, ord)
#     end
#     return nothing
# end

"Delete order with given price/tick_id from book"
function pop_order!(
    sb::OneSidedBook{Sz,Px,Oid,Aid}, price::Px, orderid::Oid
) where {Oid,Aid,Sz<:Real,Px<:Real}
    # Get price queue and delete order from it
    order_queue = _get_price_queue(sb, price)
    if !isnothing(order_queue)
        Δvolm = order_queue.total_volume[] # get stats before deletion
        ord = popat_orderid!(order_queue, orderid)
        Δvolm -= order_queue.total_volume[] # get stats after deletion

        # If order deletion depleted price queue, delete the whole queue
        if isempty(order_queue)
            _popat_queue!(sb, price) # note: this function will update price
        end

        # Update Onesidedbook info
        sb.num_orders -= 1
        sb.total_volume -= Δvolm
        sb.total_volume_funds -= Float64(Δvolm*price)

        return ord # return popped order, is nothing if no order found
    else
        return nothing
    end
end


text:
using LinearAlgebra

# wall constraints
function c_wall(x, wall)
    return dot( x - wall.pos, wall.dir )
end
# friction constraint
function c_wall_friction(x_old,x,wall)
    return dot( x_old - x, [0 -1; 1 0] * wall.dir)
end

max_stars_repo_name: SteffenPL/ContactAlgorithms

max_issues_repo_name: SteffenPL/ContactAlgorithms


text:
# fin_simulate.jl
# Function to run the provider simulation.
#
# Inputs:
# day_of_entry: array containing, for each patient, when the patient is admitted.
# rday_of_surgery_completion: array containing, for each patient, how many days after the day_of_entry did the patient complete surgery.
# rday_of_discharge: array containing, for each patient, how many days after the day_of_entry did the patient get discharged.
# predict_0: array containing, for each patient, the risk prediction (probability between 0 and 1) for the patient available on the day of surgery completion.
# predict_1: array containing, for each patient, the risk prediction (probability between 0 and 1) for the patient available 1 day after the day of surgery completion.
# predict_2: array containing, for each patient, the risk prediction (probability between 0 and 1) for the patient available 2 days or more after the day of surgery completion.
# endpoint: array containing, for each patient, a 0/1 to indicate whether the patient is readmitted or not.
# HCUP_cost: array containing, for each patient, the cost of that patient's readmission based on HCUP data from Bailey et al. (see paper for more details).
# dayofweek: array containing, for each day in the simulation, the day of the week, numbered from 1 (Monday) to 7 (Sunday).
# day_list: array containing the days in the simulation. (This array must be formatted as n, n+1, ..., n + k, where n is an integer and k is the duration of the simulation. For example,
# day_list could be [45, 46, 47, 48, 49, 50, 51], so 45 would be the first day in the simulation.)
# visit_days: array containing a subset of the values [1,2,3,4,5,6,7], indicating the provider's schedule (e.g., visit_days = [1,3] corresponds to a Monday/Wednesday schedule.)
# capacity_per_day: integer that represents the maximum number of patients the provider can see on any day in the schedule encoded by visit_days.
# atdischarge: boolean, indicating whether the provider can only see patients who are being discharged on the current day or not.
#
# Outputs: (see paper for definitions)
# PS: patients seen metric
# RA: readmissions anticipated
# ERC: expected readmission cost savings
# This function takes the above data and runs the provider simulation.

function fin_simulate(day_of_entry, rday_of_surgery_completion, rday_of_discharge, predict_0, predict_1, predict_2, endpoint, HCUP_cost, dayofweek, day_list,
						visit_days, capacity_per_day, atdischarge)

	num_patients = length(endpoint) # Determine number of patients
	is_seen = zeros(Bool, num_patients) # For each patient, track whether they have been seen by the provider yet, or not.
	num_days = length(day_list) # Determine the number of days in the simulation from day_list.

	# Iterate through each day in the simulation
	for i in 1:num_days
		day = day_list[i] # Retrieve the day ID

		if (dayofweek[i] in visit_days) # Check if the day of the week of this day falls in the provider's schedule

			# Has the patient completed surgery by this day?
			cond1 = (day_of_entry + rday_of_surgery_completion .- 1) .<= day
			if (atdischarge)
				# If atdischarge is true, has the patient been discharged yet?
				cond1 = day .>= (day_of_entry + rday_of_discharge .- 1)
			end

			# Is this day on or before the patient's day of discharge?
			cond2 = day .<= (day_of_entry + rday_of_discharge .- 1)

			# Has the patient not yet been selected by the provider?
			cond3 = .!is_seen

			# Find all patients satisfying these three criteria; these are the eligible patients.
			eligible_patients = findall(cond1 .& cond2 .& cond3)

			if (!isempty(eligible_patients))

				# For these eligible patients, compute what their current risk prediction is.
				# The conditional in the loop below ensures that we use the most recent risk prediction.
				risk = zeros(length(eligible_patients))

				for j in 1:length(eligible_patients)
					p = eligible_patients[j]
					if (day == (day_of_entry[p] + rday_of_surgery_completion[p] - 1) )
						risk[j] = predict_0[p]
					elseif (day == (day_of_entry[p] + rday_of_surgery_completion[p] - 1 + 1))
						risk[j] = predict_1[p]
					elseif (day >= (day_of_entry[p] + rday_of_surgery_completion[p] - 1 + 2) )
						risk[j] = predict_2[p]
					else
						error("Weird issue with day; check!")
					end
				end

				# Order the patients from highest to lowest risk.
				priority_order = sortperm(risk, rev = true )

				# Determine how many patients will be seen; this is the lower of the capacity_per_day,
				# or however many patients are eligible
				min_capacity_available = min(length(eligible_patients), capacity_per_day)

				# Flag the top patients
				for j in 1:min_capacity_available
					p = eligible_patients[ priority_order[j] ]
					is_seen[p] = true
				end
			end
		end
	end

	# Add up how many patients have been flagged
	PS = sum(is_seen)

	# Add up the 0/1 binary variable for the endpoint, in order to obtain RA
	RA = sum(endpoint[ is_seen ])

	# Add up the HCUP cost for those patients that go on to be readmitted, and
	# our provider had selected them.
	ERC = sum(  HCUP_cost[ (endpoint .== 1) .& is_seen ])

	return PS, RA, ERC
end

max_stars_repo_name: vvmisic/finsim-code

max_issues_repo_name: vvmisic/finsim-code

max_forks_repo_name: vvmisic/finsim-code


text:
using Cliffords, Formatting

"""
Clifford circuit randomization and compilation with recovery operation.

Since the recovery operation is included, there is no mask to output
(it essentially corresponds to the recovery operation.)
"""
function compare(pseqs,rseqs)

    maxlen = maximum(map(length,pseqs))
    rmaxlen = maximum(map(length,rseqs))

    seq   = Vector{Clifford}(maxlen)
    rseq   = Vector{Clifford}(rmaxlen)

    # now on to experiments with non-zero length (that is why we start at index 2)
    for n in 1:length(pseqs)
        # accumulate sequences
        inter = 1 # identity index
        for i in 1:length(pseqs[n])
            inter   = multiply(pseqs[n][i],inter)
        end
        rinter = 1
        for i in 1:length(rseqs[n])
            rinter  = multiply(rseqs[n][i],rinter)
        end
        # check that sequences are same up to final random Z90 rot
        index = multiply(inter,Cliffords.localcliffordindex(inv(localclifford(rinter))))
        if !(index in [1,9])
            error("Mismatch in line $(n)")
        else
            index == 1 ? (@printf "I") : (@printf "Z")
        end

    end
end

begin
    global multiply

    table = Dict{Tuple{Int,Int},Int}()

    function multiply(c1,c2)
        c1_ = round(Int,c1==0 ? 1 : c1)
        c2_ = round(Int,c1==0 ? 1 : c2)
        if !((c1_,c2_) in keys(table))
            table[(c1_,c2_)] = Cliffords.localcliffordindex(localclifford(c1_)*localclifford(c2_))
        end
        return table[(c1_,c2_)]
    end
end

# begin
#     println("Comparing gst-1q-1024.csv and $(ARGS[1])")
#
#     seqs = readlines(open("gst-1q-1024.csv"))
#     pseqs = [map(x->parse(Int,x), split(s,',')) for s in seqs];
#     seqs = readlines(open(ARGS[1]))
#     rseqs = [map(x->parse(Int,x), split(s,',')) for s in seqs];
#
#     compare(pseqs,rseqs)
#     @printf "\n"
#     println("Done!")
# end

begin
    infile = open(ARGS[1])
    seqs = readlines(infile)
    close(infile)

    #range = parse(Int,ARGS[2])
    #count = parse(Int,ARGS[3])

    pseqs = [map(x->parse(Int,x), split(s,',')) for s in seqs];

    #L = ceil(Int,log(10,range)+1)
    for k in 1:10
        # i = rand(1:range)
        label = format("{1:02d}",k)
        println("Comparing $(ARGS[1]) and $(replace(ARGS[1],r".csv$" => "-r$(label).csv"))")
        rfile = open(replace(ARGS[1],r".csv$" => "-r$(label).csv"))
        rseqs = [map(x->parse(Int,x), split(s,',')) for s in readlines(rfile)];
        compare(pseqs, rseqs)
        @printf "\n"
        close(rfile)
    end
end

max_stars_repo_name: marcusps/RGST-analysis

max_issues_repo_name: marcusps/RGST-analysis

max_forks_repo_name: marcusps/RGST-analysis


text:
# MicroscopyLabels isn't published yet so we need to add the source directly
push!(LOAD_PATH, normpath(joinpath(@__DIR__, "..", "src")))

using Documenter, MicroscopyLabels

makedocs(sitename="MicroscopyLabels.jl")

deploydocs(
    repo = "github.com/tlnagy/MicroscopyLabels.jl.git",
)


text:
import Base.==
import Base.convert

export Goal, PackedGoal, BasicPackedGoal, InterleavedPackedGoal, print_goal, compose, convert, ==, read_plu

# A goal describes a boolean function with one or more inputs and one or more
# outputs. It defines the number of inputs and the values of each output for
# all possible input combinations.
immutable Goal{N}
    num_inputs::Integer
    truth_table::NTuple{N, BitString}
end

Goal{N}(n::Integer, t::NTuple{N, Integer}) = Goal(n, convert(NTuple{N, BitString}, t))
Goal(n::Integer, b::Vector{BitString}) = Goal(n, ntuple(i -> b[i], length(b)))
Goal(c::Chromosome, b::Vector{BitString}) = Goal(c.num_inputs, b)

# The equality operator for Goals is true if the number of inputs, number of
# outputs, and the output truth tables are identical and false otherwise.
=={N}(g::Goal{N}, h::Goal{N}) = (g.num_inputs == h.num_inputs) && (g.truth_table == h.truth_table)
=={G, H}(g::Goal{G}, h::Goal{H}) = false

# Prints a goal in hex format. The truth table entries are printed from highest
# to lowest index. In other words: `print_goal(Goal(1, (0x8, 0x6)))` will
# print "(1, 2, [0x6, 0x8])".
function print_goal{N}(g::Goal{N})
    @printf("(%d, %d, [", g.num_inputs, N)
    for i in N:-1:2
        @printf("%#x ",g.truth_table[i])
    end
    @printf("%#x])\n",g.truth_table[1])
end

# Goals can have multiple outputs.
# There are three formats for goals:
# unpacked:  each output is stored in a separate BitString (this is the output format of execute chromosome)
# packed (also called non-interleved):  All outputs are stored in one BitString.
#      All bits of an output are consecutive.
# interleaved:  All outputs are stored in one BitString.
#      The bits of outputs are interleaved.  This format is used by the g_compose f function.
# Packed goals combine all outputs into a single bitstring and, consequently,
# can only handle functions with an appropriately small number of outputs.
abstract PackedGoal

immutable BasicPackedGoal <: PackedGoal
    num_inputs::Integer
    num_outputs::Integer
    truth_table::BitString
end

immutable InterleavedPackedGoal <: PackedGoal
    num_inputs::Integer
    num_outputs::Integer
    truth_table::BitString
end

function =={T <: PackedGoal}(g::T, h::T)
    (g.num_inputs == h.num_inputs) &&
    (g.num_outputs == h.num_outputs) &&
    (g.truth_table == h.truth_table)
end
==(g::BasicPackedGoal, h::InterleavedPackedGoal) = error("Cannot compare basic and interleaved goals")
==(g::InterleavedPackedGoal, h::BasicPackedGoal) = h == g
=={T <: PackedGoal}(g::Goal, h::T) = error("Cannot compare unpacked and packed goals")
=={T <: PackedGoal}(g::T, h::Goal) = h == g

# Prints a packed goal in hex by converting it to an unpacked goal.
function print_goal{T <: PackedGoal}(g::T)
    print_goal(convert(Goal, g))
end

function convert(::Type{Goal}, g::BasicPackedGoal)
    temp_ttable = g.truth_table
    ttable = Array(BitString, g.num_outputs)
    mask = output_mask(g.num_inputs)
    for i = 1:g.num_outputs
        ttable[i] = mask & temp_ttable
        temp_ttable >>= (2 ^ g.num_inputs)
    end
    return Goal(g.num_inputs, ttable)
end

function convert{N}(::Type{BasicPackedGoal}, g::Goal{N})
    if (2 ^ g.num_inputs) * N > 8*sizeof(BitString)  # 8 is the number of bits in a byte
        error("The goal is too large to be packed.")
    end

    ttable = g.truth_table[N]
    for i = (N-1):-1:1
        ttable <<= (2 ^ g.num_inputs)
        ttable $= g.truth_table[i]
    end
    return BasicPackedGoal(g.num_inputs, N, ttable)
end

function convert(::Type{InterleavedPackedGoal}, g::BasicPackedGoal)
    ttable = convert(BitString, 0)
    one = convert(BitString, 1)
    p = 2 ^ g.num_inputs
    for j = (p-1):-1:0
        for k = (g.num_outputs-1):-1:0
            r = g.truth_table >> (k * (p - 1) + j + k)
            ttable = (ttable << 1) $ (r & one)
        end
    end
    return InterleavedPackedGoal(g.num_inputs, g.num_outputs, ttable)
end

function convert(::Type{BasicPackedGoal}, g::InterleavedPackedGoal)
    ttable = convert(BitString, 0)
    one = convert(BitString, 1)
    t_g = g.truth_table
    p = 2 ^ g.num_inputs
    for j = (g.num_outputs-1):-1:0
        for k = p-1:-1:0
            r = t_g >> (k * g.num_outputs + j)
            ttable = (ttable << 1) $ (r & one)
        end
    end
    return BasicPackedGoal(g.num_inputs, g.num_outputs, ttable)
end

function convert(::Type{Goal}, g::InterleavedPackedGoal)
    return convert(Goal, convert(BasicPackedGoal, g))
end

function convert{N}(::Type{InterleavedPackedGoal}, g::Goal{N})
    return convert(InterleavedPackedGoal, convert(BasicPackedGoal, g))
end

# Compose two goals g, and h. The number of outputs from h must be
# equal to the number of inputs to g.
function compose(g::InterleavedPackedGoal, h::InterleavedPackedGoal)
    if g.num_inputs != h.num_outputs
        error("Cannot compose g and h, input-output mismatch.")
    end

    ones = typemax(BitString)

    mask_g = (ones << g.num_outputs) $ ones
    mask_h = (ones << h.num_outputs) $ ones

    ttable = convert(BitString, 0)

    t_g = g.truth_table
    t_h = h.truth_table

    for i = 0:(2 ^ h.num_outputs - 1)
        g_i = t_g & mask_g
        h_g_i = (t_h >> (h.num_outputs * g_i)) & mask_h
        ttable = ttable $ h_g_i << (i * h.num_outputs)
        t_g = t_g >> g.num_outputs
    end
    return InterleavedPackedGoal(g.num_inputs, h.num_outputs, ttable)
end

# Reads the ".plu" file whose name is fname.
# "*.plu" files define goals in Julian Miller's version 1.1 of GGP.
# This version is limited to goals with at most 4 inputs and at most 4 outputs.
function read_plu(fname)
    num_inputs = 0
    num_outputs = 0
    outputs = BitString[]
    f = open(fname,"r")
    for line in eachline(f)
        fields = split(line,' ')
        if fields[1][1] == '.'
            if fields[1][2] == 'i'
                num_inputs = parse(Int, fields[2])
                if num_inputs > 4
                    error("the number of inputs in this version of read_plu is limited to 4")
                end
            elseif fields[1][2] == 'o'
                num_outputs = parse(Int, fields[2])
                if num_outputs > 4
                    error("the number of outputs in this version of read_plu is limited to 4")
                end
            end
        else
            input_flag = true
            for f in fields
                if length(f) == 0
                    input_flag = false
                    #println()
                elseif f != "\r\n"
                    #@printf(" %#X",int(f))
                    if !input_flag
                        push!(outputs,convert(BitString, parse(Int, f)))
                    end
                end
            end
            #println()
        end
    end
    word_size = 2^num_inputs
    truth_table = convert(BitString,0)
    for o in outputs
        #@printf(" %#X",o)
        truth_table = truth_table << word_size
        truth_table = truth_table $ o
    end
    return BasicPackedGoal(num_inputs, num_outputs, truth_table)
end


text:
@doc Markdown.doc"""
    norm(m::T, I::NfOrdIdl) where T <: Map{AnticNumberField, AnticNumberField} -> NfOrdIdl

Given an embedding $m:k\to K$ of number fields and an integral ideal in $K$, find the norm
$N_{K/k}(I)$.
"""
function norm(m::T, I::NfOrdIdl; order = maximal_order(domain(m))) where T <: Map{AnticNumberField, AnticNumberField}
  K = codomain(m)
  @assert K == nf(Hecke.order(I))
  k = domain(m)
  zk = order
  @assert nf(zk) === k
  if degree(k) == 1
    return ideal(zk, norm(I))
  end
  if I.is_principal == 1
    if isdefined(I, :princ_gen)
      return ideal(zk, zk(norm(m, (I.princ_gen).elem_in_nf)))
    elseif isdefined(J,:princ_gen_special)
      el = J.princ_gen_special[2] + J.princ_gen_special[3]
      return ideal(zk, zk(norm(m, el)))
    end
  end
  assure_2_normal(I)
  J = ideal(zk, I.gen_one^div(degree(K), degree(k)), zk(norm(m, I.gen_two.elem_in_nf)))
  J.gens_normal = I.gens_normal
  return J
end

function norm(m::T, I::NfOrdFracIdl) where T <: Map{AnticNumberField, AnticNumberField}
  return norm(m, numerator(I))//denominator(I)^div(degree(codomain(m)), degree(domain(m)))
end

#TODO: intersect_nonindex uses a worse algo in a more special case. Combine.
#  for prime ideals, the gcd's can be done in F_p/ F_q hence might be faster
@doc Markdown.doc"""
    minimum(m::T, I::NfOrdIdl) where T <: Map{AnticNumberField, AnticNumberField} -> NfOrdIdl

Given an embedding $m:k\to K$ of number fields and an integral ideal in $K$, find the
intersect $I \cap \mathbf{Z}_k$.
"""
function minimum(m::T, I::NfOrdIdl) where T <: Map{AnticNumberField, AnticNumberField}
  K = codomain(m)
  @assert K == nf(order(I))
  k = domain(m)
  zk = maximal_order(k)
  if degree(k) == 1
    return ideal(zk, minimum(I))
  end
  assure_2_normal(I) # basically implies order(I) is maximal

  if !isone(gcd(minimum(I), index(order(I))))
    bk = map(m, basis(maximal_order(k), k))
    bK = map(K, basis(I))
    d = lcm(lcm(map(denominator, bk)), lcm(map(denominator, bK)))
    F = FreeModule(FlintZZ, degree(K))

    hsk = ModuleHomomorphism(FreeModule(FlintZZ, degree(k)), F, elem_type(F)[F(matrix(FlintZZ, 1, degree(K), coefficients(d*x))) for x = bk])
    hsK = ModuleHomomorphism(F, F, elem_type(F)[F(matrix(FlintZZ, 1, degree(K), coefficients(d*x))) for x = bK])
    sk = image(hsk)
    imhsK = image(hsK)
    mm = intersect(sk[1], imhsK[1])
    return ideal(zk, elem_type(zk)[zk(x.v) for x = map(x->preimage(hsk, sk[2](mm[2](x))), gens(mm[1]))])
  end

  # Here is a version with abelian groups
  #if !isone(gcd(minimum(I), index(order(I))))
  #  bk = map(m, basis(maximal_order(k), k))
  #  bK = map(K, basis(I))
  #  d = lcm(lcm(map(denominator, bk)), lcm(map(denominator, bK)))
  #  F = abelian_group([0 for i in 1:degree(K)])
  #  hsk = hom(abelian_group[0 for i in 1:degree(k)], F, elem_type(F)[F(matrix(FlintZZ, 1, degree(K), coefficients(d*x))) for x = bk])
  #  hsK = hom(F, F, elem_type(F)[F(matrix(FlintZZ, 1, degree(K), coefficients(d*x))) for x = bK])
  #  sk = image(hsk)
  #  imhsK = image(hsK)
  #  mm = intersect(sk[1], imhsK[1], false)
  #  return ideal(zk, elem_type(zk)[zk(collect(x.coeff)) for x = map(x->preimage(hsk, sk[2](mm[2](x))), gens(mm[1]))])
  #end

  @assert K == nf(order(I))
  k = domain(m)
  kt, t = PolynomialRing(k, cached = false)
  Qt = parent(K.pol)
  h = gcd(gen(k) - evaluate(Qt(m(gen(k))), t), evaluate(K.pol, t))
  g, ai, _ = gcdx(evaluate(Qt(I.gen_two.elem_in_nf), t) % h, h)
  @assert g == 1
  #so ai * a = 1 in K/k
  c = content_ideal(ai, zk)
  n,dd = integral_split(c)
  J = ideal(zk, I.gen_one) + dd
  J.gens_normal = I.gens_normal
  return J
end

function minimum(m::T, I::NfOrdFracIdl) where T <: Map{AnticNumberField, AnticNumberField}
  return minimum(m, numerator(I))//denominator(I)
end

################################################################################
#
#  Prime decomposition
#
################################################################################

@doc Markdown.doc"""
    intersect_prime(f::Map, P::NfOrdIdl, O_k::NfOrd) -> NfOrdIdl

Given a prime ideal $P$ in $K$ and the inclusion map $f:k \to K$
of number fields, find the unique prime $p$ in $k$ below.
$p$ will be in the order $O_k$ which defaults to "the" maximal order of $k$.
"""
function intersect_prime(f::Map, P::NfOrdIdl, Ok::NfOrd = maximal_order(domain(f)))
  @assert is_prime(P)
  p = minimum(P)
  if isone(degree(Ok))
    res = ideal(Ok, p)
    res.is_prime = 1
    res.splitting_type = (1, 1)
    return res
  end
  k = domain(f)
  K = codomain(f)
  OK = maximal_order(K)
  if !is_index_divisor(Ok, p) && !is_index_divisor(OK, p)
    return intersect_nonindex(f, P, Ok)
  end
  d = degree(P)
  lp = prime_decomposition(Ok, p, d, 1)
  for (Q, v) in lp
    el = Q.gen_two
    if f(K(el)) in P
      return Q
    end
  end
  error("Restriction not found!")

end

function intersect_nonindex(f::Map, P::NfOrdIdl, Zk::NfOrd = maximal_order(domain(f)))
  @assert is_prime(P)
  #let g be minpoly of k, G minpoly of K and h in Qt the primitive
  #element of k in K (image of gen(k))
  #then
  #  g(h) = 0 mod G
  k = domain(f)
  K = codomain(f)
  G = K.pol
  Qx = parent(G)
  g = change_base_ring(base_ring(Qx), k.pol, parent = Qx)
  h = Qx(f(gen(k)))

  Fp, xp = PolynomialRing(GF(Int(minimum(P)), cached=false), cached=false)
  gp = factor(Fp(g))
  hp = Fp(h)
  Gp = gcd(Fp(K(P.gen_two)), Fp(G))
  for (s, e) in gp
    if iszero(s(hp) % Gp)
      p = ideal_from_poly(Zk, Int(minimum(P)), s, e)
      return p
    end
  end
end

@doc Markdown.doc"""
    prime_decomposition_nonindex(f::Map, p::NfOrdIdl, Z_K::NfOrd) -> Vector{Tuple{NfOrdIdl, Int}}

Given a map $f: k\to K$ of number fields defined over $\mathbb Q$ and
a prime ideal in the maximal order of $k$, find all prime ideals in
the maximal order of $K$ above.
The ideals will belong to $Z_K$ which defaults to "the" maximal order of $K$.
"""
function prime_decomposition(f::Map, p::NfOrdIdl, ZK::NfOrd = maximal_order(codomain(f)))
  @assert p.is_prime == 1
  k = domain(f)
  K = codomain(f)
  if degree(k) == 1
    return prime_decomposition(ZK, minimum(p))
  end
  if !divisible(index(ZK), minimum(p))
    return prime_decomposition_nonindex(f, p, ZK)
  end
  # TODO: Implement for nonindex divisors seriously,
  # splitting the algebra.
  lp = prime_decomposition(ZK, minimum(p))
  res = Tuple{NfOrdIdl, Int}[]
  el = f(p.gen_two.elem_in_nf)
  for (P,_) in lp
    v = valuation(el, P)
    # p has a two-normal presentation, so to test the ramification
    # I only need to test the second element.
    if v > 0
      push!(res, (P, v))
    end
  end
  return res

end

function prime_decomposition_nonindex(f::Map, p::NfOrdIdl, ZK = maximal_order(codomain(f)))

  k = domain(f)
  K = codomain(f)
  G = K.pol
  Qx = parent(G)
  res = Tuple{NfOrdIdl, Int}[]
  if fits(Int, minimum(p))
    Fp = PolynomialRing(GF(Int(minimum(p)), cached = false), cached = false)[1]
    Gp = factor(ppio(Fp(G), Fp(f(p.gen_two.elem_in_nf)))[1])
    for (ke, e) in Gp
      P = ideal_from_poly(ZK, Int(minimum(p)), ke, e)
      push!(res, (P, divexact(e, ramification_index(p))))
    end
  else
    Fp1 = PolynomialRing(GF(minimum(p), cached = false), cached = false)[1]
    Gp1 = factor(ppio(Fp1(G), Fp1(Qx(f(K(p.gen_two)))))[1])
    for (ke, e) in Gp1
      P = ideal_from_poly(ZK, minimum(p), ke, e)
      push!(res, (P, divexact(e, ramification_index(p))))
    end
  end
  return res
end

function prime_decomposition_type(f::Map, p::NfOrdIdl, ZK = maximal_order(codomain(f)))

  if !is_index_divisor(ZK, minimum(p)) && !is_ramified(ZK, minimum(p))
    return prime_decomposition_type_nonindex(f, p, ZK)
  end
  lp = prime_decomposition(f, p, ZK)
  res = Vector{Tuple{Int, Int}}(undef, length(lp))
  for i = 1:length(lp)
    res[i] = (divexact(degree(lp[i][1]), degree(p)), lp[i][2])
  end
  return res

end

function prime_decomposition_type_nonindex(f::Map, p::NfOrdIdl, ZK = maximal_order(codomain(f)))
  k = domain(f)
  K = codomain(f)
  G = K.pol
  Qx = parent(G)

  if fits(Int, minimum(p, copy = false))
    Fp = PolynomialRing(GF(Int(minimum(p)), cached = false), cached = false)[1]
    Gp = factor_shape(gcd(Fp(f(K(p.gen_two))), Fp(G)))
  else
    Fpp = PolynomialRing(GF(minimum(p), cached = false), cached = false)[1]
    Gp = factor_shape(gcd(Fpp(f(K(p.gen_two))), Fpp(G)))
  end
  res = Vector{Tuple{Int, Int}}(undef, sum(values(Gp)))
  ind = 1
  for (d, e) in Gp
    for i = 1:e
      res[ind] = (d, 1)
      ind += 1
    end
  end
  return res
end


text:
module KeyedFrames

using DataFrames

import DataFrames: deletecols!, deleterows!, first, index, last, ncol, nonunique, nrow,
    permutecols!, rename, rename!, SubDataFrame, unique!

struct KeyedFrame <: AbstractDataFrame
    frame::DataFrame
    key::Vector{Symbol}

    function KeyedFrame(df::DataFrame, key::Vector{<:Symbol})
        key = unique(key)
        df_names = names(df)

        if !issubset(key, df_names)
            throw(
                ArgumentError(
                    string(
                        "The columns provided for the key ($key) must all be",
                        "present in the DataFrame ($df_names)."
                    )
                )
            )
        end

        return new(df, key)
    end
end

function KeyedFrame(df::DataFrame, key::Vector{<:AbstractString})
    return KeyedFrame(df, map(Symbol, key))
end

KeyedFrame(df::DataFrame, key::Symbol) = KeyedFrame(df, [key])

"""
    KeyedFrame(df::DataFrame, key::Vector{Symbol})

Create an `KeyedFrame` using the provided `DataFrame`; `key` specifies the columns
to use by default when performing a `join` on `KeyedFrame`s when `on` is not provided.

When performing a `join`, if only one of the arguments is an `KeyedFrame` and `on` is not
specified, the frames will be joined on the `key` of the `KeyedFrame`. If both
arguments are `KeyedFrame`s, `on` will default to the intersection of their respective
indices. In all cases, the result of the `join` will share a type with the first argument.

When calling `unique` (or `unique!`) on a KeyedFrame without providing a `cols` argument,
`cols` will default to the `key` of the `KeyedFrame` instead of all columns. If you wish to
remove only rows that are duplicates across all columns (rather than just across the key),
you can call `unique!(kf, names(kf))`.

When `sort`ing, if no `cols` keyword is supplied, the `key` is used to determine precedence.

When testing for equality, `key` ordering is ignored, which means that it's possible to have
two `KeyedFrame`s that are considered equal but whose default sort order will be different
by virtue of having the columns listed in a different order in their `key`s.
"""
KeyedFrame

DataFrame(kf::KeyedFrame) = frame(kf)
Base.copy(kf::KeyedFrame) = KeyedFrame(copy(DataFrame(kf)), copy(keys(kf)))
Base.deepcopy(kf::KeyedFrame) = KeyedFrame(deepcopy(DataFrame(kf)), deepcopy(keys(kf)))

Base.convert(::Type{DataFrame}, kf::KeyedFrame) = frame(kf)

SubDataFrame(kf::KeyedFrame, args...) = SubDataFrame(frame(kf), args...)

##### EQUALITY #####

Base.:(==)(a::KeyedFrame, b::KeyedFrame) = frame(a) == frame(b) && sort(keys(a)) == sort(keys(b))

Base.isequal(a::KeyedFrame, b::KeyedFrame) = isequal(frame(a), frame(b))&&isequal(keys(a), keys(b))
Base.isequal(a::KeyedFrame, b::AbstractDataFrame) = false
Base.isequal(a::AbstractDataFrame, b::KeyedFrame) = false

Base.hash(kf::KeyedFrame, h::UInt) = hash(keys(kf), hash(frame(kf), h))

##### SIZE #####

nrow(kf::KeyedFrame) = nrow(frame(kf))
ncol(kf::KeyedFrame) = ncol(frame(kf))

##### ACCESSORS #####

index(kf::KeyedFrame) = index(frame(kf))
Base.names(kf::KeyedFrame) = names(frame(kf))

##### INDEXING #####

const ColumnIndex = Union{Real, Symbol}

frame(kf::KeyedFrame) = getfield(kf, :frame)
Base.keys(kf::KeyedFrame) = getfield(kf, :key)
Base.setindex!(kf::KeyedFrame, value, ind...) = setindex!(frame(kf), value, ind...)

# I don't want to have to write the same function body several times, so...
function _kf_getindex(kf::KeyedFrame, index...)
    # If indexing by column, some keys might be removed.
    df = getindex(frame(kf), index...)
    return KeyedFrame(DataFrame(df), intersect(names(df), keys(kf)))
end

# Returns a KeyedFrame
Base.getindex(kf::KeyedFrame, ::Colon) = copy(kf)
Base.getindex(kf::KeyedFrame, ::Colon, ::Colon) = copy(kf)

# Returns a KeyedFrame
Base.getindex(kf::KeyedFrame, col::AbstractVector) = _kf_getindex(kf, col)

# Returns a column
Base.getindex(kf::KeyedFrame, col::ColumnIndex) = frame(kf)[col]

# Returns a KeyedFrame or a column (depending on the type of col)
Base.getindex(kf::KeyedFrame, ::Colon, col) = kf[col]

# Returns a scalar
Base.getindex(kf::KeyedFrame, row::Integer, col::ColumnIndex) = frame(kf)[row, col]

# Returns a KeyedFrame
Base.getindex(kf::KeyedFrame, row::Integer, col::AbstractVector) = _kf_getindex(kf, row, col)

# Returns a column
Base.getindex(kf::KeyedFrame, row::AbstractVector, col::ColumnIndex) = frame(kf)[row, col]

# Returns a KeyedFrame
function Base.getindex(kf::KeyedFrame, row::AbstractVector, col::AbstractVector)
    return _kf_getindex(kf, row, col)
end

# Returns a KeyedFrame
function Base.getindex(kf::KeyedFrame, row::AbstractVector, col::Colon)
    return _kf_getindex(kf, row, col)
end

# Returns a KeyedFrame
Base.getindex(kf::KeyedFrame, row::Integer, col::Colon) = kf[[row], col]

##### SORTING #####

function Base.sort(kf::KeyedFrame, cols=nothing; kwargs...)
    return KeyedFrame(sort(frame(kf), cols === nothing ? keys(kf) : cols; kwargs...), keys(kf))
end

function Base.sort!(kf::KeyedFrame, cols=nothing; kwargs...)
    sort!(frame(kf), cols === nothing ? keys(kf) : cols; kwargs...)
    return kf
end

function Base.issorted(kf::KeyedFrame, cols=nothing; kwargs...)
    return issorted(frame(kf), cols === nothing ? keys(kf) : cols; kwargs...)
end

##### PUSH/APPEND/DELETE #####

function Base.push!(kf::KeyedFrame, data)
    push!(frame(kf), data)
    return kf
end

function Base.append!(kf::KeyedFrame, data)
    append!(frame(kf), data)
    return kf
end

function deleterows!(kf::KeyedFrame, ind)
    deleterows!(frame(kf), ind)
    return kf
end

deletecols!(kf::KeyedFrame, ind::Union{Integer, Symbol}) = deletecols!(kf, [ind])
deletecols!(kf::KeyedFrame, ind::Vector{<:Integer}) = deletecols!(kf, names(kf)[ind])

function deletecols!(kf::KeyedFrame, ind::Vector{<:Symbol})
    deletecols!(frame(kf), ind)
    filter!(x -> !in(x, ind), keys(kf))
    return kf
end

##### RENAME #####

function rename!(kf::KeyedFrame, nms)
    rename!(frame(kf), nms)

    for (from, to) in nms
        i = findfirst(isequal(from), keys(kf))
        if i !== nothing
            keys(kf)[i] = to
        end
    end

    return kf
end

rename!(kf::KeyedFrame, nms::Pair{Symbol, Symbol}...) = rename!(kf, collect(nms))
rename!(f::Function, kf::KeyedFrame) = rename!(kf, [(nm => f(nm)) for nm in names(kf)])

rename(kf::KeyedFrame, args...) = rename!(copy(kf), args...)
rename(f::Function, kf::KeyedFrame) = rename!(f, copy(kf))

##### UNIQUE #####

_unique(kf::KeyedFrame, cols) = KeyedFrame(unique(frame(kf), cols), keys(kf))
function _unique!(kf::KeyedFrame, cols)
    unique!(frame(kf), cols)
    return kf
end

Base.unique(kf::KeyedFrame, cols::AbstractVector) = _unique(kf, cols)
Base.unique(kf::KeyedFrame, cols::Union{Integer, Symbol, Colon}) = _unique(kf, cols)
Base.unique(kf::KeyedFrame) = _unique(kf, keys(kf))
unique!(kf::KeyedFrame, cols::Union{Integer, Symbol, Colon}) = _unique!(kf, cols)
unique!(kf::KeyedFrame, cols::AbstractVector) = _unique!(kf, cols)
unique!(kf::KeyedFrame) = _unique!(kf, keys(kf))

nonunique(kf::KeyedFrame) = nonunique(frame(kf), keys(kf))

##### JOIN #####

# Returns a KeyedFrame
function Base.join(a::KeyedFrame, b::KeyedFrame; on=nothing, kind=:inner, kwargs...)
    df = join(
        frame(a),
        frame(b);
        on=on === nothing ? intersect(keys(a), keys(b)) : on,
        kind=kind,
        kwargs...,
    )

    if kind in (:semi, :anti)
        key = intersect(keys(a), names(df))
    else
        # A join can sometimes rename columns, meaning some of the key columns "disappear"
        key = intersect(union(keys(a), keys(b)), names(df))
    end

    return KeyedFrame(df, key)
end

# Returns a KeyedFrame
function Base.join(a::KeyedFrame, b::AbstractDataFrame; on=nothing, kwargs...)
    df = join(frame(a), b; on=on === nothing ? intersect(keys(a), names(b)) : on, kwargs...)

    # A join can sometimes rename columns, meaning some of the key columns "disappear"
    return KeyedFrame(df, intersect(keys(a), names(df)))
end

# Does NOT return a KeyedFrame
function Base.join(a::AbstractDataFrame, b::KeyedFrame; on=nothing, kwargs...)
    return join(a, frame(b); on=on === nothing ? intersect(keys(b), names(a)) : on, kwargs...)
end

##### FIRST/LAST #####

first(kf::KeyedFrame, r::Int) = KeyedFrame(first(frame(kf), r), keys(kf))
last(kf::KeyedFrame, r::Int) = KeyedFrame(last(frame(kf), r), keys(kf))

##### PERMUTE #####

function permutecols!(kf::KeyedFrame, index::AbstractVector)
    permutecols!(frame(kf), index)
    return kf
end

export KeyedFrame

end


text:
using Pkg

Pkg.add(PackageSpec(url="https://github.com/tolcz/ABCDeGraphGenerator.jl"))


text:
struct TextBuffer{N}
    positions::GPUVector{Point{N, Float32}}
    offsets::GPUVector{Vec2f0}
    rotations::GPUVector{Vec4f0}
    colors::GPUVector{RGBAf0}
    uv_offset_width::GPUVector{Vec4f0}
    scale::GPUVector{Vec2f0}
    text::Vector{Char}
    robj::RenderObject
    range::Signal{Int}
    cursors
end

get_atlas(x::TextBuffer) = x.robj[:atlas]

function TextBuffer(pos::Point{N, <: AbstractFloat} = Point3f0(0)) where N
    positions = gpuvec(Point{N, Float32}[])
    offset = gpuvec(Vec2f0[])
    rotations = gpuvec(Vec4f0[])
    colors = gpuvec(RGBAf0[])
    uv_offset_width = gpuvec(Vec4f0[])
    scale = gpuvec(Vec2f0[])
    atlas = GLVisualize.get_texture_atlas()
    range = Signal(0)
    robj = visualize(
        (DISTANCEFIELD, positions.buffer),
        atlas = atlas,
        offset = offset.buffer,
        rotation = rotations.buffer,
        color = colors.buffer,
        uv_offset_width = uv_offset_width.buffer,
        scale = scale.buffer,
        distancefield = atlas.images,
        indices = range,
        boundingbox = Signal(AABB{Float32}())
    )
    TextBuffer{N}(
        positions,
        offset,
        rotations,
        colors,
        uv_offset_width,
        scale,
        Char[],
        robj.children[],
        range,
        nothing
    )
end

function Base.empty!(tb::TextBuffer{N}) where N
    resize!(tb.positions, 0)
    resize!(tb.offsets, 0)
    resize!(tb.rotations, 0)
    resize!(tb.colors, 0)
    resize!(tb.uv_offset_width, 0)
    resize!(tb.scale, 0)
    push!(tb.range, 0)
    Reactive.set_value!(tb.robj.boundingbox, AABB{Float32}())
    push!(tb.robj.boundingbox, AABB{Float32}())
    return
end

function Base.append!(tb::TextBuffer, startpos::StaticVector{N}, str::String, scale, color, rot, align, font = GLVisualize.defaultfont()) where N
    atlas = get_atlas(tb)
    pos = Point{N, Float32}(startpos)
    rscale = Float32(scale)
    position = GLVisualize.calc_position(str, Point2f0(0), rscale, font, atlas)
    toffset = GLVisualize.calc_offset(str, rscale, font, atlas)
    aoffsetvec = Vec2f0(alignment2num.(align))
    aoffset = align_offset(rot, Point2f0(0), position[end], atlas, rscale, font, aoffsetvec)
    aoffsetn = Point{N, Float32}(to_nd(aoffset, Val{N}, 0f0))
    uv_offset_width = Vec4f0[GLVisualize.glyph_uv_width!(atlas, c, font) for c = str]
    scale = Vec2f0[GLVisualize.glyph_scale!(atlas, c, font, rscale) for c = str]
    position = map(position) do p
        pn = qmul(rot, Point{N, Float32}(to_nd(p, Val{N}, 0f0)))
        pn .+ (pos .+ aoffsetn)
    end

    append!(tb.positions, position)
    append!(tb.offsets, Vec2f0.(toffset))
    append!(tb.rotations, fill(rot, length(position)))
    append!(tb.colors, fill(to_color(color), length(position)))
    append!(tb.uv_offset_width, uv_offset_width)
    append!(tb.scale, scale)

    bb = value(tb.robj.boundingbox)
    for (s, pos) in zip(scale, position)
        pos3d = Vec{3, Float32}(to_nd(pos, Val{3}, 0))
        bb = update(bb, pos3d)
        bb = update(bb, pos3d .+ Vec{3, Float32}(to_nd(s, Val{3}, 0)))
    end
    push!(tb.robj.boundingbox, bb)
    Reactive.set_value!(tb.robj.boundingbox, bb)
    push!(tb.range, length(tb.positions))
    return
end

function align_offset(rot, startpos, lastpos, atlas, rscale, font, align)
    xscale, yscale = GLVisualize.glyph_scale!('X', rscale)
    xmove = (lastpos-startpos)[1] + xscale
    xmove, yscale, z = qmul(rot, Vec3f0(xmove, yscale, 0))
    if isa(align, GeometryTypes.Vec)
        return -Vec2f0(xmove, yscale) .* align
    elseif align == :top
        return -Vec2f0(xmove/2f0, yscale)
    elseif align == :right
        return -Vec2f0(xmove, yscale/2f0)
    else
        error("Align $align not known")
    end
end

function alignment2num(x::Symbol)
    (x in (:hcenter, :vcenter)) && return 0.5f0
    (x in (:left, :bottom)) && return 0.0f0
    (x in (:right, :top)) && return 1.0f0
    0.0f0 # 0 default, or better to error?
end

struct RichText{C, R, F, P, A}
    text::String
    color::C
    rotation::R
    font::F
    position::P
    align::A
end


text:
# initialization

# CUDA packages require complex initialization (discover CUDA, download artifacts, etc)
# that can't happen at module load time, so defer that to run time upon actual use.

"""
    functional(show_reason=false)

Check if the package has been configured successfully and is ready to use.

This call is intended for packages that support conditionally using an available GPU. If you
fail to check whether CUDA is functional, actual use of functionality might warn and error.
"""
function functional(show_reason::Bool=false)
    try
        version()
        toolkit()
        return true
    catch
        show_reason && rethrow()
        return false
    end
end

## deferred initialization API

function __init__()
    # register device overrides
    precompiling = ccall(:jl_generating_output, Cint, ()) != 0
    if !precompiling
        eval(Expr(:block, overrides...))
        empty!(overrides)

        @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" begin
            include("device/intrinsics/special_math.jl")
            eval(Expr(:block, overrides...))
            empty!(overrides)
        end
    end

    # ensure that operations executed by the REPL back-end finish before returning,
    # because displaying values happens on a different task (CUDA.jl#831)
    if isdefined(Base, :active_repl_backend)
        push!(Base.active_repl_backend.ast_transforms, synchronize_cuda_tasks)
    end
end

function synchronize_cuda_tasks(ex)
    quote
        try
            $(ex)
        finally
            $task_local_state() !== nothing && $device_synchronize()
        end
    end
end

@noinline function __init_driver__()
    if haskey(ENV, "_") && basename(ENV["_"]) == "rr"
        @error("Running under rr, which is incompatible with CUDA")
        return
    end

    if version() < v"10.1"
        @warn "This version of CUDA.jl only supports NVIDIA drivers for CUDA 10.1 or higher (yours is for CUDA $(version()))"
    end

    if version() < v"11.2"
        @warn """The NVIDIA driver on this system only supports up to CUDA $(version()).
                 For performance reasons, it is recommended to upgrade to a driver that supports CUDA 11.2 or higher."""
    end

    haskey(ENV, "JULIA_CUDA_MEMORY_LIMIT") &&
        @warn "Support for GPU memory limits (JULIA_CUDA_MEMORY_LIMIT) has been removed."

    haskey(ENV, "JULIA_CUDA_MEMORY_POOL") &&
    ENV["JULIA_CUDA_MEMORY_POOL"] != "none" && ENV["JULIA_CUDA_MEMORY_POOL"] != "cuda" &&
        @warn "Support for memory pools (JULIA_CUDA_MEMORY_POOL) other than 'cuda' and 'none' has been removed."

    # enable generation of FMA instructions to mimic behavior of nvcc
    LLVM.clopts("-nvptx-fma-level=1")

    cuInit(0)

    return
end

function __init_toolkit__()
    if toolkit_release() < v"10.1"
        @warn "This version of CUDA.jl only supports CUDA 10.1 or higher (your toolkit provides CUDA $(toolkit_release()))"
    elseif release() < v"11" && toolkit_release() > release()
        @warn """You are using CUDA toolkit $(toolkit_release()) with a driver that only supports up to $(release()).
                 It is recommended to upgrade your driver, or switch to automatic installation of CUDA."""
    elseif release() >= v"11" && toolkit_release().major > release().major
        @warn """You are using CUDA toolkit $(toolkit_release()) with a driver that only supports up to $(release().major).x.
                 It is recommended to upgrade your driver, or switch to automatic installation of CUDA."""
    end

    return
end

## convenience functions

# TODO: update docstrings

export has_cuda, has_cuda_gpu

# backwards compatibility
export has_cusolvermg, has_cudnn, has_cutensor, has_cupti, has_nvtx

"""
    has_cuda()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit.
Use this function if your code loads packages that require CUDA.jl.

Note that CUDA-dependent packages might still fail to load if the installation is broken,
so it's recommended to guard against that and print a warning to inform the user:

```
using CUDA
if has_cuda()
    try
        using CuArrays
    catch ex
        @warn "CUDA is installed, but CuArrays.jl fails to load" exception=(ex,catch_backtrace())
    end
end
```
"""
has_cuda(show_reason::Bool=false) = functional(show_reason)

"""
    has_cuda_gpu()::Bool

Check whether the local system provides an installation of the CUDA driver and toolkit, and
if it contains a CUDA-capable GPU. See [`has_cuda`](@ref) for more details.

Note that this function initializes the CUDA API in order to check for the number of GPUs.
"""
has_cuda_gpu(show_reason::Bool=false) = has_cuda(show_reason) && length(devices()) > 0


text:
using DelimitedFiles
using Test

"""
    count_chars(data)

Count the characters according to part 1's rules.
"""
function count_chars(data)::Int
    total = 0
    for (r, t, s) in eachrow(data)
        a, b = [parse(Int, ri) for ri in split(r, "-")]
        n = length(collect(m.match for m in eachmatch(Regex("$(t[1])"), s)))
        if a <= n <= b
            total += 1
        end
    end
    return total
end

test = readdlm("./data/day02_test.txt")
@test count_chars(test) == 2

data = readdlm("./data/day02.txt")
println("Part 1 count: $(count_chars(data))")

"""
    count_chars_again(data)

Count according to the new rules.
"""
function count_chars_again(data)::Int
    total = 0
    for (r, t, s) in eachrow(data)
        a, b = [parse(Int, ri) for ri in split(r, "-")]
        n = (s[a] == t[1]) + (s[b] == t[1])
        if n == 1
            total += 1
        end
    end
    return total
end

@test count_chars_again(test) == 1

println("Part 2 count: $(count_chars_again(data))")


text:
using BenchmarkTools

using Rectangle
using Rectangle: _search

include("../bstutils.jl")

begin
    t = BinarySearchTree{Int, Int}()
    t.unique = false
    n = 1000
    a = []
    for i = 1:n
        j = rand(1:10000000000)
        push!(a, j, i)
    end
    for i = 1:n
        insert!(t, a[i], i)
    end
    print("Base.findfirst call to find a number in array: \t\t")
    @btime begin
        findfirst(x -> x == 15, a)
    end

    print("_search a number in a BinarySearchTree : \t\t")
    @btime begin
        n, d = _search(t, t.root, 15)
    end

    @assert parentvalidity(t)
    @assert bstvalidity(t)
    print("delete 50% from BinarySearchTree and add again: \t")
    @btime begin
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
    end
end

begin
    t = RBTree{Int, Int}()
    t.unique = false
    n = 10000
    a = []
    for i = 1:n
        j = rand(1:10000000000)
        push!(a, j, i)
    end
    for i = 1:n
        insert!(t, a[i], i)
    end
    print("Base.find call to find a number in array: \t\t")
    @btime findfirst(x -> x == 15, a)
    print("_search a number in a RBTree : \t\t\t\t")
    @btime begin
        n, d = _search(t, t.root, 15)
    end
    @assert parentvalidity(t)
    @assert bstvalidity(t)
    print("delete 50% from RBTree and add again: \t\t\t")
    @btime begin
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
    end
end


text:
# Autogenerated wrapper script for libass_jll for x86_64-apple-darwin
export libass

using FreeType2_jll
using FriBidi_jll
using HarfBuzz_jll
using Bzip2_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("libass")
JLLWrappers.@declare_library_product(libass, "@rpath/libass.9.dylib")
function __init__()
    JLLWrappers.@generate_init_header(FreeType2_jll, FriBidi_jll, HarfBuzz_jll, Bzip2_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libass,
        "lib/libass.9.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using StoppingTutorials
using Documenter
using Printf
using Literate

assets_src = joinpath(@__DIR__,"..","assets")
assets_dst = joinpath(@__DIR__,"src","assets")

Sys.rm(assets_dst;recursive=true,force=true)

Sys.cp(assets_src,assets_dst)

repo_src = joinpath(@__DIR__,"..","src")
pages_dir = joinpath(@__DIR__,"src","pages")
notebooks_dir = joinpath(@__DIR__,"src","notebooks")

# Add index.md file as introduction to navigation menu
pages = ["Introduction"=> "index.md"]

binder_logo = "https://mybinder.org/badge_logo.svg"
nbviwer_logo = "https://img.shields.io/badge/show-nbviewer-579ACA.svg"

for (i,(title,filename)) in enumerate(StoppingTutorials.files)
  # Generate strings
  tutorial_prefix = string("t",@sprintf "%03d_" i)
  tutorial_title = string("# # Tutorial ", i, ": ", title)
  tutorial_file = string(tutorial_prefix,splitext(filename)[1])
  notebook_filename = string(tutorial_file, ".ipynb")
  binder_url = joinpath("@__BINDER_ROOT_URL__","notebooks", notebook_filename)
  nbviwer_url = joinpath("@__NBVIEWER_ROOT_URL__","notebooks", notebook_filename)
  binder_badge = string("# [![](",binder_logo,")](",binder_url,")")
  nbviwer_badge = string("# [![](",nbviwer_logo,")](",nbviwer_url,")")

  # Generate notebooks
  function preprocess_notebook(content)
    return string(tutorial_title, "\n\n", content)
  end
  Literate.notebook(joinpath(repo_src,filename), notebooks_dir; name=tutorial_file, preprocess=preprocess_notebook, documenter=false, execute=false)

  # Generate markdown
  function preprocess_docs(content)
    return string(tutorial_title, "\n", binder_badge, "\n", nbviwer_badge, "\n\n", content)
  end
  Literate.markdown(joinpath(repo_src,filename), pages_dir; name=tutorial_file, preprocess=preprocess_docs, codefence="```julia" => "```")

  # Generate navigation menu entries
  ordered_title = string(i, " ", title)
  path_to_markdown_file = joinpath("pages",string(tutorial_file,".md"))
  push!(pages, (ordered_title=>path_to_markdown_file))
end

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

makedocs(;
    modules=[StoppingTutorials],
    authors="Tangi Migot tangi.migot@gmail.com",
    repo="https://github.com/SolverStoppingJulia/StoppingTutorials.jl/blob/{commit}{path}#{line}",
    sitename="StoppingTutorials.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://SolverStoppingJulia.github.io/StoppingTutorials.jl",
        assets=String[],
    ),
    pages=pages
)

deploydocs(;
    repo="github.com/SolverStoppingJulia/StoppingTutorials.jl",
    devbranch="main",
)


text:
"""
    AVLTree

struct
"""
mutable struct AVLTree{K,D}
    root::Union{Node{K,D},Nothing}
end

AVLTree() = AVLTree{Any,Any}(nothing)
AVLTree{K,D}() where {K,D} = AVLTree{K,D}(nothing)

Base.eltype(::Type{AVLTree{K,D}}) where {K,D} = Tuple{K,D}
Base.getindex(tr::AVLTree{K,D}, k::K) where {K,D} = Base.getkey(tr, k)
Base.setindex!(tr::AVLTree{K,D}, k::K, d::D) where {K,D} = AVLTrees.insert!(tr, k, d)
Base.haskey(tr::AVLTree{K,D}, k::K) where {K,D} = !(find_node(tr, k) === nothing)
Base.length(tr::AVLTree{K,D}) where {K,D} = AVLTrees.size(tr)
Base.isempty(tr::AVLTree{K,D}) where {K,D} = tr.root === nothing
Base.in(x::K, tr::AVLTree{K,D}) where {K,D} = find_node(tr, x) !== nothing

function Base.getkey(tr::AVLTree{K,D}, k::K) where {K,D}
    d = findkey(tr, k)
    if d === nothing throw(KeyError(k)) else d end
end

function Base.size(tree::AVLTree)
    return __size(tree.root)
end # function

@inline function __size(node::Union{Nothing,Node})
    if node === nothing
        return 0
    end
    return __size(node.left) + __size(node.right) + 1
end

"""
    insert!(args)

documentation
"""
function Base.insert!(tree::AVLTree{K,D}, key, data) where {K,D}
    parent = nothing
    node = tree.root

    while node !== nothing
        parent = node
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            node.data = data
            return
        end
    end

    if parent === nothing
        tree.root = Node{K,D}(key, data)
    elseif key < parent.key
        parent.left = Node{K,D}(key, data, parent)
        balance_insertion(tree, parent, true)
    elseif key > parent.key
        parent.right = Node{K,D}(key, data, parent)
        balance_insertion(tree, parent, false)
    end

    return
end # function

macro rebalance!(_tree, _node, _height_changed)
    tree = esc(_tree)
    node = esc(_node)
    height_changed = esc(_height_changed)

    return :(
        if $(node).bf == 2
            $(node), $(height_changed) = _rebalance_barrier_p2($(tree), $(node), $(node).right)
        elseif $(node).bf == -2
            $(node), $(height_changed) = _rebalance_barrier_m2($(tree), $(node), $(node).left)
        else
            $(height_changed) = $(node).bf == zero(Int8)
        end
    )
end

@inline function _rebalance_barrier_p2(tree::AVLTree{K,D}, node::Node{K,D}, node_right::Node{K,D}) where {K,D}
    height_changed = node_right.bf != zero(Int8)
    if node_right.bf == -one(Int8)
        rotate_right(tree, node_right, node_right.left)
    end
    rotate_left(tree, node, node.right), height_changed
end

@inline function _rebalance_barrier_m2(tree::AVLTree{K,D}, node::Node{K,D}, node_left::Node{K,D}) where {K,D}
    height_changed = node_left.bf != zero(Int8)
    if node_left.bf == one(Int8)
        rotate_left(tree, node_left, node_left.right)
    end
    rotate_right(tree, node, node.left), height_changed
end

"""
    balance_insertion(tree::AVLTree{K,D},node::Node{K,D},left_insertion::Bool) where {K,D}

documentation
"""
@inline function balance_insertion(
    tree::AVLTree{K,D},
    node::Node{K,D},
    left_insertion::Bool,
) where {K,D}
    while true
        node.bf += ifelse(left_insertion, -one(Int8), one(Int8))
        height_changed = false
        @rebalance!(tree, node, height_changed)

        height_changed && break
        node_parent = node.parent
        if node_parent !== nothing
            left_insertion = node_parent.left == node
            node = node_parent
        else
            break
        end
    end
end # function

@inline function rotate_left(t::AVLTree{K,D}, x::Node{K,D}, x_right::Node{K,D}) where {K,D}
    y = x_right

    if y.left !== nothing
        x.right = y.left
        y.left.parent = x
    else
        x.right = nothing
    end
    y.left = x

    xp = x.parent
    if xp === nothing
        t.root = y
    else
        if xp.left == x
            xp.left = y
        else
            xp.right = y
        end
    end

    y.parent = xp
    x.parent = y

    x.bf -= y.bf * (y.bf >= zero(Int8)) + one(Int8)
    y.bf += x.bf * (x.bf < zero(Int8)) - one(Int8)

    return y
end

@inline function rotate_right(t::AVLTree{K,D}, x::Node{K,D}, x_left::Node{K,D}) where {K,D}
    y = x_left

    if y.right !== nothing
        x.left = y.right
        y.right.parent = x
    else
        x.left = nothing
    end
    y.right = x

    xp = x.parent
    if xp === nothing
        t.root = y
    else
        if xp.left == x
            xp.left = y
        else
            xp.right = y
        end
    end

    y.parent = xp
    x.parent = y

    x.bf -= y.bf * (y.bf < zero(Int8)) - one(Int8)
    y.bf += x.bf * (x.bf >= zero(Int8)) + one(Int8)

    return y
end

"""
    delete!(tree::AVLTree{K,D}, node::Node{K,D}) where {K,D}

documentation
"""
function Base.delete!(tree::AVLTree{K,D}, node::Node{K,D}) where {K,D}
    if node.left !== nothing
        node_right = node.right
        if node_right !== nothing
            # left != nothing && right != nothing
            temp = node_right
            temp_left = temp.left
            while temp_left !== nothing
                temp = temp_left
                temp_left = temp.left
            end
            # switch spots completely
            node.key = temp.key
            node.data = temp.data
            delete!(tree, temp)
        else
            # left != nothing && right == nothing
            dir = __parent_replace(tree, node, node.left)
            balance_deletion(tree, node.parent, dir)
        end
    else
        node_right = node.right
        if node_right !== nothing
            # left == nothing && right != nothing
            dir = __parent_replace(tree, node, node_right)
            balance_deletion(tree, node.parent, dir)
        else
            # left == nothing && right == nothing
            dir = __parent_replace(tree, node, nothing)
            balance_deletion(tree, node.parent, dir)
        end
    end
    return
end # function

function Base.delete!(tree::AVLTree{K,D}, key::K) where {K,D}
    node = find_node(tree, key)
    if node !== nothing
        delete!(tree, node)
    end
end # function

@inline balance_deletion(tree::AVLTree, node::Nothing, left_delete::Bool) where {K,D} = return

@inline function balance_deletion(
    tree::AVLTree{K,D},
    node::Node{K,D},
    left_delete::Bool,
) where {K,D}
    while node !== nothing
        node.bf += ifelse(left_delete, one(Int8), -one(Int8))
        height_changed = false
        @rebalance!(tree, node, height_changed)

        !height_changed && break
        node_parent = node.parent
        if node_parent !== nothing
            left_delete = node_parent.left == node
            node = node_parent
        else
            break
        end
    end
end # function

#    __parent_replace(tree::AVLTree{K,D}, node::Node{K,D}, replacement::Node{K,D})
#
# Replaces node with its only child. Used on nodes with a single child when erasing a node.
@inline function __parent_replace(
    tree::AVLTree{K,D},
    node::Node{K,D},
    replacement::Node{K,D},
) where {K,D}
    node_parent = node.parent
    if node_parent !== nothing
        replacement.parent = node_parent
        if node_parent.right == node
            node_parent.right = replacement
            return false
        else
            node_parent.left = replacement
            return true
        end
    else
        replacement.parent = nothing
        tree.root = replacement
        return false
    end
end # function

#    __parent_replace(tree::AVLTree{K,D}, node::Node{K,D}, replacement::Nothing)
# Replaces node with nothing. Used on leaf nodes when erasing a node.
@inline function __parent_replace(
    tree::AVLTree{K,D},
    node::Node{K,D},
    replacement::Nothing,
) where {K,D}
    node_parent = node.parent
    if node_parent !== nothing
        if node_parent.right == node
            node_parent.right = replacement
            return false
        else
            node_parent.left = replacement
            return true
        end
    else
        tree.root = replacement
        return false
    end
end # function

"""
    find(tree::AVLTree{K,D}, key::K) where {K,D}

    Warning: do not use it to check whether `key` is in the `tree`.
    It returns the node.data if found which can be `nothing`.
"""
@inline function findkey(tree::AVLTree{K,D}, key::K) where {K,D}
    node = tree.root
    while node !== nothing
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            return node.data
        end
    end
    return nothing
end # function

"""
    find_node(args)
"""
@inline function find_node(tree::AVLTree{K,D}, key::K) where {K,D}
    node = tree.root
    while node !== nothing
        if key < node.key
            node = node.left
        elseif key > node.key
            node = node.right
        else
            return node
        end
    end
    return nothing
end # function

# Iteration interface

function Base.iterate(tree::AVLTree)
    if tree.root === nothing
        return nothing
    end
    node = tree.root
    while node.left !== nothing
        node = node.left
    end
    return (node.key, node.data), node
end

function Base.iterate(tree::AVLTree, node::Node)
    if node.right !== nothing
        node = node.right
        while node.left !== nothing
            node = node.left
        end
    else
        prev = node
        while node !== nothing && node.left != prev
            prev = node
            node = node.parent
        end
    end

    if node === nothing
        return nothing
    end

    return (node.key, node.data), node
end

# Pop and get methods

function Base.popfirst!(tree::AVLTree)
    # traverse to left-most node
    if tree.root === nothing
        return
    end
    node = tree.root
    while node.left !== nothing
        node = node.left
    end
    # delete node and return data
    node_data = node.data
    delete!(tree, node)
    return node_data
end

function Base.pop!(tree::AVLTree{K,D}, key::K) where {K,D}
    node = AVLTrees.find_node(tree, key)
    if node !== nothing
        node_dat = node.data
        delete!(tree, node)
        return node_dat
    else
        return
    end
end

function Base.firstindex(tree::AVLTree)
    # traverse to left-most node
    if tree.root === nothing
        return
    end
    node = tree.root
    while node.left !== nothing
        node = node.left
    end
    # return node key
    return node.key
end

## Print and Show methods

function Base.print(io::IO, tree::AVLTree{K,D}) where {K,D}
    str_lst = Vector{String}()
    for (k, v) in Base.Iterators.take(tree, 10)
        push!(str_lst, "$k => $v")
    end
    print(io, "AVLTree{$K,$D}(")
    print(io, join(str_lst, ", "))
    length(str_lst) == 10 && print(io, ", ⋯ ")
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", tree::AVLTree{K,D}) where {K,D}
    str_lst = Vector{String}()
    indent_str = "  "
    for (k, v) in Base.Iterators.take(tree, 10)
        push!(str_lst, indent_str * "$k => $v")
    end
    if length(str_lst) > 0
        print(io, "AVLTree{$K,$D} with $(length(tree)) entries:\n")
        print(io, join(str_lst, "\n"))
    else
        print(io, "AVLTree{$K,$D}()")
    end
    length(str_lst) == 10 && print(io, "\n", indent_str * "⋮ => ⋮ \n")
end


text:
function findchart(charts, chart_tree, point)

    for box in boxes(chart_tree, (c,s)->boxesoverlap(c,s,point,0))
        for i in box
            isinclosure(charts[i], point) && return i
        end
    end

    return nothing
end


text:
# geo field countries

# GeoMakie

#Pkg.add("GeoMakie")
#Pkg.add("Proj4")

using CairoMakie
using GeoMakie
using GeoMakie.GeoInterface
using GeoMakie.GeoJSON
using Proj4
CairoMakie.activate!()
let
    worldCountries = GeoJSON.read(read("countries.geojson"))
    print(worldCountries)
    n = length(GeoInterface.features(worldCountries))
    lons = -180:180
    lats = -90:90
    field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]

    trans = Proj4.Transformation("+proj=longlat +datum=WGS84", "+proj=wintri",
        always_xy=true)
    #+proj=wintri, natearth2
    ptrans = Makie.PointTrans{2}(trans)
    fig = Figure(resolution = (1200,800), fontsize = 22)
    ax = Axis(fig[1,1], aspect = DataAspect(), title = "World Countries")
    # all input data coordinates are projected using this function
    ax.scene.transformation.transform_func[] = ptrans
    # add some limits, still it needs to be manual
    points = [Point2f0(lon, lat) for lon in lons, lat in lats]
    rectLimits = FRect2D(Makie.apply_transform(ptrans, points))
    limits!(ax, rectLimits)

    hm1 = GeoMakie.surface!(ax, lons, lats, field, shading = false, overdraw = false)

    hm2 = CairoMakie.poly!(ax, worldCountries, color= 1:n, colormap = Reverse(:plasma),
        strokecolor = :black, strokewidth = 0.25, overdraw = true,
        transformation=(:xy, 10))

    Colorbar(fig[1,2], hm1, label = "variable, color code",
        height = Relative(0.65))

#     drawing the correct lines and ticks is still a little bit cumbersome and manual
#     you could ignore the rest and just hide everything
    lonrange = -180:60:180
    latrange = -90:30:90

    lonlines = [Point2f0(j,i) for i in lats, j in lonrange]
    latlines = [Point2f0(j,i) for j in lons, i in latrange]

    [lines!(ax, lonlines[:,i], color = (:black,0.25),
        linestyle = :dash, overdraw = true) for i in 1:size(lonlines)[2]]
    [lines!(ax, latlines[:,i], color = (:black,0.25), linestyle = :dash,
        overdraw = true) for i in 1:size(latlines)[2]]

    xticks = first.(trans.(Point2f0.(lonrange, -90)))
    yticks = last.(trans.(Point2f0.(-180,latrange)))
    ax.xticks = (xticks, string.(lonrange, 'ᵒ'))
    ax.yticks = (yticks, string.(latrange, 'ᵒ'))
    hidedecorations!(ax, ticks = false, label = false, ticklabels=false)
    hidespines!(ax)
    fig

end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using MIToS.PDB
using BenchmarkTools
pdb_filepath = "1ake.pdb"
struc = read(pdb_filepath, PDBFile)
struc2 = read(pdb_filepath, PDBFile)

@show rmsd(struc,struc2,superimposed=true)

@show @benchmark rmsd(struc,struc2,superimposed=true)