text:
using Minibatch
using Revtok
using DataStructures
using GPUArrays
using NNlib
using Flux

english = """This is an account of how magical thinking made us modern.
When people talk about magical thinking, it is usually as a cognitive feature of children, uneducated people, the mushy-minded, or the mentally ill.
If we notice magical thinking in ourselves, it is with a pang of shame: literate adults are supposed to be more sophisticated than that.
At the same time, magical thinking is obviously rampant in the world.
It’s hard not to be fascinated, even if it’s a horrified fascination.
I think that the reason that it has been so difficult to precisely define “magical thinking” is that what we call “magical thinking” is a collection of stigmatized examples of a more general, and generally useful, cognitive capacity.
This is the ability to think in “as if” mode: “as if” inanimate objects had minds, “as if” thoughts could affect reality, “as if” symbols had power over their referents."""

spanish = """Este es un relato de cómo el pensamiento mágico nos hizo modernos.
Cuando la gente habla de pensamiento mágico, generalmente es como un rasgo cognitivo de los niños, de las personas sin educación, de los que tienen la mente blanda o de los enfermos mentales.
Si nos damos cuenta del pensamiento mágico en nosotros mismos, es con una espiral de vergüenza: se supone que los adultos alfabetizados son más sofisticados que eso.
Al mismo tiempo, el pensamiento mágico es obviamente desenfrenado en el mundo.
Es difícil no estar fascinado, aunque sea una fascinación horrorizada.
Creo que la razón por la que ha sido tan difícil definir con precisión el “pensamiento mágico” es que lo que llamamos "pensamiento mágico" es una colección de ejemplos estigmatizados de una capacidad cognitiva más general, y en general útil.
Esta es la capacidad de pensar en modo “como si”: “como si” los objetos inanimados tuvieran mente, “como si” los pensamientos pudieran afectar a la realidad, “como si” los símbolos tuvieran poder sobre sus referentes."""

struct Vocab
    itos::Vector
    stoi::Associative
end
Base.length(vocab::Vocab) = length(vocab.itos)

function process(corpus, vocabsize=200, batchsize=4)
    corpus = tokenize.(split(corpus, '\n'))
    segments = buildvocab(counter(Iterators.flatten(corpus)), vocabsize)
    corpus = segment.(corpus, segments)
    itos = keys(segments)
    stoi = Dict((reverse(p) for p in enumerate(itos)))
    vocab = Vocab(collect(itos), stoi)
    corpus = [[vocab.stoi[x] for x in ex] for ex in corpus]
    corpus = Iterators.partition(corpus, batchsize)
    corpus = map(b -> MaskedBatch(b, (true,)), corpus)
    return corpus, vocab
end

# first step to patching issues that come up when we try to use Flux.Tracker
function Base.getindex(a::Flux.TrackedArray{T, 2} where T, ::Colon, b::MaskedBatch{<:AbstractVector})
    return invoke(getindex, Tuple{AbstractMatrix, Colon, MaskedBatch{<:AbstractVector}}, a, :, b)
end

function posenc(timestep::Integer, channel::Integer, nchannels::Integer)
    if iseven(channel)
        return sin(timestep/(10000^(channel/nchannels)))
    else
        return cos(timestep/(10000^((channel-1)/nchannels)))
    end
end

function posenc(idx::CartesianIndex, nchannels::Integer)
    return posenc(idx[2], idx[1], nchannels)
end

function posenc!(A::GPUArray, state, nchannels::Integer)
    idx = @cartesianidx A state
    @inbounds A[idx] += posenc(idx, nchannels)
    return A
end

function posenc!(A::GPUArray)
    nchannels = size(A, 1)
    gpucall(posenc, A, (nchannels,))
    return A
end

function posenc!(A::AbstractArray)
    nchannels = size(A, 1)
    for idx in CartesianRange(size(A))
        # @inbounds
        A[idx] += posenc(idx, nchannels)
    end
    return A
end

function posenc!(B::MaskedBatch)
    posenc!(B.data)
    B.data = B.data .* B.mask # TODO in-place
    return B
end

function posenc!(A::Flux.TrackedArray)
    return posenc!(Flux.Tracker.data(A))
end

struct Embedding
    W
end
Embedding(vocabsize, embedsize) = Embedding(param(randn(embedsize, vocabsize)))
(embedding::Embedding)(x) = embedding.W[:, x]
Flux.treelike(Embedding)

# struct LayerNorm
#     γ
#     β
# end
# LayerNorm(nchannels) = LayerNorm(param(ones(nchannels)), param(zeros(nchannels)))
# (l::LayerNorm)(x) = l.γ .* (x .- mean(x, 1)) ./ (std(x, 1) .+ ϵ) .+ l.β

FeedForward(d_model, d_hidden) = Chain(Dense(d_model, d_hidden, relu), Dense(d_hidden, d_model))

struct ResidualBlock
    l
    norm
    # TODO dropout
end
ResidualBlock(l, d_model::Integer) = ResidualBlock(l, LayerNorm(d_model))
(l::ResidualBlock)(x...) = x[1] .+ l.norm(l.l(x...))
Flux.treelike(ResidualBlock)

struct Attention
    scale
    # TODO dropout
    # TODO causal
end
Attention(d_key::Integer, causal=false) = Attention(sqrt(d_key))
function (l::Attention)(q, k, v)
    alpha = k'*q
    # TODO causal mask
    return v*softmax(alpha, 1)
end
Flux.mapchildren(f, l::Attention) = Attention(l.scale)

# TODO multihead

struct EncoderLayer
    selfattn
    feedforward
end
EncoderLayer(d_model::Integer, d_hidden::Integer) = EncoderLayer(
    ResidualBlock(Attention(d_model), d_model),
    ResidualBlock(FeedForward(d_model, d_hidden), d_model))
(l::EncoderLayer)(x) = l.feedforward(l.selfattn(x, x, x))
Flux.treelike(EncoderLayer)

struct Encoder
    embed
    layers
end
Encoder(d_model, d_hidden, n_layers, vocabsize) = Encoder(
    Embedding(vocabsize, d_model),
    Chain((EncoderLayer(d_model, d_hidden) for _ in 1:n_layers)...))
function (l::Encoder)(x)
    x = posenc!(l.embed(x))
    encoding = [x = l(x) for l in l.layers.layers]
    return encoding
end
Flux.treelike(Encoder)

#corpus = zip(english, spanish)
en, vocab_en = process(english)
es, vocab_es = process(spanish)

# d = 4
# embed = Embedding(length(vocab_en), d)
# layernorm = LayerNorm(d)
# linear = Dense(d, d)
# feedforward = FeedForward(d, 2d)
# residual = ResidualBlock(FeedForward(d, 2d), d)
#
# x = en[1]
# x = embed(x)
# x = posenc!(x)
# x = layernorm(x)
# x = linear(x)
# x = feedforward(x)
# x = residual(x)

encoder = Encoder(4, 8, 3, length(vocab_en))
encoder = Flux.mapleaves(Flux.Tracker.data, encoder)
encoder(en[1])


text:
export get_current_user,
    get_user,
    modify_current_user,
    get_current_user_guilds,
    leave_guild,
    create_dm,
    create_group_dm

"""
    get_current_user(c::Client) -> User

Get the [`Client`](@ref) [`User`](@ref).
"""
function get_current_user(c::Client)
    return Response{User}(c, :GET, "/users/@me")
end

"""
    get_user(c::Client, user::Integer) -> User

Get a [`User`](@ref).
"""
function get_user(c::Client, user::Integer)
    return Response{User}(c, :GET, "/users/$user")
end

"""
    modify_current_user(c::Client; kwargs...) -> User

Modify the [`Client`](@ref) [`User`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/user#modify-current-user).
"""
function modify_current_user(c::Client; kwargs...)
    return Response{User}(c, :PATCH, "/users/@me"; body=kwargs)
end

"""
    get_user_guilds(c::Client; kwargs...) -> Vector{Guild}

Get a list of [`Guild`](@ref)s the [`Client`](@ref) [`User`](@ref) is a member of.
More details [here](https://discordapp.com/developers/docs/resources/user#get-current-user-guilds).
"""
function get_current_user_guilds(c::Client; kwargs...)
    return Response{Vector{Guild}}(c, :GET, "/users/@me/connections"; kwargs...)
end

"""
    leave_guild(c::Client, guild::Integer)

Leave a [`Guild`](@ref).
"""
function leave_guild(c::Client, guild::Integer)
    return Response(c, :DELETE, "/users/@me/guilds/$guild")
end

"""
    create_dm(c::Client; kwargs...) -> DiscordChannel

Create a DM [`DiscordChannel`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/user#create-dm).
"""
function create_dm(c::Client; kwargs...)
    return Response{DiscordChannel}(c, :POST, "/users/@me/channels"; body=kwargs)
end

"""
    create_group_dm(c::Client; kwargs...) -> DiscordChannel

Create a group DM [`DiscordChannel`](@ref).
More details [here](https://discordapp.com/developers/docs/resources/user#create-group-dm).
"""
function create_group_dm(c::Client; kwargs...)
    return Response{DiscordChannel}(c, :POST, "/users/@me/channels"; body=kwargs)
end

ext: jl


text:
using ImageView, FixedPointNumbers, Reactive
using Base.Test

@testset "contrast GUI" begin
    # test for the fix in #119
    img = rand(N0f16, 100, 100)
    clim = Signal(CLim(extrema(img)...))
    imgsig = Signal(img)
    enabled, histsig, imgc = ImageView.prep_contrast(imgsig, clim)
    push!(enabled, true)
    yield()
    ret = ImageView.contrast_gui(enabled, histsig, clim)
    @test isa(ret, Dict)
    Gtk.destroy(ret["window"])
end


text:
_parse(::Type{T}, x::AbstractString) where {T} = parse(T, x)
_parse(::Type{String}, x::AbstractString) = String(x)

_read(io::IO, ::Type{T}) where {T} = _parse(T, strip(readline(io)))

function _read(io::IO, args...)
    x = split(strip(readline(io)))
    @assert length(x) == length(args)
    return _parse.(args, x)
end

mutable struct _CBFReadData
    scalar_vars::Vector{MOI.VariableIndex}
    psd_vars::Vector{Vector{MOI.VariableIndex}}
    psd_side_dims::Vector{Int}
    psd_row_starts::Vector{Int}
    psd_row_terms::Vector{Vector{MOI.ScalarAffineTerm{Float64}}}
    psd_row_constants::Vector{Float64}
    obj_terms::Vector{MOI.ScalarAffineTerm{Float64}}
    obj_constant::Float64
    power_cone_alpha::Vector{Vector{Float64}}
    dual_power_cone_alpha::Vector{Vector{Float64}}
    con_cones::Vector{Tuple{String,Int}}
    row_terms::Vector{Vector{MOI.ScalarAffineTerm{Float64}}}
    row_constants::Vector{Float64}

    function _CBFReadData()
        return new(
            MOI.VariableIndex[],
            Vector{MOI.VariableIndex}[],
            Int[],
            Int[],
            Vector{MOI.ScalarAffineTerm{Float64}}[],
            Float64[],
            MOI.ScalarAffineTerm{Float64}[],
            0.0,
            Vector{Float64}[],
            Vector{Float64}[],
            Tuple{String,Int}[],
            Vector{MOI.ScalarAffineTerm{Float64}}[],
            Float64[],
        )
    end
end

# Convert a pair of row and column indices of a symmetric matrix into a vector
# index for the row-wise lower triangle
function _mat_to_vec_idx(i::Int, j::Int)
    if i < j
        return div((j - 1) * j, 2) + i
    else
        return div((i - 1) * i, 2) + j
    end
end

function _cbf_to_moi_cone(
    data::_CBFReadData,
    cone_str::AbstractString,
    cone_dim::Int,
)
    if cone_str == "F"
        return MOI.Reals(cone_dim)
    elseif cone_str == "L="
        return MOI.Zeros(cone_dim)
    elseif cone_str == "L-"
        return MOI.Nonpositives(cone_dim)
    elseif cone_str == "L+"
        return MOI.Nonnegatives(cone_dim)
    elseif cone_str == "Q"
        @assert cone_dim >= 2
        return MOI.SecondOrderCone(cone_dim)
    elseif cone_str == "QR"
        @assert cone_dim >= 3
        return MOI.RotatedSecondOrderCone(cone_dim)
    elseif cone_str == "EXP"
        return MOI.ExponentialCone()
    elseif cone_str == "EXP*"
        return MOI.DualExponentialCone()
    elseif startswith(cone_str, "@")
        raw_powcone_info = split(cone_str[2:end], ":")
        powcone_idx = parse(Int, raw_powcone_info[1]) + 1
        if raw_powcone_info[2] == "POW"
            alpha = data.power_cone_alpha[powcone_idx]
            return MOI.PowerCone{Float64}(first(alpha) / sum(alpha))
        elseif raw_powcone_info[2] == "POW*"
            alpha = data.dual_power_cone_alpha[powcone_idx]
            return MOI.DualPowerCone{Float64}(first(alpha) / sum(alpha))
        end
    end
    return error("CBF cone name $cone_str is not recognized or supported.")
end

function _read_VER(io::IO)
    version = _read(io, Int)
    if !(1 <= version <= 3)
        error("CBF version number $(version) is not yet supported.")
    end
    return
end

function _read_POWCONES(io::IO, model::Model, alpha::Vector{Vector{Float64}})
    num_powcone, num_lines = _read(io, Int, Int)
    alpha_idx = 0
    for _ in 1:num_powcone
        num_alpha = _read(io, Int)
        if num_alpha != 2
            error("Only 3-dimensional power cones are supported.")
        end
        push!(alpha, [_read(io, Float64) for _ in 1:num_alpha])
        alpha_idx += num_alpha
    end
    @assert num_lines == alpha_idx
    return
end

function _read_OBJSENSE(io::IO, model::Model)
    obj_sense = strip(readline(io))
    if obj_sense == "MIN"
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    else
        @assert obj_sense == "MAX"
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    end
    return
end

function _read_PSDVAR(io::IO, model::Model, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        side_dim = _read(io, Int)
        cone_dim = div(side_dim * (side_dim + 1), 2)
        psd_vars_k = MOI.add_variables(model, cone_dim)
        push!(data.psd_vars, psd_vars_k)
        MOI.add_constraint(
            model,
            MOI.VectorOfVariables(psd_vars_k),
            MOI.PositiveSemidefiniteConeTriangle(side_dim),
        )
    end
    return
end

function _read_VAR(io::IO, model::Model, data::_CBFReadData)
    num_var, num_lines = _read(io, Int, Int)
    append!(data.scalar_vars, MOI.add_variables(model, num_var))
    var_idx = 0
    for _ in 1:num_lines
        cone_str, cone_dim = _read(io, String, Int)
        if cone_str == "F"
            var_idx += cone_dim
            continue # Free cones (no constraint).
        end
        indices = 1:cone_dim
        if cone_str == "EXP" || cone_str == "EXP*"
            indices = (3:-1:1)
        end
        MOI.add_constraint(
            model,
            MOI.VectorOfVariables(data.scalar_vars[var_idx.+indices]),
            _cbf_to_moi_cone(data, cone_str, cone_dim),
        )
        var_idx += cone_dim
    end
    @assert var_idx == num_var
    return
end

function _read_INT(io::IO, model::Model, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        MOI.add_constraint(
            model,
            data.scalar_vars[_read(io, Int)+1],
            MOI.Integer(),
        )
    end
    return
end

function _read_PSDCON(io::IO, data::_CBFReadData)
    idx = 0
    for _ in 1:_read(io, Int)
        side_dim = _read(io, Int)
        push!(data.psd_side_dims, side_dim)
        push!(data.psd_row_starts, idx)
        idx += div(side_dim * (side_dim + 1), 2)
    end
    for _ in 1:idx
        push!(data.psd_row_terms, MOI.ScalarAffineTerm{Float64}[])
        push!(data.psd_row_constants, 0.0)
    end
    return
end

function _read_CON(io::IO, data::_CBFReadData)
    num_rows, num_lines = _read(io, Int, Int)
    row_idx = 0
    for _ in 1:num_lines
        cone_str, cone_dim = _read(io, String, Int)
        push!(data.con_cones, (cone_str, cone_dim))
        row_idx += cone_dim
    end
    for _ in 1:num_rows
        push!(data.row_terms, MOI.ScalarAffineTerm{Float64}[])
        push!(data.row_constants, 0.0)
    end
    return
end

function _read_OBJFCOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        psd_var_idx, i, j, coef = _read(io, Int, Int, Int, Float64)
        if i != j
            coef += coef # scale off-diagonals
        end
        push!(
            data.obj_terms,
            MOI.ScalarAffineTerm{Float64}(
                coef,
                data.psd_vars[psd_var_idx+1][_mat_to_vec_idx(i + 1, j + 1)],
            ),
        )
    end
    return
end

function _read_OBJACOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        var_idx, coef = _read(io, Int, Float64)
        push!(
            data.obj_terms,
            MOI.ScalarAffineTerm{Float64}(coef, data.scalar_vars[var_idx+1]),
        )
    end
    return
end

function _read_OBJBCOORD(io::IO, data::_CBFReadData)
    data.obj_constant += _read(io, Float64)
    return
end

function _read_FCOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        row_idx, psd_var_idx, i, j, coef =
            _read(io, Int, Int, Int, Int, Float64)
        if i != j
            coef += coef # scale off-diagonals
        end
        push!(
            data.row_terms[row_idx+1],
            MOI.ScalarAffineTerm{Float64}(
                coef,
                data.psd_vars[psd_var_idx+1][_mat_to_vec_idx(i + 1, j + 1)],
            ),
        )
    end
    return
end

function _read_ACOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        row, col, coef = _read(io, Int, Int, Float64)
        push!(
            data.row_terms[row+1],
            MOI.ScalarAffineTerm{Float64}(coef, data.scalar_vars[col+1]),
        )
    end
    return
end

function _read_BCOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        row, coef = _read(io, Int, Float64)
        data.row_constants[row+1] = coef
    end
    return
end

function _read_HCOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        psd_idx, var_idx, i, j, coef = _read(io, Int, Int, Int, Int, Float64)
        row_idx = data.psd_row_starts[psd_idx+1] + _mat_to_vec_idx(i + 1, j + 1)
        push!(
            data.psd_row_terms[row_idx],
            MOI.ScalarAffineTerm{Float64}(coef, data.scalar_vars[var_idx+1]),
        )
    end
    return
end

function _read_DCOORD(io::IO, data::_CBFReadData)
    for _ in 1:_read(io, Int)
        psd_idx, i, j, coef = _read(io, Int, Int, Int, Float64)
        row_idx = data.psd_row_starts[psd_idx+1] + _mat_to_vec_idx(i + 1, j + 1)
        data.psd_row_constants[row_idx] += coef
    end
    return
end

"""
    Base.read!(io::IO, model::FileFormats.CBF.Model)

Read `io` in the Conic Benchmark Format and store the result in `model`.
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end
    data = _CBFReadData()
    while !eof(io)
        line = strip(readline(io))
        if isempty(line) || startswith(line, "#")
            # Skip blank lines and comments.
        elseif line == "VER"
            _read_VER(io)
        elseif line == "POWCONES"
            _read_POWCONES(io, model, data.power_cone_alpha)
        elseif line == "POW*CONES"
            _read_POWCONES(io, model, data.dual_power_cone_alpha)
        elseif line == "OBJSENSE"
            _read_OBJSENSE(io, model)
        elseif line == "PSDVAR"
            _read_PSDVAR(io, model, data)
        elseif line == "VAR"
            _read_VAR(io, model, data)
        elseif line == "INT"
            _read_INT(io, model, data)
        elseif line == "PSDCON"
            _read_PSDCON(io, data)
        elseif line == "CON"
            _read_CON(io, data)
        elseif line == "OBJFCOORD"
            _read_OBJFCOORD(io, data)
        elseif line == "OBJACOORD"
            _read_OBJACOORD(io, data)
        elseif line == "OBJBCOORD"
            _read_OBJBCOORD(io, data)
        elseif line == "FCOORD"
            _read_FCOORD(io, data)
        elseif line == "ACOORD"
            _read_ACOORD(io, data)
        elseif line == "BCOORD"
            _read_BCOORD(io, data)
        elseif line == "HCOORD"
            _read_HCOORD(io, data)
        elseif line == "DCOORD"
            _read_DCOORD(io, data)
        else
            error("Failed to parse CBF file due to corrupted line: $line")
        end
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(data.obj_terms, data.obj_constant),
    )

    # Non-PSD constraints.
    row_idx = 0
    for (cone_str, cone_dim) in data.con_cones
        con_func = if cone_str == "EXP" || cone_str == "EXP*"
            # Reverse order of indices.
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm{Float64}(4 - l, t) for
                    l in 1:cone_dim for t in data.row_terms[row_idx+l]
                ],
                data.row_constants[row_idx.+(3:-1:1)],
            )
        else
            MOI.VectorAffineFunction(
                [
                    MOI.VectorAffineTerm{Float64}(l, t) for l in 1:cone_dim
                    for t in data.row_terms[row_idx+l]
                ],
                data.row_constants[row_idx.+(1:cone_dim)],
            )
        end
        con_set = _cbf_to_moi_cone(data, cone_str, cone_dim)
        MOI.add_constraint(model, con_func, con_set)
        row_idx += cone_dim
    end

    # PSD constraints.
    for psd_idx in eachindex(data.psd_side_dims)
        row_start = data.psd_row_starts[psd_idx]
        side_dim = data.psd_side_dims[psd_idx]
        cone_dim = div(side_dim * (side_dim + 1), 2)
        con_func = MOI.VectorAffineFunction(
            [
                MOI.VectorAffineTerm{Float64}(l, t) for l in 1:cone_dim for
                t in data.psd_row_terms[row_start+l]
            ],
            data.psd_row_constants[row_start.+(1:cone_dim)],
        )
        MOI.add_constraint(
            model,
            con_func,
            MOI.PositiveSemidefiniteConeTriangle(side_dim),
        )
    end

    return
end


text:
# * Key type

mutable struct Key
    legend::String
    width::Float64
    height::Float64
end

function gen_key(args)
    args = split(args, ":")
    legend = get(args, 1, "")
    width = parse(Float64, get(args, 2, "1"))
    height = parse(Float64, get(args, 3, "1"))
    Key(legend, width, height)
end

Base.show(io::IO, key::Key) = write(io, " ", key.legend, " ")

const tikz_legends = Dict("ctrl" => "\$\ue0fb\$",
                          "ret" => "\$\ue0fc\$",
                          "esc" => "\$\ue0fd\$",
                          "cmd" => "\$\ue0fe\$",
                          "tab" => "\$\ue0ff\$",
                          "space" => "\$\ue100\$",
                          "del" => "\$\ue101\$",
                          "alt" => "\$\ue102\$",
                          "option" => "\$\ue103\$",
                          "enter" => "\$\ue105\$",
                          "shift" => "\$\ue106\$",
                          "mod1" => "\$\ue107\$",
                          "mod2" => "\$\ue108\$")

function tikz_legend(legend::String)
    l = get(tikz_legends, legend, legend)
    isascii(l) ? "\\textsf{$l}" : l
end

Base.convert(::Type{TikZnode}, key::Key, U::Length) =
    TikZnode(tikz_legend(key.legend),
             "draw", "rectangle", "rounded corners",
             "minimum width" => key.width*U,
             "minimum height" => key.height*U)

# * Space type

mutable struct Space
    width::Float64
    height::Float64
end

Base.show(io::IO, space::Space) = write(io, "   ")

function gen_space(args)
    args = split(args, ":")
    a,b = (length(args) == 3 || length(args) == 2 && args[1] == "") ? (2,3) : (1,2)
    args[a] == "" && (args[a] = "1")
    width = parse(Float64, get(args, a, "1"))
    height = parse(Float64, get(args, b, "1"))
    Space(width, height)
end

Base.convert(::Type{TikZnode}, space::Space, U::Length) =
    TikZnode("",
             "rectangle", "rounded corners",
             "minimum width" => space.width*U,
             "minimum height" => space.height*U)

# * Row type

const Row = Vector{Union{Key,Space}}

function gen_row(args)
    map(split(args, " ")) do arg
        isempty(arg) || arg[1] == ':' ? gen_space(arg) : gen_key(arg)
    end |> Row
end

# * Keyboard type

mutable struct Keyboard
    rows::Vector{Row}
    U::Length
    clearance::Length
    extra_clearance::Bool
    gravity::Symbol
end
Keyboard(rows::Vector{Row} = Row[]; U::Length = 18u"mm",
         clearance::Length = 0.2u"mm", extra_clearance::Bool = false,
         gravity::Symbol = :nw) =
             Keyboard(rows, U, clearance, extra_clearance, gravity)

function keyboard(fun::Function; kwargs...)
    kbd = Keyboard()
    fun(kbd)
    kbd
end

function nextrow!(kbd::Keyboard)
    push!(kbd.rows, Row())
end

function Base.push!(kbd::Keyboard, key::Union{Key,Space})
    isempty(kbd.rows) && nextrow!(kbd)
    if kbd.extra_clearance
        key.width += kbd.clearance*floor(Int,key.width - 1.0)/kbd.U |> NoUnits
        key.height += kbd.clearance*max(0,floor(Int,key.height - 1.0))/kbd.U |> NoUnits
    end
    push!(kbd.rows[end], key)
end

function Base.push!(kbd::Keyboard, row::Row)
    nextrow!(kbd)
    for key in row
        push!(kbd, key)
    end
end

function Base.show(io::IO, ::MIME"text/plain", kbd::Keyboard)
    for row in kbd.rows
        foreach(k -> show(io, k), row)
        println()
    end
end

# * String macros

macro key_str(args)
    kbd = esc(:kbd)
    quote
        push!($kbd, gen_key($args))
    end
end

macro space_str(args)
    kbd = esc(:kbd)
    quote
        push!($kbd, gen_space($args))
    end
end

macro row_str(args)
    kbd = esc(:kbd)
    row = gen_row(args)
    quote
        push!($kbd.rows, $row)
    end
end

macro keyboard_str(spec)
    rows = map(strip, split(spec, "\n"))
    isempty(last(rows)) && deleteat!(rows, length(rows))
    quote
        keyboard() do kbd
            for row in $rows
                push!(kbd, gen_row(row))
            end
        end
    end
end

function coordinates(kbd::Keyboard, x₀::Length=0u"mm", y₀::Length=0u"mm")
    direction = kbd.gravity == :nw ? 1 : -1
    @unpack U,clearance = kbd
    l = U + clearance
    i = 0
    map(enumerate(kbd.rows)) do (i,row)
        y = y₀ - (i-1)*l
        x = x₀
        map(enumerate(row)) do (j,key)
            ox = x
            x += direction*(clearance + key.width*U)
            Dict(:x => ox, :y => y,
                 :centerx => ox + direction*key.width*U/2,
                 :centery => y - direction*key.height*U/2,
                 :key => key)
        end
    end
end

function Base.unique!(v::Vector{T}, s::T, e::T) where T
    length(v) > length(s:e) &&
        throw(ArgumentError("Impossible to shift elements of vector with only $(length(s:e)) choices available"))

    while !allunique(v)
        function shift!(i::Int, d::Int)
            di,sel = if d < 0
                d = -d
                -1,d:i
            else
                1,i:d
            end
            for j in sel
                v[j] += di
            end
        end

        i = findfirst(i -> v[i]==v[i-1], 2:length(v))
        fwd = findfirst(j -> v[j] < (j < length(v) ? v[j+1]-1 : e), i+1:length(v))
        bwd = findfirst(j -> v[j] > (j > 1 ? v[j-1]+1 : s), 1:i-1)

        if fwd != nothing || bwd != nothing
            if fwd != nothing && bwd != nothing && fwd < abs(i-bwd) || fwd != nothing
                shift!(i+1, fwd+i)
            else
                shift!(i, -(bwd+1))
            end
        else
            throw(ArgumentError("Cannot continue, no shift will improve uniqueness"))
        end
    end

    v
end

function find_columns(kbd::Keyboard)
    coords = map(coordinates(kbd)) do row
        filter(k -> k[:key] isa Keyboards.Key, row)
    end
    num_columns = 0
    maxi = 0
    centers = map(enumerate(coords)) do (i,row)
        num_columns = max(num_columns, length(row))
        length(row) == num_columns && (maxi = i)
        map(row) do key
            key[:centerx],key[:centery]
        end
    end
    @info "$(num_columns) columns required"
    @info "Row # with max columns: $(maxi)"
    column_centers = first.(centers[maxi])

    key_mappings = map(coords) do row
        mapping = map(row) do key
            argmin(abs.(key[:centerx] .- column_centers))
        end
        unique!(mapping, 1, num_columns)
        map(enumerate(mapping)) do (i,c)
            c => row[i]
        end
    end

    columns = [x => Dict[] for x in column_centers]
    for (i,row) in enumerate(key_mappings)
        for k in row
            k[2][:row] = i
            push!(columns[k[1]][2], k[2])
        end
    end

    columns
end

# * TikZ conversion

function Base.convert(::Type{TikzPicture}, kbd::Keyboard; draw_centers::Bool=false, kwargs...)
    @unpack U,clearance = kbd
    gravity,direction = if kbd.gravity == :nw
        "north west",1
    else
        "north east",-1
    end
    coords = coordinates(kbd)
    rows = map(enumerate(kbd.rows)) do (i,row)
        y = first(coords[i])[:y]
        tikz_scope("yshift"=>y) do
            map(enumerate(row)) do (j,key)
                x = coords[i][j][:x]
                tikz_scope("xshift"=>x) do
                    key_node = convert(TikZnode, key, kbd.U)
                    push!(key_node.args, "anchor" => gravity)
                    kns = convert(MIME"text/tikz", key_node)
                    kns * if draw_centers
                        id = key isa Key ? "k"*string(hash(key)) : ""
                        center_node = TikZnode("", "red", "draw", "circle",
                                               id = id,
                                               x = coords[i][j][:centerx]-x,
                                               y = coords[i][j][:centery]-y)
                        convert(MIME"text/tikz", center_node)
                    else
                        ""
                    end
                end
            end |> k -> join(k, "\n")
        end
    end |> r -> join(r, "\n") |> indent
    columns = if draw_centers
        map(find_columns(kbd)) do (column,keys)
            c = map(keys) do key
                "(k$(hash(key[:key])))"
            end |> k -> join(k, " -- ")
            "\\draw[red] $(c);"
        end |> c -> join(c, "\n") |> indent
    else
        ""
    end
    TikzPicture(rows*"\n"*columns; kwargs...)
end

TikzPictures.save(f::S, kbd::Keyboard; kwargs...) where {S<:TikzPictures.SaveType} =
    save(f, convert(TikzPicture, kbd; preamble=preamble, kwargs...))

export @keyboard_str


text:
# get the diagonal block of the jacobian

#------------------------------------------------------------------------------
# Block diagonal matrix type

"""
  A special array type for holding a block-diagonal matrix
"""
struct DiagJac{T} <: AbstractArray{T, 2}
  A::Array{T, 3}
end

"""
  Outer constructor for [`DiagJac`](@ref).

  **Inputs**

   * T: element type of the matrix
   * blocksize: the blocksize of the diagonal blocks
   * nblock: number of blocks.

   The size of the matrix is `blocksize`*`nblock`.
   Only square block are supported at this time

"""
function DiagJac(::Type{T}, blocksize::Integer, nblock::Integer) where T

  A = zeros(T, blocksize, blocksize, nblock)

  return DiagJac{T}(A)
end

const NullDiagJac = DiagJac(PetscScalar, 0, 0)

import Base: length, size, getindex, setindex!, show

"""
  Returns the length of the array *as if it were a dense matrix*
"""
function length(A::DiagJac)

  blocksize, blocksize, nblock = size(A.A)
  n = blocksize*nblock
  return n*n
end

function size(A::DiagJac)

  blocksize, blocksize, nblock = size(A.A)

  n = blocksize*nblock
  return (n, n)
end

"""
  Internal function for computing the indices in the 3D array

  **Inputs**

   * A : DaigJac
   * i: first index
   * j: second index

  **Outputs**

   * i2: first index in 3D array
   * j2: second index in 3D array
   * i2block: 3rd index in 3D array
"""
function computeIndex(A::DiagJac, i::Integer, j::Integer)
# internal function for computing indices

  blocksize, blocksize, nblock = size(A.A)

  i2 = ((i-1) % blocksize ) + 1
  i2block = div(i-1, blocksize) + 1

  j2 = ((j-1) % blocksize) + 1
  j2block = div(j-1, blocksize) + 1

  @assert i2block <= nblock
  @assert i2block > 0
  @assert j2block <= nblock
  @assert j2block > 0

  @assert i2 > 0
  @assert i2 <= blocksize
  @assert j2 > 0
  @assert j2 <= blocksize

  @assert i2block == j2block

  return i2, j2, i2block
end

"""
  Indexing using non-block indices
"""
function getindex(A::DiagJac, i::Integer, j::Integer)

  i2, j2, i2block = computeIndex(A, i, j)
  return A.A[i2, j2, i2block]
end

"""
  Indexing using non-block indices
"""
function setindex!(A::DiagJac, v, i::Integer, j::Integer)

  i2, j2, i2block = computeIndex(A, i, j)
  return A.A[i2, j2, i2block] = v
end

function show(io::IO, A::DiagJac)

  m, n = size(A)
  println("DiagJac: $m x $n block diagonal matrix with block size: $(size( A.A, 1)), number of blocks: $(size(A.A, 3))")
  for i=1:size(A.A, 3)
    println("block ", i, " = ")
    println(A.A[:, :, i])
  end

  return nothing
end

import PETSc2.MatZeroEntries
function MatZeroEntries(A::DiagJac)

  fill!(A.A, 0.0)
end

"""
  Matrix-vector multiplication function for [`DiagJac`](@ref)

  **Inputs**

   * A: a `DiagJac`
   * mesh: the mesh object for the discretiztaion where `A` was evaluated
   * x: vector to multiply against

  **Inputs/Outputs**

   * b: vector to be overwritten with the result

  **Keyword Arguments**

   * zero_output: if true, overwrite `b`, otherwise add to it, default true
   * trans: if true, multiply by the transposed matrix, default false

  Aliasing Restrictions: no aliasing allowed
"""
function diagMatVec(A::DiagJac, mesh::AbstractDGMesh, x::AbstractVector, b::AbstractVector; zero_output=true, trans::Bool=false)

  blocksize, blocksize, nblock = size(A.A)

  @assert blocksize == mesh.numDofPerNode*mesh.numNodesPerElement
  if zero_output
    fill!(b, 0)
  end

  _trans::Bool = trans  # make type-stable

  # work arrays
  xblock = zeros(eltype(x), blocksize)
  bblock = zeros(eltype(b), blocksize)
  for i=1:nblock
    getValues(mesh, x, i, xblock)
    Ablock = sview(A.A, :, :, i)
    if _trans
      smallmatTvec!(Ablock, xblock, bblock)
    else
      smallmatvec!(Ablock, xblock, bblock)
    end
    setValues(mesh, bblock, i, b)
  end  # end loop i

  return nothing
end

#------------------------------------------------------------------------------
# utility functions
"""
  Gets the values out of a `eqn.q_vec` shaped vector for a given element.

  **Inputs**

   * mesh
   * x: the vector
   * elnum: element number

  **Inputs/Outputs**

   * workvec: vector of length `mesh.numDofPerNode*mesh.numNodesPerElement`
              to overwrite with the solution.
"""
function getValues(mesh::AbstractMesh, x::AbstractVector, elnum::Integer,
                   workvec::AbstractVector)

  # get the x values
  # this could be faster if we assume the dofs on each element are numbered
  # sequentially
  pos = 1
  @simd for j=1:mesh.numNodesPerElement
    @simd for k=1:mesh.numDofPerNode
      dof_i = mesh.dofs[k, j, elnum]
      workvec[pos] = x[dof_i]
      pos += 1
    end
  end

  return nothing
end

"""
  Inverse of [`getValues`](@ref), puts values back into a vector

  **Inputs**

   * mesh
   * workvec
   * elnum

  **Inputs/Outputs**

   * b: relevent entries are added to (not overwritten)
"""
function setValues(mesh::AbstractMesh, workvec::AbstractVector, elnum::Integer,
                   b::AbstractVector)

  # put entries back into b using mesh.dofs
  pos = 1
  for j=1:mesh.numNodesPerElement
    for k=1:mesh.numDofPerNode
      dof_i = mesh.dofs[k, j, elnum]
      b[dof_i] += workvec[pos]
      pos += 1
    end
  end

  return nothing
end

#------------------------------------------------------------------------------
# New AssembleElementData type for getting the diagonal only

mutable struct AssembleDiagJacData{T} <: AssembleElementData
  A::DiagJac{T}
end

function AssembleDiagJacData(mesh, sbp, eqn, opts, jac::DiagJac{T}) where T

  return AssembleDiagJacData{T}(jac)
end

function NullAssembleDiagJacData(::Type{T}) where {T}
  return AssembleDiagJacData{T}(DiagJac(T, 0, 0))
end

"""
  Assembles the volume integral contribution.
"""
function assembleElement(helper::AssembleDiagJacData, mesh::AbstractMesh,
                         elnum::Integer, jac::AbstractArray{T, 4}) where T

  # depending on how dofs were assigned, this might not result in helper.A
  # having the same layout as the real Jacobian.  Need to make sure the
  # mat-vec function selects the correct indices to multiply against

  for q=1:mesh.numNodesPerElement
    for p=1:mesh.numNodesPerElement

      # get values
      @simd for j=1:mesh.numDofPerNode
        @simd for i=1:mesh.numDofPerNode
          i1 = i + (p-1)*mesh.numDofPerNode
          j1 = j + (q-1)*mesh.numDofPerNode
          helper.A.A[i1, j1, elnum] += real(jac[i, j, p, q])
        end
      end

    end  # end loop p
  end  # end loop q

  return nothing
end

"""
  Assembles contributions into an element-block matrix.  `jacLR` and `jacRL`
  are not used.
"""
function assembleInterface(helper::AssembleDiagJacData,
                           sbpface::DenseFace,
                           mesh::AbstractMesh, iface::Interface,
                           jacLL::AbstractArray{T, 4},
                           jacLR::AbstractArray{T, 4},
                           jacRL::AbstractArray{T, 4},
                           jacRR::AbstractArray{T, 4}) where T

  numNodesPerElement = size(jacLL, 4)
  numDofPerNode = size(jacLL, 1)

  permL = sview(sbpface.perm, :, iface.faceL)
  permR = sview(sbpface.perm, :, iface.faceR)

  for q=1:sbpface.stencilsize
    qL = sbpface.perm[q, iface.faceL]
    qR = sbpface.perm[q, iface.faceR]

    for p=1:sbpface.stencilsize
      pL = sbpface.perm[p, iface.faceL]
      pR = sbpface.perm[p, iface.faceR]

      # put values into 2 x 2 block matrix
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          i1 = i + (pL-1)*mesh.numDofPerNode
          j1 = j + (qL-1)*mesh.numDofPerNode
          i2 = i + (pR-1)*mesh.numDofPerNode
          j2 = j + (qR-1)*mesh.numDofPerNode
          helper.A.A[i1, j1, iface.elementL] +=  jacLL[i, j, pL, qL]
          helper.A.A[i2, j2, iface.elementR] +=  jacRR[i, j, pR, qR]
        end
      end

    end # end loop p
  end  # end loop q

  return nothing
end

function assembleInterface(helper::AssembleDiagJacData,
                           sbpface::SparseFace,
                           mesh::AbstractMesh, iface::Interface,
                           jacLL::AbstractArray{T, 4},
                           jacLR::AbstractArray{T, 4},
                           jacRL::AbstractArray{T, 4},
                           jacRR::AbstractArray{T, 4}) where T

  numNodesPerElement = size(jacLL, 4)
  numDofPerNode = size(jacLL, 1)

  for p=1:sbpface.numnodes
    # row indices
    iR = sbpface.nbrperm[p, iface.orient]
    pL = sbpface.perm[p, iface.faceL]
    pR = sbpface.perm[iR, iface.faceR]

    # column indices
    # this matches the way SBP puts the data into the jac arrays,
    # but its a little weird
    qL = pL
    qR = pR

    # put values into 2 x 2 block matrix
    @simd for j=1:numDofPerNode
      @simd for i=1:numDofPerNode
        i1 = i + (pL-1)*mesh.numDofPerNode
        j1 = j + (qL-1)*mesh.numDofPerNode
        i2 = i + (pR-1)*mesh.numDofPerNode
        j2 = j + (qR-1)*mesh.numDofPerNode
        helper.A.A[i1, j1, iface.elementL] +=  jacLL[i, j, pL, qL]
        helper.A.A[i2, j2, iface.elementR] +=  jacRR[i, j, pR, qR]
      end
    end

  end  # end loop p

  return nothing
end

"""
  Assembles contribution of the shared face terms into the element-block
  matrix.  `jacLR` is not used.
"""
function assembleSharedFace(helper::AssembleDiagJacData, sbpface::DenseFace,
                            mesh::AbstractMesh,
                            iface::Interface,
                            jacLL::AbstractArray{T, 4},
                            jacLR::AbstractArray{T, 4}) where T

#    for p in permR # =1:numNodesPerElement
  for q=1:sbpface.stencilsize
    qL = sbpface.perm[q, iface.faceL]

    for p=1:sbpface.stencilsize
      pL = sbpface.perm[p, iface.faceL]

      # put values into 2 x 2 block matrix
      for j=1:mesh.numDofPerNode
        for i=1:mesh.numDofPerNode
          i1 = i + (pL-1)*mesh.numDofPerNode
          j1 = j + (qL-1)*mesh.numDofPerNode

          helper.A.A[i1, j1, iface.elementL] += real(jacLL[i, j, pL, qL])
        end
      end

    end  # end loop p
  end

  return nothing
end

function assembleSharedFace(helper::AssembleDiagJacData, sbpface::SparseFace,
                            mesh::AbstractMesh,
                            iface::Interface,
                            jacLL::AbstractArray{T, 4},
                            jacLR::AbstractArray{T, 4}) where T

  numNodesPerElement = size(jacLL, 4)
  numDofPerNode = size(jacLL, 1)

  for p=1:sbpface.numnodes
    iR = sbpface.nbrperm[p, iface.orient]
    qL = sbpface.perm[p, iface.faceL]
    qR = sbpface.perm[iR, iface.faceR]

    pL = qL

    for j=1:mesh.numDofPerNode
      for i=1:mesh.numDofPerNode
        i1 = i + (pL-1)*mesh.numDofPerNode
        j1 = j + (qL-1)*mesh.numDofPerNode

        helper.A.A[i1, j1, iface.elementL] += real(jacLL[i, j, pL, qL])
      end
    end

  end  # end loop p

   return nothing
end

"""
  Assembles the boundary terms into an element-block matrix.
"""
function assembleBoundary(helper::AssembleDiagJacData, sbpface::DenseFace,
                            mesh::AbstractMesh,
                            bndry::Boundary,
                            jac::AbstractArray{T, 4}) where T
  for q=1:sbpface.stencilsize
    qL = sbpface.perm[q, bndry.face]

    for p=1:sbpface.stencilsize
      pL = sbpface.perm[p, bndry.face]

      # get dofs for node p
      #dof1 = mesh.dofs[1, pL, bndry.element]
      #blockidx = div(dof1 - 1, mesh.numDofPerNode) + 1

      # get values
      for j=1:mesh.numDofPerNode
        for i=1:mesh.numDofPerNode
          i1 = i + (pL-1)*mesh.numDofPerNode
          j1 = j + (qL-1)*mesh.numDofPerNode

          helper.A.A[i1, j1, bndry.element] += jac[i, j, pL, qL]
        end
      end

    end  # end loop p
  end  # end loop q

  return nothing
end

function assembleBoundary(helper::AssembleDiagJacData, sbpface::SparseFace,
                            mesh::AbstractMesh,
                            bndry::Boundary,
                            jac::AbstractArray{T, 4}) where T

  elnum = bndry.element
  numNodesPerElement = size(jac, 4)
  numDofPerNode = size(jac, 1)

  for p=1:sbpface.numnodes
    pL = sbpface.perm[p, bndry.face]
    qL = pL

    # get dofs for node p
    #dof1 = mesh.dofs[1, pL, bndry.element]
    #blockidx = div(dof1 - 1, mesh.numDofPerNode) + 1

    # get values
    for j=1:mesh.numDofPerNode
      for i=1:mesh.numDofPerNode
        i1 = i + (pL-1)*mesh.numDofPerNode
        j1 = j + (qL-1)*mesh.numDofPerNode

        helper.A.A[i1, j1, bndry.element] += jac[i, j, pL, qL]
      end
    end

  end  # end loop p

  return nothing
end

#------------------------------------------------------------------------------
# apply the stabilization to the diagonal Jacobian

"""
  This function takes the block-diagonal Jacobian `A` and a solution
  vector `q_vec` and removes the unstable modes using
  [`findStablePerturbation!`](@ref).

  **Inputs**

   * mesh: the mesh that was used in the computation of `A`
   * q_vec: solution vector (the local part)

  **Inputs/Outputs**

   * A: a DiagJac containing the block diagonal Jacobian.  On exit, it will
        have the unstable modes removed.
"""
function filterDiagJac(mesh::AbstractDGMesh, q_vec::AbstractVector{T2}, A::DiagJac{T}) where {T, T2}

  blocksize, blocksize, nblock = size(A.A)
  T3 = promote_type(T, T2)

  @assert blocksize == mesh.numDofPerNode*mesh.numNodesPerElement
  nwork = div(blocksize*(blocksize+1),2)
  workvec = zeros(T3, nwork)  # work array for inner function
  Ablock = zeros(T, blocksize, blocksize)  # copy array into this to negate it
  Ablock2 = zeros(T, blocksize, blocksize)
  ublock = zeros(T2, blocksize)
  for k=1:nblock
    # because the inner function assumes the residual is defined as
    # du/dt + R(u) = 0, but Ticon writes it as du/dt = R(u), we have to
    # multiply the Jacobian by -1 to make the unstable modes the negative
    # eigenvalues
    for j=1:blocksize
      for i=1:blocksize
        Ablock[i, j] = -A.A[i, j, k]
      end
    end

    # get the entries of q_vec
    idx = 1
    for j=1:mesh.numNodesPerElement
      for i=1:mesh.numDofPerNode
        ublock[idx] = q_vec[mesh.dofs[i, j, k]]
        idx += 1
      end
    end

    findStablePerturbation!(Ablock, ublock, workvec, Ablock2)
#    removeUnstableModes!(Ablock, u_k)

    # copy back
    for j=1:blocksize
      for i=1:blocksize
        A.A[i, j, k] = Ablock2[i, j]
      end
    end

  end  # end loop k

  return nothing
end

"""
  modifies matrix `Jac` such that `u.'*sym(Jac)*u` is strictly positive

  **Inputs**

   * `u`: an given vector whose dimensions are consistent with `Jac`
   * `A`: a work vector needed by this function (overwritten).  The
          element type should be the "maximum" type of the element types
          of `u` and `Jac`

  **Inputs/Outputs**

   * `Jac`: matrix being modified
"""
function removeUnstableModes!(Jac::AbstractMatrix,
                              u::AbstractVector,
                              A::AbstractVector{T}) where T
  @assert( size(Jac,1) == size(Jac,2) == length(u) )
  n = size(Jac,1)
  # compute baseline product, 0.5*u.'*(Jac^T + Jac)*u
  prod = zero(T)
  for i = 1:n
    for j = 1:n
      prod += 0.5*(Jac[i,j] + Jac[j,i])*u[i]*u[j]
    end
  end
  if prod > 0
    # nothing to do
    return
  end
  # array A stores the entries in the contraint Jacobian
#  A = zeros(div(n*(n+1),2))
  for i = 1:n
    A[div(i*(i-1),2)+i] = u[i]*u[i]
    for j = 1:(i-1)
      A[div(i*(i-1),2)+j] = 2.0*u[i]*u[j]
    end
  end
  scale!(A, -prod/dot(A,A))
#  A *= -prod/dot(A,A)
  for i = 1:n
    Jac[i,i] += A[div(i*(i-1),2)+i]
    for j = 1:(i-1)
      Jac[i,j] += A[div(i*(i-1),2)+j]
      Jac[j,i] += A[div(i*(i-1),2)+j]
    end
  end

  return nothing
end

"""
  returns matrix `Jacpert` such that `u.'*sym(Jac + Jacpert)*u` is strictly
  positive

  **Inputs**

   * `u`: an given vector whose dimensions are consistent with `Jac`
   * `Jac`: matrix that needs to be perturbed
   * `A`: a work vector needed by this function (overwritten).  The
          element type should be the "maximum" type of the element types
          of `u` and `Jac`

  **Inputs/Outputs**

   * `Jacpert`: matrix perturbation
"""
function findStablePerturbation!(Jac::AbstractMatrix,
                                 u::AbstractVector,
                                 A::AbstractVector{T},
                                 Jacpert::AbstractMatrix) where T
  @assert( size(Jac,1) == size(Jac,2) == size(Jacpert,1) == size(Jacpert,2)
           == length(u) )
  n = size(Jac,1)
  # compute baseline product, 0.5*u.'*(Jac^T + Jac)*u
  prod = zero(T)
  for i = 1:n
    for j = 1:n
      prod += 0.5*(Jac[i,j] + Jac[j,i])*u[i]*u[j]
    end
  end
  if prod > 0
    # nothing to do
    return
  end
  # array A stores the entries in the contraint Jacobian
#  A = zeros(div(n*(n+1),2))
  for i = 1:n
    A[div(i*(i-1),2)+i] = u[i]*u[i]
    for j = 1:(i-1)
      A[div(i*(i-1),2)+j] = 2.0*u[i]*u[j]
    end
  end
#  A *= -prod/dot(A,A)
  scale!(A, -prod/dot(A, A))
  fill!(Jacpert, 0.0)
  for i = 1:n
    Jacpert[i,i] += A[div(i*(i-1),2)+i]
    for j = 1:(i-1)
      Jacpert[i,j] += A[div(i*(i-1),2)+j]
      Jacpert[j,i] += A[div(i*(i-1),2)+j]
    end
  end
end


text:
using jInv.Utils
using Base.Test
using OED
using jInv.ForwardShare
using jInv.Mesh

n = 100
nrhs = 2
# A = abs(sprandn(N,n,.2))
A = sdiag(1./linspace(1,n,n).^4)
A = A[1:end-2,:]
N = size(A,1)
x = rand(n,nrhs); x[4:7,:] = 1;
D = A*x  + 0* randn(N,nrhs)/1e3
OEDparamA = OEDTypeAParam(A,D,zeros(0,0))

# setup OED type B param
r = 20;
l = 8;
Mimg  = getRegularMesh([-1 1 -1 1.],[10,10])
Mgrid = getRegularMesh([-1.2 2.4 -1.6 1.4],[20,20])
OEDparamB = getOEDTypeBParam(r,l,x,Mimg,Mgrid,noiseLevel=0.)

OEDs = (OEDparamA, OEDparamB)

for k=1:2
	display(OEDs[k])

	L = 1*speye(n)
	Ce = ones(1,n)
	ce = Ce*x
	pFor = EqconstrainedOEDParam(OEDs[k], L,Ce,ce, zeros(0,0), [])

	println("--- solve LS problem ---")
	display(pFor)
	np  =  getNumberOfDesignParameters(OEDs[k])
	p0      = rand(np)
	xr, = getData(p0,pFor)
	@test norm(Ce*xr-ce) < 1e-12
	println("rel.error: $(norm(xr-x)/norm(x))")

	println("--- test derivative of LS solution w.r.t. design ---")
	pass, = checkDerivative(rand(np),pFor)
	@test pass

	println("--- adjoint test ---")
	nw,nv = getSensMatSize(pFor)
	v = randn(nv)
	w = randn(nw)
	t1 = dot(w,getSensMatVec(v,p0,pFor))
	t2 = dot(v,getSensTMatVec(w,p0,pFor))
	@test abs(t1-t2)/abs(t1) < 1e-10
end


text:
# Note: This is an experimental feature and may be changed in future releases without notice.
mutable struct ElementContainerP4est{NDIMS, RealT<:Real, uEltype<:Real, NDIMSP1, NDIMSP2, NDIMSP3} <: AbstractContainer
  # Physical coordinates at each node
  node_coordinates      ::Array{RealT, NDIMSP2}   # [orientation, node_i, node_j, node_k, element]
  # Jacobian matrix of the transformation
  # [jacobian_i, jacobian_j, node_i, node_j, node_k, element] where jacobian_i is the first index of the Jacobian matrix,...
  jacobian_matrix       ::Array{RealT, NDIMSP3}
  # Contravariant vectors, scaled by J, in Kopriva's blue book called Ja^i_n (i index, n dimension)
  contravariant_vectors ::Array{RealT, NDIMSP3}   # [dimension, index, node_i, node_j, node_k, element]
  # 1/J where J is the Jacobian determinant (determinant of Jacobian matrix)
  inverse_jacobian      ::Array{RealT, NDIMSP1}   # [node_i, node_j, node_k, element]
  # Buffer for calculated surface flux
  surface_flux_values   ::Array{uEltype, NDIMSP2} # [variable, i, j, direction, element]

  # internal `resize!`able storage
  _node_coordinates     ::Vector{RealT}
  _jacobian_matrix      ::Vector{RealT}
  _contravariant_vectors::Vector{RealT}
  _inverse_jacobian     ::Vector{RealT}
  _surface_flux_values  ::Vector{uEltype}
end

@inline nelements(elements::ElementContainerP4est) = size(elements.node_coordinates, ndims(elements) + 2)
@inline Base.ndims(::ElementContainerP4est{NDIMS}) where NDIMS = NDIMS
@inline Base.eltype(::ElementContainerP4est{NDIMS, RealT, uEltype}) where {NDIMS, RealT, uEltype} = uEltype

# Only one-dimensional `Array`s are `resize!`able in Julia.
# Hence, we use `Vector`s as internal storage and `resize!`
# them whenever needed. Then, we reuse the same memory by
# `unsafe_wrap`ping multi-dimensional `Array`s around the
# internal storage.
function Base.resize!(elements::ElementContainerP4est, capacity)
  @unpack _node_coordinates, _jacobian_matrix, _contravariant_vectors,
    _inverse_jacobian, _surface_flux_values = elements

  n_dims = ndims(elements)
  n_nodes = size(elements.node_coordinates, 2)
  n_variables = size(elements.surface_flux_values, 1)

  resize!(_node_coordinates, n_dims * n_nodes^n_dims * capacity)
  elements.node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
    (n_dims, ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_jacobian_matrix, n_dims^2 * n_nodes^n_dims * capacity)
  elements.jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix),
    (n_dims, n_dims, ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_contravariant_vectors, length(_jacobian_matrix))
  elements.contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors),
    size(elements.jacobian_matrix))

  resize!(_inverse_jacobian, n_nodes^n_dims * capacity)
  elements.inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian),
    (ntuple(_ -> n_nodes, n_dims)..., capacity))

  resize!(_surface_flux_values,
    n_variables * n_nodes^(n_dims-1) * (n_dims*2) * capacity)
  elements.surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
    (n_variables, ntuple(_ -> n_nodes, n_dims-1)..., n_dims*2, capacity))

  return nothing
end

# Create element container and initialize element data
function init_elements(mesh::P4estMesh{NDIMS, RealT}, equations,
                       basis, ::Type{uEltype}) where {NDIMS, RealT<:Real, uEltype<:Real}
  nelements = ncells(mesh)

  _node_coordinates = Vector{RealT}(undef, NDIMS * nnodes(basis)^NDIMS * nelements)
  node_coordinates = unsafe_wrap(Array, pointer(_node_coordinates),
    (NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _jacobian_matrix = Vector{RealT}(undef, NDIMS^2 * nnodes(basis)^NDIMS * nelements)
  jacobian_matrix = unsafe_wrap(Array, pointer(_jacobian_matrix),
    (NDIMS, NDIMS, ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _contravariant_vectors = similar(_jacobian_matrix)
  contravariant_vectors = unsafe_wrap(Array, pointer(_contravariant_vectors),
    size(jacobian_matrix))

  _inverse_jacobian = Vector{RealT}(undef, nnodes(basis)^NDIMS * nelements)
  inverse_jacobian = unsafe_wrap(Array, pointer(_inverse_jacobian),
    (ntuple(_ -> nnodes(basis), NDIMS)..., nelements))

  _surface_flux_values = Vector{uEltype}(undef,
    nvariables(equations) * nnodes(basis)^(NDIMS-1) * (NDIMS*2) * nelements)
  surface_flux_values = unsafe_wrap(Array, pointer(_surface_flux_values),
    (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., NDIMS*2, nelements))

  elements = ElementContainerP4est{NDIMS, RealT, uEltype, NDIMS+1, NDIMS+2, NDIMS+3}(
    node_coordinates, jacobian_matrix, contravariant_vectors,
    inverse_jacobian, surface_flux_values,
    _node_coordinates, _jacobian_matrix, _contravariant_vectors,
    _inverse_jacobian, _surface_flux_values)

  init_elements!(elements, mesh, basis)
  return elements
end

mutable struct InterfaceContainerP4est{NDIMS, uEltype<:Real, NDIMSP2} <: AbstractContainer
  u             ::Array{uEltype, NDIMSP2}       # [primary/secondary, variable, i, j, interface]
  element_ids   ::Matrix{Int}                   # [primary/secondary, interface]
  node_indices  ::Matrix{NTuple{NDIMS, Symbol}} # [primary/secondary, interface]

  # internal `resize!`able storage
  _u            ::Vector{uEltype}
  _element_ids  ::Vector{Int}
  _node_indices ::Vector{NTuple{NDIMS, Symbol}}
end

@inline ninterfaces(interfaces::InterfaceContainerP4est) = size(interfaces.element_ids, 2)
@inline Base.ndims(::InterfaceContainerP4est{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(interfaces::InterfaceContainerP4est, capacity)
  @unpack _u, _element_ids, _node_indices = interfaces

  n_dims = ndims(interfaces)
  n_nodes = size(interfaces.u, 3)
  n_variables = size(interfaces.u, 2)

  resize!(_u, 2 * n_variables * n_nodes^(n_dims-1) * capacity)
  interfaces.u = unsafe_wrap(Array, pointer(_u),
    (2, n_variables, ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(_element_ids, 2 * capacity)
  interfaces.element_ids = unsafe_wrap(Array, pointer(_element_ids), (2, capacity))

  resize!(_node_indices, 2 * capacity)
  interfaces.node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, capacity))

  return nothing
end

# Create interface container and initialize interface data.
function init_interfaces(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_interfaces = count_required_surfaces(mesh).interfaces

  _u = Vector{uEltype}(undef, 2 * nvariables(equations) * nnodes(basis)^(NDIMS-1) * n_interfaces)
  u = unsafe_wrap(Array, pointer(_u),
    (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_interfaces))

  _element_ids = Vector{Int}(undef, 2 * n_interfaces)
  element_ids = unsafe_wrap(Array, pointer(_element_ids), (2, n_interfaces))

  _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces)
  node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces))

  interfaces = InterfaceContainerP4est{NDIMS, uEltype, NDIMS+2}(u, element_ids, node_indices,
                                                                _u, _element_ids, _node_indices)

  init_interfaces!(interfaces, mesh)

  return interfaces
end

mutable struct BoundaryContainerP4est{NDIMS, uEltype<:Real, NDIMSP1} <: AbstractContainer
  u           ::Array{uEltype, NDIMSP1}       # [variables, i, j, boundary]
  element_ids ::Vector{Int}                   # [boundary]
  node_indices::Vector{NTuple{NDIMS, Symbol}} # [boundary]
  name        ::Vector{Symbol}                # [boundary]

  # internal `resize!`able storage
  _u          ::Vector{uEltype}
end

@inline nboundaries(boundaries::BoundaryContainerP4est) = length(boundaries.element_ids)
@inline Base.ndims(::BoundaryContainerP4est{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(boundaries::BoundaryContainerP4est, capacity)
  @unpack _u, element_ids, node_indices, name = boundaries

  n_dims = ndims(boundaries)
  n_nodes = size(boundaries.u, 2)
  n_variables = size(boundaries.u, 1)

  resize!(_u, n_variables * n_nodes^(n_dims-1) * capacity)
  boundaries.u = unsafe_wrap(Array, pointer(_u),
    (n_variables, ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(element_ids, capacity)

  resize!(node_indices, capacity)

  resize!(name, capacity)

  return nothing
end

# Create interface container and initialize interface data in `elements`.
function init_boundaries(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_boundaries = count_required_surfaces(mesh).boundaries

  _u = Vector{uEltype}(undef, nvariables(equations) * nnodes(basis)^(NDIMS-1) * n_boundaries)
  u = unsafe_wrap(Array, pointer(_u),
    (nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_boundaries))

  element_ids  = Vector{Int}(undef, n_boundaries)
  node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, n_boundaries)
  names        = Vector{Symbol}(undef, n_boundaries)

  boundaries = BoundaryContainerP4est{NDIMS, uEltype, NDIMS+1}(u, element_ids,
                                                               node_indices, names, _u)

  init_boundaries!(boundaries, mesh)

  return boundaries
end

# Container data structure (structure-of-arrays style) for DG L2 mortars
#
# The positions used in `element_ids` are 1:3 (in 2D) or 1:5 (in 3D), where 1:2 (in 2D)
# or 1:4 (in 3D) are the small elements numbered in z-order and 3 or 5 is the large element.
# The solution values on the mortar element are saved in `u`, where `position` is the number
# of the small element that corresponds to the respective part of the mortar element.
# The first dimension `small/large side` takes 1 for small side and 2 for large side.
#
# Illustration of the positions in `element_ids` in 3D, where ξ and η are the local coordinates
# of the mortar element, which are precisely the local coordinates that span
# the surface of the smaller side.
# Note that the orientation in the physical space is completely irrelevant here.
#   /---------------------------\  /----------------------------\
#   |             |             |  |                            |
#   |    small    |    small    |  |                            |
#   |      3      |      4      |  |                            |
#   |             |             |  |           large            |
#   |-------------|-------------|  |             5              |
# η |             |             |  |                            |
#   |    small    |    small    |  |                            |
# ^ |      1      |      2      |  |                            |
# | |             |             |  |                            |
# | \---------------------------/  \----------------------------/
# |
# ⋅----> ξ
mutable struct MortarContainerP4est{NDIMS, uEltype<:Real, NDIMSP1, NDIMSP3} <: AbstractContainer
  u             ::Array{uEltype, NDIMSP3} # [small/large side, variable, position, i, j, mortar]
  element_ids   ::Matrix{Int}             # [position, mortar]
  node_indices  ::Matrix{NTuple{NDIMS, Symbol}} # [small/large, mortar]

  # internal `resize!`able storage
  _u            ::Vector{uEltype}
  _element_ids  ::Vector{Int}
  _node_indices ::Vector{NTuple{NDIMS, Symbol}}
end

@inline nmortars(mortars::MortarContainerP4est) = size(mortars.element_ids, 2)
@inline Base.ndims(::MortarContainerP4est{NDIMS}) where NDIMS = NDIMS

# See explanation of Base.resize! for the element container
function Base.resize!(mortars::MortarContainerP4est, capacity)
  @unpack _u, _element_ids, _node_indices = mortars

  n_dims = ndims(mortars)
  n_nodes = size(mortars.u, 4)
  n_variables = size(mortars.u, 2)

  resize!(_u, 2 * n_variables * 2^(n_dims-1) * n_nodes^(n_dims-1) * capacity)
  mortars.u = unsafe_wrap(Array, pointer(_u),
    (2, n_variables, 2^(n_dims-1), ntuple(_ -> n_nodes, n_dims-1)..., capacity))

  resize!(_element_ids, (2^(n_dims-1) + 1) * capacity)
  mortars.element_ids = unsafe_wrap(Array, pointer(_element_ids),
    (2^(n_dims-1) + 1, capacity))

  resize!(_node_indices, 2 * capacity)
  mortars.node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, capacity))

  return nothing
end

# Create mortar container and initialize mortar data.
function init_mortars(mesh::P4estMesh, equations, basis, elements)
  NDIMS = ndims(elements)
  uEltype = eltype(elements)

  # Initialize container
  n_mortars = count_required_surfaces(mesh).mortars

  _u = Vector{uEltype}(undef,
    2 * nvariables(equations) * 2^(NDIMS-1) * nnodes(basis)^(NDIMS-1) * n_mortars)
  u = unsafe_wrap(Array, pointer(_u),
    (2, nvariables(equations), 2^(NDIMS-1), ntuple(_ -> nnodes(basis), NDIMS-1)..., n_mortars))

  _element_ids = Vector{Int}(undef, (2^(NDIMS-1) + 1) * n_mortars)
  element_ids = unsafe_wrap(Array, pointer(_element_ids), (2^(NDIMS-1) + 1, n_mortars))

  _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_mortars)
  node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_mortars))

  mortars = MortarContainerP4est{NDIMS, uEltype, NDIMS+1, NDIMS+3}(u, element_ids, node_indices,
                                                                   _u, _element_ids, _node_indices)

  init_mortars!(mortars, mesh)

  return mortars
end

function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache)
  # Re-initialize elements container
  @unpack elements = cache
  resize!(elements, ncells(mesh))
  init_elements!(elements, mesh, dg.basis)

  required = count_required_surfaces(mesh)

  # resize interfaces container
  @unpack interfaces = cache
  resize!(interfaces, required.interfaces)

  # resize boundaries container
  @unpack boundaries = cache
  resize!(boundaries, required.boundaries)

  # resize mortars container
  @unpack mortars = cache
  resize!(mortars, required.mortars)

  # re-initialize containers together to reduce
  # the number of iterations over the mesh in p4est
  init_surfaces!(interfaces, mortars, boundaries, mesh)
end

# Iterate over all interfaces and count
# - (inner) interfaces
# - mortars
# - boundaries
# and collect the numbers in `user_data` in this order.
function count_surfaces_iter_face(info, user_data)
  if info.sides.elem_count == 2
    # Two neighboring elements => Interface or mortar

    # Extract surface data
    sides = (unsafe_load_sc(p4est_iter_face_side_t, info.sides, 1),
             unsafe_load_sc(p4est_iter_face_side_t, info.sides, 2))

    if sides[1].is_hanging == 0 && sides[2].is_hanging == 0
      # No hanging nodes => normal interface
      # Unpack user_data = [interface_count] and increment interface_count
      ptr = Ptr{Int}(user_data)
      id = unsafe_load(ptr, 1)
      unsafe_store!(ptr, id + 1, 1)
    else
      # Hanging nodes => mortar
      # Unpack user_data = [mortar_count] and increment mortar_count
      ptr = Ptr{Int}(user_data)
      id = unsafe_load(ptr, 2)
      unsafe_store!(ptr, id + 1, 2)
    end
  elseif info.sides.elem_count == 1
    # One neighboring elements => boundary

    # Unpack user_data = [boundary_count] and increment boundary_count
    ptr = Ptr{Int}(user_data)
    id = unsafe_load(ptr, 3)
    unsafe_store!(ptr, id + 1, 3)
  end

  return nothing
end

function count_required_surfaces(mesh::P4estMesh)
  # Let p4est iterate over all interfaces and call count_surfaces_iter_face
  iter_face_c = @cfunction(count_surfaces_iter_face, Cvoid, (Ptr{p4est_iter_face_info_t}, Ptr{Cvoid}))

  # interfaces, mortars, boundaries
  user_data = [0, 0, 0]

  iterate_faces(mesh, iter_face_c, user_data)

  # Return counters
  return (interfaces = user_data[1],
          mortars    = user_data[2],
          boundaries = user_data[3])
end

# Let p4est iterate over all interfaces and execute the C function iter_face_c
function iterate_faces(mesh::P4estMesh, iter_face_c, user_data)
  if user_data isa AbstractArray
    user_data_ptr = pointer(user_data)
  else
    user_data_ptr = pointer_from_objref(user_data)
  end
  GC.@preserve user_data begin
    p4est_iterate(mesh.p4est,
                  C_NULL, # ghost layer
                  user_data_ptr,
                  C_NULL, # iter_volume
                  iter_face_c, # iter_face
                  C_NULL) # iter_corner
  end

  return nothing
end

# Convert sc_array of type T to Julia array
function unsafe_wrap_sc(::Type{T}, sc_array) where T
  element_count = sc_array.elem_count
  element_size = sc_array.elem_size

  @assert element_size == sizeof(T)

  return [unsafe_wrap(T, sc_array.array + element_size * i) for i in 0:element_count-1]
end

# Load the ith element (1-indexed) of an sc array of type T
function unsafe_load_sc(::Type{T}, sc_array, i=1) where T
  element_size = sc_array.elem_size

  @assert element_size == sizeof(T)

  return unsafe_wrap(T, sc_array.array + element_size * (i - 1))
end

# Convert Tuple node_indices to actual indices.
# E.g., (:one, :i, :j) will be (1, i, j) for some i and j,
# (:i, :end) will be (i, size[2]),
# (:one, :i_backwards) will be (1, size[1] - i + 1).
function evaluate_index(indices, size, dim, i, j=0)
  if indices[dim] === :i
    return i
  elseif indices[dim] === :i_backwards
    return size[dim] - i + 1
  elseif indices[dim] === :j
    return j
  elseif indices[dim] === :j_backwards
    return size[dim] - j + 1
  elseif indices[dim] === :one
    return 1
  elseif indices[dim] === :end
    return size[dim]
  end

  error("Invalid identifier: Only :one, :end, :i, :j, :i_backwards, :j_backwards are valid index identifiers")
end

# Remove :one and :end to index surface_flux_values properly (2D version).
#
# Suppose some element is indexed relative to some interface as `(:end, :i_backwards)` (in 2D).
# This interface will be at the right face of the element, and it will be reversely oriented.
# When copying data to the interface, the value `u[v, end, end - i + 1, element]` will be copied
# to `u[2, v, i, interface]` in the interface container (assuming this is a secondary element).
# Now, the calculated flux at interface node i needs to be copied back to
# `surface_flux_values[v, end - i + 1, 2, element]`.
# This is the same index as in `evaluate_index` but without the `:one` or `:end`.
#
# Dispatch by dimension to ensure type stability
@inline function evaluate_index_surface(indices::NTuple{2, Symbol}, size, dim, i)
  if indices[1] in (:one, :end)
    indices_surface = indices[2:2]
  else # indices[2] in (:one, :end)
    indices_surface = indices[1:1]
  end

  return evaluate_index(indices_surface, size, dim, i)
end

# 3D version
@inline function evaluate_index_surface(indices::NTuple{3, Symbol}, size, dim, i, j=0)
  if indices[1] in (:one, :end)
    indices_surface = indices[2:3]
  elseif indices[2] in (:one, :end)
    indices_surface = (indices[1], indices[3])
  else # indices[3] in (:one, :end)
    indices_surface = indices[1:2]
  end

  return evaluate_index(indices_surface, size, dim, i)
end

# Return direction of the face, which is indexed by node_indices
@inline function indices2direction(indices)
  if indices[1] in (:one, :end)
    orientation = 1
  elseif indices[2] in (:one, :end)
    orientation = 2
  else # indices[3] in (:one, :end)
    orientation = 3
  end
  negative_direction = orientation * 2 - 1

  if indices[orientation] === :one
    return negative_direction
  else # indices[orientation] === :end
    return negative_direction + 1
  end
end

include("containers_2d.jl")


text:
@testset "XTbML" begin
    @testset "dict parse" begin
        g = MortalityTables.get_and_parse
        d = Dict(:a => "1.2", :c => "a")

        @test g(d, :a) == 1.2
        @test ismissing(g(d, :b))
        @test ismissing(g(d, :b))
        @test_throws Exception g(d, :c)

    end

    @testset "XTbML loading" begin
        pth = joinpath(soa_tbl_dir,"t1076.xml")
        file = MortalityTables.open_and_read(pth) |> MortalityTables.getXML
        xtbl = MortalityTables.parseXTbMLTable(file, pth)
        @test isa(xtbl, MortalityTables.XTbMLTable)
    end

    @testset "Ultimate Only" begin
        pth = joinpath(soa_tbl_dir,"t17.xml")
        file = MortalityTables.open_and_read(pth) |> MortalityTables.getXML
        xtbl = MortalityTables.parseXTbMLTable(file, pth)

        mt = MortalityTables.XTbML_Table_To_MortalityTable(xtbl)
        @test isa(mt, MortalityTable)

        @test mt.ultimate[0] ≈ 0.00245
        @test mt.ultimate[100] ≈ 1.0
        @test_throws BoundsError mt.ultimate[101]
    end

    @testset "XTbML to MortalityTable" begin
        @testset "Select and Ultimate" begin
            pth = joinpath(soa_tbl_dir,"t1076.xml")
            file = MortalityTables.open_and_read(pth) |> MortalityTables.getXML
            xtbl = MortalityTables.parseXTbMLTable(file, pth)

            mt = MortalityTables.XTbML_Table_To_MortalityTable(xtbl)
            @test isa(mt, MortalityTable)

            @test mt.select[35][35] ≈ 0.00037
            @test mt.ultimate[16] ≈ 0.00041
            @test mt.select[35][59] ≈ 0.00508
            @test mt.select[35][60] ≈ 0.00621
        end

        @testset "Ultimate Only, not begin at age 0" begin
            pth = joinpath(soa_tbl_dir,"t18.xml")
            file = MortalityTables.open_and_read(pth) |> MortalityTables.getXML
            xtbl = MortalityTables.parseXTbMLTable(file, pth)
            mt = MortalityTables.XTbML_Table_To_MortalityTable(xtbl)
            @test isa(mt, MortalityTable)
        end
    end
end


text:
"_sample_bar_points(samples, maxwidth; first_sample_no = 1)"
function _sample_bar_points(samples, maxwidth; first_sample_no = 1)
    pts = Vector{Point}()
    for (i, val) in enumerate(samples)
        x = maxwidth * (i + first_sample_no - 1) / length(samples)
        y = - val
        push!(pts,  Point(x, 0*x))
        push!(pts,  Point(x, 0*x) + Point(0y, y))
    end
    pts
end

"""
    _draw_bars(origo, samples, maxwidth;
                   rotatehue_degrees_total = 0°, first_sample_no = 1,
                   top_circle_radius = 0.04EM)))
"""
function _draw_bars(origo, samples, maxwidth;
                   rotatehue_degrees_total = 0°, first_sample_no = 1, top_circle_radius = 0.04EM)
    startcolo = get_current_RGB()
    samplepoints = _sample_bar_points(samples,  maxwidth; first_sample_no)
    function foovertex(n, pt)
        rotatedeg = (n - 1) * rotatehue_degrees_total / length(samplepoints)
        sethue(rotate_hue(startcolo, rotatedeg))
        circle(pt, top_circle_radius, :stroke)
    end
    @layer begin
        for i in range(1, length(samplepoints) - 1, step = 2)
            pt = origo + samplepoints[i]
            npt = origo + samplepoints[i + 1]
            line(pt, npt, :stroke)
            foovertex(i, npt)
        end
    end
end

"""
    draw_barplot(origo, samples, width;
                first_sample_no = 1,
                rotatehue_degrees_total = 270°, firstcolor = PALETTE[3],
                height = 2EM, top_circle_radius = 0.04EM)

    - origo is a Point.
    - samples is the vector to plot as bars
    - width is the horizontal distance over which to distribute bar

Draw a bar plot without decorations; arrows and the first sample in the same colour.

"""
function draw_barplot(origo, samples, width;
                      first_sample_no = 1,
                      rotatehue_degrees_total = 270°, firstcolor = PALETTE[3],
                      height = 2EM, top_circle_radius = 0.04EM)
    @layer begin
        sethue(firstcolor)
        arrow(origo, origo + (0height,  -height))
        wplus = width * (length(samples) + first_sample_no - 1) / length(samples)
        arrow(origo, origo + (wplus, 0wplus) + (0.5EM, 0.0))
        sethue(color_from_palette("red"))
        _draw_bars(origo, samples, width; rotatehue_degrees_total, first_sample_no, top_circle_radius)
    end
end

"""
    color_matrix(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
                 legend = missing, kwargs...)
    -> (color_matrix, legend)

Return a matrix of colors suitable for display with 'place_image'. Colors
are found by calling 'legend' with single values out of function f_xy.

'legend' can be input as a keyword argument. If not, a legend is generated
internally, as modified by keyword arguments relevant to the proposed legend type.
See: BinLegend or BinLegendVector.

The legend is also returned. Internally generated legends are of type 'AbstractColorLegend',
and can be displayed with 'draw_legend'.
"""
function color_matrix(f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
    legend = missing, kwargs...)

    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
    valmat = [f_xy(x, y) for y in ys, x in xs];
    _color_matrix(valmat, f_xy;  legend, kwargs...)
end
function _color_matrix(valmat::Matrix{T}, f_xy;  legend = missing, centered = true, kwargs...) where T<:Tuple
    magnitude_matrix = hypot.(valmat)
    if ismissing(legend)
        legend = propose_legend(magnitude_matrix; kwargs...)
    end
    streamlinepixels = falses(size(valmat))
    streamlines_matrix!(streamlinepixels, f_xy)
    complex_convolution_matrix = convolution_matrix(f_xy, streamlinepixels)
    curmat = lic_matrix_current(complex_convolution_matrix, 0, zeroval = -0.0)
    colmat = map(legend.(magnitude_matrix), curmat) do col, lu
        color_with_lumin(col, 50 + 50 * lu)
    end
    colmat, legend
end
function _color_matrix(valmat::Matrix, f_xy;  legend = missing, centered = true, kwargs...)
    @error "Not implemented"
end

"""
    propose_legend(scalmat::Matrix{T}, kwargs...) where T <: Number
    -> AbstractColorLegend

Return a proposed legend based on the unit and values of magmat.
Keyword arguments set here overrule proposed arguments. See:
BinLegend or BinLegendVector. Possible values include 'missing',
which will also overrule proposals from this function.
"""
function propose_legend(mat::Matrix{T}; kwargs...) where T <: Number
    lum = 50
    mi, ma = lenient_min_max(mat)
    maxlegend = ma  # kwargs overrule if present
    noofbins = if !haskey(kwargs, :binwidth) && !haskey(kwargs, :noofbins) && !haskey(kwargs, :binbounds)
        # None of three present
        7
    else
        # Sufficient arguments are present.
        # IF noofbins is actually contained in kwargs, it will be placed after in the call, and thus overrule:
        missing
    end
    if T <: Velocity
        name = :Speed # kwargs overrule if present
        colorscheme = ColorSchemes.isoluminant_cgo_80_c38_n256.colors .|> co -> color_with_lumin(co, lum)
        BinLegend(;maxlegend, noofbins, colorscheme, name, kwargs...)
    else
        @error "$T scalar quantity not implemented."
    end
end


text:
"""
    FDEsolver(F, tSpan, y0, β, ::Nothing, par...; h = 0.01, nc = 3, StopIt = "Standard", tol = 10e-10, itmax = 10)

Solves fractional differential equations with a predictor-corrector approach.

    FDEsolver(F, tSpan, y0, β, J, par...; h = 0.01, nc = 3, StopIt = "Standard", tol = 10e-10, itmax = 10)

Additionally takes the Jacobian of the system J to evaluate the solution.

# Arguments
- `F`: the right side of the system of differential equations. It must be expressed
   in the form of a function and return a vector function with the same number of
   entries of order of derivatives. This function can also include a vector of
   parameters: `par...` .
- `tSpan::Vector{<:Real}`: the time span along which computation is performed.
- `y0`: the initial values in the form of a row vector (`Vector{<:Real}`) for ``β <= 1``
   and a matrix (`Matrix{<:Real}`) for ``β > 1``, where each column corresponds to the
   initial values of one differential equation and each row to the order of derivation.
- `β::Vector{<:Real}`: the orders of derivation in the form of a row vector, where
   each element corresponds to the order of one differential equation. It can take
   decimal as well as integer values.
- `J::Matrix{<:Real}`: the Jacobian of F. If not provided, its position within the
   function definition has to be replaced with `nothing`.
- `par...`: additional parameters for the function F.
- `h::Number`: the step size for correction.
- `nc::Int64`: the desired number of corrections.
- `StopIt::String`: the method to stop correction. It can take either "Standard"
   (by default) or "Convergence". In the former case, the function will repeat
   correction as many times as specified in nc; in the latter case, correction will
   stop only when tolerance (tol) or the iteration max (itmax) is reached.
- `tol::Float64`: the tolerance.
- `ìtmax::Int64`: the maximal number of iterations.
"""
module FdeSolver

using SpecialFunctions
using LinearAlgebra

"""
    greet()
Motivates to get started with FdeSolver.
"""
greet() = print("Hey, let's solve some FDEs!")

include("main.jl")
include("main_Jacob.jl")
include("SupFuns.jl")

export(FDEsolver)

end


text:
# todo: find a meaningful name for this file
facts("testing more mpb interface methods") do
    m = MathProgBase.LinearQuadraticModel(SCIPSolver("display/verblevel", 0))

    # take a bit more time integer problem
    A = [314 867 860; 87 875 -695]
    l = [-10, -1000, -1000]
    u = [Inf Inf Inf]
    c = [-1,1,1]
    lb = [40600 -92375]
    ub = lb
    MathProgBase.loadproblem!(m, A, l, u, c, lb, ub, :Min)
    @fact MathProgBase.numconstr(m) --> 2
    @fact MathProgBase.numvar(m) --> 3

    MathProgBase.setvartype!(m, [:Int, :Int, :Int])
    MathProgBase.optimize!(m)
    @fact MathProgBase.status(m) --> :Optimal
    @fact MathProgBase.getsolvetime(m) --> greater_than(0.1)
    @fact MathProgBase.getsolution(m) --> roughly([750, -200, -25])
end


text:
const MUTATE = true
const COSMIC_RAYS_PROBABILITY = 1e-9
const FLAW_PROBABILITY = 1e-4

const SLICE_SIZE = 10

const MAX_SEARCH_DISTANCE = 256
const MAX_TEMPLATE_LENGTH = 16

const MIN_DAUGHTER_SIZE = 8
const MAX_DAUGHTER_GROWTH = 4

const MAX_MEMORY_CAPACITY = 0.5

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
using Distributions

include("./src/smote_exs.jl")
include("./src/ub_smote.jl")

m = 150
X = rand(m, 10)
y = ones(m)
X = hcat(X, y)

smote_exs(X, 11)

w1 = [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]

cases_needed(w1)
pct_needed(w1)

X = randn(150, 10);
y = vcat(zeros(140), ones(10));
X2, y2 = smote(X, y, 5);

length(y2)
countmap(y2)

X = randn(150, 10);
y = vcat(zeros(130), ones(20));
X2, y2 = smote(X, y, 5)

length(y2)
countmap(y2)


text:
"""
# Mineos

The Julia `Mineos` module provides an interface to computing normal
mode eigenfrequencies using the Mineos software.

# List of functions
- `eigengfrequencies`: Compute the normal mode eigenfrequencies for
  a range of radial and angular orders for a 1D model of the Earth.
"""
module Mineos

export
    eigenfrequencies

using DataStructures: OrderedDict

using SeisModels: LinearLayeredModel, write_mineos
using Mineos_jll: minos_bran

"""
    eigenfrequencies(m::LinearLayeredModel, freq=1.0; kwargs...) -> freqs

Compute the eigenfrequencies of normal modes for the model `m`.
`m` must be a `LinearLayeredModel` from `SeisModels`.

Eigenfrequencies are returned as an `OrderedDict` where the keys are
a tuple of the radial order ``n``, mode type (`:S` for spheroidal, `:T`
for toroidal, `:C` for inner core toroidal), and angular order ℓ.

For example, to retrieve ₀S₂, do `freqs[(0,:S,2)]`, and for ₃T₈ write
`freqs[(3,:T,8)]`.  The syntax `freqs[3,:T,8]` is also allowed.

### Keyword arguments
#### Mode selection
- `radial=true`: Calculate (or not) radial modes.  In this case, `nmin`
  and `nmax` are ignored.
- `spheroidal=true`: Calculate (or not) spheroidal modes.
- `toroidal=true`: Calculate (or not) toroidal modes.
- `ic_toroidal=true`: Calculate (or not) inner core toroidal modes.
#### Calculation parameters
- `eps=1e-10`: Accuracy of the integration scheme.  Eigenfrequencies are
  accurate relatively to about `2eps` to `3eps`.  `eps` can be `1e-7` for
  frequencies less than 100 mHz (periods over 10 s); for frequencies of
  100 to 200 mHz (periods 5 to 10 s), `eps` should be between `1e-12` and
  `1e-10`.
- `wgrav=10`: Frequency in mHz above which gravitational terms are neglected.
- `lmin=1`, `lmax=20`: Minimum and maximum angular order of modes.
- `wmin=0.0`, `wmax=166.0`: Minimum and maximum mode eigenfrequencies
  to consider.
- `nmin=0`, `nmax=10`: Minimum and maximum dispersion branch numbers to
  compute.
"""
function eigenfrequencies(m::LinearLayeredModel, freq=1.0;
        # Which modes to calculate
        radial=true, toroidal=true, spheroidal=true, ic_toroidal=true,
        # Mineos calculation parameters
        eps=1e-10, wgrav=10, lmin=1, lmax=20, wmin=0.0, wmax=166.0,
        nmin=0, nmax=10)
    out = OrderedDict{Tuple{Int,Symbol,Int}, Float64}()
    for (jcom, mode_type) in enumerate((:radial, :toroidal, :spheroidal, :ic_toroidal))
        mode_type === :radial && !radial && continue
        mode_type === :toroidal && !toroidal && continue
        mode_type === :spheroidal && !spheroidal && continue
        mode_type === :ic_toroidal && !ic_toroidal && continue
        mktempdir() do dir
            cd(dir) do
                model_in = "model.in"
                model_out = "model.out"
                eigfuncs = "eigenfunctions.out"
                control = "control"
                open(control, "w") do f
                    println(f, """
                        $model_in
                        $model_out
                        $eigfuncs
                        $eps $wgrav
                        $jcom
                        $lmin $lmax $wmin $wmax $nmin $nmax
                        """)
                end
                write_mineos(model_in, m, freq)
                output = minos_bran() do minos_bran_path
                    String(read(pipeline(control, `$minos_bran_path`)))
                end
                _check_minos_bran_stdout(output)
                freqs = read_eigenfrequencies(model_out)
                _check_minos_bran_rayleigh_quotient(freqs, eps, mode_type, scale=10)
                sym = mode_type in (:radial, :spheroidal) ? :S :
                      mode_type === :toroidal ? :T :
                      :C
                for (k, v) in freqs
                    out[(k[1],sym,k[2])] = v.frequency
                end
            end
        end
    end
    out
end

"""
    read_eigenfrequencies(file)

Return the set of eigenfrequencies for a single Mineos model
output file, as produced by the `minos_bran` program.

These are contained in a Dictionary whose keys are tuples
of the angular and radial order of the mode and whose values
are named tuples containing the calculation output.

For example, to access the ₀S₂ output, ask for
`freqs[(0,2)]`.  To get the frequency in mHz, do `freqs[(0,2)].frequency`.

The named tuple names are:
- `phase_vel`: Phase velocity
- `frequency`: Eigenfrequency in mHz
- `period`: Mode period in s
- `group_vel`: Group velocity in km/s
- `Q`: Quality factor
- `rayleigh_quotient`: Ratio of kinetic energy to potential energy
  minus one, which according to the Mineos manual should be small
  (of order `eps` used in the calculation) if the eigenfunction is
  accurate and if the number of radial knots is sufficient.  Strongly
  exponential modes such as Stonely modes are expected to show a
  larger value of `rayleigh_quotient`.

"""
function read_eigenfrequencies(file, type=nothing)
    freqs = OrderedDict{Tuple{Int,Int},Any}()
    open(file, "r") do f
        lines = readlines(f)
        iline_modes = findfirst(x->occursin(r"^\s*mode", x), lines)
        iline_modes === nothing &&
            error("cannot find mode header line in output file")
        iline_modes += 2
        while iline_modes <= length(lines)
            mode_line = lines[iline_modes]
            if isempty(mode_line)
                iline_modes += 1
                continue
            end
            !isempty(mode_line) && length(mode_line) < 108 &&
                error("line $iline_modes is too short")
            radial_order = parse(Int, mode_line[1:5])
            # Check on output mode type if requested
            if type !== nothing
                mode_type = strip(mode_line[6:7])
                mode_type == type ||
                    error("mode type in file is $mode_type, but expecting $type")
            end
            angular_order = parse(Int, mode_line[8:12])
            phase_vel = parse(Float64, mode_line[13:28])
            frequency = parse(Float64, mode_line[29:44])
            period = parse(Float64, mode_line[45:60])
            group_vel = parse(Float64, mode_line[61:76])
            Q = parse(Float64, mode_line[77:92])
            rayleigh_quotient = parse(Float64, mode_line[93:108])
            freqs[(radial_order, angular_order)] = (phase_vel=phase_vel,
                frequency=frequency, period=period, group_vel=group_vel,
                Q=Q, rayleigh_quotient=rayleigh_quotient)
            iline_modes += 1
        end
    end
    freqs
end

function _check_minos_bran_stdout(out)
    lines = split(chomp(out), '\n')
    length(lines) == 12 ||
        error("unexpected output from `minos_bran`:\n$out")
end

"Ensure that the Rayleigh quotient is similar to eps,
and warn about modes for which this is not true."
function _check_minos_bran_rayleigh_quotient(freqs, eps, mode_type; scale=5)
    for (k, v) in freqs
        if v.rayleigh_quotient > scale*eps
            mode_name = if mode_type in (:radial, :spheroidal)
                "S"
            elseif mode_type === :toroidal
                "T"
            else
                "C"
            end
            mode_name = "$(k[1])$(mode_name)$(k[2])"
            @warn("Rayleigh quotient greater than $(scale)*eps for $mode_name")
        end
    end
end

end # module


text:
@testset "grid" begin
    max_x = 3e-6
    max_y = 10e-6
    max_t = 0.1e-12

    nx = 60
    ny = 200

    Δx = max_x / nx
    Δy = max_y / ny

    Δt = 1 / C / √(1/Δx^2 + 1/Δy^2)
    nt = round(Int, max_t/Δt)

    axes_x = Δx * (Base.OneTo(nx) .- nx/2)
    axes_y = Δy * (Base.OneTo(ny) .- Δy/2)

    grid = Grid(max_x, max_y, max_t, nx, ny)

    @test grid.Δx == Δx
    @test grid.Δy == Δy
    @test grid.Δt == Δt
    @test grid.nt == nt

    @test size(grid) == (nx, ny)
    @test size(grid, 1) == nx
    @test size(grid, 2) == ny
    @test axes(grid) == (axes_x, axes_y)
    @test axes(grid, 1) == axes_x
    @test axes(grid, 2) == axes_y
    @test boundary(grid) == (max_x, max_y)
    @test boundary(grid, 1) == max_x
    @test boundary(grid, 2) == max_y

    @test build(grid) == cat(
        repeat(axes_x, 1, ny, 1),
        repeat(axes_y', nx, 1, 1),
        dims=3
    )
end

@testset "Light" begin
    λ = 2.04e-6

    light = Light(λ)

    @test light.λ == λ
    @test light.k == 2π/λ
end

@testset "permittivity" begin
    max_x = 3e-6
    max_y = 10e-6
    max_t = 0.1e-12

    nx = 60
    ny = 200

    ϵ = 9.

    grid = Grid(max_x, max_y, max_t, nx, ny)
    permittivity = Permittivity(ϵ, grid)

    @test permittivity.ϵ[1] == ϵ
    @test permittivity.ϵx[1] == C * grid.Δt/grid.Δx / ϵ
    @test permittivity.ϵy[1] == C * grid.Δt/grid.Δy / ϵ
    @test all_eq(permittivity.ϵ)
    @test all_eq(permittivity.ϵx)
    @test all_eq(permittivity.ϵy)

    ϵ_defect = 1.
    xs_defect = [0, 1e-6, -1e-6]
    ys_defect = [1e-6, 2e-6, 3e-6]
    rs_defect = [0.5e-6, 0.1e-6, 0.2e-6]
    implant!(permittivity, ϵ_defect, xs_defect, ys_defect, rs_defect, grid)

    in_circle(i, j) = true in [
        √(
            (axes(grid, 1)[i] - xs_defect[c])^2 +
            (axes(grid, 2)[j] - ys_defect[c])^2
        ) < rs_defect[c]
        for c in 1:length(rs_defect)
    ]

    @test all(
        in_circle(i, j) ? permittivity.ϵ[i, j] == ϵ_defect : permittivity.ϵ[i, j] == ϵ
        for i in 1:size(grid, 1), j in 1:size(grid, 2)
    )
end

@testset "permeability" begin
    max_x = 3e-6
    max_y = 10e-6
    max_t = 0.1e-12

    nx = 60
    ny = 200

    μ = 1.

    grid = Grid(max_x, max_y, max_t, nx, ny)
    permeability = Permeability(μ, grid)

    @test permeability.μ[1] == μ
    @test permeability.μx[1] == C * grid.Δt/grid.Δx / μ
    @test permeability.μy[1] == C * grid.Δt/grid.Δy / μ
    @test all_eq(permeability.μ)
    @test all_eq(permeability.μx)
    @test all_eq(permeability.μy)
end

@testset "simulation" begin
    # ##########
    # # const. #
    # ##########

    max_x = 3e-6
    max_y = 10e-6
    max_t = 0.1e-12

    nx = 60
    ny = 200

    λ = 2.04e-6

    ϵ = 9.
    μ = 1.

    # ##############
    # # components #
    # ##############

    grid = Grid(max_x, max_y, max_t, nx, ny)
    light = Light(λ)
    permittivity = Permittivity(ϵ, grid)
    permeability = Permeability(μ, grid)

    # #############
    # # simulator #
    # #############

    s = Simulator(nx=nx, ny=ny)
    simulate!(s)

    s = Simulator(grid, light, permittivity, permeability)
    simulate!(s)
end


text:
@doc """
	clip(A::AbstractPoint2D, B::AbstractPoint2D) -> C, D

Clip the line segment with endpoints `A` and `B` to the bounding box.
The returned points `C` and `D` are the endpoints of the intersection.

If the line segment is not intersecting the bounding box, both `C` and `D` are `Void`.
"""->
function clip(A::AbstractPoint2D, B::AbstractPoint2D)
	# TODO: Use Line2D from GeometricalPredicates as input?
	if isinside(A) && isinside(B)
		return A, B
	end

	t0 = 0.0
	t1 = 1.0
	D = B - A
	# left, right, bottom, top
	# Parametrization of line segment from A to B
	p = [-getx(D) ; getx(D) ; -gety(D) ; gety(D)]
	q = [getx(A)-LEFT ; RIGHT-getx(A) ; gety(A)-LOWER ; UPPER-gety(A)]
	for k in 1:4
		if p[k] == 0.0
			# Line parallel with k'th edge
			if q[k] < 0.0
				return nothing, nothing
			end
		elseif p[k] < 0.0
			# Outside to inside
			u = q[k] / p[k]
			if u > t1
				return nothing, nothing
			end
			t0 = max( u, t0 )
		else # p[k] > 0.0
			# Inside to outside
			u = q[k] / p[k]
			if u < t0
				return nothing, nothing
			end
			t1 = min( u, t1 )
		end
	end

	return A + t0*D, A + t1*D
end

@doc """
	isinside(p::Point2D) -> Bool

Test if the point `p` is inside the bounding box.
"""->
function isinside(p::AbstractPoint2D)
	LEFT <= getx(p) <= RIGHT && LOWER <= gety(p) <= UPPER
end

@doc """
	isoutside(edge::VoronoiEdge) -> Bool

Test if `edge` intersects the bounding box.
"""->
function isoutside(edge::VoronoiDelaunay.VoronoiEdge{IndexablePoint2D})
	A, B = clip( geta(edge), getb(edge) )

	if isa(A,Void) && isa(B,Void)
		return true
	else
		return false
	end
end


text:
using Pkg

Pkg.add("https://github.com/csimal/NetworkEpidemics.jl")
Pkg.add("LightGraphs")
Pkg.add("Statistics")
Pkg.add("Plots")
Pkg.add("Colorschemes")
Pkg.add("IJulia") # Jupyter Notebook kernel for Julia

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
# Parent file to run all scripts which may generate
# some output that you want to display on the website.
# this can be used as a tester to check that all the code
# on your website runs properly.

dir = @__DIR__

"""
    genplain(s)

Small helper function to run some code and redirect the output (stdout) to a file.
"""
function genplain(s::String)
    open(joinpath(dir, "scripts/output", "$(splitext(s)[1]).out"), "w") do outf
        redirect_stdout(outf) do
            include(joinpath(dir, "scripts", s))
        end
    end
end
#genplain("gravities.jl")
genplain("flying.jl")
#genplain("volumeTransparent.jl")

#genplain("submarineCables.jl")
#genplain("submarineCables3D.jl")

#genplain("proj4_to_Sphere.jl")
#genplain("delaunayTriangulation.jl")
#genplain("RRGraph3D.jl")
#genplain("RRGraph.jl")
#genplain("volumeScatters.jl")
#genplain("volume.jl")
#genplain("strippedBar.jl")
#genplain("tesseralSphericalHSliders.jl")
#genplain("animScatters.jl")
#genplain("animLinePoint.jl")
#genplain("heatmapScaleAdd.jl")
#genplain("heatmapScaleDetail.jl")
#genplain("heatmapScaleCrosshair.jl")
#genplain("heatmapScaleSections.jl")
#genplain("theme_dark_astronauts.jl")
#genplain("heatmapIrregularCategories.jl")
#genplain("temCubeESDL.jl")
#genplain("theme_dark_Contour3d.jl")
#genplain("surfaceConstraints.jl")
#genplain("candlestick.jl")
#genplain("twinAxis.jl")
#genplain("ScatterHistDensitySides.jl")
#genplain("timeSeries.jl")
#genplain("heatmapText.jl")
#genplain("ContourOverHeatmap.jl")
#genplain("surfWireContour.jl")
#genplain("Fig2Lines.jl")
#genplain("irisDataSet.jl")
#genplain("revolutionSurf.jl")
#genplain("revolutionSurf2.jl")
#genplain("KleinBottle.jl")
#genplain("wireTori.jl")
#genplain("wireSurfToris.jl")
#genplain("tesseralSphericalH.jl")
#genplain("RGBcube.jl")
#genplain("line3D.jl")
#genplain("surface1Color.jl")
#genplain("surfaceComplexF.jl")
#genplain("contourQubit.jl")
#genplain("BlueMarbel.jl")
#genplain("moon.jl")
#genplain("landSea.jl")
#genplain("temPrecESDL.jl")
#genplain("contour1.jl")
#genplain("ContourComplexF.jl")
#genplain("complexPolyaField.jl")
#genplain("themeBlack.jl")
#genplain("heatmapLogIrregular.jl")
#genplain("topoLinesFuji.jl")

#genplain("gppBeer2010.jl")
#genplain("topography.jl")
#genplain("topographySphere.jl")
# output
#using Glob
#files = glob("./scripts/*.jl")
#for f in files #[30:end]
#    println("Generating output for "*f[11:end])
#    genplain(dir*f[11:end])
#end

max_stars_repo_name: lazarusA/BeautifulMakie

max_issues_repo_name: lazarusA/BeautifulMakie

max_forks_repo_name: lazarusA/BeautifulMakie


text:
abstract type BSONValidator end

struct UncheckedBSONValidator <: BSONValidator end
struct LightBSONValidator <: BSONValidator end
struct StrictBSONValidator <: BSONValidator end

@inline validate_field(
    ::BSONValidator,
    type::UInt8,
    p::Ptr{UInt8},
    field_length::Integer,
    available_length::Integer
) = nothing

@inline function validate_field(
    ::LightBSONValidator,
    type::UInt8,
    p::Ptr{UInt8},
    field_length::Integer,
    available::Integer
)
    field_length > available && throw(BSONValidationError("Field length too long, $field_length > $available"))
    nothing
end

@inline function validate_field(
    ::StrictBSONValidator,
    type::UInt8,
    p::Ptr{UInt8},
    field_length::Integer,
    available::Integer
)
    field_length > available && throw(BSONValidationError("Field length too long, $field_length > $available"))
    if type == BSON_TYPE_DOCUMENT || type == BSON_TYPE_ARRAY
        field_length < 5 && throw(BSONValidationError("Document or array under 5 bytes ($field_length)"))
        unsafe_load(p + field_length - 1) != 0 && throw(BSONValidationError("Document or array missing null terminator"))
    end
    nothing
end

@inline validate_root(::BSONValidator, src::DenseVector{UInt8}) = nothing

@inline function validate_root(validator::LightBSONValidator, src::DenseVector{UInt8})
    GC.@preserve src begin
        p = pointer(src)
        len = ltoh(unsafe_load(Ptr{Int32}(p)))
        validate_field(validator, BSON_TYPE_DOCUMENT, p, len, length(src))
    end
end

@inline function validate_root(validator::StrictBSONValidator, src::DenseVector{UInt8})
    GC.@preserve src begin
        p = pointer(src)
        len = ltoh(unsafe_load(Ptr{Int32}(p)))
        validate_field(validator, BSON_TYPE_DOCUMENT, p, len, length(src))
        len < length(src) && throw(BSONValidationError("Garbage after end of document"))
    end
end

@inline validate_bool(::BSONValidator, x::UInt8) = nothing

@inline function validate_bool(::StrictBSONValidator, x::UInt8)
    x == 0x0 || x == 0x1 || throw(BSONValidationError("Invalid bool value $x"))
    nothing
end

@inline validate_string(::BSONValidator, p::Ptr{UInt8}, len::Integer) = nothing

@inline function validate_string(::StrictBSONValidator, p::Ptr{UInt8}, len::Integer)
    unsafe_load(p + len) != 0 && throw(BSONValidationError("Code string missing null terminator"))
    s = unsafe_string(p, len)
    isvalid(s) || throw(BSONValidationError("Invalid UTF-8 string"))
    nothing
end

@inline validate_binary_subtype(::BSONValidator, p::Ptr{UInt8}, len::Integer, subtype::UInt8) = nothing

@inline function validate_binary_subtype(::StrictBSONValidator, p::Ptr{UInt8}, len::Integer, subtype::UInt8)
    if subtype == BSON_SUBTYPE_UUID || subtype == BSON_SUBTYPE_UUID_OLD || subtype == BSON_SUBTYPE_MD5
        len != 16 && throw(BSONValidationError("Subtype $subtype must be 16 bytes long"))
    elseif subtype == BSON_SUBTYPE_BINARY_OLD
        nested_len = ltoh(unsafe_load(Ptr{Int32}(p + 5)))
        nested_len != len - 4 && throw(
            BSONValidationError("Binary (Old) subtype invalid length $nested_len, expected $(len - 4)")
        )
    end
    nothing
end


text:
using Base: depwarn
DeprecatedAlgorithms = Union{Otsu, Yen, MinimumError, Moments,
                             MinimumIntermodes, Intermodes, Balanced,
                             Entropy, UnimodalRosin}
function find_threshold(algorithm::DeprecatedAlgorithms,  histogram::AbstractArray, edges::AbstractRange)
    depwarn("find_threshold(alg, histogram, edges) is deprecated, use find_threshold(histogram, edges, alg) instead", :find_threshold)
    return find_threshold(histogram, edges, algorithm)

end

max_stars_repo_name: zygmuntszpak/HistogramThresholding

max_forks_repo_name: zygmuntszpak/HistogramThresholding


text: Dict{String, Any}("string_array" => Any[Dict{String, Any}("value" => "all", "type" => "string"), Dict{String, Any}("value" => "strings", "type" => "string"), Dict{String, Any}("value" => "are the same", "type" => "string"), Dict{String, Any}("value" => "type", "type" => "string")])