text:
import LaTeXStrings
import PGFPlotsX
import StatsBase
import ValueHistories

"""
"""
function plotlearningcurve(
        inputobject::AbstractFittable,
        curvetype::Symbol = :loss_vs_iteration;
        kwargs...,
        )
    history::ValueHistories.MultivalueHistory = get_history(inputobject)
    result = plotlearningcurve(
        history,
        curvetype;
        kwargs...,
        )
    return result
end

"""
"""
function plotlearningcurve(
        history::ValueHistories.MultivalueHistory,
        curvetype::Symbol = :loss_vs_iteration;
        window::Integer = 0,
        legend_pos::AbstractString = "outer north east",
        sampleevery::Integer = 1,
        startat::Union{Integer, Symbol} = :start,
        endat::Union{Integer, Symbol} = :end,
        include_tuning::Bool = true,
        show_raw::Bool = true,
        show_smoothed::Bool = true,
        )
    legend_pos::String = convert(String, legend_pos)
    if curvetype==:loss_vs_iteration
        has_tuning = include_tuning &&
                haskey(history, :tuning_loss_at_iteration)
        xlabel = "Iteration"
        ylabel = "Loss"
        legendentry = "Loss function"
        training_xvalues, training_yvalues = ValueHistories.get(
            history,
            :training_loss_at_iteration,
            )
        if has_tuning
            tuning_xvalues, tuning_yvalues = ValueHistories.get(
                history,
                :tuning_loss_at_iteration,
                )
        end
    elseif curvetype==:loss_vs_epoch
        has_tuning = include_tuning &&
           haskey(history, :tuning_loss_at_epoch)
        xlabel = "Epoch"
        ylabel = "Loss"
        legendentry = "Loss function"
        training_xvalues, training_yvalues = ValueHistories.get(
            history,
            :training_loss_at_epoch,
            )
        if has_tuning
            tuning_xvalues, tuning_yvalues = ValueHistories.get(
                history,
                :tuning_loss_at_epoch,
                )
        end
    else
        error("\"curvetype\" must be one of: :loss_vs_iteration, :loss_vs_epoch")
    end
    if length(training_xvalues) != length(training_yvalues)
        error("length(training_xvalues) != length(training_yvalues)")
    end
    if has_tuning
        if length(training_xvalues) != length(tuning_yvalues)
            error("length(training_xvalues) != length(tuning_yvalues)")
        end
        if length(training_xvalues) != length(tuning_xvalues)
            error("length(training_xvalues) != length(tuning_xvalues)")
        end
        if !all(training_xvalues .== tuning_xvalues)
            error("!all(training_xvalues .== tuning_xvalues)")
        end
    end
    if startat == :start
        startat = 1
    elseif typeof(startat) <: Symbol
        error("$(startat) is not a valid value for startat")
    end
    if endat == :end
        endat = length(training_xvalues)
    elseif typeof(endat) <: Symbol
        error("$(endat) is not a valid value for endat")
    end
    if startat > endat
        error("startat > endat")
    end
    training_xvalues = training_xvalues[startat:endat]
    training_yvalues = training_yvalues[startat:endat]
    if has_tuning
        tuning_xvalues = tuning_xvalues[startat:endat]
        tuning_yvalues = tuning_yvalues[startat:endat]
    end
    if length(training_xvalues) != length(training_yvalues)
        error("length(training_xvalues) != length(training_yvalues)")
    end
    if has_tuning
        if length(tuning_xvalues) != length(tuning_yvalues)
            error("length(tuning_xvalues) != length(tuning_yvalues)")
        end
        if length(training_xvalues) != length(tuning_xvalues)
            error("length(training_xvalues) != length(tuning_xvalues)")
        end
        if !all(training_xvalues .== tuning_xvalues)
            error("!all(training_xvalues .== tuning_xvalues)")
        end
        result = plotlearningcurve(
            training_xvalues,
            training_yvalues,
            xlabel,
            ylabel,
            legendentry;
            window = window,
            legend_pos = legend_pos,
            sampleevery = sampleevery,
            tuning_yvalues = tuning_yvalues,
            show_raw = show_raw,
            show_smoothed = show_smoothed,
            )
    else
        result = plotlearningcurve(
            training_xvalues,
            training_yvalues,
            xlabel,
            ylabel,
            legendentry;
            window = window,
            legend_pos = legend_pos,
            sampleevery = sampleevery,
            show_raw = show_raw,
            show_smoothed = show_smoothed,
            )
    end
    return result
end

"""
"""
function plotlearningcurve(
        xvalues::AbstractVector{<:Real},
        training_yvalues::AbstractVector{<:Real},
        xlabel::AbstractString,
        ylabel::AbstractString,
        legendentry::AbstractString;
        window::Integer = 0,
        legend_pos::AbstractString = "outer north east",
        sampleevery::Integer = 1,
        tuning_yvalues::Union{Nothing, AbstractVector{<:Real}} = nothing,
        show_raw::Bool = true,
        show_smoothed::Bool = true,
        )
    legend_pos::String = convert(String, legend_pos)
    if !show_raw && !show_smoothed
        error("At least one of show_raw, show_smoothed must be true.")
    end
    if is_nothing(tuning_yvalues)
        has_tuning = false
    else
        has_tuning = true
    end
    if has_tuning
        training_legendentry = string(strip(legendentry),
                ", training set")
        tuning_legendentry = string(strip(legendentry),
                ", tuning set")
    else
        training_legendentry = string(strip(legendentry),
                ", training set")
    end
    if sampleevery < 1
        error("sampleevery must be >=1")
    end
    if length(xvalues) != length(training_yvalues)
        error("length(xvalues) != length(training_yvalues)")
    end
    if length(xvalues) == 0
        error("length(xvalues) == 0")
    end
    xvalues = xvalues[1:sampleevery:end]
    training_yvalues = training_yvalues[1:sampleevery:end]
    if has_tuning
        tuning_yvalues = tuning_yvalues[1:sampleevery:end]
    end
    all_plots_and_legends = []
    if show_raw
        training_linearplotobject_yraw = PGFPlotsX.@pgf(
            PGFPlotsX.Plot(
                PGFPlotsX.Coordinates(
                    xvalues,
                    training_yvalues,
                    )
                )
           )
        training_legendentry_yraw = PGFPlotsX.@pgf(
            PGFPlotsX.LegendEntry(
                LaTeXStrings.LaTeXString(training_legendentry)
                )
            )
        push!(all_plots_and_legends, training_linearplotobject_yraw)
        push!(all_plots_and_legends, training_legendentry_yraw)
        if has_tuning
            tuning_linearplotobject_yraw = PGFPlotsX.@pgf(
                PGFPlotsX.Plot(
                    PGFPlotsX.Coordinates(
                        xvalues,
                        tuning_yvalues,
                        )
                    )
               )
            tuning_legendentry_yraw = PGFPlotsX.@pgf(
                PGFPlotsX.LegendEntry(
                    LaTeXStrings.LaTeXString(tuning_legendentry)
                    )
                )
            push!(all_plots_and_legends, tuning_linearplotobject_yraw)
            push!(all_plots_and_legends, tuning_legendentry_yraw)
        end
    end
    if show_smoothed && window > 0
        training_yvaluessmoothed = simple_moving_average(
            training_yvalues,
            window,
            )
        training_legendentry_smoothed = string(
                strip(training_legendentry),
                " (smoothed)",
                )
        training_linearplotobject_ysmoothed = PGFPlotsX.@pgf(
            PGFPlotsX.Plot(
                PGFPlotsX.Coordinates(
                    xvalues,
                    training_yvaluessmoothed,
                    )
                )
           )
        training_legendentry_ysmoothed = PGFPlotsX.@pgf(
            PGFPlotsX.LegendEntry(
                LaTeXStrings.LaTeXString(training_legendentry_smoothed)
                )
            )
        push!(all_plots_and_legends, training_linearplotobject_ysmoothed)
        push!(all_plots_and_legends, training_legendentry_ysmoothed)
        if has_tuning
            tuning_yvaluessmoothed = simple_moving_average(
                tuning_yvalues,
                window,
                )
            tuning_legendentry_smoothed = string(
                    strip(tuning_legendentry),
                    " (smoothed)",
                    )
            tuning_linearplotobject_ysmoothed = PGFPlotsX.@pgf(
                PGFPlotsX.Plot(
                    PGFPlotsX.Coordinates(
                        xvalues,
                        tuning_yvaluessmoothed,
                        )
                    )
               )
            tuning_legendentry_ysmoothed = PGFPlotsX.@pgf(
                PGFPlotsX.LegendEntry(
                    LaTeXStrings.LaTeXString(tuning_legendentry_smoothed)
                    )
                )
            push!(all_plots_and_legends, tuning_linearplotobject_ysmoothed)
            push!(all_plots_and_legends, tuning_legendentry_ysmoothed)
        end
    end
    p = PGFPlotsX.@pgf(
        PGFPlotsX.Axis(
            {
                xlabel = LaTeXStrings.LaTeXString(xlabel),
                ylabel = LaTeXStrings.LaTeXString(ylabel),
                no_markers,
                legend_pos = legend_pos,
                },
            all_plots_and_legends...,
            ),
        )
    wrapper = PGFPlotsXPlot(p)
    return wrapper
end


text:
using Oceananigans.BoundaryConditions: UVelocityBoundaryConditions, VVelocityBoundaryConditions, TracerBoundaryConditions
using Oceananigans.Fields: XFaceField, YFaceField, CenterField, TracerFields

function HorizontalVelocityFields(velocities, arch, grid)
    u = XFaceField(arch, grid, UVelocityBoundaryConditions(grid))
    v = YFaceField(arch, grid, VVelocityBoundaryConditions(grid))
    return u, v
end

function HydrostaticFreeSurfaceTendencyFields(velocities, free_surface, arch, grid, tracer_names)
    u, v = HorizontalVelocityFields(velocities, arch, grid)
    η = FreeSurfaceDisplacementField(velocities, arch, grid)
    tracers = TracerFields(tracer_names, arch, grid)
    return merge((u=u, v=v, η=η), tracers)
end

ext: jl


text:
@discord_object struct EmbedField
    name::String
    value::String
    inline::Bool
end


text:
"""
    default_equality(reference, actual) -> f

Infer a suitable equality comparison method `f` according to input types.

`f` is a function that satisfies signature `f(reference, actual)::Bool`. If `f` outputs
`true`, it indicates that `reference` and `actual` are "equal" in the sense of `f`.
"""
default_equality(reference, actual) = isequal
function default_equality(
    reference::AbstractArray{<:Colorant},
    actual::AbstractArray{<:Colorant})

    return psnr_equality()
end

# ---------------------------------
# Image

"""
    psnr_equality(threshold=25) -> f

Generates an equality comparison function in terms of Peak Signal-to-Noise Ratio (PSNR).

The return function `f` accepts two images `reference` and `actual` as its inputs;
`f(reference, actual) == true` if `psnr(reference, actual) >= threshold`.

This function is useful for image comparison, for example:

```julia
@test_reference "references/camera.png" testimage("cameraman") by=psnr_equality(20)
```
"""
function psnr_equality(threshold=25)
    function (ref, x)
        if size(ref) != size(x)
            # trigger a test fail instead of throwing error
            @warn("test fails because size(ref) $(size(ref)) != size(x) $(size(x))")
            return false
        end
        rst = _psnr(ref, x)
        if rst >= threshold
            return true
        else
            @warn("test fails because PSNR $rst < $threshold")
            return false
        end
    end
end

# a simplified PSNR is sufficient since we only use it to approximately compare two images
_psnr(ref::AbstractArray{<:Color3}, x::AbstractArray{<:Color3}) =
    _psnr(channelview(RGB.(ref)), channelview(RGB.(x)), 1.0)

_psnr(ref::AbstractArray{<:ColorTypes.Transparent3}, x::AbstractArray{<:ColorTypes.Transparent3}) =
    _psnr(channelview(ARGB.(ref)), channelview(ARGB.(x)), 1.0)

_psnr(ref::AbstractArray{<:AbstractGray}, x::AbstractArray{<:AbstractGray}) =
    _psnr(channelview(ref), channelview(x), 1.0)

_psnr(ref::AbstractArray{<:Real}, x::AbstractArray{<:Real}, peakval::Real) =
    20log10(peakval) - 10log10(euclidean(float.(ref), float.(x)))


text:
using Stochastics
using Test

@testset "Regression.jl" begin
    x₁ = [56.7, 66.6, 55.7, 63.5, 63.1, 57.6, 57.8, 60.4, 57.4, 59.2, 59.3]
    y₁ = [7.5, -0.6,8.5, 4.0, 4.2, 7.6, 7.7, 5.8, 6.0, 7.0, 7.0]
    d = Stochastics.regression_linear_simple(x₁, y₁)
    @test round(get(d, "α", nothing), digits=4) == 48.9646
    @test round(get(d, "β", nothing), digits=4) == -0.7210
    @test round(√(get(d, "s²", nothing)), digits=4) == 0.9431
end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
using Test, ComputedGeodesicEquations

@testset "BoyerLindquist Coordinates" begin
    M = 1.0
    a = 1.0
    u = [0.0, 100.0, π / 2.0, 0.0]
    v = [0.0, -1.0, 0.0, 1e-4]

    @test BoyerLindquistCoords.null_constrain(u, v, M, a) ≈ 1.020404071219846

    v[1] = 1.020404071219846
    @test all(
        BoyerLindquistCoords.geodesic_eq(u, v, M, a) .≈ (
            0.00020818408358881096,
            -4.0208173160182054e-8,
            -2.4612832272924962e-26,
            2.020616346724209e-6,
        ),
    )

    m = BoyerLindquist(M = M, a = a)
    @test BoyerLindquistCoords.null_constrain(u, v, M, a) ≈ null_constrain(u, v, m)
end

@testset "EddingtonFinkelstein Coordinates" begin
    M = 1.0
    u = [0.0, 100.0, π / 2.0, 0.0]
    v = [0.0, -1.0, 0.0, 1e-4]

    @test EddingtonFinkelsteinCoords.null_constrain(u, v, M) ≈ -1.0408663253056722

    v[1] = -1.0408663253056722
    @test all(
        EddingtonFinkelsteinCoords.geodesic_eq(u, v, M) .≈
        (-0.00041648353577666787, 9.699999999999988e-7, 6.123233995736766e-25, 2.0e-6),
    )
end


text:
"""
region_min_dist is the minimum horizontal distance between regions, in ppm.
"""
function prepareoptim(config_path::String,
    molecule_names,
    hz2ppmfunc,
    U_cost0, y_cost0::Vector{Complex{T}},
    As;
    region_min_dist = 0.1) where T <: Real

    ## parse from config file.
    config_dict = Dict()
    if ispath(config_path)

        # TODO add error-handling if name is not found in the dictionary or filename does not exist.
        config_dict = JSON.parsefile(config_path)
    end

    N_compounds = length(molecule_names)

    λ_lbs = Vector{Vector{T}}(undef, N_compounds)
    λ_ubs = Vector{Vector{T}}(undef, N_compounds)
    Δsys_cs = Vector{Vector{T}}(undef, N_compounds)
    κs_β_orderings = Vector{Vector{Vector{Int}}}(undef, N_compounds)
    κs_β_DOFs = Vector{Vector{Int}}(undef, N_compounds)

    for n = 1:N_compounds

        dict = config_dict[molecule_names[n]] # TODO graceful error-handle.

        λ_lbs[n] = convert(Vector{T}, dict["λ_lb"])
        λ_ubs[n] = convert(Vector{T}, dict["λ_ub"])
        Δsys_cs[n] = convert(Vector{T}, dict["maximum chemical shift"])
        κs_β_orderings[n] = convert(Vector{Vector{Int}}, dict["κs_β ordering"])
        κs_β_DOFs[n] = convert(Vector{Int}, dict["κs_β degrees of freedom"])

    end

    # cs_delta_group = NMRSpecifyRegions.extractinfofromconfig( cs_config_path, molecule_names)
    # Δsys_cs = NMRSpecifyRegions.condenseΔcsconfig(cs_delta_group)

    ## get regions.
    ΩS0 = NMRSpectraSimulator.getΩS(As)
    ΩS0_ppm = NMRSpectraSimulator.getPs(ΩS0, hz2ppmfunc)

    exp_info = NMRSpecifyRegions.setupexperimentresults(molecule_names, ΩS0_ppm, Δsys_cs;
    min_dist = region_min_dist)

    P_cost0 = hz2ppmfunc.(U_cost0)
    cost_inds, cost_inds_set = NMRSpecifyRegions.getcostinds(exp_info, P_cost0)

    U_cost = U_cost0[cost_inds]
    P_cost = P_cost0[cost_inds]
    y_cost = y_cost0[cost_inds]

    return Δsys_cs, y_cost, U_cost, P_cost, exp_info, cost_inds, cost_inds_set,
    λ_lbs,
    λ_ubs,
    κs_β_orderings,
    κs_β_DOFs
end

function fitβLSκ(U_rad_cost,
    y_cost::Vector{Complex{T}},
    LS_inds,
    Es::Vector{NMRSpectraSimulator.κCompoundFIDType{T, SST}},
    As::Vector{NMRSpectraSimulator.CompoundFIDType{T, SST}},
    β_initial::Vector{T},
    β_lb,
    β_ub;
    optim_algorithm::Symbol = :GN_ESCH,
    max_iters = 5000,
    xtol_rel = 1e-7,
    ftol_rel = 1e-12,
    maxtime = Inf,
    w = ones(T, length(Es)),
    κ_lb_default = 0.2,
    κ_ub_default = 50.0) where {T <: Real, SST}

    @assert length(U_rad_cost) == length(y_cost) == length(LS_inds)

    q, updateβfunc, updateκfunc,
    E_BLS, κ_BLS, b_BLS, getβfunc = setupcostβLS(Es, As, LS_inds, U_rad_cost, y_cost)

    # updateκfunc(0.0)
    # NMRCalibrate.parseκ!(Es, κ_BLS)
    obj_func = pp->costβLS(U_rad_cost, y_cost, updateβfunc, updateκfunc, pp, Es, E_BLS, κ_BLS, b_BLS, q)
    grad_func = xx->FiniteDiff.finite_difference_gradient(obj_func, xx)

    # force eval for debugging.
    println("Timing: obj_func")
    @time cost_initial = obj_func(β_initial)
    cost_initial = obj_func(β_initial)
    println("obj_func(β_initial) = ", cost_initial)
    println()

    opt = NLopt.Opt(optim_algorithm, length(β_initial)) # global evolutionary.
println(β_initial)
    minf, minx, ret, numevals = runNLopt!(  opt,
        β_initial,
        obj_func,
        grad_func,
        β_lb,
        β_ub;
        max_iters = max_iters,
        xtol_rel = xtol_rel,
        ftol_rel = ftol_rel,
        maxtime = maxtime)

    println("Timing: obj_func")
    @time cost_final = obj_func(minx)
    println("obj_func(minx) = ", cost_final)
    println()

    return minx, q, κ_BLS, getβfunc, obj_func
end


text:
"""
    run_unit_benchmark(
        model::JuMP.Model;
        direct_model::Bool,
        atol::Float64 = 1e-10,
    )

Given a JuMP Model `model`, run a series of unit tests. This can be helpful to
test differences between JuMP and SymbolicAD (using an atol of `atol`), and also
to check relative performance.

If `direct_model`, build the SymbolicAD backend directly from JuMP's NLP
datastructures. If `direct_model == false`, build from `MOI.objective_expr` and
`MOI.constraint_expr`.
"""
function run_unit_benchmark(
    model::JuMP.Model;
    direct_model::Bool,
    atol::Float64 = 1e-10,
    rtol::Float64 = 0.0,
)
    @info "Constructing oracles"

    @time d = JuMP.NLPEvaluator(model)
    @time nlp_block = if direct_model
        _nlp_block_data(model; backend = DefaultBackend())
    else
        _nlp_block_data(d; backend = DefaultBackend())
    end
    oracle = nlp_block.evaluator

    @info "MOI.initialize"

    @time MOI.initialize(d, [:Grad, :Jac, :Hess])
    @time MOI.initialize(oracle, [:Grad, :Jac, :Hess])

    # !!! note
    #     When trying to micro-benchmark this, be careful, because JuMP caches
    #     the last x and avoids redundant computation. Therefore calling
    #     MOI.eval_constraint multiple times with the same x will be optimistic
    #     compared with calling it with different x

    x = map(JuMP.all_variables(model)) do xi
        lb = JuMP.has_lower_bound(xi) ? JuMP.lower_bound(xi) : -1e6
        ub = JuMP.has_upper_bound(xi) ? JuMP.upper_bound(xi) : 1e6
        return lb + (ub - lb) * rand()
    end

    if d.has_nlobj
        # MOI.eval_objective

        @info "MOI.eval_objective"

        @time MOI.eval_objective(d, x)
        @time MOI.eval_objective(d, x)
        @time MOI.eval_objective(oracle, x)
        @time MOI.eval_objective(oracle, x)

        d_eval_objective = MOI.eval_objective(d, x)
        oracle_eval_objective = MOI.eval_objective(oracle, x)
        Test.@test isapprox(
            d_eval_objective,
            oracle_eval_objective;
            atol = atol,
            rtol = rtol,
        )
        println("Error = ", d_eval_objective - oracle_eval_objective)

        # MOI.eval_objective_gradient

        @info "MOI.eval_objective_gradient"

        f = zeros(length(x))
        f1 = zeros(length(x))

        @time MOI.eval_objective_gradient(d, f1, x)
        @time MOI.eval_objective_gradient(d, f1, x)
        @time MOI.eval_objective_gradient(oracle, f, x)
        @time MOI.eval_objective_gradient(oracle, f, x)

        Test.@test isapprox(f, f1; atol = atol, rtol = rtol)
        println("Error = ", maximum(abs.(f .- f1)))
    end

    # MOI.eval_constraint

    @info "MOI.eval_constraint"

    g = zeros(length(oracle.constraints))
    g1 = zeros(length(oracle.constraints))

    @time MOI.eval_constraint(d, g1, x)
    @time MOI.eval_constraint(d, g1, x)
    @time MOI.eval_constraint(oracle, g, x)
    @time MOI.eval_constraint(oracle, g, x)

    Test.@test isapprox(g, g1; atol = atol, rtol = rtol)
    println("Error = ", maximum(abs.(g .- g1)))

    # MOI.jacobian_structure

    @info "MOI.jacobian_structure"

    @time J1 = MOI.jacobian_structure(d)
    @time J = MOI.jacobian_structure(oracle)

    p, p1 = sortperm(J), sortperm(J1)
    Test.@test J[p] == J1[p1]

    # MOI.eval_constraint_jacobian

    @info "MOI.eval_constraint_jacobian"

    J_nz, J_nz_1 = zeros(length(J)), zeros(length(J1))

    @time MOI.eval_constraint_jacobian(d, J_nz_1, x)
    @time MOI.eval_constraint_jacobian(d, J_nz_1, x)
    @time MOI.eval_constraint_jacobian(oracle, J_nz, x)
    @time MOI.eval_constraint_jacobian(oracle, J_nz, x)

    Test.@test isapprox(J_nz[p], J_nz_1[p1]; atol = atol, rtol = rtol)
    println("Error = ", LinearAlgebra.norm(J_nz[p] .- J_nz_1[p1]))

    # MOI.hessian_lagrangian_structure

    @info "MOI.hessian_lagrangian_structure"

    @time H1 = MOI.hessian_lagrangian_structure(d)
    @time H = MOI.hessian_lagrangian_structure(oracle)

    # MOI.eval_hessian_lagrangian

    @info "MOI.eval_hessian_lagrangian"

    μ, H_nz, H_nz_1 = rand(length(g)), zeros(length(H)), zeros(length(H1))
    σ = rand()
    @time MOI.eval_hessian_lagrangian(d, H_nz_1, x, σ, μ)
    @time MOI.eval_hessian_lagrangian(d, H_nz_1, x, σ, μ)
    @time MOI.eval_hessian_lagrangian(oracle, H_nz, x, σ, μ)
    @time MOI.eval_hessian_lagrangian(oracle, H_nz, x, σ, μ)

    # For the Hessian, we won't compute identical sparsity patterns due to
    # duplicates, etc. What we need to do is check the full sparse matrix.

    function _to_sparse(IJ, V)
        I, J = [i for (i, _) in IJ], [j for (_, j) in IJ]
        n = max(maximum(I), maximum(J))
        return SparseArrays.sparse(I, J, V, n, n)
    end

    d_sparse = _to_sparse(H, H_nz)
    oracle_sparse = _to_sparse(H1, H_nz_1)
    Test.@test isapprox(d_sparse, oracle_sparse; atol = atol, rtol = rtol)
    println("Error = ", LinearAlgebra.norm(d_sparse - oracle_sparse))
    return
end

"""
    run_solution_benchmark(model::JuMP.Model, optimizer; atol = 1e-6)

Given a JuMP Model `model`, solve the problem using JuMP's built-in AD, as well
as the `DefaultBackend` and `ThreadedBackend` of SymbolicAD and check that they
find the same solutions.

`optimizer` must be something that can be passed to `MOI.instantiate`.
"""
function run_solution_benchmark(
    model::JuMP.Model,
    optimizer;
    atol::Float64 = 1e-6,
    rtol::Float64 = 0.0,
)
    @info "Solving with JuMP"

    jump_model = MOI.instantiate(optimizer)
    MOI.copy_to(jump_model, model)
    MOI.optimize!(jump_model)
    jump_solution = MOI.get(
        jump_model,
        MOI.VariablePrimal(),
        MOI.get(jump_model, MOI.ListOfVariableIndices()),
    )
    jump_nlp_block = MOI.get(model, MOI.NLPBlock())

    @info "Solving with serial SymbolicAD"

    serial_model = MOI.instantiate(optimizer)
    MOI.copy_to(serial_model, model)
    serial_nlp_block =
        _nlp_block_data(JuMP.NLPEvaluator(model); backend = DefaultBackend())
    MOI.set(serial_model, MOI.NLPBlock(), serial_nlp_block)
    MOI.optimize!(serial_model)
    serial_solution = MOI.get(
        serial_model,
        MOI.VariablePrimal(),
        MOI.get(serial_model, MOI.ListOfVariableIndices()),
    )

    @info "Solving with threaded SymbolicAD"

    threaded_model = MOI.instantiate(optimizer)
    MOI.copy_to(threaded_model, model)
    threaded_nlp_block =
        _nlp_block_data(JuMP.NLPEvaluator(model); backend = ThreadedBackend())
    MOI.set(threaded_model, MOI.NLPBlock(), threaded_nlp_block)
    MOI.optimize!(threaded_model)
    threaded_solution = MOI.get(
        threaded_model,
        MOI.VariablePrimal(),
        MOI.get(threaded_model, MOI.ListOfVariableIndices()),
    )

    @info "Validating solutions"
    println("Errors = ", extrema(jump_solution .- serial_solution))
    Test.@test ≈(jump_solution, serial_solution; atol = atol, rtol = rtol)
    Test.@test ≈(jump_solution, threaded_solution; atol = atol, rtol = rtol)

    @info "Timing statistics"

    println("JuMP timing statistics")
    println(" - ", jump_nlp_block.evaluator.eval_constraint_timer)
    println(" - ", jump_nlp_block.evaluator.eval_constraint_jacobian_timer)
    println(" - ", jump_nlp_block.evaluator.eval_hessian_lagrangian_timer)

    println("Serial SymbolicAD timing statistics")
    println(" - ", serial_nlp_block.evaluator.eval_constraint_timer)
    println(" - ", serial_nlp_block.evaluator.eval_constraint_jacobian_timer)
    println(" - ", serial_nlp_block.evaluator.eval_hessian_lagrangian_timer)

    println("Threaded SymbolicAD timing statistics")
    println(" - ", threaded_nlp_block.evaluator.eval_constraint_timer)
    println(" - ", threaded_nlp_block.evaluator.eval_constraint_jacobian_timer)
    println(" - ", threaded_nlp_block.evaluator.eval_hessian_lagrangian_timer)
    return
end


text:
"""
    allele_avg(data::PopObj, rounding::Bool = true)
Returns a NamedTuple of the average number of alleles ('avg') and standard
deviation (`stdev`) of a `PopObj`. Use `false` as second argument (no keyword)
to not round results. Default (`true`) rounds to 4 digits.
"""
function allele_avg(data::PopObj, rounding::Bool = true)
    num_alleles = map(i -> Iterators.flatten(i |> skipmissing) |> collect |> unique |> length, eachcol(data.loci, false))
    avg = mean(num_alleles)
    sd = std(num_alleles)
    if rounding == true
        return (avg = round(avg, digits = 4), stdev = round(sd, digits = 4))
    else
        return (avg = avg, stdev = sd)
    end
end

"""
    richness(data::PopObj)
Calculates various allelic richness and returns vector of per-locus allelic richness.
To be called internally by functions calculating overall or per-population richness
either rarefied or not.
"""
function richness(data::PopObj)
    # collapse genotypes into array of alleles |> count number of unique alleles
    rich = map(eachcol(data.loci, false)) do i
        Iterators.flatten(i |> skipmissing) |> collect |> unique |> length
    end
    return DataFrame(locus = names(data.loci), richness = rich)
end

"""
    summary(data::PopObj)
Prints a summary of the information contained in a PopObj
"""
function Base.summary(data::PopObj)

    println(" Object of type PopObj")
    if typeof(skipmissing(data.loci[!, :1])[1][1]) == Int16
        marker = "Microsatellite"
    else
        marker = "SNP"
    end
    println(" Marker type: ", marker)
    ploidy = unique(data.samples.ploidy) |> sort
    length(ploidy) == 1 && println(" Ploidy: ", ploidy |> join)
    if length(ploidy) != 1
        print(" Ploidy (varies): ")
        print(ploidy[1]), [print(", $i") for i in ploidy[2:end]]
    end
    println("\n Number of individuals: ", length(data.samples.name))
    println(" Number of loci: ", size(data.loci,2))
    if ismissing.(data.samples.longitude) |> all == true
        print(" Longitude:") ; printstyled(" absent\n", color = :yellow)
    else
        println(" Longitude: present with ", count(i -> i === missing, data.samples.longitude), " missing")
    end
    if ismissing.(data.samples.longitude) |> all == true
        print(" Latitude:") ; printstyled(" absent\n", color = :yellow)
    else
        println(" Latitude: present with ", count(i -> i === missing, data.samples.latitude), " missing")
    end
    print("\n Population names and counts:")
    DataFrames.show(populations(data), summary = false, rowlabel = :num, allrows = true)
end


text:
module AssemblersTests

using Test
using Gridap
using SparseArrays
using SparseMatricesCSR

for T in [SparseMatrixCSC, SparseMatrixCSR, SymSparseMatrixCSR]
    model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))

    tags = [1,2,3,4,6,5]

    order = 1
    fespace = H1ConformingFESpace(Float64,model,order,tags)

    assem = SparseMatrixAssembler(T,fespace,fespace)

    trian = Triangulation(model)

    quad = CellQuadrature(trian,degree=2)

    basis = CellBasis(fespace)

    a(v,u) = varinner(∇(v), ∇(u))

    bfun(x) = x[2]

    b(v) = varinner(v,CellField(trian,bfun))

    mmat = integrate(a(basis,basis),trian,quad)

    bvec = integrate(b(basis),trian,quad)

    vec = assemble(assem, bvec)

    mat = assemble(assem, mmat)

    x = mat \ vec

    assemble!(vec,assem, bvec)

    assemble!(mat,assem, mmat)

    x2 = mat \ vec

    @test x ≈ x2

    @test vec ≈ [0.0625, 0.125, 0.0625]

    @test mat[1, 1]  ≈  1.333333333333333
    @test mat[2, 1]  ≈ -0.33333333333333
    @test mat[1, 2]  ≈ -0.33333333333333
    @test mat[2, 2]  ≈ 2.666666666666666
    @test mat[3, 2]  ≈ -0.33333333333333
    @test mat[2, 3]  ≈ -0.33333333333333
    @test mat[3, 3]  ≈ 1.333333333333333

    cellids = IdentityCellNumber(Int,length(bvec))

    vec2 = assemble(assem, (bvec,cellids))

    @test vec2 == vec

    vec3 = assemble(assem, (bvec,cellids), (bvec,cellids))
    @test vec3 ≈ 2*vec

    assemble!(vec3, assem, (bvec,cellids), (bvec,cellids))
    @test vec3 ≈ 2*vec

    mat2 = assemble(assem, (mmat,cellids,cellids) )

    @test mat2 == mat

    mat3 = assemble(assem, (mmat,cellids,cellids), (mmat,cellids,cellids))
    @test mat3 ≈ 2*mat

    assemble!(mat3, assem, (mmat,cellids,cellids), (mmat,cellids,cellids))
    @test mat3 ≈ 2*mat
end

end # module AssemblersTests


text:
# Using Convex.jl to test if the probability distributions in Homework problem 1
# Sets up the optimization problem for the Cut Inflation for the triangle causal
# model with the distribution:
#
# P(a, b, c) = if (a + b + c == 1) => 1/3
#              else                => 0

using Convex
using GLPK
using Inflation

# the probability matrix is set up as
# P_{A2,B1} = sum_{C1}P_{A2,B1,C1}   --> 'yyx'  (3rd term is integrated)
# P_{A2,C1} = sum_{B1}P_{A2,B1,C1}   --> 'yxy'  (2nd term is integrated)
# P_{B1,C1} = sum_{A2}P_{A2,B1,C1}   --> 'xyy'  (1st term is integrated)
infl_strings = ["yyx", "yxy", "xyy"]

# creating the M-matrix and b-array to be used in the Convex.jl optimization
M = get_M_matrix(infl_strings)
b = get_b_array(infl_strings, pfunc2, cut)

# the variables that get optimized by Convex.jl
v = get_v_array_Variable(infl_strings)

# set up the objective and constraints
objective   = 0
constraints = [v >= 0, M*v == b]

# create the problem and solve
problem = Convex.maximize(objective, constraints)
solve!(problem, GLPK.Optimizer)

# return the results
if problem.status == Convex.MathOptInterface.INFEASIBLE
	println("Solution is infeasible.")
else
	println("Solution is feasible.")
	println("objective = $(problem.optval)")
	println("v = $(v.value)")
end


text:
"""
# MRMP: Multi-Robot Motion Planning Library

## minimum examples

```jl
# 1. generate instance
ins = MRMP.gen_random_instance_StatePoint2D(;
    N = 5,
    rad = 0.1,
    num_obs = 3,
    rad_obs = 0.1,
)
config_init, config_goal, obstacles, ins_params... = ins         # example of ins_params: radius, base positions of arms

# 2. define functions
connect = gen_connect(config_init[1], obstacles, ins_params...)  # connection checker
collide = gen_collide(config_init[1], ins_params...)             # collision checker
check_goal = gen_check_goal(config_goal)                         # goal judge

# 3. solve
@time solution, roadmaps = MRMP.Solvers.SSSP(
    config_init,
    config_goal,
    connect,
    collide,
    check_goal;
    TIME_LIMIT = 10,
)
validate(config_init, connect, collide, check_goal, solution)    # check validity of solution

# 4. refine
(TPG, solution, cost) = smoothing(solution, connect, collide)
println(cost)

# 5. visualize
plot_anim!(config_init, config_goal, obstacles, ins_params...; solution=solution, interpolate_depth=2)
```

## model definition

A model X (in `./models`) is defined as follows

- `StateX <: AbsState` defines configuration space
- `get_mid_status(q1::StateX, q2::StateX)` returns intermediate state
- `dist(q1::StateX, q2::StateX)` returns distance between two states
- `gen_connect(q::StateX,  obstacles, other_params...)` generates connect check function
- `gen_collide(q::StateX, other_params...)` generates inter-agent collision check function
- `gen_uniform_sampling(q::StateX)` generates uniform sampling function for StateX
- `plot_motion!(q_from::StateX, q_to::StateX, other_params...)` plots motion of one agent
- `plot_start_goal!(q_init::StateX, q_goal::StateX, other_params...)` plots start and goal statues
- `plot_agent!(q::StateX, other_params...)` plots an agent
- `gen_random_instance_StateX(; params...)` generates random instance, see `./utils.jl`

"""
module MRMP

import DataStructures: PriorityQueue, enqueue!, dequeue!
import Printf: @printf, @sprintf
import Base: @kwdef
import Random: randperm
using LinearAlgebra: norm, dot, normalize

# parameters
const STEP_DIST = 0.01
const SAFETY_DIST_LINE = 0.01
const DUBINS_TURN_RADIUS = 1.0

# state (i.e., sampled point), defined in ./models
abstract type AbsState end

# vertex in roadmaps
mutable struct Node{State<:AbsState}
    q::State   # corresponding state
    id::Int64  # unique id
    neighbors::Vector{Int64}  # indexes
end

include("obstacle.jl")
include("smoother.jl")
include("viz.jl")
include("utils.jl")

# models
include("models/point.jl")
include("models/point2d.jl")
include("models/point3d.jl")
include("models/line2d.jl")
include("models/capsule3d.jl")
include("models/arm22.jl")
include("models/arm33.jl")
include("models/snake2d.jl")
include("models/dubins2d.jl")

# solvers
include("solvers/solvers.jl")

export gen_connect, gen_collide, gen_check_goal
export plot_res!, plot_anim!, plot_instance!, plot_tpg!
export smoothing, get_solution_cost, get_tpg_cost, is_valid_instance, validate
export gen_random_instance_StateArm22,
    gen_random_instance_StateArm33,
    gen_random_instance_StateDubins,
    gen_random_instance_StateLine2D,
    gen_random_instance_StatePoint2D,
    gen_random_instance_StatePoint3D,
    gen_random_instance_StateSnake2D

using .Solvers

end


text:
include("Grover.jl")
using Test, BitBasis

"""traditional grover search algorithm."""
function grover_search(oracle::AbstractBlock{N}, gen::AbstractBlock{N}=repeat(N,H,1:N)) where N
    reg = zero_state(N) |> gen
    for i = 1:num_grover_step(oracle, gen)
        grover_step!(reg, oracle, gen)
    end
    return reg
end

function grover_circuit(oracle::AbstractBlock{N}, gen::AbstractBlock{N}, niter::Int=num_grover_step(oracle, gen)) where {N}
    chain(N, chain(oracle, reflect_circuit(gen)) for i = 1:niter)
end

#################### Tests ##################

@testset "oracle" begin
    oracle = inference_oracle(3, [2,-1,3])
    # ≈ method for Identity/PermuteMultiply/Sparse
    # add and mimus between sparse matrices.
    # alway use sorted CSC format.
    v = ones(1<<3)
    v[Int(0b110)+1] *= -1
    @test mat(oracle) ≈ Diagonal(v)
    @test mat(chain(phase(π), Z)) ≈ -mat(Z)
end

@testset "Grover Search" begin
    ####### Construct Grover Search Using Reflection Block
    num_bit = 12
    oracle = inference_oracle(num_bit, push!(collect(Int, 1:num_bit-1), num_bit))

    psi = grover_search(oracle)
    target_state = zeros(1<<num_bit); target_state[end] = 1
    @test isapprox(abs(statevec(psi)'*target_state), 1, atol=1e-3)
end

@testset "groverblock" begin
    gen = repeat(5, H, 1:5)
    or = inference_oracle(5, [-1,2,5,4,3])
    gb = grover_circuit(or, gen)
    @test apply!(zero_state(5) |> gen, gb) == grover_search(or)
end

@testset "test inference" begin
    Random.seed!(2)
    num_bit = 12
    gen = dispatch!(variational_circuit(num_bit), :random)
    psi0 = zero_state(num_bit) |> gen
    #psi0 = uniform_state(num_bit)
    evidense = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    #evidense = collect(1:num_bit)

    # the desired subspace
    basis = collect(UInt, 0:1<<num_bit-1)
    subinds = [itercontrol(num_bit, abs.(evidense), Int.(evidense.>0))...]

    v_desired = statevec(psi0)[subinds .+ 1]
    p = norm(v_desired)^2
    v_desired[:] ./= sqrt(p)

    # search the subspace
    psi = grover_search(inference_oracle(num_bit, evidense), gen)
    @test isapprox((psi.state[subinds .+ 1]'*v_desired) |> abs2, 1, atol=3e-2)
end


text:
@testset "disparity" begin
    ft = @load_toyfairtensor
    M = [true_positive_rate, false_positive_rate, ppv]
    @test_throws ArgumentError disparity(M, ft)
    d = disparity(M, ft; refGrp="Education")
    @test names(d) == ["labels", "true_positive_rate_disparity", "false_positive_rate_disparity",
                        "positive_predictive_value_disparity"]
    @test size(d)==(3, 4)
end

@testset "Disparity" begin
    X, y = @load_toydata()
    model = ConstantClassifier()
    result = evaluate(model, X, y, measure=Disparity(tpr, grp=:Sex, refGrp="M", func=-))
    @test length(keys(result.per_fold[1][1])) == 3
    result = evaluate(model, X, y, measures=Disparities([tpr, fpr], grp=:Sex, refGrp="M", func=/))
    @test length(keys(result.per_fold[1][1])) == 3
    @test length(result.per_fold) == 2
end

@testset "Parity" begin
    ft = @load_toyfairtensor
    M = [true_positive_rate]
    df = disparity(M, ft; refGrp="Board")
    df = parity(df, 0.2)
    @test names(df)[3] == "true_positive_rate_parity"
    arr = df[:, :true_positive_rate_parity]
    @test arr[1] && !arr[2] && !arr[3]
end


text:
import Kuber, MbedTLS, Base, Swagger, YAML, Base64, HTTP

"""
    K8sClient(uri::String, sslconfig::MbedTLS.SSLConfig, headers::Dict{String, String})

Holds the configurations for making HTTP k8s connections.
"""
mutable struct K8sClient
    server::String
    sslconfig::MbedTLS.SSLConfig
    headers::Dict{String, String}
    function K8sClient(server::String, sslconfig::MbedTLS.SSLConfig; headers::Dict{String, String}=Dict{String, String}())
        new(server, sslconfig, headers)
    end
end

"""
    request(client::K8sClient, method::String, path::String)::HTTP.Messages.Response

Makes a request to the Kubernetes server using the given client on the method and path.
"""
function request(client::K8sClient, method::String, path::String)::HTTP.Messages.Response
    uri = HTTP.URI(client.server)
    uri = merge(uri; path=path)
    return HTTP.request(method, string(uri), headers=client.headers, sslconfig=client.sslconfig, status_exception=false)
end

function request(client::K8sClient, method::String, path::String, body::String)::HTTP.Messages.Response
    uri = HTTP.URI(client.server)
    uri = merge(uri; path=path)
    return HTTP.request(method, string(uri), client.headers, body, sslconfig=client.sslconfig, status_exception=false)
end

"""
    client()

Creates a K8sClient based on whether it detects running in a pod or locally.
"""
function client()::K8sClient
    envs = Base.EnvDict()
    host = get(envs, "KUBERNETES_SERVICE_HOST", "")
    if host == ""
        @info "detected running locally, attempting kubeconfig connection"
        return client("")
    end

    @info "detected running in pod, attempting in-cluster connection"
    return incluster_client()
end

"""
    client(kubeconfigpath::string)

Creates a K8sClient configured from a kubeconfig file. If kubeconfig path is
an empty string it will look for the KUBCONFIG envvar or default to ~/.kube/config
"""
function client(kubeconfigpath::String)::K8sClient
    envs = Base.EnvDict()

    if kubeconfigpath == ""
        default = Base.Filesystem.joinpath(Base.Filesystem.homedir(), ".kube/config")
        kubeconfigpath = get(envs, "KUBECONFIG", default)
    end
    @info "kubeconfig filepath: " kubeconfigpath
    kubecfg = YAML.load(open(kubeconfigpath))

    currentcontext = get(kubecfg, "current-context", "")
    if currentcontext == ""
        error("current kubeconfig context is empty")
    end

    @info "using current context: " currentcontext

    clustername = ""
    username = ""
    for context in kubecfg["contexts"]
        if context["name"] == currentcontext
            clustername = context["context"]["cluster"]
            username = context["context"]["user"]
        end
    end

    if clustername == "" || username == ""
        error(string("could not find user or cluster for current context: ", currentcontext))
    end

    # TODO: handle other auth scenarios
    cadata = ""
    server = ""
    for cluster in kubecfg["clusters"]
        if cluster["name"] == clustername
            cadata = cluster["cluster"]["certificate-authority-data"]
            server = cluster["cluster"]["server"]
        end
    end

    if server == ""
        error(string("could not find cluster for current context: ", currentcontext))
    end

    @info "using server: " server

    clientcert = ""
    clientkey = ""
    for user in kubecfg["users"]
        if user["name"] == username
            clientcert = user["user"]["client-certificate-data"]
            clientkey = user["user"]["client-key-data"]
        end
    end
    if clientcert == "" || clientkey == ""
        error("could not find client cert or key from current context: ", currentcontext)
    end

    @info "using user: " username

    # save ssl info as temp files
    # TODO: would prefer to find a way of doing this without temp files
    capath, io = Base.Filesystem.mktemp()
    write(io, Base64.base64decode(cadata))
    close(io)
    @debug "capath: " capath

    certpath, io = Base.Filesystem.mktemp()
    write(io, Base64.base64decode(clientcert))
    close(io)
    @debug "certpath: " certpath

    keypath, io = Base.Filesystem.mktemp()
    write(io, Base64.base64decode(clientkey))
    close(io)
    @debug "keypath: " keypath

    conf = MbedTLS.SSLConfig(certpath, keypath)
    MbedTLS.config_defaults!(conf)

    cacert = MbedTLS.crt_parse_file(capath)
    @debug "capath: " cacert
    MbedTLS.ca_chain!(conf, cacert)

    return K8sClient(server, conf)
end

"""
    incluster_client()

Creates a K8sClient configured from a pods enviornment.
"""
function incluster_client()::K8sClient
    tokenFile  = "/var/run/secrets/kubernetes.io/serviceaccount/token"
    rootCAFile = "/var/run/secrets/kubernetes.io/serviceaccount/ca.crt"

    envs = Base.EnvDict()

    host = get(envs, "KUBERNETES_SERVICE_HOST", "")
    port = get(envs, "KUBERNETES_SERVICE_PORT", "")

    if host == "" || port == ""
        error("could not find host and/or port from environment,
         must not be running in a pod")
    end

    token = read(tokenFile, String)
    conf = MbedTLS.SSLConfig()
    MbedTLS.ca_chain!(conf, MbedTLS.crt_parse_file(rootCAFile))

    uri = string("https://", host, ":", port)
    @info "using server: " uri

    # TODO: this needs fixed upstream to handle token refreshes
    headers = Dict("Authorization" => string("bearer ", token))

    return K8sClient(uri, conf; headers=headers)
end

"""
    kubercontext(client::K8sClient)::Kuber.KuberContext

Creates a kubercontext from a client
"""
function kubercontext(client::K8sClient)::Kuber.KuberContext
    kctx = Kuber.KuberContext()
    Kuber.set_server(kctx, client.uri, true, headers=client.headers, sslconfig=client.sslconfig, require_ssl_verification=true)
    return kctx
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
#=
hypograph of root-determinant of real symmetric or complex Hermitian
positive semidefinite matrices (svec scaled triangle)
(u, w) : W ⪰ 0, u ≤ rtdet(W)
where rtdet(W) = (∏ᵢ λᵢ(W))^(1/d)

dual cone is
(u, w) : u ≤ 0, W ⪰ 0, u ≥ -d * rtdet(W)
=#

mutable struct HypoRootdetTri{C <: RealCompF} <: Cache
    oa_s::Vector{AE}
    d::Int
    w_temp::Vector{RealF}
    W_temp::Matrix{C}
    HypoRootdetTri{C}() where {C <: RealCompF} = new{C}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.HypoRootdetTriCone{RealF, C},
    ::Optimizer,
) where {C <: RealCompF}
    cache = HypoRootdetTri{C}()
    cache.oa_s = oa_s
    dim = MOI.dimension(cone)
    d = cache.d = svec_side(C, dim - 1)
    cache.w_temp = zeros(RealF, dim - 1)
    cache.W_temp = zeros(C, d, d)
    return cache
end

function MOIPajarito.Cones.add_init_cuts(cache::HypoRootdetTri{C}, opt::Optimizer) where {C}
    # add variable bounds W_ii ≥ 0
    d = cache.d
    @views w = cache.oa_s[2:end]
    W_diag = [w[svec_idx(C, i, i)] for i in 1:d]
    JuMP.@constraint(opt.oa_model, W_diag .>= 0)
    opt.use_init_fixed_oa || return

    # add cut (-d, I) on (u, W), a linearization at W = I
    u = cache.oa_s[1]
    JuMP.@constraint(opt.oa_model, -d * u + sum(W_diag) >= 0)
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::HypoRootdetTri,
    opt::Optimizer,
)
    p = z[1]
    R = cache.W_temp
    @views svec_to_smat!(R, z[2:end], rt2)
    F = eigen(Hermitian(R, :U))
    V = F.vectors
    ω = F.values

    cuts = AE[]
    if p > -1e-7
        # add eigenvector cuts
        num_pos = count(>(1e-7), ω)
        @views V_neg = V[:, 1:num_pos] * Diagonal(sqrt.(ω[1:num_pos]))
        @views w = cache.oa_s[2:end]
        cuts = _get_psd_cuts(V_neg, w, cache, opt)
    end

    # add rootdet cut
    ω_pos = max.(ω, 1e-7)
    cut = _get_cut(ω_pos, V, cache, opt)
    push!(cuts, cut)
    return cuts
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::HypoRootdetTri,
    opt::Optimizer,
)
    us = s[1]
    Ws = cache.W_temp
    @views svec_to_smat!(Ws, s[2:end], rt2)
    F = eigen(Hermitian(Ws, :U))
    V = F.vectors
    ω = F.values
    num_neg = count(<(-opt.tol_feas), ω)
    if iszero(num_neg) && us < opt.tol_feas
        return AE[]
    end
    @assert issorted(ω)

    cuts = AE[]
    if !iszero(num_neg)
        # add eigenvector cuts
        V_neg = V[:, 1:num_neg]
        @views w = cache.oa_s[2:end]
        cuts = _get_psd_cuts(V_neg, w, cache, opt)
    end

    geomean(ω) - us > -opt.tol_feas && return AE[]
    # gradient cut is (-1, V * Diagonal(geom(ω) / d ./ ω) * V')
    ω_pos = max.(ω, 1e-7)
    rω = (geomean(ω_pos) / cache.d) ./ ω_pos
    cut = _get_cut(rω, V, cache, opt)
    push!(cuts, cut)
    return cuts
end

function _get_cut(
    rω::Vector{RealF},
    V::Matrix{C},
    cache::HypoRootdetTri{C},
    opt::Optimizer,
) where {C}
    # strengthened cut is (-d * geom(rω), V * Diagonal(rω) * V')
    p = -cache.d * geomean(rω)
    R = V * Diagonal(rω) * V'
    r = smat_to_svec!(cache.w_temp, R, rt2)
    u = cache.oa_s[1]
    @views w = cache.oa_s[2:end]
    return JuMP.@expression(opt.oa_model, p * u + JuMP.dot(r, w))
end


text:
using LightGraphs

"""
    get_features_instance(graph, path_matpower)

Return the dict of the instances features.
"""
function get_features_instance(graph, path_matpower)
    features = Dict("graph" => readf.getfeatures_graph(graph))
    merge!(features, Dict("OPF" => readf.getfeatures_OPF(path_matpower)))
    merge!(features, Dict("kernel" => readf.getfeatures_kernel(graph)))
    return features
end

"""
    get_features_decomposition(chordal_graph, nb_added_edges, cliques, cliquetree)

Return the dict of the decompositions features.
"""
function get_features_decomposition(chordal_graph, nb_added_edges, cliques, cliquetree)
    features = Mongoc.BSON("nb_added_edges_chordal_extension" => nb_added_edges)
    chordal_graph_features = Mongoc.BSON("chordal_graph" => readf.getfeatures_graph(chordal_graph))
    cliques_features = Mongoc.BSON("cliques" => readf.getfeatures_cliques(chordal_graph, cliques))
    kernel_graph_features = Mongoc.BSON("kernel" => readf.getfeatures_kernel(chordal_graph))
    merge!(features, chordal_graph_features, cliques_features, kernel_graph_features)
    return features
end

function recompute_chordal_graph_features(manager::ExperimentManager)
    MPI.Init()
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)
    instance_graphs = Dict()
    for instance in manager.instances
        instance_graphs[instance["_id"]] = construct_network_graph(instance["paths"]["matpower"])
    end
    if rank == 0
        @info "Size: $size"
        indexes = get_indexes(manager, :solve, size,
                              [true,false],
                              [true,false],
                              [true,false])
        for (i, index) in enumerate(indexes)
            @info "[$rank] Rank $(i - 1): $(index[begin]) -> $(index[end]) | $(length(index))"
            MPI.Isend(index, i - 1, 0, comm)
        end
        if length(indexes) < size
            for i in length(indexes) + 1:size
                MPI.Isend(Vector{Int}(undef,0), i - 1, 0, comm)
            end
        end
    end
    MPI.Barrier(comm)
    @info "[$rank] Retrieving the number of instances/decompositions to process"
    status = MPI.Probe(0, 0, comm)
    count = MPI.Get_count(status, Int)
    if count == 0
        @info "[$rank] Exciting because no instances/decompositions to process"
        MPI.Finalize()
        return
    end
    indexes = Array{Int}(undef, count)
    MPI.Irecv!(indexes, 0, 0, comm)
    @info "[$rank] Processing"
    decompositions = collect(manager.decompositions)
    for (i, id) in enumerate(indexes)
        d = decompositions[id]
        @info "[$rank] Processing $(d["_id"]["instance_name"]): $(d["_id"]["oid"]): $(i) / $(length(indexes))"
        graph = construct_graph_
        features = Mongoc.BSON("instance_graph" => graph_features,
                               "kernel" => kernel_features)
        db.setfeatures_decomposition!(manager.decompositions,
                                      d["_id"]["instance_name"],
                                      d["_id"]["oid"],
                                      features)
    end
    @info "[$rank] Process done."
    MPI.Finalize()
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
# This script generates a `.hpp` file (C++ header) containing
# an LDPC code stored in compressed sparse column (CSC) format.
# See command line help for how to use it.

using SparseArrays
using LinearAlgebra
using DelimitedFiles

using ArgParse # Must be installed (use the Pkg package manager)

header = raw"""
// This is an automatically generated file.
// It contains row indices that should be combined to reduce the number of syndrome bits
// of an LDPC code in a rate-adaptive manner.
// These row indices are obtained via optimization of the particular LDPC matrix.

#include <cstdint>
#include <array>

namespace AutogenRateAdapt {

"""

struct ValidPath
    val::String
    ValidPath(s::AbstractString) = isfile(s) ? new(s) : error("No such file: '$s'")
end

function parse_my_args(args)
    s = ArgParseSettings(
        "Create a C++ code-file containing line indices for combination,
        to be applied to a given LDPC code in order to achieve rate adaption."
        ; allow_ambiguous_opts=false)

    @add_arg_table! s begin
        "rate_adaption_file_path"
            required = true
            arg_type = ValidPath
            help = "Path to file containing the line indices for rate adaption."
        "--output_path"
            default = "autogen_rate_adaption.hpp"
            arg_type = String
            help = "Path to put the automatically generated C++ file containing the rate adaption."
        "--debug_use_only_part"
            help = "If true, only a small part of the provided data processed."
            action = :store_true
    end

    parsed_args = parse_args(args, s)
    println("Parsed args:")
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end

    return abspath(parsed_args["rate_adaption_file_path"].val),
            abspath(parsed_args["output_path"]),
            parsed_args["debug_use_only_part"]
end

function write_cpp_constexpr_rate_adaption(
    rate_adaption_rows::AbstractArray{T, 2} where T <: Integer,
    output_cpp_path::AbstractString
    ; only_debug_mode = false
    )
    size(rate_adaption_rows, 2) == 2 || throw(
        ArgumentError("Expecting input of shape N x 2, where N is the number of time that two rows are combined.
                        Received shape $(size(rate_adaption_rows))"))

    if only_debug_mode  # make input data smaller...
        cut_at = minimum([10, size(rate_adaption_rows)...])
        rate_adaption_rows = rate_adaption_rows[1:cut_at, :]
    end

    max_row_idx = maximum(rate_adaption_rows)
    if log2(max_row_idx) < 16
        row_idx_type = "std::uint16_t"
    elseif log2(max_row_idx) < 32
        row_idx_type = "std::uint32_t"
    elseif log2(max_row_idx) < 64
        row_idx_type = "std::uint64_t"
    else
        throw(ArgumentError("Input insanely large? Has $max_row_idx entries..."))
    end

    open(output_cpp_path, "w+") do f
        print(f, header)

        println(f, """
        constexpr std::size_t num_combined_rows = $(length(rate_adaption_rows));\n
        constexpr std::array<$row_idx_type, num_combined_rows> rows = {""")

        for (i, idx) in enumerate(transpose(rate_adaption_rows))
            print(f, "0x$(string(idx, base=16))")
            if i != length(rate_adaption_rows)
                print(f, ",")
            end
            if mod(i, 100) == 0
                println(f, "")  # for formatting.
            end
        end
        println(f, "\n};\n")

        println(f, "} // namespace RALDPC")
    end
end

function read_rate_adaption_file(file_path::AbstractString)
    return DelimitedFiles.readdlm(file_path, ',', Int)
end

function write_array_to_csv(arr::AbstractArray, out_path::AbstractString)
    open(out_path, "w") do io
       writedlm(io, arr, ',')
    end
end

function main(args)
    rate_adaption_file_path, output_path, only_debug_mode = parse_my_args(args)
    rate_adaption_rows = read_rate_adaption_file(rate_adaption_file_path)
    write_cpp_constexpr_rate_adaption(rate_adaption_rows, output_path; only_debug_mode)

    @info "Saved C++ file ($(Base.stat(output_path).size) bytes) at '$(output_path)'."
end

main(ARGS)


text:
puzzleinput = readlines("inputs/day17.txt")

initstate2D = map(reshape(vcat(collect.(puzzleinput)...), 8, 8)) do x
    x == '#'
end

const maxsz = 26

# First Half

const sum_states3D = tuple(false, false, false, true, zeros(Bool, 23)...),
tuple(false, false, true, true, zeros(Bool, 23)...)

initstate3D = reshape(fill(false, maxsz^3), maxsz, maxsz, maxsz)
view(initstate3D, maxsz ÷ 2 - 4:maxsz ÷ 2 + 3, maxsz ÷ 2 - 4:maxsz ÷ 2 + 3, maxsz ÷ 2) .= initstate2D
# Moves the initial state 8 x 8, as a 1-dim slice,
# to the center of the 3D initial state

function sim1(state, cycles)
    @inline function offsets3D()
        # Returns the Moore neighborhood in 3D as an iterable of offsets

        return ((i, j, k) for k in -1:1, j in -1:1, i in -1:1 if (i, j, k) != (0, 0, 0))
    end

    function hood(grid, i)
        # Return the state of the Moore neighbors of cell i in grid

        c = Tuple(i)
        return [grid[CartesianIndex(c .+ i)] for i in offsets3D()]
    end

    life(hood, cstate) = sum_states3D[cstate + 1][sum(hood) + 1]

    for cycle = 1:cycles
        next = reshape(fill(false, maxsz^3), maxsz, maxsz, maxsz)
        for i in CartesianIndices((2:(maxsz - 1), 2:(maxsz - 1), 2:(maxsz - 1)))
            # For every cell in the state calculates the life of next generation

            next[i] = life(hood(state, i), state[i])
        end
        state = next
    end
    return state
end

@time count(sim1(initstate3D, 6))

# Second Half

const sum_states4D = tuple(false, false, false, true, zeros(Bool, 77)...),
tuple(false, false, true, true, zeros(Bool, 77)...)

initstate4D = reshape(fill(false, maxsz^4), maxsz, maxsz, maxsz, maxsz)
view(initstate4D, maxsz ÷ 2 - 4:maxsz ÷ 2 + 3, maxsz ÷ 2 - 4:maxsz ÷ 2 + 3, maxsz ÷ 2, maxsz ÷ 3) .=
    initstate2D
# Moves the initial state 8 x 8, as a 1-dim slice,
# to the center of the 4D initial state

function sim2(state, cycles)
    @inline function offsets4D()
        # Returns the Moore neighborhood in 4D as an iterable of offsets

        return (
            (i, j, k, l) for l in -1:1, k in -1:1, j in -1:1, i in -1:1 if (i, j, k, l) != (0, 0, 0, 0)
        )
    end

    function hood(grid, cell)
        # Return the state of the Moore neighbors of cell in grid

        c = Tuple(cell)
        return [grid[CartesianIndex(c .+ i)] for i in offsets4D()]
    end

    life(hood, cstate) = sum_states4D[cstate + 1][sum(hood) + 1]

    for cycle = 1:cycles
        next = reshape(fill(false, maxsz^4), maxsz, maxsz, maxsz, maxsz)
        for i in CartesianIndices((2:(maxsz - 1), 2:(maxsz - 1), 2:(maxsz - 1), 2:(maxsz - 1)))
            # For every cell in the state calculates the life of next generation

            next[i] = life(hood(state, i), state[i])
        end
        state = next
    end
    return state
end

@time count(sim2(initstate4D, 6))


text:
module SmallFixedSizeSets

export Set0, Set1, Set2, Set3, Set4, Set5, SmallFixedSizeSet
export SFSSFromOrderedVec, SFSSFromVec, SFSSFromSet
export SFSS_issubset, SFSS_intersect, SFSS_setdiff, SFSS_union

const Set0 = NTuple{0,Int32}
const Set1 = NTuple{1,Int32}
const Set2 = NTuple{2,Int32}
const Set3 = NTuple{3,Int32}
const Set4 = NTuple{4,Int32}
const Set5 = NTuple{5,Int32}
const SmallFixedSizeSet = Union{Set0, Set1, Set2, Set3, Set4, Set5}

"""
SFSS_FromOrderedVec
----------

SFSSFromOrderedVec(v::Vector{Int32})

Create a small fixed-size set from a vector v

v must contain be of length at most 5
"""
function SFSSFromOrderedVec(v::Vector{Int32})
    if length(v) == 0; return Set0(); end
    if length(v) == 1; return Set1(v); end
    if length(v) == 2; return Set2(v); end
    if length(v) == 3; return Set3(v); end
    if length(v) == 4; return Set4(v); end
    if length(v) == 5; return Set5(v); end
    error("Set size is too large (maximum size is 5)")
end

SFSSFromOrderedVec(v::Vector{Int64}) =
    SFSSFromOrderedVec(convert(Vector{Int32}, v))
SFSSFromVec(v::Vector{Int64}) =
    SFSSFromOrderedVec(sort(v))
SFSSFromVec(v::Vector{Int32}) =
    SFSSFromOrderedVec(sort(v))
SFSSFromSet(S::Set{Int64}) =
    SFSSFromOrderedVec(sort(collect(S)))
SFSSFromSet(S::Set{Int32}) =
    SFSSFromOrderedVec(sort(collect(S)))

# Assumption: sets are sorted
function SFSS_issubset(A::SmallFixedSizeSet, B::SmallFixedSizeSet)
    nA, nB = length(A), length(B)
    Aind, Bind = 1, 1
    while Aind <= nA && Bind <= nB
        aval, bval = A[Aind], B[Bind]
        if   aval < bval;  return false
        else Bind += 1;
        end
        Aind += (aval == bval)
    end
    return Aind > nA
end

"""
SFSS_intersect
-----------------

SFSS_intersect(A::SmallFixedSizeSet, B::SmallFixedSizeSet)

Returns the intersect of small fixed-size sets A and B as a small fixed-size set.
"""
function SFSS_intersect(A::SmallFixedSizeSet, B::SmallFixedSizeSet)
    intersect_vec = Int32[]
    Aind, Bind = 1, 1
    nA, nB = length(A), length(B)
    while Aind <= nA && Bind <= nB
        aval, bval = A[Aind], B[Bind]
        if aval == bval
            push!(intersect_vec, aval)
            Aind += 1
            Bind += 1
        elseif aval < bval; Aind += 1;
        else                Bind += 1;
        end
    end
    return SFSSFromOrderedVec(intersect_vec)
end

"""
SFSS_setdiff
------------

SFSS_setdiff(A::SmallFixedSizeSet, B::SmallFixedSizeSet)

Returns the elements in A that are not also in B as a small fixed-size set.
"""
function SFSS_setdiff(A::SmallFixedSizeSet, B::SmallFixedSizeSet)
    setdiff_vec = Int32[]
    Aind, Bind = 1, 1
    nA, nB = length(A), length(B)
    while Aind <= nA && Bind <= nB
        aval, bval = A[Aind], B[Bind]
        if aval == bval
            Aind += 1
            Bind += 1
        elseif aval < bval
            push!(setdiff_vec, aval)
            Aind += 1
        else
            Bind += 1
        end
    end
    while Aind <= nA
        push!(setdiff_vec, A[Aind])
        Aind += 1
    end
    return SFSSFromOrderedVec(setdiff_vec)
end

"""
SFSS_union
--------------

SFSS_union(A::SmallFixedSizeSet, B::SmallFixedSizeSet)

Returns the union of small fixed-size sets A and B.
"""
function SFSS_union(A::SmallFixedSizeSet, B::SmallFixedSizeSet)
    union_vec = Int32[]
    Aind, Bind = 1, 1
    nA, nB = length(A), length(B)
    while Aind <= nA && Bind <= nB
        aval, bval = A[Aind], B[Bind]
        if aval == bval
            push!(union_vec, aval)
            Aind += 1
            Bind += 1
        elseif aval < bval
            push!(union_vec, aval)
            Aind += 1
        else
            push!(union_vec, bval)
            Bind += 1
        end
    end
    while Aind <= nA
        push!(union_vec, A[Aind])
        Aind += 1
    end
    while Bind <= nB
        push!(union_vec, B[Bind])
        Bind += 1
    end
    return SFSSFromOrderedVec(union_vec)
end

end # module


text:
module Classification
    using Polyester

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss = zeros(Int64, size(output_matrix)[2])
        @batch for b in axes(output_matrix, 2)
            if findmax(output_matrix[:,b])[2]!=findmax(sample_matrix[:,b])[2]
                loss[b] += 1
            end
        end
        return sum(loss)
    end
end


text:
module SimpleAssoc

import Base: copy, delete!, done, eachindex, eltype, empty!, get!,
    get, getkey, haskey, isempty, keys, length, next, pop!, setindex!,
    start, values

export AssocArray, AssocTuple
# export pop, push, setindex

################################################################################

"""
    type AssocArray{K,V}

A simple associative collection based on (unsorted) arrays

An `AssocArray` is a simple associative collection. Functionally it is
equivalent to a `Dict`. The major difference is that an `AssocArray`
uses a simpler internal representation, based on an unsorted array.
This can be faster than a `Dict` if the collection is small, or if the
collection is searched only infrequently.
"""
type AssocArray{K,V} <: Associative{K,V}
    elems::Vector{Pair{K,V}}
    AssocArray() = new(Pair{K,V}[])
    AssocArray(elem::Pair{K,V}) = new([elem])
    AssocArray(elems::Dict{K,V}) = new(collect(elems))
    AssocArray(a::AssocArray) = new(copy(a.elems))
end

function (::Type{AssocArray{K,V}}){K,V}(elem::Pair)
    AssocArray{K,V}(Pair{K,V}(K(elem.first) => V(elem.second)))
end
(::Type{AssocArray{K,V}}){K,V}(iter) = AssocArray{K,V}(Dict{K,V}(iter))
(::Type{AssocArray{K,V}}){K,V}(elems::Pair...) = AssocArray{K,V}(elems)

(::Type{AssocArray}){K,V}(elem::Pair{K,V}) = AssocArray{K,V}(elem)
function (::Type{AssocArray})(iter)
    dict = Dict(iter)
    K = keytype(dict)
    V = valtype(dict)
    AssocArray{K,V}(dict)
end
(::Type{AssocArray})(elems::Pair...) = AssocArray(elems)
(::Type{AssocArray}){K,V}(a::AssocArray{K,V}) = AssocArray{K,V}(a)

copy(a::AssocArray) = AssocArray(a)

eltype{K,V}(::AssocArray{K,V}) = Pair{K,V}

isempty(a::AssocArray) = isempty(a.elems)
length(a::AssocArray) = length(a.elems)

start(a::AssocArray) = start(a.elems)
done(a::AssocArray, st) = done(a.elems, st)
next(a::AssocArray, st) = next(a.elems, st)

immutable AssocArrayIter{I,K,V}
    arr::AssocArray{K,V}
end
eltype{I}(ai::AssocArrayIter{I}) = fieldtype(eltype(ai.arr), I)
isempty(ai::AssocArrayIter) = isempty(ai.arr)
length(ai::AssocArrayIter) = length(ai.arr)
start(ai::AssocArrayIter) = start(ai.arr)
done(ai::AssocArrayIter, st) = done(ai.arr, st)
function next{I}(ai::AssocArrayIter{I}, st)
    it, st = next(ai.arr, st)
    getfield(it, I), st
end
keys{K,V}(a::AssocArray{K,V}) = AssocArrayIter{:first,K,V}(a)
values{K,V}(a::AssocArray{K,V}) = AssocArrayIter{:second,K,V}(a)
eachindex(a::AssocArray) = keys(a)

empty!(a::AssocArray) = empty!(a.elems)

function pop!(a::AssocArray)
    if isempty(a)
        throw(BoundsError(a))
    end
    pop!(a.elems)
end

function getkey(a::AssocArray, k, d)
    for el in a.elems
        if el.first === k
            return el.first
        end
    end
    d
end

function get(a::AssocArray, k, d)
    for el in a.elems
        if el.first === k
            return el.second
        end
    end
    d
end

function get!(a::AssocArray, k, d)
    for el in a.elems
        if el.first === k
            return el.second
        end
    end
    a[k] = d
end

function setindex!{K,V}(a::AssocArray{K,V}, v, k)
    for i in eachindex(a.elems)
        if a.elems[i].first === k
            a.elems[i] = Pair{K,V}(K(k), V(v))
            return v
        end
    end
    push!(a.elems, Pair{K,V}(K(k), V(v)))
    v
end

function delete!(a::AssocArray, k)
    for i in eachindex(a.elems)
        if a.elems[i].first === k
            if i < length(a.elems)
                a.elems[i] = pop!(a.elems)
            else
                pop!(a.elems)
            end
            return a
        end
    end
    a
end

################################################################################

# immutable UnsafeAssertIsUnique end

"""
    type AssocTuple{K,V,N}

A simple associative collection based on (unsorted) tuples

An `AssocTuple` is a simple associative collection. Functionally it is
equivalent to a immutable `Dict`. The major difference is that an
`AssocTuple` uses a simpler internal representation, based on an
unsorted tuple. This can be faster than a `Dict` if the collection is
small, or if the collection is searched only infrequently.
"""
immutable AssocTuple{K,V,N} <: Associative{K,V}
    elems::NTuple{N, Pair{K,V}}
    function AssocTuple()
        @assert N == 0
        new(())
    end
    function AssocTuple(elem::Pair{K,V})
        @assert N == 1
        new((elem,))
    end
    function AssocTuple(elems::Dict{K,V})
        @assert N == length(elems)
        new((elems...))
    end
    # function AssocTuple(::Type{UnsafeAssertIsUnique},
    #                     elems::NTuple{N, Pair{K,V}})
    #     new(elems)
    # end
end

function (::Type{AssocTuple{K,V,N}}){K,V,N}(elem::Pair)
    AssocTuple{K,V,N}(K(elem.first) => V(elem.second))
end

(::Type{AssocTuple{K,V}}){K,V}() = AssocTuple{K,V,0}()
(::Type{AssocTuple{K,V}}){K,V}(elem::Pair) = AssocTuple{K,V,1}(elem)
function (::Type{AssocTuple{K,V}}){K,V}(iter)
    dict = Dict{K,V}(iter)
    N = length(dict)
    AssocTuple{K,V,N}(dict)
end
(::Type{AssocTuple{K,V}}){K,V}(elems::Pair...) = AssocTuple{K,V}(elems)

(::Type{AssocTuple})() = AssocTuple{Any,Any}()
(::Type{AssocTuple}){K,V}(elem::Pair{K,V}) = AssocTuple{K,V,1}(elem)
function (::Type{AssocTuple})(iter)
    dict = Dict(iter)
    K = keytype(dict)
    V = valtype(dict)
    N = length(dict)
    AssocTuple{K,V,N}(dict)
end
(::Type{AssocTuple})(elems::Pair...) = AssocTuple(elems)

eltype{K,V}(::AssocTuple{K,V}) = Pair{K,V}

isempty(a::AssocTuple) = isempty(a.elems)
length(a::AssocTuple) = length(a.elems)

start(a::AssocTuple) = start(a.elems)
done(a::AssocTuple, st) = done(a.elems, st)
next(a::AssocTuple, st) = next(a.elems, st)

immutable AssocTupleIter{I,K,V}
    tup::AssocTuple{K,V}
end
eltype{I}(ai::AssocTupleIter{I}) = fieldtype(eltype(ai.tup), I)
length(ai::AssocTupleIter) = length(ai.tup)
start(ai::AssocTupleIter) = start(ai.tup)
done(ai::AssocTupleIter, st) = done(ai.tup, st)
function next{I}(ai::AssocTupleIter{I}, st)
    it, st = next(ai.tup, st)
    getfield(it, I), st
end
keys{K,V}(a::AssocTuple{K,V}) = AssocTupleIter{:first,K,V}(a)
values{K,V}(a::AssocTuple{K,V}) = AssocTupleIter{:second,K,V}(a)
eachindex(a::AssocTuple) = keys(a)

# function pop{K,V,N}(a::AssocTuple{K,V,N})
#     if isempty(a)
#         throw(BoundsError(a))
#     end
#     a.elems[end], AssocTuple{K,V,N-1}(UnsafeAssertIsUnique, a.elems[1:end-1])
# end

function getkey(a::AssocTuple, k, d)
    for el in a.elems
        if el.first === k
            return el.first
        end
    end
    d
end

function get(a::AssocTuple, k, d)
    for el in a.elems
        if el.first === k
            return el.second
        end
    end
    d
end

# function get!(a::AssocTuple, k, d)
#     for el in a.elems
#         if el.first === k
#             return el.second
#         end
#     end
#     a[k] = d
# end

# function setindex!(a::AssocTuple, v, k)
#     for i in eachindex(a.elems)
#         if a.elems[i].first === k
#             a.elems[i] = k => v
#             return v
#         end
#     end
#     push!(a.elems, k => v)
#     v
# end

# push!(a, elem) = a[elem.first] = elem.second

# function delete!(a::AssocTuple, k)
#     for i in eachindex(a.elems)
#         if a.elems[i].first === k
#             if i < length(a.elems)
#                 a.elems[i] = pop!(a.elems)
#             else
#                 pop!(a.elems)
#             end
#             return a
#         end
#     end
#     a
# end

end


text:
"""
Public holidays for the Australian Stock Exchange (ASX).
"""
struct AustraliaASX <: HolidayCalendar end

"""
Public holidays for the Australian states and territories.
Although some holidays are common to all states and territories, such as Christmas Day, each state and territory also has its own additional holidays.
Therefore the set of relevant holidays depends on which state/territory you are concerned with.

The Australian states and territories are:
- Australian Capital Territory (ACT)
- New South Wales (NSW)
- Northern Territory (NT)
- Queensland (QLD)
- South Australia (SA)
- Tasmania (TAS)
- Western Australia (WA)
- Victoria (VIC)

For example: cal = Australia(:QLD)
"""
struct Australia <: HolidayCalendar
    state::Symbol

    function Australia(state::Symbol)
        states = Set([:ACT, :NSW, :NT, :QLD, :SA, :TAS, :WA, :VIC])
        @assert state ∈ states "$(state) is not a valid Australian state or territory. Must be one of :ACT, :NSW, :NT, :QLD, :SA, :TAS, :WA, :VIC."
        new(state)
    end
end

Base.:(==)(a1::Australia, a2::Australia) = a1.state == a2.state
Base.hash(a::Australia) = hash(a.state)

function isholiday(::AustraliaASX, dt::Dates.Date)
    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)
    easter_sunday = BusinessDays.easter_date(Dates.Year(yy))
    is_australian_national_holiday(dt, yy, mm, dd, easter_sunday) && return true
    mm == 6 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    false
end

function isholiday(cal::Australia, dt::Dates.Date)
    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)
    easter_sunday = BusinessDays.easter_date(Dates.Year(yy))
    is_australian_national_holiday(dt, yy, mm, dd, easter_sunday) && return true
    is_australian_state_holiday(Val{cal.state}, dt, yy, mm, dd, easter_sunday)
end

################################################################################
function is_australian_national_holiday(dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    mm == 1  && dd == 1  && return true  # New year's day
    mm == 1  && dd == 26 && return true  # Australia Day
    mm == 4  && dd == 25 && return true  # ANZAC Day
    mm == 12 && dd == 25 && return true  # Christmas Day
    mm == 12 && dd == 26 && return true  # Boxing Day
    dt == easter_sunday - Dates.Day(2) && return true  # Good Friday
    dt == easter_sunday + Dates.Day(1) && return true  # Easter Monday
    false
end

function is_australian_state_holiday(::Type{Val{:ACT}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday
    dt == easter_sunday          && return true  # Easter Sunday
    mm == 3  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Canberra Day (2nd Monday of March)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 10 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Labour Day (1st Monday of October)
    if yy >= 2018 && mm == 5  # Reconciliation Day (Last Monday of May)
        first_june   = Dates.Date(yy, 6, 1)
        last_mon_may = Dates.toprev(first_june, Dates.Mon)
        dt == last_mon_may && return true
    end
    false
end

function is_australian_state_holiday(::Type{Val{:NSW}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday
    dt == easter_sunday          && return true  # Easter Sunday
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 8  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Bank holiday (1st Monday of August)
    mm == 10 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Labour Day (1st Monday of October)
    false
end

function is_australian_state_holiday(::Type{Val{:NT}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday (Easter Sunday is not a public holiday in the Northern Territory)
    mm == 5  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # May Day (1st Monday of May)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 8  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Picnic Day (1st Monday of August)
    false
end

function is_australian_state_holiday(::Type{Val{:QLD}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday
    dt == easter_sunday          && return true  # Easter Sunday
    mm == 5  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Labour Day (1st Monday of May)

    # Royal Brisbane Show (Brisbane area only).
    if mm == 8 && Dates.dayofweek(dt) == Dates.Wed
        end_aug      = Dates.Date(yy, 8, 31)
        last_wed_aug = Dates.dayofweek(end_aug) == Dates.Wed ? end_aug : Dates.toprev(end_aug, Dates.Wed)
        n_wed_in_aug = Dates.day(last_wed_aug) >= 29 ? 5 : 4
        n_wed_in_aug == 4 && Dates.dayofweekofmonth(dt) == 2 && return true  # 2nd Wednesday of August if there are 4 Wednesdays in August
        n_wed_in_aug == 5 && Dates.dayofweekofmonth(dt) == 3 && return true  # 3rd Wednesday of August if there are 5 Wednesdays in August
    end

    # Queen's Birthday holiday
    if yy == 2012 || yy >= 2016
        mm == 10 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # 1st Monday of October
    else
        mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # 2nd Monday of June
    end
    false
end

function is_australian_state_holiday(::Type{Val{:SA}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday (Easter Sunday is not a public holiday in South Australia)
    mm == 3  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # March public holiday (2nd Monday of March)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 10 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Labour Day (1st Monday of October)
    false
end

function is_australian_state_holiday(::Type{Val{:TAS}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    # Neither Easter Saturday nor Easter Sunday are public holidays in Tasmania
    mm == 2  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Royal Hobart Regatta (2nd Monday of February)
    mm == 3  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Labour Day (2nd Monday of March)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 11 && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Recreation Day (1st Monday of November)
    false
end

function is_australian_state_holiday(::Type{Val{:WA}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    # Neither Easter Saturday nor Easter Sunday are public holidays in Western Australia
    mm == 3  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Labour Day (1st Monday of March)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 1 && return true # Western Australia Day (1st Monday of June)

    # Queen's Birthday holiday (proclaimed each year by the Governor, usually the last Monday of September)
    if yy == 2011
        mm == 10 && dd == 28 && return true
    elseif yy == 2012
        mm == 10 && dd == 1 && return true
    else
        end_sep      = Dates.Date(yy, 9, 30)
        last_mon_sep = Dates.dayofweek(end_sep) == Dates.Mon ? end_sep : Dates.toprev(end_sep, Dates.Mon)
        dt == last_mon_sep && return true
    end
    false
end

function is_australian_state_holiday(::Type{Val{:VIC}}, dt::Dates.Date, yy::Int, mm::Int, dd::Int, easter_sunday::Dates.Date)
    dt == easter_sunday - Dates.Day(1) && return true  # Easter Saturday
    dt == easter_sunday          && return true  # Easter Sunday
    mm == 3  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Labour Day (2nd Monday of March)
    mm == 6  && Dates.dayofweek(dt) == Dates.Mon && Dates.dayofweekofmonth(dt) == 2 && return true # Queen's Birthday holiday (2nd Monday of June)
    mm == 11 && Dates.dayofweek(dt) == Dates.Tue && Dates.dayofweekofmonth(dt) == 1 && return true # Melbourne Cup (1st Tuesday of November)
    if yy >= 2015 && Dates.dayofweek(dt) == Dates.Fri && (mm == 9 || mm == 10)               # Friday before AFL Grand Final (Friday closest to 30th September)
        first_oct     = Dates.Date(yy, 10, 1)  # 1st October
        first_sat_oct = Dates.dayofweek(first_oct) == Dates.Sat ? first_oct     : Dates.tonext(first_oct, Dates.Sat)  # 1st Saturday in October
        afl_final     = Dates.day(first_sat_oct)   <= 4   ? first_sat_oct : first_sat_oct - Dates.Day(7)
        dt == afl_final - Dates.Day(1) && return true
    end
    false
end


text:
#=
Simple support for FFTs using:
https://github.com/JuliaMath/AbstractFFTs.jl

Does not (yet) cover plan_fft & friends,
because extracting the dimensions from those is tricky
=#

using AbstractFFTs

for fun in [:fft, :ifft, :bfft, :rfft]
    @eval function AbstractFFTs.$fun(A::Union{KeyedArray,NdaKa}, dims = ntuple(+,ndims(A)))
        numerical_dims = NamedDims.dim(A, dims)
        data = $fun(keyless(A), numerical_dims)
        keys = fft_keys($fun, A, numerical_dims, data)
        KeyedArray(data, keys)
    end
end

function AbstractFFTs.irfft(A::Union{KeyedArray,NdaKa}, len::Integer, dims = ntuple(+,ndims(A)))
    numerical_dims = NamedDims.dim(A, dims)
    data = irfft(keyless(A), len, numerical_dims)
    keys = fft_keys(irfft, A, numerical_dims, data)
    KeyedArray(data, keys)
end

for shift in [:fftshift, :ifftshift]
    @eval function AbstractFFTs.$shift(A::Union{KeyedArray,NdaKa}, dims = ntuple(+,ndims(A)))
        numerical_dims = NamedDims.dim(A, dims)
        data = $shift(keyless(A), numerical_dims)
        keys = ntuple(ndims(A)) do d
            k = axiskeys(A,d)
            d in numerical_dims || return k
            $shift(k)
        end
        KeyedArray(data, keys)
    end
end

# copy(fftfreq(10, 1)) isa Vector # perhaps that needs a method:
# Base.copy(x::Frequencies{<:Number}) = x  # fixed in AbstractFFts 1.0

fft_keys(f, A, dims, B) = ntuple(ndims(A)) do d
        k = axiskeys(A,d)
        d in dims || return k
        if k isa Frequencies
            return fft_un_freq(k)
        elseif k isa AbstractRange{<:Number}
            return fftfreq(length(k), inv(step(k)))
        end
        return k
    end

fft_keys(::typeof(rfft), A, dims, B) = ntuple(ndims(A)) do d
        k = axiskeys(A,d)
        d in dims || return k
        if k isa Frequencies || k isa AbstractRange{<:Number}
            return rfftfreq(length(k), inv(step(k)))
        end
        return axes(B,d)
    end

fft_keys(::typeof(irfft), A, dims, B) = ntuple(ndims(A)) do d
        k = axiskeys(A,d)
        d in dims || return k
        if k isa Frequencies || k isa AbstractRange{<:Number}
            return irfft_un_freq(k, size(B,d))
        end
        return axes(B,d)
    end

# This just chooses a nicer zero when inverting FFT:
function fft_un_freq(x::Frequencies)
    s = inv(x.multiplier * x.n)
    range(zero(s), step = s, length = x.n)
end

function irfft_un_freq(x, len)
    s = inv(step(x) * len)
    range(zero(s), step = s, length = len)
end


text:
using LightGraphs, MetaGraphs, Missings
using LightGraphs: SimpleEdge

struct SimpleNode{RT, NL} <: AbstractNode{RT, NL}
    id::Int
end
==(e1::SimpleNode, e2::SimpleNode) = (e1.id == e2.id)

struct SimpleBranch{RT, NL} <: AbstractBranch{RT, NL}
    src::Int
    dst::Int
end
SimpleBranch{RT, NL}(edge::SimpleEdge{Int}) where {RT, NL} =
    SimpleBranch{RT, NL}(src(edge), dst(edge))
==(e1::SimpleBranch, e2::SimpleBranch) = (e1.src == e2.src && e1.dst == e2.dst)

abstract type MetaTree{RT, NL, TD} <:
    AbstractTree{OneTree, RT, NL, SimpleNode{RT, NL}, SimpleBranch{RT, NL}} end

_noderecordtype(::Type{<:MetaTree}) = Dict{Symbol, Any}
_branchrecordtype(::Type{<:MetaTree}) = Dict{Symbol, Any}

mutable struct UnrootedMetaTree{NL, TD} <: MetaTree{Unrooted, NL, TD}
    mg::MetaGraph{Int, Float64}
    isvalid::Union{Bool, Missing}
    name::Union{String, Missing}
    tipdata::TD
end

mutable struct RootedMetaTree{RT <: Rooted, NL, TD} <: MetaTree{RT, NL, TD}
    mg::MetaDiGraph{Int, Float64}
    roots::Vector{SimpleNode{RT, NL}}
    isvalid::Union{Bool, Missing}
    name::Union{String, Missing}
    tipdata::TD
end

const RootedTree = RootedMetaTree{OneRoot, String, Dict{String, Any}}
const ManyRootTree = RootedMetaTree{ManyRoots, String, Dict{String, Any}}
const UnrootedTree = UnrootedMetaTree{String, Dict{String, Any}}
function RootedTree(n::Int = 0)
    tree = RootedTree(MetaDiGraph(), Vector{SimpleNode{OneRoot, String}}(),
                      missing, missing, Dict{String, Any}())
    set_indexing_prop!(tree.mg, :name)
    createnodes!(tree, n)

    return tree
end
function RootedTree(names::Vector{String})
    tree = RootedTree(MetaDiGraph(), Vector{SimpleNode{OneRoot, String}}(),
                      missing, missing, Dict{String, Any}())
    set_indexing_prop!(tree.mg, :name)
    createnodes!(tree, names)

    return tree
end

import Phylo.API._nnodes
_nnodes(tree::MetaTree) = nv(tree.mg)
import Phylo.API._getnodes
_getnodes(tree::T) where {RT, NL, TD, T <: MetaTree{RT, NL, TD}} =
    SimpleNode{RT, NL}.(vertices(tree.mg))
import Phylo.API._getnodenames
_getnodenames(tree::MetaTree) =
    String[_getnodename(tree, v) for v in _getnodes(tree)]
import Phylo.API._getnode
_getnode(tree::MetaTree{RT, NL, TD}, nodename::NL) where {RT, NL, TD} =
    SimpleNode{RT, NL}(tree[nodename, :name])
_getnode(tree::MetaTree{RT, NL, TD}, nodename::NL) where {RT, NL<:Integer, TD} =
    SimpleNode{RT, NL}(tree[:name][nodename])
import Phylo.API._getnodename
_getnodename(tree::MetaTree, node::SimpleNode) = tree.mg[node.id, :name]
import Phylo.API._hasnode
_hasnode(tree::MetaTree, node::SimpleNode) = 1 ≤ node.id ≤ _nnodes(tree)
_hasnode(tree::MetaTree, id::Int) = 1 ≤ id ≤ _nnodes(tree)

import Phylo.API._nbranches
_nbranches(tree::MetaTree) = ne(tree.mg)
import Phylo.API._getbranches
_getbranches(tree::T) where {RT, NL, TD, T <: MetaTree{RT, NL, TD}} =
    SimpleBranch{RT, NL}.(edges(tree.mg))
import Phylo.API._getbranchnames
_getbranchnames(tree::MetaTree) =
    [node + (idx - 1) * _nnodes(tree)
     for node in vertices(tree)
     for idx in Base.OneTo(length(fadj(tree, node)))]
import Phylo.API._getbranch
function _getbranch(tree::MetaTree{RT, NL, TD}, id::Int) where {RT, NL, TD}
    node = id % _nnodes(tree)
    idx = id ÷ _nnodes(tree)
    return SimpleBranch{RT, NL}(node, fadj(tree, node)[idx])
end
import Phylo.API._getbranchname
function _getbranchname(tree::MetaTree, branch::SimpleBranch)
    idx = findfirst(dest -> dest == branch.dst, fadj(tree, branch.src))
    idx > 0 || error("Branch not found")
    return branch.src + (idx - 1) * _nnodes(tree)
end
import Phylo.API._hasbranch
function _hasbranch(tree::MetaTree, branch::SimpleBranch)
    idx = findfirst(dest -> dest == branch.dst, fadj(tree, branch.src))
    return idx > 0
end
function _hasbranch(tree::MetaTree{RT, NL, TD}, id::Int) where {RT, NL, TD}
    node = id % _nnodes(tree)
    idx = id ÷ _nnodes(tree)
    return SimpleBranch{RT, NL}(node, fadj(tree, node)[idx])
end

import Phylo.API._src
_src(tree::RootedMetaTree, branch::SimpleBranch{<:Rooted}) =
     _getnode(tree, branch.src)
_src(tree::RootedMetaTree, id::Int) =
     _getnodename(tree, _getbranch(tree, id).src)

 import Phylo.API._dst
_dst(tree::RootedMetaTree, branch::SimpleBranch{<:Rooted}) =
     _getnode(tree, branch.dst)
_dst(tree::RootedMetaTree, id::Int) =
     _getnodename(tree, _getbranch(tree, id).dst)

 import Phylo.API._getroots
_getroots(tree::RootedMetaTree) = tree.roots

import Phylo.API._indegree
_indegree(tree::RootedMetaTree, node) =
    length(badj(tree.mg, _getnode(tree, node).id))
import Phylo.API._outdegree
_outdegree(tree::RootedMetaTree, node) =
    length(fadj(tree.mg, _getnode(tree, node).id))
import Phylo.API._degree
_degree(tree::UnrootedMetaTree, node) =
    length(fadj(tree.mg, _getnode(tree, node).id))

import Phylo.API._getparent
_getparent(tree::RootedMetaTree{RT, NL, TD},
           node::SimpleNode{RT, NL}) where {RT, NL, TD} =
    SimpleNode{RT, NL}(first(badj(tree.mg, node.id)))
_getparent(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    _getnodename(tree, _getparent(tree, _getnode(tree, node)))
import Phylo.API._getinbound
_getinbound(tree::RootedMetaTree{RT, NL, TD},
            node::SimpleNode{RT, NL}) where {RT, NL, TD} =
    SimpleBranch{RT, NL}(first(badj(tree.mg, node.id)), node.id)
_getinbound(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    _getbranchname(tree, _getinbound(tree, _getnode(tree, node)))

import Phylo.API._getchildren
_getchildren(tree::RootedMetaTree{RT, NL, TD},
              node::SimpleNode{RT, NL}) where {RT, NL, TD} =
    [SimpleNode{RT, NL}(out) for out in fadj(tree.mg, node.id)]
_getchildren(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    [_getnodename(tree, out)
     for out in _getchildren(tree, _getnode(tree, node))]
import Phylo.API._getoutbounds
_getoutbounds(tree::RootedMetaTree{RT, NL, TD},
              node::SimpleNode{RT, NL}) where {RT, NL, TD} =
    [SimpleBranch{RT, NL}(node, out)
     for out in fadj(tree.mg, _getnode(tree, node).id)]
_getoutbounds(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    [_getbranchname(tree, branch)
     for branch in _getoutbounds(tree, _getnode(tree, node))]

import Phylo.API._getsiblings
_getsiblings(tree::RootedMetaTree{RT, NL, TD},
             node::SimpleNode{RT, NL}) where {RT, NL, TD} =
    [SimpleNode{RT, NL}(out) for out in fadj(tree.mg, node.id)]
_getsiblings(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    [_getnodename(tree, out)
     for out in _getchildren(tree, _getnode(tree, node))]
import Phylo.API._getconnections
_getconnections(tree::UnrootedMetaTree{NL, TD}, node) where {NL, TD} =
    SimpleNode{Unrooted, NL}.(fadj(tree.mg, _getnode(tree, node).id))
_getconnections(tree::RootedMetaTree{RT, NL, TD}, node::NL) where {RT, NL, TD} =
    [_getbranchname(tree, branch)
     for branch in _getconnections(tree, _getnode(tree, node))]

import Phylo.API._addinbound!
_addinbound!(::RootedMetaTree, node) = Nothing()
import Phylo.API._removeinbound!
_removeinbound!(::RootedMetaTree, node) = Nothing()
import Phylo.API._addoutbound!
_addoutbound!(::RootedMetaTree, node) = Nothing()
import Phylo.API._removeoutbound!
_removeoutbound!(::RootedMetaTree, node) = Nothing()
import Phylo.API._addconnection!
_addconnection!(::UnrootedMetaTree, node) = Nothing()
import Phylo.API._removeconnection!
_removeconnection!(::UnrootedMetaTree, node) = Nothing()

import Phylo.API._createbranch!
function _createbranch!(tree::MetaTree{RT, NL, TD},
                        source::SimpleNode{RT, NL},
                        dest::SimpleNode{RT, NL},
                        length::Float64, data) where {RT, NL, TD}
    if ismissing(data)
        add_edge!(tree, source.id, dest.id)
        add_prop!(tree, source.id, dest.id, :length, length)
    else
        data[:length] = length
        add_edge!(tree, source.id, dest.id, data)
    end
    branch = SimpleBranch{RT, NL}(source.id, dest.id)
    filter!(n -> n != dest.id, tree.roots)
    return branch
end
_createbranch!(tree::MetaTree{RT, NL, TD}, source::NL, dest::NL,
               length::Float64, data) where {RT, NL, TD} =
    _createbranch!(tree, _getnode(tree, source), _getnode(tree, dest),
                   length, data)

import Phylo.API._removebranch!
function _removebranch!(tree::MetaTree{RT, NL, TD},
                        source::SimpleNode{RT, NL},
                        dest::SimpleNode{RT, NL},
                        length::Float64, data) where {RT, NL, TD}
    rem_edge!(tree, source.id, dest.id, data)
    push!(tree.roots, dest)
end
_removebranch!(tree::MetaTree{RT, NL, TD}, source::NL, dest::NL) where
    {RT, NL, TD} =
    _removebranch!(tree, _getnode(tree, source), _getnode(tree, dest))

import Phylo.API._createnode!
function _createnode!(tree::T, nodename::NL, data) where
    {RT, NL, TD, T <: MetaTree{RT, NL, TD}}
    if ismissing(data)
        add_vertex!(tree.mg, noderecordtype(T)())
    else
        add_vertex!(tree.mg, data)
    end
    node = SimpleNode{RT, NL}(nv(tree.mg))
    set_indexing_prop!(tree.mg, node.id, :name, nodename)
    return node
end

import Phylo.API._deletenode!
function _deletenode!(tree::RootedMetaTree, node::SimpleNode)
    rem_vertex!(tree.mg, node.id)
    filter!(n -> n != node.id, tree.roots)
end
function _deletenode!(tree::UnrootedMetaTree, node::SimpleNode)
    rem_vertex!(tree.mg, node.id)
end
_deletenode!(tree::MetaTree{RT, NL, TD}, nodename::NL) where {RT, NL, TD} =
    _deletenode!(tree, _getnode(tree, nodename))


text:
using CommonRLInterface
using Documenter

makedocs(;
    modules=[CommonRLInterface],
    authors="Zachary Sunberg <sunbergzach@gmail.com> and contributors",
    repo="https://github.com/JuliaReinforcementLearning/CommonRLInterface.jl/blob/{commit}{path}#L{line}",
    sitename="CommonRLInterface.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaReinforcementLearning.github.io/CommonRLInterface.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "required.md",
        "multiplayer.md",
        "optional.md",
        "wrappers.md",
        "faqs.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaReinforcementLearning/CommonRLInterface.jl",
)