text:
""" Contains data for a constraint that is repeated. Add to GlobalModel using `add_linked_constraint`."""
@with_kw mutable struct LinkedClassifier
    vars::Array{JuMP.VariableRef,1}                    # JuMP variables (flat)
    relax_var::Union{Real, JuMP.VariableRef} = 0.      # slack variable
    mi_constraints::Dict = Dict{Int64, Array{JuMP.ConstraintRef}}()
    leaf_variables::Dict = Dict{Int64, Tuple{JuMP.VariableRef, Array}}()
    active_leaves::Array = []                          # leaf of last solution
    feas_gap::Array = []                               # Feasibility gaps of solutions
    equality::Bool = false
end

function Base.show(io::IO, lc::LinkedClassifier)
    println(io, "LinkedClassifier with: ")
    println(io, "variables: $(lc.vars) ") # TODO: improve printing
end

""" Contains data for a constraint that is repeated. Add to GlobalModel using `add_linked_constraint`."""
@with_kw mutable struct LinkedRegressor
    vars::Array{JuMP.VariableRef,1}                    # JuMP variables (flat)
    dependent_var::JuMP.VariableRef                    # Dependent variable
    relax_var::Union{Real, JuMP.VariableRef} = 0.      # slack variable
    mi_constraints::Dict = Dict{Int64, Array{JuMP.ConstraintRef}}() # and their corresponding MI constraints,
    leaf_variables::Dict = Dict{Int64, Tuple{JuMP.VariableRef, Array, JuMP.VariableRef}}() # and leaf variables.
    active_leaves::Array = []                          # leaf of last solution
    optima::Array = []
    actuals::Array = []
    feas_gap::Array = []                               # Feasibility gaps of solutions
    equality::Bool = false
end

function Base.show(io::IO, lr::LinkedRegressor)
    println(io, "LinkedRegressor with: ")
    println(io, "variables: $(lr.vars) ") # TODO: improve printing
    println(io, "and dependent variable: $(lr.dependent_var)")
end

""" Union of LinkedClassifier and LinkedRegressor."""
LinkedLearner = Union{LinkedClassifier, LinkedRegressor}

"""
    $(TYPEDSIGNATURES)

Allows for approximation of continuous functions or objectives using ORTs.
Can be added to GlobalModel using `add_nonlinear_constraint` or `add_nonlinear_or_compatible`.
"""
@with_kw mutable struct BlackBoxRegressor
    constraint::Union{Nothing, JuMP.ConstraintRef, Expr}        # The "raw" constraint
    vars::Array{JuMP.VariableRef,1}                    # JuMP variables (flat)
    dependent_var::JuMP.VariableRef                    # Dependent variable
    name::String = ""                                  # Function name
    expr_vars::Array                                   # Function inputs (nonflat JuMP variables)
    varmap::Union{Nothing,Array} = get_varmap(expr_vars, vars)     # ... with the required varmapping.
    datamap::Union{Nothing,Array} = get_datamap(expr_vars, vars)     # ... with the required datamapping.
    f::Union{Nothing, Function} = functionify(constraint)         # ... and actually evaluated f'n
    g::Union{Nothing, Function} = gradientify(constraint, expr_vars)   # ... and its gradient f'n
    X::DataFrame = DataFrame(string.(vars) .=> [Float64[]
    for i=1:length(vars)]) # Function samples
    Y::Array = []                                      # Function values
    gradients::Union{Nothing, DataFrame} = nothing     # Gradients
    curvatures::Union{Nothing, Array} = nothing        # Curvature around the points
    infeas_X::DataFrame = DataFrame(string.(vars) .=> [Float64[]
    for i=1:length(vars)]) # Infeasible samples, if any
    equality::Bool = false                             # Equality check
    learners::Array{Union{IAI.OptimalTreeRegressor, IAI.OptimalTreeClassifier,
                          IAI.Heuristics.RandomForestRegressor}} = []     # Learners...
    learner_kwargs = []                                # and their kwargs...
    thresholds::Array{Pair} = []                       # For thresholding.
    ul_data::Array{Dict} = Dict[]                      # Upper/lower bounding data
    active_trees::Dict{Int64, Union{Nothing, Pair}} = Dict() # Currently active tree indices
    mi_constraints::Dict = Dict{Int64, Array{JuMP.ConstraintRef}}() # and their corresponding MI constraints,
    leaf_variables::Dict = Dict{Int64, Tuple{JuMP.VariableRef, Array, JuMP.VariableRef}}() # and leaf variables.
    active_leaves::Array = []                          # leaf of last solution
    optima::Array = []
    actuals::Array = []
    feas_gap::Array = []                               # Feasibility gaps of solutions
    relax_var::Union{Real, JuMP.VariableRef} = 0.      # slack variable
    accuracies::Array{Float64} = []                    # and the tree MSE scores.
    lls::Array{LinkedRegressor} = []                   # Linked regressor mi_constraints and leaf_variables
    bigM::Bool = false
    convex::Bool = false
    local_convexity::Float64 = 0.
    vexity::Dict = Dict{Int64, Tuple}()                # Size and convexity of leaves
    knn_tree::Union{KDTree, Nothing} = nothing         # KNN tree
    params::Dict = bbr_defaults(length(vars))          # Relevant settings
    max_Y::Union{Nothing, Real} = nothing
    min_Y::Union{Nothing, Real} = nothing
end

function Base.show(io::IO, bbr::BlackBoxRegressor)
    if bbr.equality
        println(io, "BlackBoxRegressor EQUALITY " * bbr.name * " with $(length(bbr.vars)) independent variables, ")
    else
        println(io, "BlackBoxRegressor inequality " * bbr.name * " with $(length(bbr.vars)) independent variables, ")
    end
    println(io, "and dependent variable $(bbr.dependent_var).")
    println(io, "Sampled $(length(bbr.Y)) times, and has $(length(bbr.learners)) trained ORTs.")
    if !isempty(bbr.lls)
        println(io, "Has $(length(bbr.lls)) linked constraints.")
    end
end

"""
    $(TYPEDSIGNATURES)

Allows for approximation of constraints using OCTs.
Can be added to GlobalModel using `add_nonlinear_constraint` or `add_nonlinear_or_compatible`.
"""
@with_kw mutable struct BlackBoxClassifier
    constraint::Union{Nothing, JuMP.ConstraintRef, Expr} # The "raw" constraint
    vars::Array{JuMP.VariableRef,1}                    # JuMP variables (flat)
    name::String = ""                                  # Function name
    expr_vars::Array                                   # Function inputs (nonflat JuMP variables)
    varmap::Union{Nothing,Array} = get_varmap(expr_vars, vars)     # ... with the required varmapping.
    datamap::Union{Nothing,Array} = get_datamap(expr_vars, vars)     # ... with the required datamapping.
    f::Union{Nothing, Function} = functionify(constraint)         # ... and actually evaluated f'n
    g::Union{Nothing, Function} = gradientify(constraint, expr_vars)   # ... and its gradient f'n
    X::DataFrame = DataFrame(string.(vars) .=> [Float64[] for i=1:length(vars)])
                                                       # Function samples
    Y::Array = []                                      # Function values
    gradients::Union{Nothing, DataFrame} = nothing     # Gradients
    curvatures::Union{Nothing, Array} = nothing        # Curvature around the points
    feas_ratio::Float64 = 0.                           # Feasible sample proportion
    bigM::Bool = false                                 # Whether approximations should use big-M formulations
    convex::Bool = false
    local_convexity::Float64 = 0.
    vexity::Dict = Dict{Int64, Tuple}()                # Size and convexity of leaves
    equality::Bool = false                             # Equality check
    learners::Array{Union{IAI.OptimalTreeClassifier,
                          IAI.Heuristics.RandomForestClassifier}} = []    # Learners...
    learner_kwargs = []                                # And their kwargs...
    mi_constraints::Dict = Dict{Int64, Array{JuMP.ConstraintRef}}() # and their corresponding MI constraints,
    leaf_variables::Dict = Dict{Int64, Tuple{JuMP.VariableRef, Array}}() # and their leaf variables
    active_leaves::Array = []                          # Leaf of last solution
    feas_gap::Array = []                               # Feasibility gaps of solutions
    lls::Array{LinkedClassifier} = []                  # LinkedClassifiers
    relax_var::Union{Real, JuMP.VariableRef} = 0.    # slack variable
    accuracies::Array{Float64} = []                    # and the tree misclassification scores.
    knn_tree::Union{KDTree, Nothing} = nothing         # KNN tree
    params::Dict = bbc_defaults(length(vars))          # Relevant settings
    max_Y::Union{Nothing, Real} = nothing
    min_Y::Union{Nothing, Real} = nothing
end

function Base.show(io::IO, bbc::BlackBoxClassifier)
    if bbc.equality
        println(io, "BlackBoxClassifier EQUALITY " * bbc.name * " with $(length(bbc.vars)) variables: ")
    else
        println(io, "BlackBoxClassifier inequality " * bbc.name * " with $(length(bbc.vars)) variables: ")
    end
    println(io, "Sampled $(length(bbc.Y)) times, and has $(length(bbc.learners)) trained OCTs.")
    if !isempty(bbc.lls)
        println(io, "Has $(length(bbc.lls)) linked constraints.")
    end
    if get_param(bbc, :ignore_feasibility)
        if get_param(bbc, :ignore_accuracy)
            println(io, "Ignores training accuracy and data_feasibility thresholds.")
        else
            println(io, "Ignores data feasibility thresholds.")
        end
    else
        if get_param(bbc, :ignore_accuracy)
            println("Ignores training accuracy thresholds.")
        end
    end
end

""" Union of BlackBoxClassifier and BlackBoxRegressor types. """
BlackBoxLearner = Union{BlackBoxClassifier, BlackBoxRegressor}

set_param(bbl::BlackBoxLearner, key::Symbol, val) = set_param(bbl.params, key, val)
set_param(bbls::Array{BlackBoxLearner}, key::Symbol, val) = foreach(bbl -> set_param(bbl, key, val), bbls)
get_param(bbl::BlackBoxLearner, key::Symbol) = get_param(bbl.params, key)
get_param(bbls::Array{BlackBoxLearner}, key::Symbol) = [get_param(bbl, key) for bbl in bbls]

"""
    $(TYPEDSIGNATURES)

Adds data to BlackBoxLearner.
"""
function add_data!(bbc::BlackBoxClassifier, X::DataFrame, Y::Array)
    @assert length(Y) == size(X, 1)
    if size(bbc.X,1) == 0
        bbc.feas_ratio = sum(Y .>= 0)/length(Y)
    else
        bbc.feas_ratio = (bbc.feas_ratio*size(bbc.X,1) + sum(Y .>= 0))/(size(bbc.X, 1) + length(Y))
    end
    if isnothing(bbc.gradients) && get_param(bbc, :gradients) # TODO: improve the gradient DF init.
        bbc.gradients = DataFrame(string.(bbc.vars) .=> [Union{Missing, Float64}[] for i=1:length(bbc.vars)])
        bbc.curvatures = Union{Missing, Float64}[]
    end
    if !isnothing(bbc.gradients)
        append!(bbc.gradients, DataFrame(string.(bbc.vars) .=>
            [missings(size(X, 1)) for i=1: length(bbc.vars)]), cols=:intersect)
        append!(bbc.curvatures, missings(size(X,1)))
    end
    append!(bbc.X, X, cols=:intersect)
    append!(bbc.Y, Y)
    return
end

function add_data!(bbr::BlackBoxRegressor, X::DataFrame, Y::Array)
    @assert length(Y) == size(X, 1)
    infeas_idxs = findall(x -> isinf(x), Y)
    length(infeas_idxs) >= 0.5*length(Y) &&
    throw(OCTHaGOnException("A majority of samples on BBR $(bbr.name) evaluated as infeasible. " *
        "Please check your expressions. "))
    if isnothing(bbr.gradients) && get_param(bbr, :gradients)
        bbr.gradients = DataFrame(string.(bbr.vars) .=> [Union{Missing, Float64}[] for i=1:length(bbr.vars)])
        bbr.curvatures = Union{Missing, Float64}[]
    end
    if !isempty(infeas_idxs)
        append!(bbr.infeas_X, X[infeas_idxs, :], cols=:intersect)
        clean_X = delete!(X, infeas_idxs)
        if !isnothing(bbr.gradients)
            append!(bbr.gradients, DataFrame(string.(bbr.vars) .=>
                        [missings(size(clean_X, 1)) for i=1:length(bbr.vars)]), cols=:intersect)
            append!(bbr.curvatures, missings(size(clean_X,1)))
        end
        append!(bbr.X, clean_X, cols=:intersect)
        append!(bbr.Y, deleteat!(Y, infeas_idxs))
    else
        if !isnothing(bbr.gradients)
            append!(bbr.gradients, DataFrame(string.(bbr.vars) .=>
                        [missings(size(X, 1)) for i=1:length(bbr.vars)]), cols=:intersect)
            append!(bbr.curvatures, missings(size(X,1)))
        end
        append!(bbr.X, X, cols=:intersect)
        append!(bbr.Y, Y)
    end
    return
end

"""
    $(TYPEDSIGNATURES)

Evaluates constraint violation on data in the variables, and returns distance from set.
        Note that the keys of the Dict have to be of the same type.
"""
function evaluate(bbl::BlackBoxLearner, data::Union{Dict, DataFrame})
    @assert !isnothing(bbl.constraint)
    clean_data = data_to_DataFrame(data);
    if isnothing(bbl.f)
        constr_obj = constraint_object(bbl.constraint)
        rhs_const = get_constant(constr_obj.set)
        vals = [JuMP.value(constr_obj.func, i -> clean_data[!,string(i)][j]) for j=1:size(clean_data,1)]
        if any(isinf.(vals)) || any(isnan.(vals))
            throw(OCTHaGOnException(string("Constraint ", bbl.constraint, " returned an infinite or NaN value.",
                                      "Please check your variable definitions in bbl ", bbl, " . ")))
        end
        if isnothing(rhs_const)
            vals = [-1*distance_to_set(vals[i], constr_obj.set) for i=1:length(vals)]
        else
            coeff=1
            if constr_obj.set isa MOI.LessThan
                coeff=-1
            end
            vals = coeff.*(vals .- rhs_const)
        end
        length(vals) == 1 && return vals[1]
        return vals
    else
        arrs = deconstruct(clean_data, bbl.vars, bbl.expr_vars, bbl.varmap)
        vals = []
        for arr in arrs
            try
                push!(vals, Base.invokelatest(bbl.f, (arr...)))
            catch DomainError
                push!(vals, -Inf)
            end
        end
        return vals
    end
end

"""
    evaluate_gradient(bbl::BlackBoxLearner, data::Union{Array, DataFrame})

Evaluates gradient of function through ForwardDiff.
TODO: speed-ups!.
"""
function evaluate_gradient(bbl::BlackBoxLearner, data::DataFrame)
    @assert get_param(bbl, :gradients)
    @assert !isnothing(bbl.constraint)
    @assert(size(data, 2) == length(bbl.vars))
    if bbl.constraint isa JuMP.ConstraintRef
        return DataFrame(hcat([bbl.g(Array(row))
                for row in eachrow(data[:, string.(bbl.vars)])]...)', string.(bbl.vars))
    elseif bbl.constraint isa Expr
        arrs = deconstruct(data, bbl.vars, bbl.expr_vars, bbl.varmap)
        grads = hcat([bbl.g(Float64[flat(arr)...])
                            for arr in arrs]...)'
        return DataFrame(grads[:, bbl.datamap], string.(bbl.vars))
    end
end

"""
    update_gradients(bbl::BlackBoxLearner, idxs = collect(1:size(bbl.X,1)))

Updates gradient information of selected points.
"""
function update_gradients(bbl::BlackBoxLearner, idxs::Array = collect(1:size(bbl.X,1)))
    @assert get_param(bbl, :gradients)
    isempty(idxs) && return
    empties = idxs[findall(idx -> any(ismissing.(values(bbl.gradients[idx,:]))), idxs)]
    isempty(empties) && return
    bbl.gradients[empties, :] = evaluate_gradient(bbl, bbl.X[empties, :])
    return
end

"""
    function (bbl::BlackBoxLearner)(x::Union{DataFrame,Dict,DataFrameRow})

Makes the BBC or BBR callable as a function.
"""
function (bbl::BlackBoxLearner)(X::Union{DataFrame,Dict,DataFrameRow})
    return evaluate(bbl, X)
end

""" Evaluates BlackBoxLearner at all X and stores the resulting data. """
function eval!(bbl::BlackBoxLearner, X::DataFrame)
    Y = bbl(X);
    add_data!(bbl, X, Y)
end

"""
    clear_tree_data!(bbc::BlackBoxClassifier)
    clear_tree_data!(bbr::BlackBoxRegressor)

Deletes tree data associated with object.
"""
function clear_tree_data!(bbc::BlackBoxClassifier)
    bbc.learners = [];
    bbc.learner_kwargs = []
    bbc.accuracies = []
end

function clear_tree_data!(bbr::BlackBoxRegressor)
    bbr.learners = [];
    bbr.learner_kwargs = []
    bbr.thresholds = []
    bbr.ul_data = Dict[]
end

"""
    $(TYPEDSIGNATURES)

Deletes all data (NOT constraints) associated with object.
Please clear constraints using the `clear_tree_constraints!` function.
"""
function clear_data!(bbc::BlackBoxClassifier)
    bbc.X = DataFrame(string.(bbc.vars) .=> [Float64[] for i=1:length(bbc.vars)])
    bbc.Y = [];
    bbc.gradients = nothing
    bbc.curvatures = nothing
    bbc.feas_ratio = 0
    clear_tree_data!(bbc)
end

function clear_data!(bbr::BlackBoxRegressor)
    bbr.X = DataFrame(string.(bbr.vars) .=> [Float64[] for i=1:length(bbr.vars)])
    bbr.Y = []
    bbr.gradients = nothing
    bbr.curvatures = nothing
    bbr.infeas_X = DataFrame(string.(bbr.vars) .=> [Float64[] for i=1:length(bbr.vars)]);
    clear_tree_data!(bbr)

end

""" Helper function to find which binary leaf variables are one. """
function active_leaves(bbl::Union{BlackBoxLearner, LinkedLearner})
    leaf_in = []
    for (leaf, var_tuple) in bbl.leaf_variables
        if isapprox(getvalue(var_tuple[1]), 1; atol=1e-5)
            push!(leaf_in, leaf)
        end
    end
    bbl.active_leaves = leaf_in
end

"""
    $(TYPEDSIGNATURES)

Returns all JuMP.ConstraintRefs associated with BlackBox/LinkedLearner.
Hypertype input optional for upper/lower bounding learners.
"""
function all_mi_constraints(bbl::Union{BlackBoxLearner, LinkedLearner}, hypertype::Union{String, Nothing} = nothing)
    all_constraints = []
    if hypertype == nothing
        for (leaf, constraints) in bbl.mi_constraints
            push!(all_constraints, constraints...)
        end
    elseif hypertype in valid_lowers
        for (leaf, constraints) in bbl.mi_constraints
            if leaf >= 0
                push!(all_constraints, constraints...)
            end
        end
    end
    if hypertype in valid_uppers
        for (leaf, constraints) in bbl.mi_constraints
            if leaf <= 0
                push!(all_constraints, constraints...)
            end
        end
    end
    return all_constraints
end

"""
    active_lower_tree(bbr::BlackBoxRegressor)

Returns the index of currently active lower bounding tree of BBR.
"""
function active_lower_tree(bbr::BlackBoxRegressor)
    if length(bbr.active_trees) > 2
        throw(OCTHaGOnException("Regressor $(bbr.name) has too many active trees. There is a bug in the update. "))
    elseif length(bbr.active_trees) == 1
        if collect(values(bbr.active_trees))[1].first in valid_lowers
            return collect(keys(bbr.active_trees))[1]
        else
            return # no active lower trees
        end
    elseif length(bbr.active_trees) == 2
        tree_keys = collect(keys(bbr.active_trees))
        tree_hypertypes = collect(values(bbr.active_trees))
        if tree_hypertypes[1].first == "lower"
            return tree_keys[1]
        elseif tree_hypertypes[2].first == "lower"
            return tree_keys[2]
        elseif tree_hypertypes[1].first in valid_lowers
            return tree_keys[1]
        else
            return tree_keys[2]
        end
    else
        return # no active trees
    end
end

"""
    active_upper_tree(bbr::BlackBoxRegressor)

Returns the index of currently active active_upper_tree bounding tree of BBR.
"""
function active_upper_tree(bbr::BlackBoxRegressor)
    if length(bbr.active_trees) == 1
        if collect(values(bbr.active_trees))[1].first in valid_uppers
            return collect(keys(bbr.active_trees))[1]
        else
            return # no active upper trees
        end
    elseif length(bbr.active_trees) == 2
        tree_keys = collect(keys(bbr.active_trees))
        tree_hypertypes = collect(values(bbr.active_trees))
        if tree_hypertypes[1].first == "upper"
            return tree_keys[1]
        elseif tree_hypertypes[2].first == "upper"
            return tree_keys[2]
        elseif tree_hypertypes[1].first in valid_uppers
            return tree_keys[1]
        else
            return tree_keys[2]
        end
    else
        return # no active trees
    end
end

"""
    update_local_convexity(bbl::BlackBoxLearner)

Checks proportion of "neighbor convex" sample points of BBL.
"""
function update_local_convexity(bbl::BlackBoxLearner)
    classify_curvature(bbl)
    bbl.local_convexity = sum(bbl.curvatures .== 1) / size(bbl.X, 1)
    return
end

"""
    update_vexity(bbl::BlackBoxLearner)

Checks whether a function is perhaps locally or globally convex.
Threshold sets the border of being considered for convex regression.
"""
function update_vexity(bbl::BlackBoxLearner, t::Int64 = 5)
    update_local_convexity(bbl)
    thresh = 1e-10 * (bbl.max_Y - bbl.min_Y)
    if bbl.local_convexity == 1.0
        if bbl isa BlackBoxRegressor
            # Checking against quasi_convexity with 5 random points
            cvx = true
            test_idxs = Int64.(ceil.(rand(t) .* size(bbl.X, 1)))
            diffs = [[Array(bbl.X[j, :]) - Array(bbl.X[i, :]) for i in test_idxs] for j in test_idxs]
            for i=1:t, j=1:t
                if i != j && !(bbl.Y[test_idxs[j]] >= bbl.Y[test_idxs[i]] -
                                sum(Array(bbl.gradients[test_idxs[i],:]) .* diffs[i][j]) - thresh)
                    cvx = false
                    @info "Points $(test_idxs[j]) and $(test_idxs[i]) broke the convexity of $(bbl.name)."
                    break
                end
            end
            if cvx
                bbl.convex = true
            end
        else
            # Checking against quasi_convexity with 5 random points
            cvx = true
            test_idxs = Int64.(ceil.(rand(t) .* size(bbl.X, 1)))
            diffs = [[Array(bbl.X[j, :]) - Array(bbl.X[i, :]) for i in test_idxs] for j in test_idxs]
            for i=1:t, j=1:t
                if i != j && !(bbl.Y[test_idxs[i]] >= bbl.Y[test_idxs[j]] +
                                sum(Array(bbl.gradients[test_idxs[i],:]) .* diffs[i][j]) - thresh)
                    cvx = false
                    @info "Points $(test_idxs[i]) and $(test_idxs[j]) broke the convexity of $(bbl.name)."
                    break
                end
            end
            if cvx
                bbl.convex = true
            end
        end
    end
    return
end

"""
    update_leaf_vexity(bbl::BlackBoxLearner)

Finds the local convexity of leaves (bbl.vexity) of the active lower bounding tree of BBL.
"""
function update_leaf_vexity(bbl::BlackBoxLearner)
    if bbl.convex
        bbl.vexity[1] = (size(bbl.X, 1), 1.0) # We only have a "root".
        return
    end
    lnr = bbl.learners[end] # last tree is always the active tree...
    leaf_idxs = IAI.apply(lnr, bbl.X)
    all_leaves = find_leaves(lnr)
    leaf_vexity = Dict()
    if lnr isa BlackBoxClassifier
        all_leaves = [i for i in all_leaves if Bool(IAI.get_classification_label(lnr, i))];
    end
    if any(ismissing.(bbl.curvatures))
        classify_curvature(bbl)
    end
    for leaf in all_leaves
        in_leaf_idxs = findall(x -> x == leaf, leaf_idxs)
        leaf_vexity[leaf] = (length(in_leaf_idxs), sum(bbl.curvatures[in_leaf_idxs] .> 0) / length(in_leaf_idxs))
    end
    bbl.vexity = leaf_vexity
    return
end


text:
"""Write a tea header to the given IO stream."""
Base.write(io::IO, metadata::TeaFileMetadata) = _write_tea(io, metadata)

function _write_tea(io::IO, metadata::TeaFileMetadata)::Int
    bytes_written = write(io, MAGIC_VALUE)
    bytes_written += write(io, metadata.item_start)
    bytes_written += write(io, metadata.item_end)
    bytes_written += write(io, Int64(length(metadata.sections)))  # section count
    for section in metadata.sections
        bytes_written += _write_section(io, section)
    end

    # We now pad this out with zeros to become a multiple of 8 bytes, so that the item area
    # starts with nice alignment for subsequent memory mapping.
    n_bytes_padding = metadata.item_start - bytes_written
    bytes_written += write(io, repeat([0x00], n_bytes_padding))

    return bytes_written
end

"""Write a section, including its metadata."""
function _write_section(io::IO, section::AbstractSection)::Int
    return (
        write(io, section_id(section)) # section_id
        + write(io, _tea_size(section))  # next_section_offset
        + _write_tea(io, section)
    )
end

"""
Get the size of this object in bytes when written.
"""
_tea_size(x::Real)::Int32 = sizeof(x)
_tea_size(x::AbstractString)::Int32 = sizeof(Int32) + sizeof(x)
_tea_size(x::Base.UUID)::Int32 = sizeof(x)
_tea_size(x::Vector)::Int32 = sizeof(Int32) + (isempty(x) ? 0 : sum(_tea_size, x))
_tea_size(x::Field)::Int32 = _tea_size(x.type_id) + _tea_size(x.offset) + _tea_size(x.name)
_tea_size(x::NameValue)::Int32 = _tea_size(x.name) + _tea_size(kind(x)) + _tea_size(x.value)

function _tea_size(x::AbstractSection)::Int32
    # For all sections, we just write every field in that section.
    return sum(_tea_size, (getproperty(x, name) for name in propertynames(x)))
end

"""Write x in the way that a TeaFile expects."""
_write_tea(io::IO, x::Real) = write(io, x)
_write_tea(io::IO, x::Base.UUID) = write(io, x.value)

function _write_tea(io::IO, string::AbstractString)::Int
    return (
        write(io, Int32(sizeof(string)))
        + write(io, string)
    )
end

function _write_tea(io::IO, xs::Vector)
    # Write the number of items as an int32, followed by the items.
    bytes_written = write(io, Int32(length(xs)))
    for x in xs
        bytes_written += _write_tea(io, x)
    end
    return bytes_written
end

function _write_tea(io::IO, field::Field)::Int
    return (
        write(io, field.type_id)
        + write(io, field.offset)
        + _write_tea(io, field.name)
    )
end

function _write_tea(io::IO, name_value::NameValue)::Int
    return (
        _write_tea(io, name_value.name)
        + write(io, kind(name_value))
        + _write_tea(io, name_value.value)
    )
end

function _write_tea(io::IO, section::ItemSection)::Int
    return (
        write(io, section.item_size)
        + _write_tea(io, section.item_name)
        + _write_tea(io, section.fields)
    )
end

function _write_tea(io::IO, section::TimeSection)::Int
    return (
        write(io, section.epoch)
        + write(io, section.ticks_per_day)
        + _write_tea(io, section.time_field_offsets)
    )
end

_write_tea(io::IO, section::ContentDescriptionSection) = _write_tea(io, section.description)
_write_tea(io::IO, section::NameValueSection) = _write_tea(io, section.name_values)

ext: jl


text:
export sslice
using PyCall
"""
	function sslice(a,b)
Usefull in Gridspec calls
"""
function sslice(a,b)
    pycall(pybuiltin("slice"), PyObject, a,b+1)
end

function significance(pval; ns="ns")
	if pval > 0.1
        sigsymbol = ns
        elseif pval > 0.05
        sigsymbol = "#"
        elseif pval > 0.01
        sigsymbol = "*"
        elseif pval > 0.001
        sigsymbol = "**"
        elseif pval < 0.001
        sigsymbol = "***"
    else
        sigsymbol = "sig unavailable"
    end
end


text:
using ComputedArrays
using Test

ca1 = ComputedArray(cospi, range(0, stop=1, length=1000));
@test size(ca1) == (1000,)
@test eltype(ca1) == Float64

pot(x,i) = ceil(Int32, i*x^i);
ca2 = ComputedArray(pot, range(0, stop=1, length=11), 1:5);
@test size(ca2) == (11,5)
@test eltype(ca2) == Int32

v(x1,y1,x2,y2) = 1/hypot(x1-x2, y1-y2);
x1s = 0 .+ range(-1, stop=1, length=10) |> collect
y1s = 0 .+ range(-1, stop=1, length=11) |> collect
x2s = 3 .+ range(-1, stop=1, length=12) |> collect
y2s = 4 .+ range(-1, stop=1, length=13) |> collect
ca4 = ComputedArray(v, x1s, y1s, x2s, y2s)
@test size(ca4) == (10,11,12,13)
@test eltype(ca4) == Float64

ca4p = ComputedArray(v, x1s, y1s, x2s, y2s, order=(2,3,4,1))
@test size(ca4p) == (11,12,13,10)
@test ca4[3,5,6,8] == ca4p[5,6,8,3]
pca4 = PermutedDimsArray(ca4, (2,3,4,1))
@test ca4p == pca4


text:
# Question 1
import Pkg; Pkg.add("Distributions"); Pkg.add("QuantEcon")
using Distributions, QuantEcon

function profit_max_q(a, c, mu, sigma, method, n)
        b = rand(LogNormal(mu, sigma), n)
        if method == "mc"
            q = Float64((a-c)/(2*mean(b)))
            return q
        elseif method == "quad"
            nodes_weights = qnwlogn(n, mu, sigma) #http://quantecon.github.io/QuantEcon.jl/latest/api/QuantEcon.html
            nodes = nodes_weights[1]
            weights = nodes_weights[2]
            mean_b = sum(nodes .* weights)
            q = Float64((a-c)/(2*mean_b))
        end
        return q
    end

profit_max_q(3,1,1,1,"mc",10)

# we dont have to worry about code_warntype
@code_warntype profit_max_q(3,1,1,1,"mc",10)

# Question 2
using Distributions
function estimate_pi(lower,upper,num_draws)
    n_land_in_circle = 0
    x = rand(Uniform(lower, upper), num_draws)
    y = rand(Uniform(lower, upper), num_draws)
        for i in 1:num_draws
        r_sq = x[i]^2 + y[i]^2
            if r_sq < 1
                n_land_in_circle += 1
            end
        end
    return n_land_in_circle / num_draws * 4
 end

estimate_pi(0,1,1000000)


text:
abstract type MatrixConeTriangle <: MOI.AbstractVectorSet end

struct EmptyCone <: MatrixConeTriangle end
side_dimension(::EmptyCone) = 0

struct PositiveSemidefinite2x2ConeTriangle <: MatrixConeTriangle end
side_dimension(::PositiveSemidefinite2x2ConeTriangle) = 2


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

using Markdown
using InteractiveUtils

# ╔═╡ b32bec97-359c-4a83-ba0b-2503e5c50122
using Pkg, DrWatson

# ╔═╡ 3a39a0db-54de-45ed-a8d5-e948e7c14f89
begin
	@quickactivate "NumericalMethodsforEngineers"
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ a68242b5-0995-428f-9b73-37cf943e9b87
md" ### f-01-02.jl"

# ╔═╡ 895bad2c-0a85-4584-9cb5-d1579385bfe3
begin
	A = [10. 1. -5.; -20. 3. 20.; 5. 3. 5.]
	b = [1., 2., 6.]
end

# ╔═╡ e15153d7-bdee-459c-b236-6c037e701e29
md" ##### Julia way: select method based on characterictics of A."

# ╔═╡ dc76295f-88ad-4919-b5e7-d428d27288fa
c = A\b

# ╔═╡ e310cdd3-2d8b-4e17-b4ef-ae8a74740eec
A * c

# ╔═╡ 144f98fa-58ce-4518-b7dc-b98459b50832
b

# ╔═╡ 532a43a1-eb7b-4fc6-861d-421a48eeee5b
md" ##### NMFE implementation:"

# ╔═╡ 00dedc2b-42c1-44ed-a998-c7f49a06c11b
# Convert to upper triang;e
function to_upper(A::Matrix, b::Vector)
	At = deepcopy(A)
	bt = deepcopy(b)
	r, c = size(At)

	for k in 1:r-1
	  if abs(At[k, k]) > 1.0e-6
		for i in k+1:r
		  x = At[i, k] / At[k, k]
		  At[i, k] = 0.0
		  for j in k+1:r
			At[i, j] -= At[k, j] * x
		  end
		  bt[i] -= bt[k] * x
		end
	  else
		md"Zero pivot found in row $k"
	  end
	end
	(A=At, b=bt)
end

# ╔═╡ 420d0142-790d-4565-b22f-accb0c049023
function direct_solve(A::Matrix, b::Vector)
	At = deepcopy(A)
	bt = deepcopy(b)
	r, c = size(At)

	for i in r:-1:1
	  x = bt[i]
	  if i < r
		for j in i+1:r
		  x -= At[i, j] * bt[j]
		end
	  end
	  bt[i] = x / At[i, i]
	end
	bt
end

# ╔═╡ 546655ea-e19b-450e-a298-3a420c8685d8
lu(A)

# ╔═╡ b9d85721-0442-43fe-8936-b3c30246e51b
md" ##### Modified Matrix:"

# ╔═╡ 843e2d08-c226-4d1d-9998-b359a9dffa10
At, bt = to_upper(A, b)

# ╔═╡ 617129ac-1321-476d-93cb-c1696ddbd40c
md" ##### NMFE solution vector:"

# ╔═╡ 119074dc-3372-4cc3-84cd-ab2abdd1a7b2
c_nmfe = direct_solve(At, bt)

# ╔═╡ f6fec610-96aa-40e2-8646-ff8fe3c2e804
A * c_nmfe

# ╔═╡ 585d3b15-949d-4aca-828d-2477a9c88599
md" ### End of f-01-02.jl"

# ╔═╡ Cell order:
# ╟─a68242b5-0995-428f-9b73-37cf943e9b87
# ╠═b32bec97-359c-4a83-ba0b-2503e5c50122
# ╠═3a39a0db-54de-45ed-a8d5-e948e7c14f89
# ╠═895bad2c-0a85-4584-9cb5-d1579385bfe3
# ╟─e15153d7-bdee-459c-b236-6c037e701e29
# ╠═dc76295f-88ad-4919-b5e7-d428d27288fa
# ╠═e310cdd3-2d8b-4e17-b4ef-ae8a74740eec
# ╠═144f98fa-58ce-4518-b7dc-b98459b50832
# ╟─532a43a1-eb7b-4fc6-861d-421a48eeee5b
# ╠═00dedc2b-42c1-44ed-a998-c7f49a06c11b
# ╠═420d0142-790d-4565-b22f-accb0c049023
# ╠═546655ea-e19b-450e-a298-3a420c8685d8
# ╟─b9d85721-0442-43fe-8936-b3c30246e51b
# ╠═843e2d08-c226-4d1d-9998-b359a9dffa10
# ╟─617129ac-1321-476d-93cb-c1696ddbd40c
# ╠═119074dc-3372-4cc3-84cd-ab2abdd1a7b2
# ╠═f6fec610-96aa-40e2-8646-ff8fe3c2e804
# ╟─585d3b15-949d-4aca-828d-2477a9c88599


text:
# test_df_scan.jl

@testset "Test standard DataFrame scan" begin
    con = DBInterface.connect(DuckDB.DB)
    df = DataFrame(a = [1, 2, 3], b = [42, 84, 42])

    DuckDB.register_data_frame(con, df, "my_df")
    GC.gc()

    results = DBInterface.execute(con, "SELECT * FROM my_df")
    GC.gc()
    df = DataFrame(results)
    @test names(df) == ["a", "b"]
    @test size(df, 1) == 3
    @test df.a == [1, 2, 3]
    @test df.b == [42, 84, 42]

    DBInterface.close!(con)
end

@testset "Test DataFrame scan with NULL values" begin
    con = DBInterface.connect(DuckDB.DB)
    df = DataFrame(a = [1, missing, 3], b = [missing, 84, missing])

    DuckDB.register_data_frame(con, df, "my_df")

    results = DBInterface.execute(con, "SELECT * FROM my_df")
    df = DataFrame(results)
    @test names(df) == ["a", "b"]
    @test size(df, 1) == 3
    @test isequal(df.a, [1, missing, 3])
    @test isequal(df.b, [missing, 84, missing])

    DBInterface.close!(con)
end

@testset "Test DataFrame scan with numerics" begin
    con = DBInterface.connect(DuckDB.DB)
    numeric_types = [Int8, Int16, Int32, Int64, UInt8, UInt16, UInt32, UInt64, Float32, Float64]
    for type in numeric_types
        my_df = DataFrame(a = [1, missing, 3], b = [missing, 84, missing])
        my_df[!, :a] = convert.(Union{type, Missing}, my_df[!, :a])
        my_df[!, :b] = convert.(Union{type, Missing}, my_df[!, :b])

        DuckDB.register_data_frame(con, my_df, "my_df")

        results = DBInterface.execute(con, "SELECT * FROM my_df")
        df = DataFrame(results)
        @test isequal(df, my_df)
    end

    DBInterface.close!(con)
end

@testset "Test DataFrame scan with various types" begin
    con = DBInterface.connect(DuckDB.DB)

    # boolean
    my_df = DataFrame(a = [true, false, missing])

    DuckDB.register_data_frame(con, my_df, "my_df")

    results = DBInterface.execute(con, "SELECT * FROM my_df")
    df = DataFrame(results)
    @test isequal(df, my_df)

    # date/time/timestamp
    my_df = DataFrame(
        date = [Date(1992, 9, 20), missing, Date(1950, 2, 3)],
        time = [Time(23, 3, 1), Time(11, 49, 33), missing],
        timestamp = [DateTime(1992, 9, 20, 23, 3, 1), DateTime(1950, 2, 3, 11, 49, 3), missing]
    )

    DuckDB.register_data_frame(con, my_df, "my_df")

    results = DBInterface.execute(con, "SELECT * FROM my_df")
    df = DataFrame(results)
    @test isequal(df, my_df)

    DBInterface.close!(con)
end

@testset "Test DataFrame scan with strings" begin
    con = DBInterface.connect(DuckDB.DB)

    # date/time/timestamp
    my_df = DataFrame(str = ["hello", "this is a very long string", missing, "obligatory mühleisen"])

    DuckDB.register_data_frame(con, my_df, "my_df")

    results = DBInterface.execute(con, "SELECT * FROM my_df")
    df = DataFrame(results)
    @test isequal(df, my_df)

    DBInterface.close!(con)
end

@testset "Test DataFrame scan projection pushdown" begin
    con = DBInterface.connect(DuckDB.DB)
    df = DataFrame(a = [1, 2, 3], b = [42, 84, 42], c = [3, 7, 18])

    DuckDB.register_data_frame(con, df, "my_df")
    GC.gc()

    results = DBInterface.execute(con, "SELECT b FROM my_df")
    GC.gc()
    df = DataFrame(results)
    @test names(df) == ["b"]
    @test size(df, 1) == 3
    @test df.b == [42, 84, 42]

    results = DBInterface.execute(con, "SELECT c, b FROM my_df")
    GC.gc()
    df = DataFrame(results)
    @test names(df) == ["c", "b"]
    @test size(df, 1) == 3
    @test df.b == [42, 84, 42]
    @test df.c == [3, 7, 18]

    results = DBInterface.execute(con, "SELECT c, a, a FROM my_df")
    GC.gc()
    df = DataFrame(results)
    @test names(df) == ["c", "a", "a_1"]
    @test size(df, 1) == 3
    @test df.c == [3, 7, 18]
    @test df.a == [1, 2, 3]
    @test df.a_1 == [1, 2, 3]

    GC.gc()

    DBInterface.close!(con)
end

max_stars_repo_name: Dtenwolde/duckdb

max_issues_repo_name: Dtenwolde/duckdb

max_forks_repo_name: Dtenwolde/duckdb


text:
#=
Implement locking in a binary tree. A binary tree node can be locked or unlocked only if all of its descendants or ancestors are not locked.

Design a binary tree node class with the following methods:

    is_locked, which returns whether the node is locked
    lock, which attempts to lock the node. If it cannot be locked, then it should return false. Otherwise, it should lock it and return true.
    unlock, which unlocks the node. If it cannot be unlocked, then it should return false. Otherwise, it should unlock it and return true.

You may augment the node to add parent pointers or any other property you would like. You may assume the class is used in a single-threaded program, so there is no need for actual locks or mutexes. Each method should run in O(h), where h is the height of the tree.
=#

using Test
include("Solutions/problem24_locking_binary_tree.jl")

#= binary tree
        0(f)
      /      \
     1(f)    2(t)
    / \      /  \
 3(t) 4(t) 5(t) 6(t)
=#
node3 = Node(3, nothing, nothing, true)
node4 = Node(4, nothing, nothing, true)
node5 = Node(5, nothing, nothing, true)
node2 = Node(2, node3, nothing, true)
node1 = Node(1, node4, node5, false)
node0 = Node(0, node1, node2, false)

@test is_locked(node3) == true
@test is_locked(node1) == false
@test lock(node0) == false
@test lock(node1) == true
@test node1.locked == true

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#=======================================================

Tests for the Node type

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

using Random, SkipLists, Test
using SkipLists: Node, LeftSentinel, RightSentinel

@testset "Node tests" begin
    Random.seed!(0)

    @testset "Construct Node" begin
        node = Node{:List}([1, 2, 3]; capacity=10)
        @test isa(node, Node{Int64,:List})
        @test height(node) > 0
        @test SkipLists.key(node) == 1
        @test SkipLists.capacity(node) == 10
        @test 1 ∈ node && 2 ∈ node && 3 ∈ node

        # If the node doesn't have any items in it, isempty(node) should be
        # true
        node = Node{:List}([]; capacity=10)
        @test isempty(node)

        # Construct left and right sentinel nodes
        lsentinel = SkipLists.LeftSentinel{Int64,:List}()
        rsentinel = SkipLists.RightSentinel{Int64,:List}()
        @test SkipLists.is_left_sentinel(lsentinel)
        @test SkipLists.is_right_sentinel(rsentinel)
    end

    @testset "Insert and delete from Node" begin
        node = Node{:List}(1:2:50; capacity=100)

        success = true
        for ii = shuffle(2:2:50)
            insert!(node, ii)
            success = success && ii ∈ node
        end
        @test success

        success = true
        for ii = shuffle(1:2:50)
            delete!(node, ii)
            success = success && ii ∉ node
        end
        @test success

        # If we attempt to insert into a node past its capacity, we should
        # get an error
        node = Node{:List}(1:10; capacity=10)
        @test SkipLists.isfull(node)
        @test_throws ErrorException insert!(node, 11)
    end

    @testset "Split Node" begin
        node = Node{:List}(1:10)

        node, right_node = SkipLists.split!(node; capacity=25)
        @test node.vals == 1:5
        @test right_node.vals == 6:10
        @test SkipLists.capacity(right_node) == 25
    end

    @testset "Compare nodes" begin
        lsentinel = LeftSentinel{Int64,:List}()
        rsentinel = RightSentinel{Int64,:List}()
        node_1 = Node{Int64}([0])
        node_2 = Node{Int64}(1:10)

        @test lsentinel < node_1     && lsentinel ≤ node_1
        @test node_1 < node_2        && node_1 ≤ node_2
        @test node_2 < rsentinel     && node_2 ≤ rsentinel
        @test lsentinel < rsentinel  && lsentinel ≤ rsentinel
        @test lsentinel ≤ lsentinel && !(lsentinel < lsentinel)
        @test rsentinel ≤ rsentinel && !(rsentinel < rsentinel)

        @test 0 == node_1 == Node{Int64}(0:5:10) == 0
        @test 1 == node_2 == Node{Int64}([1]) == 1
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using Documenter, FEMBase

DEVELOPER_GUIDE = [
          "mesh.md",
          "fields.md",
          "basis.md",
          "integration.md",
          "elements.md",
          "problems.md",
          "solvers.md",
          "postprocessing.md",
          "results.md",
          "materials.md",
          "testing.md"]

LIBRARY = ["api.md"]

PAGES = [
         "Home" => "index.md",
         "Developing JuliaFEM" => DEVELOPER_GUIDE,
         "Library" => LIBRARY
        ]

makedocs(modules=[FEMBase],
         format = :html,
         checkdocs = :all,
         sitename = "FEMBase.jl",
         authors = "Jukka Aho",
         analytics = "UA-83590644-1",
         pages = PAGES)


text:
# Show methods for the intermediate representation.
# References:
#  - base/show.jl
#  - ?show
#
#
# TODO: use  MIME"application/vnd.julia-vscode.diagnostics" for runtime diagnostics
#
import PrettyPrinting: quoteof

"Indention increment."
const indent_width = 4

"Return string of current indent size in `io`."
indent(io::IO) = repeat(' ', get(io, :indent, 0))

"Increment the `:indent` value by `indent_width`."
inc_indent(io::IO) = IOContext(io, :indent => get(io, :indent, 0) + indent_width)

#-------------------
function Base.show(io::IO, fill::Fill)
    pprint(io, fill)
end

quoteof(f::Fill) = :(Fill($(quoteof(f.color)),
        $(quoteof(f.image)),
        $(quoteof(f.gradient_color)),
        $(quoteof(f.gradient_rotation))))

#-------------------
function Base.show(io::IO, font::Font)
    pprint(io, font)
end

quoteof(f::Font) = :(Font($(quoteof(f.family)),
            $(quoteof(f.style)),
            $(quoteof(f.weight)),
            $(quoteof(f.size)),
            $(quoteof(f.align)),
            $(quoteof(f.rotation)),
            $(quoteof(f.decoration)) ))

#-------------------
function Base.show(io::IO, line::Line)
    pprint(io, line)
end

quoteof(l::Line) = :(Line($(quoteof(l.color)), $(quoteof(l.style)),
            $(quoteof(l.shape)), $(quoteof(l.width))))

#-------------------
function Base.show(io::IO, c::Coordinate)
    compact = get(io, :compact, false)
    print(io, "(", c.x, ",", c.y, ")")
end

function Base.show(io::IO, ::MIME"text/plain", c::Coordinate)
    print(io, c)
end

#-------------------
function shownames(io::IO, g::Graphics)
    print(io, "Graphics(",
        "dimension=", g.dimension,
        ", fill=",      g.fill,
        ", font=",      g.font,
        ", line=",      g.line,
        ", offset=",    g.offset,
        ", position=",  g.position, ")")
end

function Base.show(io::IO, g::Graphics)
    print(io, "Graphics(",
            g.dimension, ", ",
            g.fill, ", ",
            g.font, ", ",
            g.line, ", ",
            g.offset, ", ",
            g.position, ")")
end

#-------------------
function Base.show(io::IO, labelvector::Vector{PnmlLabel})
    show(io, MIME"text/plain"(), labelvector)
end
function Base.show(io::IO, mime::MIME"text/plain", labelvector::Vector{PnmlLabel})
    print(io, indent(io), typeof(labelvector), "[")
    io = inc_indent(io)
    for (i,label) in enumerate(labelvector)
        i > 1 && print(io, indent(io))
        pprint(io, label)
        i < length(labelvector) && print(io, "\n")
    end
    print(io, "]")

end

function Base.show(io::IO, label::PnmlLabel)
    pprint(io, label)
end

function Base.show(io::IO, mime::MIME"text/plain", label::PnmlLabel)
    pprint(io, label)
end

quoteof(l::PnmlLabel) = :(PnmlLabel($(quoteof(l.tag)), $(quoteof(l.dict))))

#-------------------
function Base.show(io::IO, elvector::Vector{AnyElement})
    show(io, MIME"text/plain"(), elvector)
end
function Base.show(io::IO, mime::MIME"text/plain", elvector::Vector{AnyElement})
    print(io, typeof(elvector), "[")
    io = inc_indent(io)
    for (i,el) in enumerate(elvector)
        print(io, "\n", indent(io), "$i: ")
        pprint(io, el)
    end
    print(io, "]")
end

function Base.show(io::IO, el::AnyElement) #TODO Make parametric.
    pprint(io, el)
end

function Base.show(io::IO, mime::MIME"text/plain", el::AnyElement) #TODO Make parametric.
    pprint(io, el)
end

quoteof(a::AnyElement) = :(AnyElement($(quoteof(a.tag)), $(quoteof(a.dict))))

#-------------------
Base.summary(io::IO, ti::ToolInfo) = print(io, summary(ti))
function Base.summary(ti::ToolInfo)
    string(typeof(ti), " name ", ti.toolname, ", version ", ti.version,
           ", ", length(ti.infos), " infos")
end

function Base.show(io::IO, toolvector::Vector{ToolInfo})
    for (i, ti) in enumerate(toolvector)
        print(io, "\n", indent(io), "$i: ")
        show(io, ti)
    end
end

function Base.show(io::IO, ti::ToolInfo)
    pprint(io, ti)
end

quoteof(ti::ToolInfo) = :($(quoteof(ti.toolname)), $(quoteof(ti.version)), $(quoteof(ti.infos)))

#-------------------
function Base.show(io::IO, tg::TokenGraphics)
    print(io, "positions: ", tg.positions)
end

function Base.show(io::IO, ::MIME"text/plain", tg::TokenGraphics)
    print(io, tg)
end

#-------------------
function Base.show(io::IO, name::Name)
    print(io, typeof(name), " '", name.text, "'")
    !isnothing(name.graphics) && print(io, ", has graphics")
    !isnothing(name.tools)    && print(io, ", ", length(name.tools), " tool info")
end
function Base.show(io::IO, ::MIME"text/plain", name::Name)
    print(io, name)
end

#-------------------
Base.summary(io::IO, oc::ObjectCommon) = print(io, summary(oc))
function Base.summary(oc::ObjectCommon)
    string(isnothing(oc.graphics) ? ", no graphics, " : ", has graphics, ",
           isnothing(oc.tools)  ? 0 : length(oc.tools),  " tools, ",
           isnothing(oc.labels) ? 0 : length(oc.labels), " labels")
end

function Base.show(io::IO, oc::ObjectCommon)
    io = inc_indent(io)
    #if !isnothing(oc.graphics) ||
    if !isnothing(oc.tools) || !isnothing(oc.labels)
        print(io, ", ")
        #!isnothing(oc.graphics) && pprint(io, oc.graphics)
        if !isnothing(oc.tools)
            println(io, "\n", indent(io), "tools:")
            show(inc_indent(io), oc.tools)
        end
        if !isnothing(oc.labels)
            println(io, "\n", indent(io), "labels:")
            show(inc_indent(io), oc.labels)
        end
    end
end

function show_common(io::IO, x::Union{PnmlNet, PnmlObject, AbstractLabel})
    isempty(x.com) && return
    #    print(io, ", ")
    show(io, MIME"text/plain"(), x.com )
end

#---------------------------------------------------------------------------------
Base.summary(io::IO, place::Place)  = summary(io, summary(place))
function Base.summary(place::Place)
    string(typeof(place))
end

function Base.show(io::IO, place::Place)
    print(io, summary(place),
          " id ", place.id,
          ", name '" , has_name(place) ? name(place) : "", "'",
          ", type ", place.sorttype,
          ", marking ", place.marking)
#          ", name ", place.name)
    show_common(io, place)
end

function Base.show(io::IO, ::MIME"text/plain", place::Place)
    show(io, place)
end
function Base.show(io::IO, ::MIME"text/plain", placevector::Vector{Place})
    show(io, placevector)
end

function Base.show(io::IO, placevector::Vector{Place})
    isempty(placevector) && return
    for (i,place) in enumerate(placevector)
        print(io, indent(io), place)
        i < length(placevector) && print(io, "\n")
    end
end

#-------------------
function Base.show(io::IO, trans::Transition)
    print(io, typeof(trans),
          " id ", trans.id,
          ", name '", has_name(trans) ? name(trans) : "", "'",
          ", condition ", trans.condition)
    show_common(io, trans)
end

function Base.show(io::IO, ::MIME"text/plain", transvector::Vector{Transition})
    show(io, transvector)
end

function Base.show(io::IO, transvector::Vector{Transition})
    for (i,trans) in enumerate(transvector)
        print(io, indent(io), trans)
        i < length(transvector) && print(io, "\n")
    end
end

function Base.show(io::IO, ::MIME"text/plain", trans::Transition)
    show(io, trans)
end

#-------------------
#TODO Make RefPlace, RefTransition an Abstract Type
function Base.show(io::IO, ::MIME"text/plain", r::ReferenceNode)
    show(io, r)
end
function Base.show(io::IO, r::ReferenceNode)
    print(io, typeof(r), " (id ", pid(r), ", ref ", refid(r))
    show_common(io, r)
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", rvector::Vector{<:ReferenceNode})
    show(io, rvector)
end
function Base.show(io::IO, rvector::Vector{<:ReferenceNode})
    for (i,r) in enumerate(rvector)
        print(io, indent(io), r)
        i < length(rvector) && print(io, "\n")
    end
end

#-------------------
function Base.show(io::IO, arc::Arc)
    print(io, typeof(arc), " id ", arc.id,
          ", name '", has_name(arc) ? name(arc) : "", "'",
          ", source: ", arc.source,
          ", target: ", arc.target,
          ", inscription: ", arc.inscription)
    show_common(io, arc)
end
function Base.show(io::IO, ::MIME"text/plain", arcvector::Vector{Arc})
    for (i,arc) in enumerate(arcvector)
        print(io, indent(io), arc)
        i < length(arcvector) && print(io, "\n")
    end
end
function Base.show(io::IO, ::MIME"text/plain", arc::Arc)
    show(io, arc)
end

#-------------------
Base.summary(io::IO, page::Page) = print(io, summary(page))
function Base.summary( page::Page)
    string(typeof(page)," id ", page.id, ", ",
           " name '", name(page), "', ",
           length(page.places), " places, ",
           length(page.transitions), " transitions, ",
           length(page.arcs), " arcs, ",
           isnothing(declarations(page)) ? 0 : length(declarations(page)), " declarations, ",
           length(page.refPlaces), " refP, ",
           length(page.refTransitions), " refT, ",
           length(page.subpages), " subpages",
           summary(page.com)
           )
end

function show_page_field(io::IO, label::AbstractString, x)
    println(io, indent(io), label)
    if !isnothing(x) && length(x) > 0
        show(inc_indent(io), MIME"text/plain"(), x)
        print(io, "\n")
    end
end

function Base.show(io::IO, ::MIME"text/plain", p::Page)
    show(io, p)
end
function Base.show(io::IO, page::Page)
    #TODO Add support for :trim and :compact
    println(io, indent(io), summary(page))
    # Start indent here. Will indent subpages.
    inc_io = inc_indent(io)

    show_page_field(inc_io, "places:",         places(page))
    show_page_field(inc_io, "transitions:",    transitions(page))
    show_page_field(inc_io, "arcs:",           arcs(page))
    show_page_field(inc_io, "declaration:",    declarations(page))
    show_page_field(inc_io, "refPlaces:",      refplaces(page))
    show_page_field(inc_io, "refTransitions:", reftransitions(page))
    show_common(io, page)
    show_page_field(inc_io, "subpages:",       pages(page))
end

function Base.show(io::IO, ::MIME"text/plain", pages::Vector{Page})
    show(io, pages)
end
function Base.show(io::IO, pages::Vector{Page})
    isempty(pages) && return
    for (i,page) in enumerate(pages)
        show(io, MIME"text/plain"(), page)
        i < length(pages) && print(io, "\n")
    end
end

#-------------------
Base.summary(io::IO, net::PnmlNet) = print(io, summary(net))
function Base.summary(net::PnmlNet)
    string( typeof(net), " id ", pid(net),
            " name '", has_name(net) ? name(net) : "", "', ",
            " type ", nettype(net), ", ",
            length(pages(net)), " pages ",
            length(declarations(net)), " declarations",
            summary(net.com))
end

# No indent here.
function Base.show(io::IO, net::PnmlNet)
    println(io, summary(net))
    iio = inc_indent(io) # Indent any declarations.
    foreach(declarations(net)) do decl
        print(iio, indent(io))
        show(iio, MIME"text/plain"(), decl)
        println(iio, "\n")
    end
    show_common(io, net)
    show(io, pages(net))
end

function Base.show(io::IO, ::MIME"text/plain", net::PnmlNet)
    show(io, net)
end

#-------------------
Base.summary(io::IO, pnml::PnmlModel) = print(io, summary(pnml))
function Base.summary(pnml::PnmlModel)
    string(typeof(pnml), " model with ",  length(nets(pnml)), " nets" )
end

# No indent done here.
function Base.show(io::IO, pnml::PnmlModel)
    println(io, summary(pnml))
    println(io, "namespace = ", namespace(pnml))
    for (i, net) in enumerate(nets(pnml))
        show(io, MIME"text/plain"(), net)
        if i < length(nets(pnml))
            print(io, "\n")
        end
    end
    # Omit display of any xml
end

function Base.show(io::IO, ::MIME"text/plain", pnml::PnmlModel)
    return show(io, pnml)
end

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

#-------------------
Base.summary(io::IO, ptm::PTMarking)  = summary(io, summary(ptm))
function Base.summary(ptm::PTMarking)
    string(typeof(ptm))
end

function Base.show(io::IO, ptm::PTMarking)
    pprint(io, ptm)
end

quoteof(m::PTMarking) = :(PTMarking($(quoteof(m.value)), $(quoteof(m.com))))

#-------------------
Base.summary(io::IO, hlm::HLMarking) = summary(hlm)
function
    Base.summary(hlm::HLMarking)
    string(typeof(hlm))
end

function Base.show(io::IO, hlm::HLMarking)
    pprint(io, hlm)
end

quoteof(m::HLMarking) = :(HLMarking($(quoteof(m.text)), $(quoteof(m.term)), $(quoteof(m.com))))

#-------------------
function Base.show(io::IO, cond::Condition)
    pprint(io, cond)
end

function Base.show(io::IO, ::MIME"text/plain", cond::Condition)
    show(io, cond)
end
quoteof(c::Condition) = :(Condition($(quoteof(c.text)), $(quoteof(c.term)), $(quoteof(c.com))))

#-------------------
function Base.show(io::IO, inscription::PTInscription)
    pprint(io, inscription)
end
function Base.show(io::IO, ::MIME"text/plain", inscription::PTInscription)
    show(io, inscription)
end
quoteof(i::PTInscription) = :(PTInscription($(quoteof(i.value)), $(quoteof(i.com))))

#-------------------
function Base.show(io::IO, inscription::HLInscription)
    pprint(io, inscription)
end
function Base.show(io::IO, ::MIME"text/plain", inscription::HLInscription)
    show(io, inscription)
end

quoteof(i::HLInscription) =
    :(HLInscription($(quoteof(i.text)), $(quoteof(i.term)), $(quoteof(i.com))))

#-------------------
function Base.show(io::IO, declarations::Vector{AbstractDeclaration})
    iio = inc_indent(io)
    print(io, indent(io), typeof(declarations), "[")

    for (i,dec) in enumerate(declarations)
        print(iio, "\n", indent(iio))
        show(inc_indent(io), MIME"text/plain"(), dec)
    end
    print(io, "]")
end
function Base.show(io::IO, declare::AbstractDeclaration)
    pprint(io, declare)
end
quoteof(i::AbstractDeclaration) = :(AbstractDeclaration($(quoteof(i.id)), $(quoteof(i.name)),
         $(quoteof(i.com))))
#-------------------
function Base.show(io::IO, terms::Vector{AbstractTerm})
    iio = inc_indent(io)
    print(io, indent(io), typeof(terms), "[")

    for (i,term) in enumerate(terms)
        print(iio, "\n", indent(iio))
        show(inc_indent(io), term)
    end
    print(io, "]")
end
function Base.show(io::IO, term::AbstractTerm)
    pprint(io, term)
end
quoteof(t::AbstractTerm) = :(AbstractTerm($(quoteof(t.tag)), $(quoteof(t.dict))))
#         $(quoteof(t.com))))

function Base.show(io::IO, term::Term)
    pprint(io, term)
end
quoteof(t::Term) = :(Term($(quoteof(t.tag)), $(quoteof(t.dict))))

#-------------------
function Base.show(io::IO, nsorts::Vector{NamedSort})

    print(io, typeof(nsorts), "[")
    for (i,dec) in enumerate(nsorts)
        print(io, "\n", indent(io))
        show(inc_indent(io), dec)
        i < length(nsorts) && print(io, "\n")
    end
    print(io, "]")
end

function Base.show(io::IO, nsort::NamedSort)
    pprint(IOContext(io, :displaysize => (24, 180)), nsort)
end

quoteof(n::NamedSort) = :(NamedSort($(quoteof(n.id)), $(quoteof(n.name)), $(quoteof(n.def))))


text:
function binarysearch(arr,x,l,r)
    if r>=l
        mid = (l + (r-l)÷2)
        if (arr[mid]==x)
            return mid
        elseif (arr[mid]>x)
            return binarysearch(arr,x,l,mid-1)
        else
            return binarysearch(arr,x,mid+1,r)
        end
    end
    return -1
end

v = [2 3 4 5 9 11 15 78 98 101 169 259]
len = length(v)
x = 170
t = binarysearch(v,x,1,len)
if(t==-1)
    println("Element = ",x," not found")
else
    println("Found element ",x," at = " ,t)
end

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
include("common.jl")
using Plots
import PGFPlots
pgfplots()
PGFPlots.pushPGFPlotsPreamble(raw"\usepackage{amsmath}")

# Conventions used for approx. eigenvalues and wavefunctions
const EPS = raw"\widetilde{\varepsilon}"
const PSI = raw"\widetilde{u}"

const defm = (3, :xcross, Plots.Stroke(1.0, :black, 0.8, :solid))

function gap_versus_ecut(data, data2=nothing, data3=nothing; iband=1, ik=1, log=false)
    Ecuts = data.Ecuts
    results = data.results
    n_bands = length(results[1][:λ][1])

    eigenvalues = [results[iEcut][:λ][ik][iband] for iEcut = 1:length(Ecuts)]
    eigenvalues_next = [results[iEcut][:λ][ik][iband+1]
                        for iEcut = 1:length(Ecuts)]
    estimated = [results[iEcut][:gap][ik][iband] for iEcut = 1:length(Ecuts)]
    computed = eigenvalues_next - eigenvalues

    kwargs = ()
    if log
        estimated = computed[end] .- estimated
        kwargs = (yaxis=:log, )
    end
    estlabel = (raw"$\mu^\ast_{" * "$(iband+1)" * raw"} - " *
                EPS * raw"_{" * "$iband" * raw",X}$")
    estlabel1 = isnothing(data2) ? estlabel : estlabel * " (\$M = $n_bands\$)"

    p = plot(legend=:bottomright; xlabel=raw"$E_\text{cut}$", kwargs...)
    if !log
        plot!(p, Ecuts, computed,
              label=(raw"$" * EPS * raw"_{" * "$(iband+1)" * raw",X} - " *
                     EPS * raw"_{" * "$iband" * raw",X}$"),
              m=defm, linestyle=:dot, color=2)
    end
    if !isnothing(data3)
        d3_Ecuts = data3.Ecuts
        @assert length(d3_Ecuts) ≤ length(Ecuts)
        n_bands3 = length(data3.results[1][:λ][1])
        estimated3 = [data3.results[iEcut][:gap][ik][iband]
                      for iEcut = 1:length(d3_Ecuts)]
        if log
            estimated3 = computed[end] .- estimated3
        end
        plot!(p, d3_Ecuts, estimated3, label=estlabel * raw" ($M = \text{dim}(X)$)", m=defm,
              linestyle=:dashdot, color=4)
    end
    if !isnothing(data2)
        d2_Ecuts = data2.Ecuts
        @assert length(d2_Ecuts) ≤ length(Ecuts)
        n_bands2 = length(data2.results[1][:λ][1])
        estimated2 = [data2.results[iEcut][:gap][ik][iband]
                      for iEcut = 1:length(d2_Ecuts)]
        if log
            estimated2 = computed[end] .- estimated2
        end
        plot!(p, d2_Ecuts, estimated2, label=estlabel * " (\$M = $n_bands2\$)", m=defm,
              linestyle=:dash, color=3)
    end
    plot!(p, Ecuts, estimated, label=estlabel1,
          m=defm, linestyle=:solid, color=1)

    p
end

function residual_total_versus_ecut(data; selected_bands=[1], ik=1, arithmetic=true)
    Ecuts = data.Ecuts
    results = data.results
    n_bands = length(results[1][:λ][ik])

    # Prepare arrays
    residual_estimated = [[results[iEcut][:residual_norm][ik][iband]
                           for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]

    # Note: Algorithm error is contained inside residual_computed already
    error_algorithm = [[results[iEcut][:error_algorithm][ik][iband]
                        for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]
    error_arithmetic = [[results[iEcut][:error_arithmetic][ik][iband]
                        for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]
    error_algoarith = [error_arithmetic[iband] + error_algorithm[iband]
                       for iband in 1:n_bands]

    p = plot(legend=:topright, xlabel=raw"$E_\text{cut}$", yaxis=:log)
    for iband in selected_bands
        kwargs = (colour=iband, )
        siband = " $iband"
        if length(selected_bands) == 1
            kwargs = NamedTuple()
            siband = ""
        end
        arithmetic || (kwargs = merge(kwargs, (m=defm, )))
        plot!(p, Ecuts, residual_estimated[iband];
              label=(raw"$\|V " * PSI * raw"\|$" * siband), kwargs...)
        if arithmetic
            plot!(p, Ecuts, error_algoarith[iband]; ls=:dash,
                  label="arithmetic + algorithm error" * siband, kwargs...)
        end
    end
    p
end

function residual_versus_ecut(data; selected_bands=[1], ik=1)
    Ecuts = data.Ecuts
    results = data.results
    n_bands = length(results[1][:λ][ik])

    # Prepare arrays
    residual_estimated = [[results[iEcut][:residual_norm][ik][iband]
                           for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]
    residual_computed = [[results[iEcut][:residual_norm_term_computed][ik][iband]
                          for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]
    residual_loc = [[results[iEcut][:residual_norm_term_AtomicLocal][ik][iband]
                     for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]
    residual_nloc = [[results[iEcut][:residual_norm_term_AtomicNonlocal][ik][iband]
                      for iEcut = 1:length(Ecuts)] for iband = 1:n_bands]

    # Derived arrays
    residual_beyond = [residual_loc[iband] .+ residual_nloc[iband] for iband = 1:n_bands]
    p = plot(legend=:bottomleft, xlabel=raw"$E_\text{cut}$", yaxis=:log)
    for iband in selected_bands
        getcolor(iband, i) = length(selected_bands) == 1 ? i : iband
        siband = " $iband"
        if length(selected_bands) == 1
            siband = ""
        end

        labels = [raw"$\|P_{Y} V " * PSI * raw"\|$" * siband,
                  raw"estimated $\|P_{Y^\perp} V " * PSI * raw"\|$" * siband,
                  raw"estimated $\|V " * PSI * raw"\|$" * siband]
        iband != selected_bands[1] && (labels[2:end] .= "")

        plot!(p, Ecuts, residual_estimated[iband]; label=labels[3],
              linestyle=:solid, color=getcolor(iband, 1))
        plot!(p, Ecuts, residual_computed[iband]; label=labels[1],
              linestyle=:dash, color=getcolor(iband, 2))
        plot!(p, Ecuts, residual_beyond[iband]; label=labels[2],
              linestyle=:dot, color=getcolor(iband, 3))
    end
    p
end

function residual_versus_ecut2(data::AbstractDict; iband=1, ik=1, iEcut=1)
    sEcuts = sort(collect(keys(data)), by=sEcut -> data[sEcut].results[iEcut][:Ecut2])
    @assert all(length(data[sEcut].results) == 1 for sEcut in sEcuts)

    # Prepare data arrays
    Ecut2s = [data[sEcut].results[iEcut][:Ecut2] for sEcut in sEcuts]
    residual_computed = [data[sEcut].results[iEcut][:residual_norm_term_computed][ik][iband]
                         for sEcut in sEcuts]
    residual_loc = [data[sEcut].results[iEcut][:residual_norm_term_AtomicLocal][ik][iband]
                    for sEcut in sEcuts]
    residual_nloc = [data[sEcut].results[iEcut][:residual_norm_term_AtomicNonlocal][ik][iband]
                     for sEcut in sEcuts]
    residual_estimated = [data[sEcut].results[iEcut][:residual_norm][ik][iband]
                          for sEcut in sEcuts]

    p = plot(legend=:bottomleft, xlabel=raw"$E_\text{cut}^{(2)}$", yaxis=:log)
    plot!(p, Ecut2s, residual_estimated,
          label=raw"estimated $\|V " * PSI * raw"\|$")
    plot!(p, Ecut2s, residual_computed;
          label=raw"$\| P_Y V " * PSI *  raw"\|$", ls=:dash)
    plot!(p, Ecut2s, residual_loc;
          label=raw"estimated $\|P_{Y^\perp} V_\text{loc} " * PSI *  raw"\|$", ls=:dot)
    plot!(p, Ecut2s, residual_nloc;
          label=raw"estimated $\|P_{Y^\perp} V_\text{nonloc} " * PSI * raw"\|$", ls=:dashdot)

    p
end

function energy_errors_versus_ecut(data; selected_bands=[1, 2], ik=1)
    results = data.results
    n_bands = length(results[1][:λ][ik])

    # Prepare arrays
    eigenvalues = [[results[iEcut][:λ][ik][iband] for iEcut = 1:length(data.Ecuts)]
                    for iband = 1:n_bands]
    bound_bauer_fike = [[results[iEcut][:bound_bauer_fike][ik][iband]
                         for iEcut = 1:length(data.Ecuts)] for iband = 1:n_bands]
    bound_kato_temple = [[results[iEcut][:bound_kato_temple][ik][iband]
                          for iEcut = 1:length(data.Ecuts)] for iband = 1:n_bands]

    Ecuts = data.Ecuts[1:end-1]
    p = plot(legend=:bottomleft, yaxis=:log, xlabel=raw"$E_\text{cut}$")
    for iband in selected_bands
        kwargs = (colour = iband, )
        siband = " $iband"
        if length(selected_bands) == 1
            kwargs = ()
            siband = ""
        end
        labels = ["``true'' error$siband", "Bauer-Fike$siband", "Kato-Temple$siband"]
        if iband != selected_bands[1]
            labels[2:end] .= ""
        end

        plot!(p, Ecuts, bound_bauer_fike[iband][1:end-1], label=labels[2],
              linestyle=:dash; kwargs...)

        kt_plot = copy(bound_kato_temple[iband][1:end-1])
        kt_plot[isinf.(kt_plot)] .= NaN
        plot!(p, Ecuts, kt_plot, label=labels[3], linestyle=:dot; kwargs...)

        error = abs.(eigenvalues[iband][1:end-1] .- eigenvalues[iband][end])
        error[iszero.(error)] .= NaN
        plot!(p, Ecuts, error, label=labels[1]; kwargs...)
    end
    p
end

function bands_errors(data; iEcut=1)
    results = data.results[iEcut]

    # Reconstruct a dummy basis
    if data.system == "Cohen-Bergstresser"
        model = model_cohen_bergstresser()
    elseif data.system == "Linear Silicon"
        model = model_linear_silicon()
    else
        error("Model $(data.system) not known.")
    end
    ksymops = [[(Mat3{Int}(I), Vec3(zeros(3)))] for _ in 1:results[:n_kpoints]]
    basis = PlaneWaveBasis(model, results[:Ecut], results[:kcoords], ksymops)

    # Compute the error to display on the bands
    results[:bound_best] = deepcopy(results[:bound_bauer_fike])
    results[:bound_best_label] = [fill("BF", length(results[:bound_bauer_fike][ik]))
                                  for ik in 1:results[:n_kpoints]]
    for ik in 1:results[:n_kpoints]
        for iband in 1:length(results[:bound_bauer_fike][ik])
            kt_lower = (results[:bound_kato_temple][ik][iband]
                        ≤ results[:bound_bauer_fike][ik][iband])
            if kt_lower
                results[:bound_best][ik][iband] = results[:bound_kato_temple][ik][iband]
                results[:bound_best_label][ik][iband] = "KT"
            end
        end
    end

    λerror = [(   results[:bound_best][ik]        # discretisation error
                + results[:error_algorithm][ik]   # algorithmic error
                + results[:error_arithmetic][ik]  # floating-point error
              ) for ik in 1:results[:n_kpoints]]

    shift = results[:εF]
    shift = results[:λ][1][4]  # Shift wrt. 4th eigenvalue at Gamma point
    plot_band_data_mod((basis=basis, λ=results[:λ], λerror=λerror),
                       results[:bound_best_label],
                       shift=shift, klabels=results[:klabels])
end

function plot_band_data_mod(band_data, best_label; shift=nothing,
                            klabels=Dict{String, Vector{Float64}}(), unit=:Ha)
    @assert !isnothing(shift)
    data = DFTK.prepare_band_data(band_data, klabels=klabels)

    Δ = 0.025  # length of error bars (in units of kdistance)

    kstart = 0
    # For each branch, plot all bands, spins and errors
    p = plot()
    for ibranch = 1:data.n_branches
        kdistances = data.kdistances[ibranch]
        for iband = 1:data.n_bands
            energies = (data.λ[ibranch][:up][iband, :] .- shift)
            yerror = data.λerror[ibranch][:up][iband, :]

            plot!(p, kdistances, energies, color=:black, label="")  # Bands

            redpoints_x = Float64[]
            redpoints_y = Float64[]
            bluepoints_x = Float64[]
            bluepoints_y = Float64[]
            for i in 1:length(kdistances)
                xs = [kdistances[i] - Δ, kdistances[i] + Δ, NaN,
                      kdistances[i] - Δ, kdistances[i] + Δ, NaN,
                      kdistances[i], kdistances[i], NaN]
                ys = [energies[i] - yerror[i], energies[i] - yerror[i], NaN,
                      energies[i] + yerror[i], energies[i] + yerror[i], NaN,
                      energies[i] - yerror[i], energies[i] + yerror[i], NaN]

                if best_label[kstart + i][iband] == "BF"
                    append!(bluepoints_x, xs)
                    append!(bluepoints_y, ys)
                else
                    append!(redpoints_x, xs)
                    append!(redpoints_y, ys)
                end
            end
            common = (linewidth=1, label="", linealpha=0.85)
            plot!(p, redpoints_x, redpoints_y; color=:red, common...)
            plot!(p, bluepoints_x, bluepoints_y; color=:blue, common...)
        end
        kstart += length(kdistances)
    end

    # X-range: 0 to last kdistance value
    xlims!(p, (0, data.kdistances[end][end]))
    xticks!(p, data.ticks["distance"], data.ticks["label"])
    ylims!(p, -0.2, 0.2)

    p
end

function setup()
    fac = 0.75
    default(size=tuple(Int.(ceil.(fac .* [600, 400]))...))
end

function main()
    setup()
    mticks(r) = 10. .^r, ["10^{$i}" for i in r]

    # Cohen-Bergstresser
    #
    begin
        cb_gamma = load_results("cb_gamma_Ecut_8bands_Float64.hdf5")
        p = gap_versus_ecut(cb_gamma)
        savefig(xlims!(p, -Inf, 15), "cb_gap_ecut.pdf")

        p = residual_total_versus_ecut(cb_gamma, selected_bands=[1], arithmetic=false)
        savefig(ylims!(xlims!(p, -Inf, 30), 1e-8, Inf), "cb_residual_versus_ecut.pdf")

        p = residual_total_versus_ecut(cb_gamma, selected_bands=[1], arithmetic=true)
        savefig(yticks!(p, mticks(0:-2:-13)...), "cb_residual_arithmetic_versus_ecut.pdf")
    end
    begin
        cb_gamma_Df64 = load_results("cb_gamma_Ecut_8bands_Double64_fast_gap.hdf5")
        savefig(energy_errors_versus_ecut(cb_gamma_Df64, selected_bands=[1]),
                "cb_energy_errors_versus_ecut.pdf")
    end
    begin
        cb_bands = load_results("cb_bands.hdf5")
        savefig(bands_errors(cb_bands), "cb_band_errors.pdf")
    end

    # Silicon
    #
    begin
        p = residual_versus_ecut2(load_results("si_gamma_residual_versus_Ecut2.hdf5",
                                               has_subdirs=true))
        savefig(p, "si_residual_versus_ecut2.pdf")
    end
    begin
        si_gamma25 = load_results("si_gamma_Ecut_25bands_Float64.hdf5")
        si_gamma50 = load_results("si_gamma_Ecut_50bands_Float64.hdf5")
        si_gamma_full = load_results("si_gamma_Ecut_full_bands_Float64.hdf5")
        p = gap_versus_ecut(si_gamma25, si_gamma50, si_gamma_full)
        savefig(xlims!(p, -Inf, 65), "si_gap_ecut.pdf")
        p = residual_versus_ecut(si_gamma25, selected_bands=[1])
        savefig(yticks!(p, mticks(0:-2:-14)...), "si_residual_versus_ecut.pdf")
        p = energy_errors_versus_ecut(si_gamma25, selected_bands=[1])
        savefig(yticks!(p, mticks(0:-2:-14)...), "si_energy_errors_versus_ecut.pdf")
    end
    begin
        si_bands = load_results("si_bands.hdf5")
        savefig(bands_errors(si_bands), "si_band_errors.pdf")
    end
end

(abspath(PROGRAM_FILE) == @__FILE__) && main()


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.2
#   kernelspec:
#     display_name: Julia 1.5.0-rc2
#     language: julia
#     name: julia-1.5
# ---

"""
    SOneTo(n)

Return a statically-sized `AbstractUnitRange` starting at `1`, functioning as the `axes` of a `StaticArray`
"""
struct SOneTo{n} <: AbstractUnitRange{Int} end

Base.first(::SOneTo) = 1
Base.last(::SOneTo{n}) where {n} = n::Int

# +
SOneTo(n::Int) = SOneTo{n}()

function SOneTo{n}(r::AbstractUnitRange) where {n}
    ((first(r) == 1) & (last(r) == n)) && return SOneTo{n}()

    errmsg(r) = throw(DimensionMismatch("$r is inconsistent with SOneTo{$n}"))
    errmsg(r)
end
# -

Base.Tuple(::SOneTo{n}) where {n} = ntuple(identity, Val(n))

Base.axes(s::SOneTo) = (s,)
Base.size(s::SOneTo) = (length(s),)
Base.length(s::SOneTo{n}) where {n} = n

Base.@propagate_inbounds function Base.getindex(s::SOneTo, i::Int)
    @boundscheck checkbounds(s, i)
    return i
end
Base.@propagate_inbounds function Base.getindex(s::SOneTo, s2::SOneTo)
    @boundscheck checkbounds(s, s2)
    return s2
end

Base.@pure Base.iterate(::SOneTo{n}) where {n} = ifelse(n::Int<1,nothing,(1,1))
Base.iterate(::SOneTo{n}, s::Int) where {n} = ifelse(s<n::Int, (s+1,s+1),nothing)

function Base.getproperty(::SOneTo{n}, s::Symbol) where {n}
    if s === :start
        return 1
    elseif s === :stop
        return n::Int
    else
        error("type SOneTo has no property $s")
    end
end

Base.show(io::IO, ::SOneTo{n}) where {n} = print(io, "SOneTo(", n::Int, ")")

# +
Base.@pure function Base.checkindex(::Type{Bool}, ::SOneTo{n1}, ::SOneTo{n2}) where {n1, n2}
    return n1::Int >= n2::Int
end

function Base.promote_rule(a::Type{Base.OneTo{T}}, ::Type{SOneTo{n}}) where {T,n}
    Base.OneTo{promote_type(T, Int)}
end


text:
immutable MvSkewTDist <: ContinuousMultivariateDistribution
    ξ::Vector{Float64}
    Ω::PDMat
    α::Vector{Float64}
    df::Float64
end

function MvSkewTDist(ξ::Vector{Float64}, Ω::Matrix{Float64}, α::Vector{Float64}, df::Float64)
    MvSkewTDist(ξ, PDMat(Ω), α, df)
end

length(dist::MvSkewTDist) = length(dist.α)

function _rand!{T<:Real}(dist::MvSkewTDist, x::AbstractVector{T})
    chisqd = Chisq(dist.df)
    w = rand(chisqd)/dist.df
    sndist = MvSkewNormal(dist.Ω.mat, dist.α)
    rand!(sndist, x)
    # x = ξ + z/sqrt(x)
    broadcast!(/, x, x, sqrt(w))
    broadcast!(+, x, x, dist.ξ)
    x
end

function _rand!(dist::MvSkewTDist, X::DenseMatrix)
    chisqd = Chisq(dist.df)
    w = rand(chisqd, 1, size(X,2))/dist.df
    sndist = MvSkewNormal(dist.Ω.mat, dist.α)
    rand!(sndist, X)
    broadcast!(/, X, X, sqrt(w))
    broadcast!(+, X, X, dist.ξ)
    X
end

function _logpdf!{T<:Real}(r::AbstractArray, dist::MvSkewTDist, X::DenseMatrix{T})
    k = length(dist)
    Y = broadcast(-, X, dist.ξ)
    Q = invquad(dist.Ω, Y)
    ω = Diagonal(sqrt(diag(dist.Ω)))
    logtd = - 0.5*logdet(dist.Ω) + _log_g(Q, dist.df, k)
    r[:] = log(2) + logtd + _logT₁(dot(dist.α,ω\Y) * sqrt((dist.df + k)./(Q + dist.df)), dist.df + k)
end

function _logpdf{T<:Real}(dist::MvSkewTDist, x::AbstractVector{T})
    k = length(dist)
    Q = invquad(dist.Ω, x - dist.ξ)
    ω = Diagonal(sqrt(diag(dist.Ω)))
    logtd = - 0.5*logdet(dist.Ω) + _log_g(Q, dist.df, k)
    log(2) + logtd + _logT₁(dot(dist.α,ω\(x-dist.ξ)) * sqrt((dist.df + k)/(Q + dist.df)), dist.df + k)
end

function μ(dist::MvSkewTDist)
    δ(dist.Ω, dist.α) * sqrt(dist.df/π) * exp(lgamma(0.5*(dist.df - 1)) - lgamma(0.5*dist.df))
end

function mean(dist::MvSkewTDist)
    ω = Diagonal(sqrt(diag(dist.Ω)))
    dist.ξ + ω * μ(dist)
end

function var(dist::MvSkewTDist)
    ω2 = Diagonal(diag(dist.Ω))
    ω = Diagonal(sqrt(ω2))
    diag(ω2*(dist.df/(dist.df-2))) - (ω*μ(dist)).^2
end

function cov(dist::MvSkewTDist)
    ω = Diagonal(sqrt(diag(dist.Ω)))
    mu = μ(dist)
    (dist.df/(dist.df-2)) * dist.Ω.mat - ω*mu*mu'ω
end


text:
#=
finding all the included files for a module

1. Revise-like approach
  * mostly adapted from https://github.com/timholy/Revise.jl/tree/b0c5c864ea78b93caaa820cb9cfc45eca47f43ff
  * only works for precompiled modules
2. CSTPraser-based approach
  * static parsing -- works for all the files but costly
  * TODO: excludes files in submodules when searched from the parent module
  * TODO: looks for non-toplevel `include` calls
=#

# Revise-like approach
# --------------------

using Base: PkgId, UUID

"""
    parentfile, included_files = modulefiles(mod::Module)

Return the `parentfile` in which `mod` was defined, as well as a list of any
other files that were `include`d to define `mod`. If this operation is unsuccessful,
`(nothing, nothing)` is returned.
All files are returned as absolute paths.
"""
function modulefiles(mod::Module)
  # NOTE: src_file_key stuff was removed when adapted
  parentfile = String(first(methods(getfield(mod, :eval))).file)
  id = Base.PkgId(mod)
  if id.name == "Base" || Symbol(id.name) ∈ stdlib_names
    parentfile = normpath(Base.find_source_file(parentfile))
    filedata = Base._included_files
  else
    use_compiled_modules() || return nothing, nothing   # FIXME: support non-precompiled packages
    _, filedata = pkg_fileinfo(id)
  end
  filedata === nothing && return nothing, nothing
  included_files = filter(mf -> mf[1] == mod, filedata)
  return fixpath(parentfile), [fixpath(mf[2]) for mf in included_files]
end

# Fix paths to files that define Julia (base and stdlibs)
function fixpath(
  filename::AbstractString;
  badpath = basebuilddir,
  goodpath = basepath("..")
)
  startswith(filename, badpath) || return fullpath(normpath(filename)) # NOTE: `fullpath` added when adapted
  filec = filename
  relfilename = relpath(filename, badpath)
  relfilename0 = relfilename
  for strippath in (joinpath("usr", "share", "julia"),)
    if startswith(relfilename, strippath)
      relfilename = relpath(relfilename, strippath)
      if occursin("stdlib", relfilename0) && !occursin("stdlib", relfilename)
        relfilename = joinpath("stdlib", relfilename)
      end
    end
  end
  ffilename = normpath(joinpath(goodpath, relfilename))
  if (isfile(filename) & !isfile(ffilename))
    ffilename = normpath(filename)
  end
  fullpath(ffilename) # NOTE: `fullpath` added when adapted
end

"""
    basebuilddir

Julia's top-level directory when Julia was built, as recorded by the entries in
`Base._included_files`.
"""
const basebuilddir = let # NOTE: changed from `begin` to `let` when adapted
  sysimg = filter(x -> endswith(x[2], "sysimg.jl"), Base._included_files)[1][2]
  dirname(dirname(sysimg))
end

use_compiled_modules() = Base.JLOptions().use_compiled_modules != 0

# For tracking Julia's own stdlibs
const stdlib_names = Set([
  :Base64,
  :CRC32c,
  :Dates,
  :DelimitedFiles,
  :Distributed,
  :FileWatching,
  :Future,
  :InteractiveUtils,
  :Libdl,
  :LibGit2,
  :LinearAlgebra,
  :Logging,
  :Markdown,
  :Mmap,
  :OldPkg,
  :Pkg,
  :Printf,
  :Profile,
  :Random,
  :REPL,
  :Serialization,
  :SHA,
  :SharedArrays,
  :Sockets,
  :SparseArrays,
  :Statistics,
  :SuiteSparse,
  :Test,
  :Unicode,
  :UUIDs,
])

function pkg_fileinfo(id::PkgId)
  uuid, name = id.uuid, id.name
  # Try to find the matching cache file
  paths = Base.find_all_in_cache_path(id)
  sourcepath = Base.locate_package(id)
  if length(paths) > 1
    fpaths = filter(path->Base.stale_cachefile(sourcepath, path) !== true, paths)
    if isempty(fpaths)
      # Work-around for #371 (broken dependency prevents tracking):
      # find the most recent cache file. Presumably this is the one built
      # to load the package.
      sort!(paths; by=path->mtime(path), rev=true)
      deleteat!(paths, 2:length(paths))
    else
      paths = fpaths
    end
  end
  isempty(paths) && return nothing, nothing
  @assert length(paths) == 1
  path = first(paths)
  provides, includes_requires = try
    parse_cache_header(path)
  catch
    return nothing, nothing
  end
  mods_files_mtimes, _ = includes_requires
  for (pkgid, buildid) in provides
    if pkgid.uuid === uuid && pkgid.name == name
      return path, mods_files_mtimes
    end
  end
end

# A near-copy of the same method in `base/loading.jl`. However, this retains the full module path to the file.
function parse_cache_header(f::IO)
  modules = Vector{Pair{PkgId,UInt64}}()
  while true
    n = read(f, Int32)
    n == 0 && break
    sym = String(read(f, n)) # module name
    uuid = UUID((read(f, UInt64), read(f, UInt64))) # pkg UUID
    build_id = read(f, UInt64) # build UUID (mostly just a timestamp)
    push!(modules, PkgId(uuid, sym) => build_id)
  end
  totbytes = read(f, Int64) # total bytes for file dependencies
  # read the list of requirements
  # and split the list into include and requires statements
  includes = Tuple{Module,String,Float64}[]
  requires = Pair{Module,PkgId}[]
  while true
    n2 = read(f, Int32)
    n2 == 0 && break
    depname = String(read(f, n2))
    mtime = read(f, Float64)
    n1 = read(f, Int32)
    mod = (n1 == 0) ? Main : Base.root_module(modules[n1].first)
    if n1 != 0
      # determine the complete module path
      while true
        n1 = read(f, Int32)
        totbytes -= 4
        n1 == 0 && break
        submodname = String(read(f, n1))
        mod = getfield(mod, Symbol(submodname))
        totbytes -= n1
      end
    end
    if depname[1] != '\0'
      push!(includes, (mod, depname, mtime))
    end
    totbytes -= 4 + 4 + n2 + 8
  end
  @assert totbytes == 12 "header of cache file appears to be corrupt"
  return modules, (includes, requires)
end

function parse_cache_header(cachefile::String)
  io = open(cachefile, "r")
  try
    !Base.isvalid_cache_header(io) && throw(ArgumentError("Invalid header in cache file $cachefile."))
    return parse_cache_header(io)
  finally
    close(io)
  end
end

# CSTParser-based approach
# ------------------------

"""
    included_files = modulefiles(mod::String, entrypath::String)::Vector{String}

Returns all the files in `mod` module that can be reached via [`include`](@ref)
  calls from `entrypath`.
Note this function currently only looks for static toplevel calls (i.e. miss the
  calls in non-toplevel scope).
"""
function modulefiles(mod::String, entrypath::String, files = Vector{String}(); inmod = false)
  isfile′(entrypath) || return files

  push!(files, entrypath)

  text = read(entrypath, String)
  items = toplevelitems(text; mod = mod, inmod = inmod)

  for item in items
    if item isa ToplevelCall
      expr = item.expr
      if isinclude(expr)
        nextfile = expr.args[3].val
        nextentrypath = joinpath(dirname(entrypath), nextfile)
        isfile(nextentrypath) || continue
        # `nextentrypath` is always in `mod`
        modulefiles(mod, nextentrypath, files; inmod = true)
      end
    end
  end

  return files
end

#=
find entry file of a module
=#

"""
    entrypath, line = moduledefinition(mod::Module)

Returns an entry file of `mod`, and its definition line.

!!! note
    This function works for non-precompiled packages.
"""
function moduledefinition(mod::Module) # NOTE: added when adapted
  evalmethod = first(methods(getfield(mod, :eval)))
  parentfile = String(evalmethod.file)
  line = evalmethod.line
  id = Base.PkgId(mod)
  if id.name == "Base" || id.name == "Core" || Symbol(id.name) ∈ stdlib_names  # NOTE: "Core" is added when adapted
    parentfile = normpath(Base.find_source_file(parentfile))
  end
  fixpath(parentfile), line
end

#=
auto-module detection for the current file

# TODO:
# use the same logics as finding module files for this auto-module detection
# ref: https://github.com/JunoLab/Juno.jl/issues/411
=#

using CodeTools, LNR

LNR.cursor(data::AbstractDict) = cursor(data["row"], data["column"])

function modulenames(data, pos)
  main = haskey(data, "module") ? data["module"] :
         haskey(data, "path") ? CodeTools.filemodule(data["path"]) :
         "Main"
  main == "" && (main = "Main")
  sub = CodeTools.codemodule(data["code"], pos)
  main, sub
end

# keeps the latest file that has been used for Main module scope
const MAIN_MODULE_LOCATION = Ref{Tuple{String, Int}}(moduledefinition(Main))

handle("module") do data
  main, sub = modulenames(data, cursor(data))

  mod = CodeTools.getmodule(main)
  smod = CodeTools.getmodule(mod, sub)

  if main == "Main" && sub == ""
    MAIN_MODULE_LOCATION[] = get!(data, "path", ""), data["row"]
  end

  return d(:main => main,
           :sub  => sub,
           :inactive => (mod==nothing),
           :subInactive => smod==nothing)
end

#=
find all modules
=#

handle("allmodules") do
  sort!([string(m) for m in CodeTools.allchildren(Main)])
end


text:
# Basic zero-mean Gaussian (uni/multi-variate) Metropolis step.
struct StaticRWM{D <: Union{Normal,MvNormal}} <: Metropolis
  proposal::D
end
function _update(rng::Random.AbstractRNG, rwm::StaticRWM, curr::T, logprob::Function) where T
  cand = curr + rand(rng, rwm.proposal)
  log_acceptance_prob = logprob(cand) - logprob(curr)
  accept = log_acceptance_prob > -Random.randexp(rng)
  draw = accept ? cand : curr
  return (draw, accept)
end


text:
"""
The HCubature module is a pure-Julia implementation of
multidimensional "h-adaptive" integration.  That is, given
an n-dimensional integral

``\\int_{a_1}^{b_1}\\int_{a_2}^{b_2}\\cdots\\int_{a_n}^{b_n} (\\vec{x}) d^n\\vec{x}``

then `hcubature(f, a, b)` computes the integral, adaptively subdividing
the integration volume into smaller and smaller pieces until convergence
is achieved to the desired tolerance (specified by optional `rtol` and
`atol` keyword arguments, described in more detail below.

Because `hcubature` is written purely in Julia, the integrand `f(x)`
can return any vector-like object (technically, any type supporting
`+`, `-`, `*` real, and `norm`: a Banach space).  You can integrate
real, complex, and matrix-valued integrands, for example.
"""
module HCubature

using StaticArrays, LinearAlgebra
import Combinatorics, DataStructures, QuadGK

export hcubature, hquadrature

include("genz-malik.jl")
include("gauss-kronrod.jl")

struct Box{n,T<:Real,TI,TE<:Real}
    a::SVector{n,T}
    b::SVector{n,T}
    I::TI
    E::TE
    kdiv::Int
end
Base.isless(i::Box, j::Box) = isless(i.E, j.E)

cubrule(v::Val{n}, ::Type{T}) where {n,T} = GenzMalik(v, T)
cubrule(::Val{1}, ::Type{T}) where {T} = GaussKronrod(T)

# trivial rule for 0-dimensional integrals
struct Trivial; end
function (::Trivial)(f, a::SVector{0}, b::SVector{0}, norm)
    I = f(a)
    return I, norm(I - I), 1
end
cubrule(::Val{0}, ::Type{T}) where {T} = Trivial()
countevals(::Trivial) = 1

function hcubature_(f, a::SVector{n,T}, b::SVector{n,T}, norm, rtol_, atol, maxevals, initdiv) where {n, T<:AbstractFloat}
    rtol = rtol_ == 0 == atol ? sqrt(eps(T)) : rtol_
    (rtol < 0 || atol < 0) && throw(ArgumentError("invalid negative tolerance"))
    maxevals < 0 && throw(ArgumentError("invalid negative maxevals"))
    initdiv < 1 && throw(ArgumentError("initdiv must be positive"))

    rule = cubrule(Val{n}(), T)
    numevals = evals_per_box = countevals(rule)

    Δ = (b-a) / initdiv
    b1 = initdiv == 1 ? b : a+Δ
    I, E, kdiv = rule(f, a,b1, norm)
    (n == 0 || iszero(prod(Δ))) && return I,E
    firstbox = Box(a,b1, I,E,kdiv)
    boxes = DataStructures.BinaryMaxHeap{typeof(firstbox)}()
    push!(boxes, firstbox)

    ma = MVector(a)
    mb = MVector(b)

    if initdiv > 1 # initial box divided by initdiv along each dimension
        skip = true # skip the first box, which we already added
        @inbounds for c in CartesianIndices(ntuple(i->Base.OneTo(initdiv), Val{n}())) # Val ntuple loops are unrolled
            if skip; skip=false; continue; end
            for i = 1:n
                ma[i] = a[i]+(c[i]-1)*Δ[i]
                mb[i] = c[i]==initdiv ? b[i] : a[i]+c[i]*Δ[i]
            end
            x = SVector(ma)
            y = SVector(mb)
            # this is shorter and has unrolled loops, but somehow creates a type instability:
            # x = SVector(ntuple(i -> a[i]+(c[i]-1)*Δ[i], Val{n}()))
            # y = SVector(ntuple(i -> c[i]==initdiv ? b[i] : a[i]+c[i]*Δ[i], Val{n}()))
            box = Box(x,y, rule(f, x,y, norm)...)
            I += box.I; E += box.E; numevals += evals_per_box
            push!(boxes, box)
        end
    end

    (E ≤ max(rtol*norm(I), atol) || numevals ≥ maxevals) && return I,E

    @inbounds while true
        # get box with largest error
        box = pop!(boxes)
        # split box along dimension kdiv
        w = (box.b[box.kdiv] - box.a[box.kdiv]) * T(0.5)
        ma[:] = box.a
        ma[box.kdiv] += w
        a′ = SVector(ma)
        mb[:] = box.b
        mb[box.kdiv] -= w
        b′ = SVector(mb)
        # evaluate and push the two new boxes
        box1 = Box(a′,box.b, rule(f, a′,box.b, norm)...)
        box2 = Box(box.a,b′, rule(f, box.a,b′, norm)...)
        push!(boxes, box1)
        push!(boxes, box2)
        I += box1.I + box2.I - box.I
        E += box1.E + box2.E - box.E
        numevals += 2*evals_per_box
        # convergence test:
        Inorm = norm(I)
        (E ≤ max(rtol*Inorm, atol) || numevals ≥ maxevals || !isfinite(Inorm)) && break
    end

    # roundoff paranoia: re-sum
    I = zero(I)
    E = zero(E)
    for i = 1:length(boxes.valtree)
        I += boxes.valtree[i].I
        E += boxes.valtree[i].E
    end
    return I,E
end

function hcubature_(f, a::SVector{n,T}, b::SVector{n,S},
                    norm, rtol, atol, maxevals, initdiv) where {n, T<:Real, S<:Real}
    F = float(promote_type(T, S))
    return hcubature_(f, SVector{n,F}(a), SVector{n,F}(b), norm, rtol, atol, maxevals, initdiv)
end
function hcubature_(f, a::AbstractVector{<:Real}, b::AbstractVector{<:Real},
           norm, rtol, atol, maxevals, initdiv)
    n = length(a)
    n == length(b) || throw(DimensionMismatch("endpoints $a and $b must have the same length"))
    hcubature_(f, SVector{n}(a), SVector{n}(b), norm, rtol, atol, maxevals, initdiv)
end
hcubature_(f, a::Tuple{Vararg{Real,n}}, b::Tuple{Vararg{Real,n}}, norm, rtol, atol, maxevals, initdiv) where {n} =
    hcubature_(f, SVector{n}(a), SVector{n}(b), norm, rtol, atol, maxevals, initdiv)

"""
    hcubature(f, a, b; norm=norm, rtol=sqrt(eps), atol=0, maxevals=typemax(Int), initdiv=1)

Compute the n-dimensional integral of f(x), where `n == length(a) == length(b)`,
over the hypercube whose corners are given by the vectors (or tuples) `a` and `b`.
That is, dimension `x[i]` is integrated from `a[i]` to `b[i]`.  The
return value of `hcubature` is a tuple `(I, E)` of the estimated integral
`I` and an estimated error `E`.

`f` should be a function `f(x)` that takes an n-dimensional vector `x`
and returns the integrand at `x`.   The integrand can be any type that supports
`+`, `-`, `*` real, and `norm` functions.  For example, the integrand
can be real or complex numbers, vectors, matrices, etcetera.

The integrand `f(x)` will be always be passed an `SVector{n,T}`,
where `SVector` is an efficient vector type defined in the `StaticArrays`
package and `T` is a floating-point type determined by promoting
the endpoint `a` and `b` coordinates to a floating-point type.
(Your integrand `f` should be type-stable: it should always return
a value of the same type, given this type of `x`.)

The integrand will never be evaluated exactly at the boundaries of the
integration volume.  (So, for example, it is possible to have an
integrand that blows up at the boundaries, as long as the integral
is finite, though such singularities will slow convergence.)

The integration volume is adaptively subdivided, using a cubature
rule due to Genz and Malik (1980), until the estimated error `E`
satisfies `E ≤ max(rtol*norm(I), atol)`, i.e. `rtol` and `atol` are
the relative and absolute tolerances requested, respectively.
It also stops if the number of `f` evaluations exceeds `maxevals`.
If neither `atol` nor `rtol` are specified, the
default `rtol` is the square root of the precision `eps(T)`
of the coordinate type `T` described above.
Initially, the volume is divided into `initdiv` segments along each dimension.

The error is estimated by `norm(I - I′)`, where `I′` is an alternative
estimated integral (via an "embedded" lower-order cubature rule.)
By default, the norm function used (for both this and the convergence
test above) is `norm`, but you can pass an alternative norm by
the `norm` keyword argument.  (This is especially useful when `f`
returns a vector of integrands with different scalings.)
"""
hcubature(f, a, b; norm=norm, rtol::Real=0, atol::Real=0,
                   maxevals::Integer=typemax(Int), initdiv::Integer=1) =
    hcubature_(f, a, b, norm, rtol, atol, maxevals, initdiv)

"""
    hquadrature(f, a, b; norm=norm, rtol=sqrt(eps), atol=0, maxevals=typemax(Int), initdiv=1)

Compute the integral of f(x) from `a` to `b`.  The
return value of `hcubature` is a tuple `(I, E)` of the estimated integral
`I` and an estimated error `E`.

The other parameters are the same as [`hcubature`](@ref).  `hquadrature`
is just a convenience wrapper around `hcubature` so that you can work
with scalar `x`, `a`, and `b`, rather than 1-component vectors.

Alternatively, for 1d integrals you can import the [`QuadGK`](@ref) module
and call the [`quadgk`](@ref) function, which provides additional flexibility
e.g. in choosing the order of the quadrature rule.
"""
hquadrature(f, a, b; norm=norm, rtol::Real=0, atol::Real=0,
                     maxevals::Integer=typemax(Int), initdiv::Integer=1) =
    hcubature_(x -> f(x[1]), SVector(a), SVector(b), norm, rtol, atol, maxevals, initdiv)

end # module


text:
using Pkg.Artifacts

artifact_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
data_hash = artifact_hash("data", artifact_toml)

if data_hash == nothing || !artifact_exists(data_hash)
    data_hash = create_artifact() do artifact_dir
        data_url_base = "https://coffeebeetle.s3.eu-central-1.amazonaws.com"
        download("$(data_url_base)/data", joinpath(artifact_dir, "data"))
    end
    bind_artifact!(artifact_toml, "data", data_hash)
end

datafile = joinpath(artifact"data", "data")


text:
nlp_test = Dict()

nlp_test[1] = function (optimizer;opt...)

    m = SimpleNL.Model(optimizer;opt...)

    x = variable(m;lb=2,ub=4)
    objective(m,(x-1)^2)

    optimize!(m)

    return compare(value(x), 2.)
end

nlp_test[2] = function (optimizer;opt...)

    m = SimpleNL.Model(optimizer;opt...)

    x = variable(m)
    constraint(m,x+1<= -1)
    constraint(m,x+1>= -1)
    constraint(m,-2<= x+1)
    constraint(m,0>=x+1)

    optimize!(m)

    return compare(value(x), -2)
end

nlp_test[3] = function (optimizer;opt...)

    m = SimpleNL.Model(optimizer;opt...)

    x = [variable(m,name="s[$i]",start = .1) for i=1:3]

    c = constraint(m,0 ==x[1]+sin(x[2])-x[3]/2 + 1. )
    objective(m,x[2]^2 + x[2]^2 + 1.)

    optimize!(m)

    return compare(value.(x), [-0.74,0,.52]) && compare(dual(c),0.)
end

nlp_test[4] = function (optimizer;opt...)

    m = SimpleNL.Model(optimizer;opt...)

    x = [variable(m;start=mod(i,2)==1 ? -1.2 : 1.) for i=1:10];
    p = [parameter(m,2) for i=1:10];
    setvalue.(p,2)

    objective(m,sum(100(x[i-1]^2-x[i])^2+(x[i-1]-1)^2 for i=2:10))

    c = [ constraint(m,3x[i+1]^3+p[i]*x[i+2]-5+sin(x[i+1]-x[i+2])sin(x[i+1]+x[i+2])+4x[i+1]-x[i]exp(x[i]-x[i+1])-3 == 0) for i=1:8 ];

    optimize!(m)
    return compare(value.(x), [-0.95055636, 0.91390082, 0.98909052, 0.99855924, 0.99980874, 0.99997459, 0.99999662, 0.99999955, 0.99999994, 0.99999993])
end

m = SimpleNL.Model()

x = variable(m)
y = variable(m)
setvalue(x,3)
set_lower_bound(x,-1)
set_upper_bound(x,1)
@test SimpleNL.index(x) == 1
@test value(x) == 3
@test lower_bound(x) == -1
@test upper_bound(x) == 1

c = constraint(m,x==1.)
constraint(m,x>=1.)
constraint(m,x<=1.)

constraint(m,1. ==y)
constraint(m,1. <=y)
constraint(m,-1. >=y)

constraint(m,x==y)
constraint(m,x-4. <=y)
constraint(m,x+4. >=y)

setvalue(c,2)
set_lower_bound(c,-1)
set_upper_bound(c,1)
@test value(c) == 2
@test lower_bound(c) == -1
@test upper_bound(c) == 1

p = parameter(m,1.)
setvalue(p,4)
@test value(p) == 4

objective(m,sum([p,x,y]))

m[:test] = "test"

@test objective_value(m) == 7.
@test num_variables(m) == 2
@test num_constraints(m) == 9
@test m[:test] == "test"

for (optimizer,opt) in [(Ipopt.Optimizer,[:print_level=>0]),
                        (MadNLP.Optimizer,[:print_level=>MadNLP.ERROR])]
    for f in values(nlp_test)
        @test f(optimizer;opt...)
    end
end


text:
abstract type AbstractNumberFormatter end

const float_regex = r"(?'mantissa'(?'before_dp'(?'sign'-?)(?'before_dp_nosign'\d+))(\.(?'after_dp'\d+))?)(?'e_or_E'e)(?'raw_exp'\+?0*(?'exp'(?'sign_exp'-?)\d+))"i

struct PlainNumberFormatter <: AbstractNumberFormatter end

(::PlainNumberFormatter)(x) = string(x)

struct PrintfNumberFormatter <: AbstractNumberFormatter
    fmt::String
end

(f::PrintfNumberFormatter)(x) = @eval @sprintf $(f.fmt) $x

struct StyledNumberFormatter <: AbstractNumberFormatter
    fmt::String

    StyledNumberFormatter(fmt::String="%.4g") = new(fmt)

end

function StyledNumberFormatter(significant_digits::Int)
    return StyledNumberFormatter("%.$(significant_digits)g")
end

(f::StyledNumberFormatter)(x) = string(x)
function (f::StyledNumberFormatter)(x::AbstractFloat)
    s = @eval @sprintf $(f.fmt) $x
    return replace(s, float_regex => s"\g<mantissa> \\mathrm{\g<e_or_E>} \g<exp>")
end

(f::StyledNumberFormatter)(x::Unsigned) = "\\mathtt{0x$(string(x; base=16, pad=2sizeof(x)))}"

struct FancyNumberFormatter <: AbstractNumberFormatter
    fmt::String
    exponent_format::SubstitutionString

    function FancyNumberFormatter(fmt::String="%.4g",
                                  exponent_format::SubstitutionString=s"\g<mantissa> \\cdot 10^{\g<exp>}")
        return new(fmt, exponent_format)
    end
end

function FancyNumberFormatter(fmt::String, mult_symbol)
    return FancyNumberFormatter(fmt, SubstitutionString("\\g<mantissa> $(escape_string(mult_symbol)) 10^{\\g<exp>}"))
end
function FancyNumberFormatter(significant_digits, mult_symbol="\\cdot")
    return FancyNumberFormatter("%.$(significant_digits)g", mult_symbol)
end

(f::FancyNumberFormatter)(x) = string(x)

function (f::FancyNumberFormatter)(x::AbstractFloat)
    s = @eval @sprintf $(f.fmt) $x
    return replace(s, float_regex => f.exponent_format)
end

(f::FancyNumberFormatter)(x::Unsigned) = "\\mathtt{0x$(string(x; base=16, pad=2sizeof(x)))}"


text:
Random.seed!(123)

@testset "Schelling example" begin
  model = instantiate_modelS(numagents=370, griddims=(20,20), min_to_be_happy=3)
  agent_properties = [:pos, :mood, :group]
  steps_to_collect_data = collect(1:2)
  data = step!(agent_step!, model, 2, agent_properties, steps_to_collect_data)

  @test data[1, :pos_1] == 261
  @test data[1, :pos_2] == 363

  i = 1
  visualize_2D_agent_distribution(data, model, Symbol("pos_$i"), types=Symbol("group_$i"), savename="step_$i", cc=Dict(0=>"blue", 1=>"red"))

  rm("./step_1.pdf")

  @test Agents.grid0D() == Agents.Graph(1)
end