text:
const USE_MIN_RANGE = true

struct AffineEAGO{N}
    c::Float64               # mid-point
    γ::SVector{N,Float64}    # affine terms
    Δ::Float64               # error term
end

function AffineEAGO(x::AffineEAGO{N}, p::Float64, q::Float64, δ::Float64) where N
    c = p * x.c + q
    γ = p .* x.γ
    Δ = p * x.Δ + δ
    AffineEAGO{N}(c, γ, δ)
end
mid(x::Interval{Float64}) = 0.5*(x.lo + x.hi)
AffineEAGO{N}(x::Float64, X::Interval{Float64}, i::Int) where N = AffineEAGO{N}(mid(X), radius(X)*seed_gradient(i, Val(N)), 0.0)

const UNIT_INTERVAL = Interval{Float64}(-1,1)
Interval(x::AffineEAGO{N}) where N = x.c + x.Δ*UNIT_INTERVAL + sum(y -> abs(y)*UNIT_INTERVAL, x.γ)
function bounds(x::AffineEAGO{N}) where N
    z = Interval(x)
    z.lo, z.hi
end

zero(::Type{AffineEAGO{N}}) where N = AffineEAGO{N}(0.0, zero(SVector{N,Float64}), 0.0)
zero(x::AffineEAGO{N}) where N = AffineEAGO{N}(0.0, zero(SVector{N,Float64}), 0.0)

one(::Type{AffineEAGO{N}}) where N = AffineEAGO{N}(1.0, zero(SVector{N,Float64}), 0.0)
one(x::AffineEAGO{N}) where N = AffineEAGO{N}(1.0, zero(SVector{N,Float64}), 0.0)

+(x::AffineEAGO{N}, y::AffineEAGO{N}) where N = AffineEAGO{N}(x.c + y.c, x.γ .+ y.γ, x.Δ + y.Δ)
+(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(α+x.c, x.γ, x.Δ)
+(α::Float64, x::AffineEAGO{N}) where N = x + α

function *(x::AffineEAGO{N}, y::AffineEAGO{N}) where N
    c = x.c * y.c
    γ = x.c .* y.γ + y.c .* x.γ
    pos_sum = 0.0
    neg_sum = 0.0
    for i = 1:N
        vt = y.γ[i]*x.γ[i]
        if vt > 0.0
            pos_sum += vt
        else
            neg_sum -= vt
        end
    end
    Δ = max(pos_sum, neg_sum)
    for i = 1:N
        xi = x.γ[i]
        yi = y.γ[i]
        for j = 1:(i-1)
            Δ += yi*x.γ[j] + y.γ[j]*xi
        end
    end
    return AffineEAGO{N}(c, γ, Δ)
end
*(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(α*x.c, α.*x.γ, abs(α)*x.Δ)
*(α::Float64, x::AffineEAGO{N}) where N = x*α

function -(x::AffineEAGO{N}, y::AffineEAGO{N}) where N
    AffineEAGO{N}(x.c - y.c, x.γ .- y.γ, x.Δ + y.Δ)
end
-(x::AffineEAGO{N}) where N = AffineEAGO{N}(-x.c, .-(x.γ), x.Δ)
-(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(x.c - α, x.γ, x.Δ)
-(α::Float64, x::AffineEAGO{N}) where N = α + (-x)

/(x::AffineEAGO{N}, α::Float64) where N = AffineEAGO{N}(x.c/α, x.γ/α, x.Δ/abs(α))
/(α::Float64, x::AffineEAGO{N}) where N = α*inv(x)

function log(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = log(a)
    fb = log(b)
    if USE_MIN_RANGE
        p = 1/b
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = 1/p
    fξ = log(ξ)
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function log10(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    fa = log10(a)
    fb = log10(b)
    if USE_MIN_RANGE
        p = 1/(b*log(10))
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = 1/p
    fξ = log10(ξ)
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_1d(x::AffineEAGO{N}, n::Number, p) where N
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    if USE_MIN_RANGE
        q = 0.5*(fa + fb - p*(a + b))
        Δ = abs(0.5*(fb - fa - p*(b - a)))
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = (p/n)^(1/(n - 1))
    fξ = ξ^n
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_even(x::AffineEAGO{N}, n::Int) where N
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    if USE_MIN_RANGE
        m = min(0.0, a, b)
        M = max(0.0, a, b)
        p = 0.0
        q = 0.5*(m + M)
        Δ = 0.5*(M - m)
        return AffineEAGO(x, p, q, Δ)
    end
    p = (fb - fa)/(b - a)
    ξ = (p/n)^(1/(n - 1))
    fξ = ξ^n
    q = 0.5*(fa + fξ - p*(a + ξ))
    Δ = abs(0.5*(fξ - fa - p*(ξ - a)))
    return AffineEAGO(x, p, q, Δ)
end

function pow_odd(x::AffineEAGO{N}, n::Int) where N
    a, b = bounds(x)
    fa = a^n
    fb = b^n
    p = (fb - fa)/(b - a)
    q = 0.0
    ξ = (p/n)^(1/(n-1))
    fξ = ξ^n
    Δ = abs(fξ - p*ξ)
    return AffineEAGO(x, p, q, Δ)
end

function ^(x::AffineEAGO{N}, n::Int) where N

    iszero(n) && zero(x)
    isone(n)  && one(x)

    xL, xU = bounds(x)
    if (xL > 0.0) || (xU < 0.0)
        return pow_1d(x, n, (n >= 0) ? n*xL^(n-1) : n*xU^(n-1))
    elseif iseven(n)
        return pow_even(x, n)
    end
    (xL < 0.0 < xU) && error("Undefined domain...")
    return pow_odd(x, n)
end

function ^(x::AffineEAGO{N}, n::Float64) where N
    if isinteger(n)
        return x^Int(n)
    end
    xL, xU = bounds(x)
    (xL < 0.0) && error("Invalid domain...")
    if (n > 1.0) || (n < 0.0)
        return pow_1d(x, n, n*xU^(n-1))
    end
    return pow_1d(x, n, n*xL^(n-1))
end

# DONE...
function ^(x::AffineEAGO{N}, n::Number) where N
    if iszero(n)
        return zero(x)
    elseif isone(n)
        return one(x)
    end
    return x^n
end

function inv(x::AffineEAGO{N}) where N
    a, b = bounds(x)
    (a < 0.0 < b) && error("Invalid domain...")
    if b < 0.0
        return -inv(-x)
    end
    if USE_MIN_RANGE
        p = -1/b^2
        q = -(p*(a + b)^2)/(2*a)
        Δ = -(p*(a - b)^2)/(2*a)
        return AffineEAGO(x, p, q, Δ)
    end
    p = -1/(a*b)
    q = -0.5*p*(sqrt(a) + sqrt(b))^2
    Δ = -0.5*p*(sqrt(a) - sqrt(b))^2
    return AffineEAGO(x, p, q, Δ)
end

struct MCAffPnt{N,T}
    v::MC{N,T}
    box::AffineEAGO{N}
end
MC(x::MCAffPnt{N,T}) where {N,T<:RelaxTag} = x.v
MC(x::MC{N,T}) where {N, T<:RelaxTag} = x

zero(::Type{MCAffPnt{N,T}}) where {N,T} = MCAffPnt{N,T}(zero(MC{N,T}), zero(AffineEAGO{N}))
zero(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(zero(x.v), zero(x.box))

one(::Type{MCAffPnt{N,T}}) where {N,T} = MCAffPnt{N,T}(one(MC{N,T}), one(AffineEAGO{N}))
one(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(one(MC{N,T}), one(AffineEAGO{N}))

+(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v + y.v, x.box + y.box)
+(x::MCAffPnt{N,T}, α::Number) where {N,T} = MCAffPnt{N,T}(x.v + α, x.box + α)
+(α::Number, x::MCAffPnt{N,T}) where {N,T} = x + α

*(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v*y.v, x.box*y.box)
*(x::MCAffPnt{N,T}, α::Number) where {N,T} = MCAffPnt{N,T}(x.v * α, x.box * α)
*(α::Number, x::MCAffPnt{N,T}) where {N,T} = x*α

-(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(x.v-y.v, x.box-y.box)
-(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(-x.v, -x.box)
-(x::MCAffPnt{N,T}, α::Number) where {N,T} = x + (-α)
-(α::Number, x::MCAffPnt{N,T}) where {N,T} = α + (-x)

/(x::MCAffPnt{N,T}, y::MCAffPnt{N,T}) where {N,T} = x*inv(y)
/(x::MCAffPnt{N,T}, α::T) where {N,T} = MCAffPnt{N,T}(x.v/α, x.box/α)
/(α::T, x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(α*inv(x.v), α*inv(x.box))

^(x::MCAffPnt{N,T}, n::Integer) where {N,T} = MCAffPnt{N,T}(x.v^n, x.box^n)
^(x::MCAffPnt{N,T}, n::Number) where {N,T} = MCAffPnt{N,T}(x.v^n, x.box^n)
inv(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(inv(x.v), inv(x.box))

log(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(log(x.v), log(x.box))
log10(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(log10(x.v), log10(x.box))

exp(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(exp(x.v), exp(x.box))
exp10(x::MCAffPnt{N,T}) where {N,T} = MCAffPnt{N,T}(exp10(x.v), exp10(x.box))

function extract_apriori_info(t::Union{RelaxAA,RelaxAAInfo}, v::VariableValues{Float64}, x::AffineEAGO{N}) where {N,T}
    z = 0.0
    for (k,i) in enumerate(t.v)
        l = _bd(v, i)
        u = ubd(v, i)
        z += x.γ[k]*(_val(v, i) - 0.5*(l + u))/(u - l)
    end
    Z = x.c + sum(z -> z*UNIT_INTERVAL, x.γ)
    xcv  = x.c + z - x.Δ
    xcc  = x.c + z + x.Δ
    xcvU = Z.hi - x.Δ
    xccL = Z.lo + x.Δ
    xcvg = x.γ
    xccg = x.γ
    return xcv, xcvU, xcc, xccL, xcvg, xccg
end

function _cut_info(v::VariableValues{Float64}, z::MC{N,T}, x::MCAffPnt{N,T}) where {N,T}
    xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), v, x.box)
    zaff = MC{N,T}(xcv, xcc, Interval(x.box), xcvg, xccg, false)
    return zaff ∩ z
end
function _cut_info(v::VariableValues{Float64}, z::MCAffPnt{N,T}, x::MCAffPnt{N,T}) where {N,T}
    xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), v, x.box)
    zaff = MC{N,T}(xcv, xcc, Interval(x.box), xcvg, xccg, false)
    return zaff ∩ z.v
end

function _cut(t::RelaxAAInfo, b, k, x::MCAffPnt{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, c::Bool, p::Bool) where {N,T<:RelaxTag}
    xt = p ? MCAffPnt{N,T}(set_value_post(x.v, v, s, ϵ), x.box) : x
    xtmc = _cut_info(v, xt, xt)
    xtaffp = MCAffPnt{N,T}(xtmc, xt.box)
    b._info[k] = xtaffp
    return
end

function _cut(t::RelaxAA, b, k, x::MC{N,T}, z::MCAffPnt{N,T}, v::VariableValues, ϵ::Float64, s::Vector{Int}, c::Bool, p::Bool) where {N,T<:RelaxTag}
    xMC = set_value_post(c ? x ∩ x.Intv ∩ Interval(z.box) : x, v, s, ϵ)
    xt = p ? xMC : x
    zt = _cut_info(v, xt, info(b, k))
    _store_set!(b, zt, k)
    return
end

relax_info(s::RelaxAA, n::Int, t::T) where T = MCAffPnt{n,T}
function f_init!(t::RelaxAA, g::DAT, b::RelaxCache)
    tinfo = RelaxAAInfo(t.v)
    for k = node_count(g):-1:1
        c = node_class(g, k)
        if c == EXPRESSION
            fprop!(tinfo, Expression(), g, b, k)
        elseif c == VARIABLE
            fprop!(tinfo, Variable(), g, b, k)
        end
        b._set[k] = info(b, k).v
    end
end

function fprop!(t::RelaxAAInfo, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
    zaff = AffineEAGO{N}(x, Interval(l,u), i)
    zinfo = MCAffPnt{N,T}(z, zaff)
    b._info[k] = zinfo
    return
end
function fprop!(t::RelaxAA, vt::Variable, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    i = first_index(g, k)
    x = val(b, i)
    l = lbd(b, i)
    u = ubd(b, i)
    z = varset(MC{N,T}, rev_sparsity(g, i, k), x, x, l, u)
    z = z ∩ _interval(b, k)
    _store_set!(b, z, k)
end

info_or_set(t::RelaxAAInfo, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag} = info(b, k)
info_or_set(t::RelaxAA, b::RelaxCache{V,N,T}, k) where {V,N,T<:RelaxTag} = set(b, k)
_store_out!(t::RelaxAAInfo, b::RelaxCache{V,N,T}, z, k) where {V,N,T<:RelaxTag} = (b._info[k] = z; nothing)
_store_out!(t::RelaxAA, b::RelaxCache{V,N,T}, z, k) where {V,N,T<:RelaxTag} = _store_set!(b, z, k)

function fprop!(t::Union{RelaxAAInfo,RelaxAA}, v::Val{EXP}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = info_or_set(t, b, child(g, 1, k))
    _cut(t, b, k, exp(x), info(b, k), b.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end
function fprop!(t::Union{RelaxAAInfo,RelaxAA}, v::Val{EXP10}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = info_or_set(t, b, child(g, 1, k))
    _cut(t, b, k, exp10(x), info(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end
function fprop!(t::Union{RelaxAAInfo,RelaxAA}, v::Val{LOG}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = info_or_set(t, b, child(g, 1, k))
    _cut(t, b, k, log(x), info(b, k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
end

function fprop!(t::Union{RelaxAA,RelaxAAInfo}, v::Val{POW}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    x_is_num = is_num(b, x)
    y_is_num = is_num(b, y)
    if y_is_num && isone(num(b, y))
        z = info_or_set(t, b, x)
        store_out!(t, b, z, k)
    elseif y_is_num && iszero(num(b, y))
        store_out!(t, b, zero(info_or_set(t, b, x)), k)
    else
        if !x_is_num && y_is_num
            z = info_or_set(t, b, x)^num(b, y)
        elseif x_is_num && !y_is_num
            z = num(b, x)^info_or_set(t, b, y)
        elseif !x_is_num && !y_is_num
            z = info_or_set(t, b, x)^info_or_set(t, b, y)
        end
    _cut(t, b, k, z, info(b,k), b.ic.v, zero(Float64), sparsity(g,k), b.cut, b.post)
    end
end

function fprop!(t::Union{RelaxAA,RelaxAAInfo}, v::Val{MINUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    x_is_num = is_num(b, x)
    if arity(g, k) == 2
        y = child(g, 2, k)
        y_is_num = is_num(b, y)
        if !x_is_num && y_is_num
            z = info_or_set(t, b, x) - num(b, y)
        elseif x_is_num && !y_is_num
            z = num(b, x) - info_or_set(t, b, y)
        else
            z = info_or_set(t, b, x) - info_or_set(t, b, y)
        end
    else
        z = -info_or_set(t, b, x)
    end
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop!(t::Union{RelaxAA,RelaxAAInfo}, v::Val{DIV}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    if !is_num(b, x) && is_num(b, y)
        z = info_or_set(t, b, x)/num(b, y)
    elseif is_num(b, x) && !is_num(b, y)
        z = num(b, x)/info_or_set(t, b, y)
    else
        z = info_or_set(t, b, x)/info_or_set(t, b, y)
    end
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_2!(t::RelaxAAInfo, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    x_is_num = is_num(b, x)
    y_is_num = is_num(b, y)
    if !x_is_num && y_is_num
        z = info(b, x) + num(b, y)
    elseif x_is_num && !y_is_num
        z = num(b, x) + info(b, y)
    else
        z = info(b, x) + info(b, y)
    end
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_2!(t::RelaxAA, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    x_is_num = is_num(b, x)
    y_is_num = is_num(b, y)
    if !x_is_num && y_is_num
        z = set(b, x) + num(b, y)
    elseif x_is_num && !y_is_num
        z = num(b, x) + set(b, y)
    else
        z = set(b, x) + set(b, y)
    end
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_n!(t::RelaxAAInfo, ::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = one(V)
    znum = one(Float64)
    count = 0
    for i in children(g, k)
        if is_num(b, i)
            znum = znum*num(b, i)
        else
            z = z*info(b, i)
        end
        count += 1
    end
    z = z*znum
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_n!(t::RelaxAA, ::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = one(MC{N,T})
    znum = one(Float64)
    count = 0
    for i in children(g, k)
        if is_num(b, i)
            znum = znum*num(b, i)
        else
            z = z*set(t, b, i)
        end
        count += 1
    end
    z = z*znum
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_2!(t::RelaxAAInfo, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    x_is_num = is_num(b, x)
    y_is_num = is_num(b, y)
    if !x_is_num && y_is_num
        zaff = *(info(b, x), num(b, y))
    elseif x_is_num && !y_is_num
        zaff = *(num(b, x), info(b, y))
    else
        xv = info(b, x)
        yv = info(b, y)
        xv_box = xv.box
        yv_box = yv.box
        z_box = xv_box*yv_box
        xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), b.v, xv_box)
        ycv, ycvU, ycc, yccL, ycvg, yccg = extract_apriori_info(RelaxAA(), b.v, yv_box)
        wIntv = xv.v.Intv*yv.v.Intv
        za_l = McCormick.mult_apriori_kernel(xv.v, yv.v, wIntv, xcv, ycv, xcvU, ycvU, xcvg, ycvg)
        za_u = McCormick.mult_apriori_kernel(-xv.v, -yv.v, wIntv, -xcc, -ycc, -xccL, -yccL, -xccg, -yccg)
        z = xv.v*yv.v
        zaff = MCAffPnt{N,T}(z, z_box)
    end
    _cut(t, b, k, zaff, info(b, k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_2!(t::RelaxAA, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    x = child(g, 1, k)
    y = child(g, 2, k)
    x_is_num = is_num(b, x)
    y_is_num = is_num(b, y)
    if !x_is_num && y_is_num
        z = *(set(b, x), num(b, y))
    elseif x_is_num && !y_is_num
        z = *(num(b, x), set(b, y))
    else
        xs = _set(b, x)
        ys = _set(b, y)
        xinfo = info(b, x)
        yinfo = info(b, y)
        xv = xinfo
        yv = yinfo
        xcv, xcvU, xcc, xccL, xcvg, xccg = extract_apriori_info(RelaxAA(), b.v, xv.box)
        ycv, ycvU, ycc, yccL, ycvg, yccg = extract_apriori_info(RelaxAA(), b.v, yv.box)
        wIntv = xv.v.Intv*yv.v.Intv
        za_l = McCormick.mult_apriori_kernel(xs, ys, wIntv, xcv, ycv, xcvU, ycvU, xcvg, ycvg)
        za_u = McCormick.mult_apriori_kernel(-xs, -ys, wIntv, -xcc, -ycc, -xccL, -yccL, -xccg, -yccg)
        z = (xs*ys) ∩ za_l ∩ za_u
    end
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

#TODO:
function fprop_n!(t::RelaxAAInfo, ::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    z = one(V)
    znum = one(Float64)
    count = 0
    for i in children(g, k)
        if is_num(b, i)
            znum = znum*num(b, i)
        else
            z = z*info(b, i)
        end
        count += 1
    end
    z = z*znum
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop_n!(t::RelaxAA, ::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    ys = one(MC{N,T})
    znum = one(Float64)
    b._mult_temp = one(MCAffPnt{N,T})
    count = 0
    first_set = true
    for i in children(g, k)
        if is_num(b, i)
            znum = znum*num(b, i)
        else
            xs = _set(b, i)
            xinfo = info(b, i)
            xi = xinfo.v
            if !first_set
                ys = b._mult_temp.v
                xcvU, xccL = extract_apriori_info(t, b.v, xinfo.box)
                ycvU, yccL = extract_apriori_info(t, b.v, b._mult_temp.box)
                xcv = xi.cv;       ycv = b._mult_temp.v.cv
                xcc = xi.cc;       ycc = b._mult_temp.v.cc
                xcvg = xi.cv_grad; ycvg = b._mult_temp.v.cv_grad
                xccg = xi.cc_grad; yccg = b._mult_temp.v.cc_grad
                wIntv = xs.Intv*ys.Intv ∩ Interval(xinfo.box*b._mult_temp.box)
                za_l = McCormick.mult_apriori_kernel(xs, ys, wIntv, xcv, ycv, xcvU, ycvU, xcvg, ycvg)
                za_u = McCormick.mult_apriori_kernel(-xs, -ys, wIntv, -xcc, -ycc, -xccL, -yccL, -xccg, -yccg)
                ys = (xs*ys) ∩ za_l ∩ za_u
                b._mult_temp = MCAffPnt{N,T}(ys, b._mult_temp.box*xinfo.box)
            else
                first_set = false
                b._mult_temp = MCAffPnt{N,T}(xi, xinfo.box)
                ys = xs
            end
        end
        count += 1
    end
    z = ys*znum
    _cut(t, b, k, z, info(b,k), b.ic.v, b.ϵ_sg, sparsity(g, k), b.cut, false)
end

function fprop!(t::Union{RelaxAA,RelaxAAInfo}, v::Val{PLUS}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    (arity(g, k) == 2) ? fprop_2!(t, Val(PLUS), g, b, k) : fprop_n!(t, Val(PLUS), g, b, k)
end
function fprop!(t::Union{RelaxAA,RelaxAAInfo}, v::Val{MULT}, g::DAT, b::RelaxCache{V,N,T}, k::Int) where {V,N,T<:RelaxTag}
    (arity(g, k) == 2) ? fprop_2!(t, Val(MULT), g, b, k) : fprop_n!(t, Val(MULT), g, b, k)
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#################### Posterior Statistics ####################
function autocor(c::AbstractChains;
                 lags::Vector=[1, 5, 10, 50],
                 relative::Bool=true)
    if relative
        lags *= step(c)
    elseif any(lags .% step(c) .!= 0)
        throw(ArgumentError("lags do not correspond to thinning interval"))
    end
    labels = map(x -> "Lag " * string(x), lags)

    (niter, nvar, nchain) = size(c.value)

    vals = zeros(nvar, length(lags), nchain)
    for k in 1:nchain
        for v in 1:nvar
            # skipping missing values
            # skipping should be done inside of autocor to ensure correct alignment
            # TODO: modify autocor to deal with missing values
            x = convert(Vector{Float64}, collect(skipmissing(c.value[:,v,k])))
            vals[v, :, k] = autocor(x, lags)
        end
    end
    ChainSummary(vals, c.names, labels, header(c))
end

function cor(c::AbstractChains)
    return ChainSummary(cor(combine(c)), c.names, c.names, header(c))
end

function changerate(c::AbstractChains)

    @assert !any(ismissing.(c.value)) "Change rate comp. doesn't support missing values."

    n, p, m = size(c.value)
    r = zeros(Float64, p, 1, 1)
    r_mv = 0.0
    delta = Array{Bool}(undef, p)
    for k in 1:m
        prev = c.value[1, :, k]
        for i in 2:n
            for j in 1:p
                x = c.value[i, j, k]
                dx = x != prev[j]
                r[j] += dx
                delta[j] = dx
                prev[j] = x
            end
            r_mv += any(delta)
        end
    end
    vals = round.([r; r_mv] / (m * (n - 1)), digits = 3)
    return ChainSummary(vals, [c.names; "Multivariate"], ["Change Rate"], header(c))
end

describe(c::AbstractChains; args...) = describe(stdout, c; args...)

function describe(io::IO,
                  c::AbstractChains;
                  q = [0.025, 0.25, 0.5, 0.75, 0.975],
                  etype=:bm,
                  args...
                 )
    ps_stats = summarystats(c; etype=etype, args...)
    ps_quantiles = quantile(c, q=q)
    println(io, ps_stats.header)
    print(io, "Empirical Posterior Estimates:\n")
    show(io, ps_stats)
    print(io, "Quantiles:\n")
    show(io, ps_quantiles)

end

function hpd(x::Vector{T}; alpha::Real=0.05) where {T<:Real}
    n = length(x)
    m = max(1, ceil(Int, alpha * n))

    y = sort(x)
    a = y[1:m]
    b = y[(n - m + 1):n]
    _, i = findmin(b - a)

    return [a[i], b[i]]
end

function hpd(c::AbstractChains; alpha::Real=0.05)
    pct = first(showoff([100.0 * (1.0 - alpha)]))
    labels = ["$(pct)% Lower", "$(pct)% Upper"]
    vals = permutedims(
        mapslices(x -> hpd(collect(skipmissing(x)), alpha=alpha), c.value, dims = [1, 3]),
        [2, 1, 3]
    )
    return ChainSummary(vals, c.names, labels, header(c))
end

function quantile(c::AbstractChains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975])
    labels = map(x -> string(100 * x) * "%", q)
    vals = permutedims(
        mapslices(x -> quantile(collect(skipmissing(x)), q), c.value, dims = [1, 3]),
        [2, 1, 3]
    )
    return ChainSummary(vals, c.names, labels, header(c))
end

function summarystats(c::AbstractChains; etype=:bm, args...)
    f = x -> [
            mean(x),
            std(x),
            sem(x),
            mcse(x, etype; args...)
    ]

    labels = ["Mean", "SD", "Naive SE", "MCSE", "ESS"]
    vals = permutedims(
        mapslices(x -> f(collect(skipmissing(x))), c.value, dims =  [1, 3]),
        [2, 1, 3]
    )
    stats = [vals  min.((vals[:, 2] ./ vals[:, 4]).^2, size(c.value, 1))]
    return ChainSummary(stats, c.names, labels, header(c))
end

ext: jl


text:
#=

    Seguir tutorial de instalação do Julia demonstrado no README.md
    Como rodar esse script: $ julia doubler.jl

    Objetivo desse script: realizar uma operação básica de dobrar o número.

=#

#=
    Declaração de variável 'a'
    Como visto, variáveis no Julia não são tipadas como em algumas linguagens,
    ou seja, você pode declará-las e o inicializador delas assumirá o valor
    que você colocar como atribuição.
    Ex.: Em a = 1 (a será um inteiro) / Em a = 'Olá' (a será uma string)
=#
x = 333

#   Função chamada 'doubler' que recebe um parâmetro
#   qualquer e retorna esse parâmetro multiplicado por 2
function doubler(x)
    return x * 2
end

#   Printando resultado da função 'doubler' passando o parâmetro 'x' que é 333.
#   Comentário para demonstrar o valor que deve ser apresentado.
#   Seria interessante colocar um assert.equals para ficar mais nítida a compreensão.

println(doubler(x)) # Deve imprimir: 666


text:
abstract type LesserElement <: Element end

mutable struct BabelCall <: LesserElement
    name::AbstractString
end

abstract type Block <: LesserElement end

mutable struct CommentBlock <: Block
    contents::Vector{AbstractString}
end
mutable struct ExampleBlock <: Block
    contents::Vector{AbstractString}
end
mutable struct ExportBlock <: Block
    backend::AbstractString
    contents::Vector{AbstractString}
end
mutable struct SourceBlock <: Block
    lang::Union{AbstractString, Nothing}
    arguments::Union{AbstractString, Nothing}
    contents::Vector{AbstractString}
end
mutable struct VerseBlock <: Block
    contents::Vector{LesserElement}
end

mutable struct Clock{T <: Union{TimestampInactive, TimestampInactiveRange}} <: LesserElement
    timestamp::T
    duration::Union{Nothing, Tuple{Integer, Integer}}
end

mutable struct DiarySexp <: LesserElement
    sexp::AbstractString
end

mutable struct Planning <: LesserElement
    deadline::Union{Timestamp, Nothing}
    scheduled::Union{Timestamp, Nothing}
    closed::Union{Timestamp, Nothing}
end

mutable struct Comment{S <: AbstractString} <: LesserElement
    contents::Vector{S}
end

mutable struct FixedWidth{S <: AbstractString} <: LesserElement
    contents::Vector{S}
end

struct HorizontalRule <: LesserElement end

mutable struct Keyword{V <: Union{<:AbstractString, Vector{Object}, Nothing}} <: LesserElement
    key::AbstractString
    value::V
end

mutable struct AffiliatedKeyword{V <: Union{<:AbstractString, Vector{Object}, Nothing}} <: LesserElement
    key::AbstractString
    optval::Union{<:AbstractString, Vector{Object}, Nothing}
    value::V
end

mutable struct LaTeXEnvironment <: LesserElement
    name::AbstractString
    contents::Vector{AbstractString}
end

mutable struct NodeProperty <: LesserElement
    name::AbstractString
    additive::Bool
    value::AbstractString
end

mutable struct Paragraph <: LesserElement
    contents::Vector{Object}
end

mutable struct TableRow <: LesserElement
    cells::Vector{TableCell}
end
struct TableHrule <: LesserElement end


text:
using Test
using BasicModelInterface
const BMI = BasicModelInterface

# dummy model structs to test the interface
struct UnknownModel end

struct KnownModel
    waterlevel::Matrix{Float64}
end

struct MyModel end

@testset "BasicModelInterface.jl" begin

@testset "no BMI implementation" begin
    # test expected errors for unknown model
    unknown_model = UnknownModel()
    @test_throws MethodError BMI.initialize(unknown_model)
    @test_throws MethodError BMI.initialize(unknown_model, "config_file")
    @test_throws MethodError BMI.initialize(unknown_model, "config_file", "one_too_many")
end

@testset "some BMI implementation" begin
    function BMI.initialize(model::KnownModel)
        known_model.waterlevel[:] .= 0.0
    end

    known_model = KnownModel(ones(3, 4))
    @test all(isone, known_model.waterlevel)
    BMI.initialize(known_model)
    @test all(iszero, known_model.waterlevel)
end

@testset "interface functions are defined" begin
    # not necessarily correct parameters
    # only needed to verify that the function is defined, but does not have methods yet
    model = MyModel()
    time = 1.0
    name = "myparam"
    dest = zeros(5)
    inds = [1, 3, 5]
    value = 3.0
    grid = 2
    spacing = zeros(5)
    origin = zeros(2, 5)
    x = collect(4:0.1:5)
    y = collect(6:0.1:9)
    z = collect(10:0.1:13)
    edge_nodes = zeros(Int, 2, 5)
    face_edges = zeros(Int, 5)
    face_nodes = zeros(Int, 5)
    nodes_per_face = zeros(Int, 5)

    @test_throws MethodError BMI.initialize(model)
    @test_throws MethodError BMI.update(model)
    @test_throws MethodError BMI.update_until(model, time)
    @test_throws MethodError BMI.finalize(model)
    @test_throws MethodError BMI.get_component_name(model)
    @test_throws MethodError BMI.get_input_item_count(model)
    @test_throws MethodError BMI.get_output_item_count(model)
    @test_throws MethodError BMI.get_input_var_names(model)
    @test_throws MethodError BMI.get_output_var_names(model)
    @test_throws MethodError BMI.get_var_grid(model, name)
    @test_throws MethodError BMI.get_var_type(model, name)
    @test_throws MethodError BMI.get_var_units(model, name)
    @test_throws MethodError BMI.get_var_itemsize(model, name)
    @test_throws MethodError BMI.get_var_nbytes(model, name)
    @test_throws MethodError BMI.get_var_location(model, name)
    @test_throws MethodError BMI.get_current_time(model)
    @test_throws MethodError BMI.get_start_time(model)
    @test_throws MethodError BMI.get_end_time(model)
    @test_throws MethodError BMI.get_time_units(model)
    @test_throws MethodError BMI.get_time_step(model)
    @test_throws MethodError BMI.get_value(model, name, dest)
    @test_throws MethodError BMI.get_value_ptr(model, name)
    @test_throws MethodError BMI.get_value_at_indices(model, name, dest, inds)
    @test_throws MethodError BMI.set_value(model, name, value)
    @test_throws MethodError BMI.set_value_at_indices(model, name, inds, value)
    @test_throws MethodError BMI.get_grid_rank(model, grid)
    @test_throws MethodError BMI.get_grid_size(model, grid)
    @test_throws MethodError BMI.get_grid_type(model, grid)
    @test_throws MethodError BMI.get_grid_shape(model, grid)
    @test_throws MethodError BMI.get_grid_spacing(model, grid, spacing)
    @test_throws MethodError BMI.get_grid_origin(model, grid, origin)
    @test_throws MethodError BMI.get_grid_x(model, grid, x)
    @test_throws MethodError BMI.get_grid_y(model, grid, y)
    @test_throws MethodError BMI.get_grid_z(model, grid, z)
    @test_throws MethodError BMI.get_grid_node_count(model, grid)
    @test_throws MethodError BMI.get_grid_edge_count(model, grid)
    @test_throws MethodError BMI.get_grid_face_count(model, grid)
    @test_throws MethodError BMI.get_grid_edge_nodes(model, grid, edge_nodes)
    @test_throws MethodError BMI.get_grid_face_edges(model, grid, face_edges)
    @test_throws MethodError BMI.get_grid_face_nodes(model, grid, face_nodes)
    @test_throws MethodError BMI.get_grid_nodes_per_face(model, grid, nodes_per_face)
end

end  # testset "BasicModelInterface.jl"


text:
using Primes

function generate_integer_hash_function(a,b,cardinality)
    # cardinality must be prime number
    # a,b: 1 ≤ a ≤ p−1  , 0 ≤ b ≤ p−1
    @assert a < cardinality
    @assert a > 0
    @assert b < cardinality
    @assert b >= 0
    @assert isprime(cardinality)
    m = cardinality
    x -> (x^a%m + b)%m  #TODO use memory-efficient modular exponentiation for better performance
end

"""
# function generate_vector_deterministic_pseudorandom
Generates a pseudorandom vector. For the same hash function it always generates the same vector.
If the hash function is an integer hash from a universal family, the elements of the vector will
obey a continuous uniform distribution.
# Arguments:
- `seed`: A random number chosen by the user to ensure the result is deterministic
- `hash_func`: An integer hash function sampled from a universal family
- `dimension`: number of elements in the vector
"""
function generate_vector_deterministic_pseudorandom(seed,hash_func,dimension)
    [hash_func(x+seed) for x in 1:dimension]
end

"""
# function generate_dataset
Generates a pseudorandom dataset that can be used for benchmarking SOMs.
The elements in the vectors are distributed uniformly between 0 and `max_value`.
# Arguments:
- `number`: Number of rows/entries/vectors/individual iris plants in the dataset
- `dimension`: Number of dimensions/elements in each vector
- `max_value`: The maximum value that an element in a vector can have
- `y`: an optional function to map the uniformly distributed elements to a desired non-uniform distribution
- `params`: Randomness parameters chosen by user to ensure deterministic result
"""
function generate_dataset(number, dimension, max_value, y=(x->x), params = Dict(:a => 13, :b => 17, :hash_seed => 23, :prime => 1013))
    hash_func = generate_integer_hash_function(params[:a],params[:b],params[:prime])

    seed = params[:hash_seed]
    matrix = zeros((number,dimension))
    for i in 1:number
        vector = generate_vector_deterministic_pseudorandom(seed,hash_func,dimension)
        vector = map((x->x*max_value/params[:prime]),vector)
        vector = map((x->(x+max_value)/2),vector)
        vector = map(y,vector)
        seed = hash_func(seed+i)
        matrix[i,:] = vector
    end
    convert(DataFrame,matrix)
end

# Benchmarking functions
# ======================

function benchmark_init(train, topol; toroidal = false,
                    normType = :zscore)

    xdim = 10
    ydim = 10

    initSOM(train, xdim, ydim, norm = normType,
                  topol = topol, toroidal = toroidal)
end

function benchmark_train(som,train, kernel = gaussianKernel)
    trainSOM(som, train, 100000, kernelFun = kernel)
end

# Example
# =======

iris = dataset("datasets", "iris")
train = iris[:,1:4]

train = generate_dataset(1000,100,1)
som = benchmark_init(train, :rectangular, toroidal = false)
# som = benchmark_train(som,train)


text:
module POMDPPlayback
    using POMDPs
    using POMDPPolicies
    using Distributions
    using POMDPSimulators

    export PlaybackPolicy
    include("playback_policy.jl")

end # module


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function draweasingfunction(f, pos, w, h)
    @layer begin
        translate(pos)
        setline(0.5)
        box(O, w, h, :stroke)
        for i in 0:0.005:1.0
            circle(Point(-w/2, h/2) + Point(w * i, -f(i, 0, h, 1)), 1)
        end
        text(string(f), Point(w/2 - 20, h/2 - 20), halign=:right)
    end
end

function testeasingfunctions(fname)
    Drawing(1200, 1200, fname)
    background("white")
    origin()

    t = Tiler(1200, 1200, 5, 5)
    margin=5

    for (pos, n) in t
        sethue("black")
        n > length(Luxor.easingfunctions) && continue
        draweasingfunction(Luxor.easingfunctions[n], pos, t.tilewidth-margin, t.tileheight-margin)
    end

    @test finish() == true
end

fname = "easingfunctions.pdf"
testeasingfunctions(fname)
println("...finished test: output in $(fname)")


text:
## assignments
# init methods
assignements(ex::EArgs, assigned=[]) =
    begin map(x -> assignements(x, assigned), ex.args); unique(assigned) end

assignements(ex::EHanded, assigned=[]) =
    begin map(x -> assignements(x, assigned), [ex.rhs; ex.lhs]);  unique(assigned) end

assignements(ex::ECall, assigned=[]) =
    begin map(x -> assignements(x, assigned), [ex.fun; ex.args]); unique(assigned) end

assignements(ex, assigned) = nothing
assignements(ex::Symbol, assigned) = nothing

# collect methods
function assignements(ex::EAssignment, assigned=[])
    assignements(ex, ex.lhs, assigned)
    assignements(ex.rhs, assigned) #propagate on the rhs
    unique(assigned)
end

assignements(ex::EAssignment, lhs::Symbol, assigned) = push!(assigned, lhs)
assignements(ex::EAssignment, lhs::ECall,  assigned) = push!(assigned, lhs.fun)
assignements(ex::EAssignment, lhs::ERef,   assigned) = push!(assigned, lhs.lhs)
assignements(ex::EAssignment, lhs::EDot,   assigned) = assigned #e.g. x.a = b, x is not assigned
assignements(ex::EAssignment, lhs::ETuple, assigned) = append!(assigned, lhs.args)

## variables
variables(ex::Union{Expr,AbstractExpr}, vars=[]) = begin map(x -> variables(x,vars), args(ex)); unique(vars) end
variables(ex::Vector, vars) = begin map(x -> variables(x,vars), ex); unique(vars) end
variables(ex, vars) = nothing

# when we have an assignement switch to 4 args form
function variables(ex::EAssignment, vars=[])
    variables(ex, ex.lhs, ex.rhs, vars)
    unique(vars)
end
# f(x,y) = x*y + C : ignore x and y
function variables(ex::EAssignment, lhs::ECall, rhs, vars)
    ignore = variables(lhs)
    rhs_vars = variables(rhs)
    append!(vars, setdiff(rhs_vars, ignore))
    unique(vars)
end
# otherwise keep going
function variables(ex::EAssignment, lhs, rhs, vars)
    variables(lhs, vars)
    variables(rhs, vars)
    unique(vars)
end

variables(ex::Symbol, vars) = push!(vars, ex)
variables(ex::ECall,  vars)  = variables(ex.args, vars) #skip the function name
variables(ex::EKeyword, vars) = variables(ex.rhs, vars) #skip the lhs
variables(ex::EUsing, vars) = nothing

function _unassigned_variables(ex, eval_type)
    vars = variables(ex)
    assigned = assignements(ex)

    uvars = setdiff(vars, assigned)
    filter(s -> filter_variable(s, eval_type), uvars)
end
unassigned_variables(ex::Expr, eval_type) = _unassigned_variables(econvert(ex), eval_type)
unassigned_variables(ex::AbstractExpr, eval_type) = _unassigned_variables(ex, eval_type)

function filter_variable(s::Symbol, eval_type)
    try
        t = eval_type(s)
        t <: DataType && return false
        t <: Module && return false
        t <: Type && return false
    catch err
    end
    true
end

function indent_body(body, tab)

    lines = split(body,"\n")
    ident_length = [length(l) - length(lstrip(l)) for l in lines]
    line_length  = [length(l) for l in lines]

    min_ident = minimum( ident_length[line_length.>0] )

    lines = [tab * l[(1+min_ident):end] for l in lines]

    join(lines, "\n")
end

function extract_method(body, tab = "    ", eval_type = s -> @eval Main typeof($s))

    ex  = try
        ex = Base.parse_input_line(body)
    catch err
        println(err)
        return ""
    end

    args = unassigned_variables(ex, eval_type)
    sargs = join(args, ", ")

    body = indent_body(body, tab)

"
function ($sargs)
$body
end
"
end


text:
using PandasLite
using Dates
using Test

df = pd.DataFrame(Dict(:name => ["a", "b"], :age => [27, 30]))
age = values(df.age)
age[2] = 31
@test df.loc[1, "age"] == 31

df = pd.read_csv(joinpath(dirname(@__FILE__), "test.csv"))
@test isa(df, PandasLite.DataFrame)

include("test_tabletraits.jl")

@test !isempty(df)
empty!(df)
@test isempty(df)

x = pd.Series([3, 5], index = [:a, :b])

@test x.a == 3
@test x["a"] == 3
@test x.loc["a"] == 3
@test x.b == 5
@test x.iloc[1] == 3
@test x.iloc[2] == 5
@test x.iloc[end] == 5
@test length(x) == 2
@test values(x + 1) == [4, 6]
@test x.sum() == 8
@test eltype(x) == Int
@test all(x.iloc[1:2] == x)

# Rolling
roll = Series([1,2,3,4,5]).rolling(3)
@test isequal(values(roll.mean()), [NaN, NaN, 2.0, 3.0, 4.0])

df = pd.DataFrame()
df["a"] = [1]
df["b"] = pd.to_datetime("2015-01-01")
df["c"] = pd.to_timedelta("0.5 hour")
df["d"] = "abcde"

@test Array(df["a"]) == values(df["a"]) == [1]
@test Array(df["b"]) == values(df["b"]) == [DateTime(2015, 1, 1)]
@test Array(df["c"]) == values(df["c"]) == [Millisecond(1800000)]
@test Array(df["d"]) == values(df["d"]) == ["abcde"]
@test Array(df["a"] / 2) == values(df["a"] / 2) == [0.5]
@test Array(2 / df["a"]) == values(2 / df["a"]) == [2]


text:
"""
    barycentric refinement(mesh) -> refined_mesh

Create the mesh obtained by inserting an extra vertex in the barycenters of all cells and
recusively creating fine cells by connecting the barycenter of a k-cell to the already
constructed refined (k-1)-cells on its boundary.
"""
function barycentric_refinement(mesh::Mesh{U,2}) where U

    # Get the points and the faces (segments) we will use in the refinement process
    Edges = skeleton(mesh, 1)

    # Note thier number and use it to find the new number of points and segments(faces)
    # after refienmnts
    NV, NE = numvertices(Edges), numcells(Edges)

    nv = NV + NE
    ne = 2NE

    # define the new vertices array that will hold the old points(coarse) and the
    # new points(finer mesh) of vertices
    verts = similar(mesh.vertices, nv)

    # Now we assign the coarse mesh vetrices (old points) in the upper half excatly
    # as the original coarse mesh
    for i in 1:NV
        verts[i] = mesh.vertices[i]
    end

    # for the second half we store the refienment points and we calculate them using
    # the old points. We insert a point between every pair of points in the original
    # coarse mesh by (point1+ point2)\2
    for (i, Face) in enumerate(cells(mesh))
        verts[NV+i] = cartesian(center(chart(mesh, Face)))
    end

    # Now we create the faces using the new points as well
    # example:
    # oringinal points_index (1,2,3,4) , original faces ((1,2), (2,3), (3,4),(4,5))
    # new points_index(1,2,3,4,5,6,7,8), new_faces((1,5),(5,2),(2,6),(6,3),(3,7),(7,4),(4,8),(8,1))
    #                  a,b      For       E=1 ->   (a,c),(c,b) and so on ..
    edges = similar(Edges.faces, ne)
    for E in 1:NE
        c = NV + E
        a = Edges.faces[E][1]
        b = Edges.faces[E][2]
        edges[2(E-1) + 1] = SVector(a,c)
        edges[2(E-1) + 2] = SVector(c,b)
    end

    # return the new mesh after refienments
    return Mesh(verts, edges)
end

function barycentric_refinement(mesh::Mesh{U,3}) where U

    D1 = 3

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE + NF
    verts = Array{vertextype(mesh)}(undef,nv)
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each edge centroid
    for (E,Edge) in enumerate(cells(edges))
        verts[NV+E] = cartesian(center(chart(edges, Edge)))
    end

    # add a vertex in each face centroid
    for (F,Face) in enumerate(cells(faces))
        verts[NV+NE+F] = cartesian(center(chart(faces, Face)))
    end

    D = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add six faces in each coarse face
    nf = 6NF
    fcs = zeros(celltype(mesh), nf)
    for F in 1 : numcells(faces)
        c = NV + NE + F

        for i in nzrange(D, F)
            E = rows[i]
            e = NV + E

            j = abs(vals[i]) # local index of edge E in face F
            a = mesh.faces[F][mod1(j+1,3)]
            b = mesh.faces[F][mod1(j+2,3)]

            fcs[6(F-1)+2(j-1)+1] = index(a,e,c)
            fcs[6(F-1)+2(j-1)+2] = index(b,c,e)
        end
    end

    Mesh(verts, fcs)
end

"""
    bisecting_refinement(mesh) -> refinement

Construct a refinement of `mesh` by inserting a new vertex on every existing edge.
Every face is subdived in four small faces.

Only defined for 2D meshes.
"""
function bisecting_refinement(mesh::Mesh{U,3}) where U

    D1 = 3

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE
    verts = Array{vertextype(mesh)}(undef,nv)

    # copy over the existing vertices
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each edge centroid
    for (E, Edge) in enumerate(cells(edges))
        verts[NV+E] = cartesian(center(chart(edges, Edge)))
    end

    # Build a matrix that given a coarse face gives
    # access to the coarse edges making up its boundary
    D = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add four faces in each coarse face
    nf = 4NF
    fcs = zeros(celltype(mesh), nf)

    for F in 1 : numcells(faces)

        # retrieve the indices of the fine vertices in the centroids
        # of the edges in the boundary of coarse face F
        es = zeros(Int,3)
        for i in nzrange(D, F)
            E = rows[i]
            j = abs(vals[i]) # local index of edge E in face F
            @assert 1 <= j <= 3
            e = NV + E
            es[j] = e
        end

        # build the four triangles
        r1 = mesh.faces[F][1]
        r2 = mesh.faces[F][2]
        r3 = mesh.faces[F][3]

        fcs[4(F-1)+1] = index(r1,es[3],es[2])
        fcs[4(F-1)+2] = index(r2,es[1],es[3])
        fcs[4(F-1)+3] = index(r3,es[2],es[1])
        fcs[4(F-1)+4] = index(es[1],es[2],es[3])
    end

    Mesh(verts, fcs)
end

max_stars_repo_name: shakisparki/CompScienceMeshes

max_issues_repo_name: shakisparki/CompScienceMeshes

max_forks_repo_name: shakisparki/CompScienceMeshes


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

@inline num_faces(elem::Tri) = 3
@inline num_faces(elem::Quad) = 4

#     compute_triangle_area(tri)
#
# Computes the area of a triangle given `tri`, which is a tuple of three points (vectors),
# using the [Shoelace_formula](https://en.wikipedia.org/wiki/Shoelace_formula).
function compute_triangle_area(tri)
    A, B, C = tri
    return 0.5 * (A[1] * (B[2] - C[2]) + B[1] * (C[2]-A[2]) + C[1] * (A[2] - B[2]))
end

#   reference_plotting_triangulation(reference_plotting_coordinates)
#
# Computes a triangulation of the points in `reference_plotting_coordinates`, which is a tuple containing
# vectors of plotting points on the reference element (e.g., reference_plotting_coordinates = (r,s)).
# The reference element is assumed to be [-1,1]^d.
#
# This function returns `t` which is a `3 x N_tri` Matrix{Int} containing indices of triangles in the
# triangulation of the plotting points, with zero-volume triangles removed.
#
# For example, r[t[1, i]] returns the first reference coordinate of the 1st point on the ith triangle.
function reference_plotting_triangulation(reference_plotting_coordinates, tol=50*eps())
  # on-the-fly triangulation of plotting nodes on the reference element
  tri_in = Triangulate.TriangulateIO()
  tri_in.pointlist = permutedims(hcat(reference_plotting_coordinates...))
  tri_out, _ = Triangulate.triangulate("Q", tri_in)
  triangles = tri_out.trianglelist

  # filter out sliver triangles
  has_volume = fill(true, size(triangles, 2))
  for i in axes(triangles, 2)
      ids = @view triangles[:, i]
      x_points = @view tri_out.pointlist[1, ids]
      y_points = @view tri_out.pointlist[2, ids]
      area = compute_triangle_area(zip(x_points, y_points))
      if abs(area) < tol
          has_volume[i] = false
      end
  end
  return permutedims(triangles[:, findall(has_volume)])
end

# This function is used to avoid type instabilities when calling `digest_solution_variables`.
function transform_to_solution_variables!(u, solution_variables, equations)
  for (i, u_i) in enumerate(u)
    u[i] = solution_variables(u_i, equations)
  end
end

#     global_plotting_triangulation_triplot(u_plot, rst_plot, xyz_plot)
#
# Returns (plotting_coordinates_x, plotting_coordinates_y, ..., plotting_values, plotting_triangulation).
# Output can be used with TriplotRecipes.DGTriPseudocolor(...).
#
# Inputs:
#   - xyz_plot = plotting points (tuple of matrices of size (Nplot, K))
#   - u_plot = matrix of size (Nplot, K) representing solution to plot.
#   - t = triangulation of reference plotting points
function global_plotting_triangulation_triplot(xyz_plot, u_plot, t)

  @assert size(first(xyz_plot), 1) == size(u_plot, 1) "Row dimension of u_plot does not match row dimension of xyz_plot"

  # build discontinuous data on plotting triangular mesh
  num_plotting_points, num_elements = size(u_plot)
  num_reference_plotting_triangles = size(t, 1)
  num_plotting_elements_total = num_reference_plotting_triangles * num_elements

  # each column of `tp` corresponds to a vertex of a plotting triangle
  tp = zeros(Int32, 3, num_plotting_elements_total)
  zp = similar(tp, eltype(u_plot))
  for e = 1:num_elements
    for i = 1:num_reference_plotting_triangles
      tp[:, i + (e-1)*num_reference_plotting_triangles] .= @views t[i, :] .+ (e-1) * num_plotting_points
      zp[:, i + (e-1)*num_reference_plotting_triangles] .= @views u_plot[t[i, :], e]
    end
  end
  return vec.(xyz_plot)..., zp, tp
end

#     mesh_plotting_wireframe(rd::RefElemData{2}, md::MeshData{2})
#
# Generates data for plotting a mesh wireframe given StartUpDG data types.
# Returns (plotting_coordinates_x, plotting_coordinates_y) for a 2D mesh wireframe.
function mesh_plotting_wireframe(rd::RefElemData{2}, md::MeshData{2}; num_plotting_points=25)

  # Construct 1D plotting interpolation matrix `Vp1D` for a single face
  @unpack N, Fmask = rd
  vandermonde_matrix_1D = StartUpDG.vandermonde(Line(), N, StartUpDG.nodes(Line(), N))
  rplot = LinRange(-1, 1, num_plotting_points)
  Vp1D = StartUpDG.vandermonde(Line(), N, rplot) / vandermonde_matrix_1D

  num_face_points = N+1
  num_faces_total = num_faces(rd.elementType) * md.num_elements
  xf, yf = map(x->reshape(view(x, Fmask, :), num_face_points, num_faces_total), md.xyz)

  num_face_plotting_points = size(Vp1D, 1)
  x_mesh, y_mesh = ntuple(_->zeros(num_face_plotting_points, num_faces_total), 2)
  for f in 1:num_faces_total
    mul!(view(x_mesh, :, f), Vp1D, view(xf, :, f))
    mul!(view(y_mesh, :, f), Vp1D, view(yf, :, f))
  end

  return x_mesh, y_mesh
end

# These methods are used internally to set the default value of the solution variables:
# - If a `cons2prim` for the given `equations` exists, use it
# - Otherwise, use `cons2cons`, which is defined for all systems of equations
digest_solution_variables(equations, solution_variables) = solution_variables
function digest_solution_variables(equations, solution_variables::Nothing)
  if hasmethod(cons2prim, Tuple{AbstractVector, typeof(equations)})
    return cons2prim
  else
    return cons2cons
  end
end

"""
    adapt_to_mesh_level!(u_ode, semi, level)
    adapt_to_mesh_level!(sol::Trixi.TrixiODESolution, level)

Like [`adapt_to_mesh_level`](@ref), but modifies the solution and parts of the
semidiscretization (mesh and caches) in place.
"""
function adapt_to_mesh_level!(u_ode, semi, level)
  # Create AMR callback with controller that refines everything towards a single level
  amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable=first), base_level=level)
  amr_callback = AMRCallback(semi, amr_controller, interval=0)

  # Adapt mesh until it does not change anymore
  has_changed = amr_callback.affect!(u_ode, semi, 0.0, 0)
  while has_changed
    has_changed = amr_callback.affect!(u_ode, semi, 0.0, 0)
  end

  return u_ode, semi
end

adapt_to_mesh_level!(sol::TrixiODESolution, level) = adapt_to_mesh_level!(sol.u[end], sol.prob.p, level)

"""
    adapt_to_mesh_level(u_ode, semi, level)
    adapt_to_mesh_level(sol::Trixi.TrixiODESolution, level)

Use the regular adaptive mesh refinement routines to adaptively refine/coarsen the solution `u_ode`
with semidiscretization `semi` towards a uniformly refined grid with refinement level `level`. The
solution and semidiscretization are copied such that the original objects remain *unaltered*.

A convenience method accepts an ODE solution object, from which solution and semidiscretization are
extracted as needed.

See also: [`adapt_to_mesh_level!`](@ref)
"""
function adapt_to_mesh_level(u_ode, semi, level)
  # Create new semidiscretization with copy of the current mesh
  mesh, _, _, _ = mesh_equations_solver_cache(semi)
  new_semi = remake(semi, mesh=deepcopy(mesh))

  return adapt_to_mesh_level!(deepcopy(u_ode), new_semi, level)
end

adapt_to_mesh_level(sol::TrixiODESolution, level) = adapt_to_mesh_level(sol.u[end], sol.prob.p, level)

# Extract data from a 2D/3D DG solution and prepare it for visualization as a heatmap/contour plot.
#
# Returns a tuple with
# - x coordinates
# - y coordinates
# - nodal 2D data
# - x vertices for mesh visualization
# - y vertices for mesh visualization
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_data_2d(center_level_0, length_level_0, leaf_cells, coordinates, levels, ndims,
                     unstructured_data, n_nodes,
                     grid_lines=false, max_supported_level=11, nvisnodes=nothing,
                     slice=:xy, point=(0.0, 0.0, 0.0))
  # Determine resolution for data interpolation
  max_level = maximum(levels)
  if max_level > max_supported_level
    error("Maximum refinement level $max_level is higher than " *
          "maximum supported level $max_supported_level")
  end
  max_available_nodes_per_finest_element = 2^(max_supported_level - max_level)
  if nvisnodes === nothing
    max_nvisnodes = 2 * n_nodes
  elseif nvisnodes == 0
    max_nvisnodes = n_nodes
  else
    max_nvisnodes = nvisnodes
  end
  nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes)
  resolution = nvisnodes_at_max_level * 2^max_level
  nvisnodes_per_level = [2^(max_level - level)*nvisnodes_at_max_level for level in 0:max_level]
  # nvisnodes_per_level is an array (accessed by "level + 1" to accommodate
  # level-0-cell) that contains the number of visualization nodes for any
  # refinement level to visualize on an equidistant grid

  if ndims == 3
    (unstructured_data, coordinates, levels,
        center_level_0) = unstructured_3d_to_2d(unstructured_data,
        coordinates, levels, length_level_0, center_level_0, slice,
        point)
  end

  # Normalize element coordinates: move center to (0, 0) and domain size to [-1, 1]²
  n_elements = length(levels)
  normalized_coordinates = similar(coordinates)
  for element_id in 1:n_elements
    @views normalized_coordinates[:, element_id] .= (
          (coordinates[:, element_id] .- center_level_0) ./ (length_level_0 / 2 ))
  end

  # Interpolate unstructured DG data to structured data
  (structured_data =
      unstructured2structured(unstructured_data, normalized_coordinates,
                              levels, resolution, nvisnodes_per_level))

  # Interpolate cell-centered values to node-centered values
  node_centered_data = cell2node(structured_data)

  # Determine axis coordinates for contour plot
  xs = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[1]
  ys = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[2]

  # Determine element vertices to plot grid lines
  if grid_lines
    mesh_vertices_x, mesh_vertices_y = calc_vertices(coordinates, levels, length_level_0)
  else
    mesh_vertices_x = Vector{Float64}(undef, 0)
    mesh_vertices_y = Vector{Float64}(undef, 0)
  end

  return xs, ys, node_centered_data, mesh_vertices_x, mesh_vertices_y
end

# Extract data from a 1D DG solution and prepare it for visualization as a line plot.
# This returns a tuple with
# - x coordinates
# - nodal 1D data
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_data_1d(original_nodes, unstructured_data, nvisnodes)
  # Get the dimensions of u; where n_vars is the number of variables, n_nodes the number of nodal values per element and n_elements the total number of elements.
  n_nodes, n_elements, n_vars = size(unstructured_data)

  # Set the amount of nodes visualized according to nvisnodes.
  if nvisnodes === nothing
    max_nvisnodes = 2 * n_nodes
  elseif nvisnodes == 0
    max_nvisnodes = n_nodes
  else
    @assert nvisnodes >= 2 "nvisnodes must be zero or >= 2"
    max_nvisnodes = nvisnodes
  end

  interpolated_nodes = Array{eltype(original_nodes),    2}(undef, max_nvisnodes, n_elements)
  interpolated_data  = Array{eltype(unstructured_data), 3}(undef, max_nvisnodes, n_elements, n_vars)

  for j in 1:n_elements
    # Interpolate on an equidistant grid.
    interpolated_nodes[:, j] .= range(original_nodes[1,1,j], original_nodes[1,end,j], length = max_nvisnodes)
  end

  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)
  nodes_out = collect(range(-1, 1, length = max_nvisnodes))

  # Calculate vandermonde matrix for interpolation.
  vandermonde = polynomial_interpolation_matrix(nodes_in, nodes_out)

  # Iterate over all variables.
  for v in 1:n_vars
    # Interpolate data for each element.
    for element in 1:n_elements
      multiply_scalar_dimensionwise!(@view(interpolated_data[:, element, v]),
        vandermonde, @view(unstructured_data[:, element, v]))
    end
  end
  # Return results after data is reshaped
  return vec(interpolated_nodes), reshape(interpolated_data, :, n_vars), vcat(original_nodes[1, 1, :], original_nodes[1, end, end])
end

# Change order of dimensions (variables are now last) and convert data to `solution_variables`
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function get_unstructured_data(u, solution_variables, mesh, equations, solver, cache)

  if solution_variables === cons2cons
    raw_data = u
    n_vars = size(raw_data, 1)
  else
    # FIXME: Remove this comment once the implementation following it has been verified
    # Reinterpret the solution array as an array of conservative variables,
    # compute the solution variables via broadcasting, and reinterpret the
    # result as a plain array of floating point numbers
    # raw_data = Array(reinterpret(eltype(u),
    #        solution_variables.(reinterpret(SVector{nvariables(equations),eltype(u)}, u),
    #                   Ref(equations))))
    # n_vars = size(raw_data, 1)
    n_vars_in = nvariables(equations)
    n_vars = length(solution_variables(get_node_vars(u, equations, solver), equations))
    raw_data = Array{eltype(u)}(undef, n_vars, Base.tail(size(u))...)
    reshaped_u = reshape(u, n_vars_in, :)
    reshaped_r = reshape(raw_data, n_vars, :)
    for idx in axes(reshaped_u, 2)
      reshaped_r[:, idx] = solution_variables(get_node_vars(reshaped_u, equations, solver, idx), equations)
    end
  end

  unstructured_data = Array{eltype(raw_data)}(undef,
                                              ntuple((d) -> nnodes(solver), ndims(equations))...,
                                              nelements(solver, cache), n_vars)
  for variable in 1:n_vars
    @views unstructured_data[.., :, variable] .= raw_data[variable, .., :]
  end

  return unstructured_data
end

# Convert cell-centered values to node-centered values by averaging over all
# four neighbors and making use of the periodicity of the solution
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function cell2node(cell_centered_data)
  # Create temporary data structure to make the averaging algorithm as simple
  # as possible (by using a ghost layer)
  tmp = similar(first(cell_centered_data), size(first(cell_centered_data)) .+ (2, 2))

  # Create output data structure
  resolution_in, _ = size(first(cell_centered_data))
  resolution_out = resolution_in + 1
  node_centered_data = [Matrix{Float64}(undef, resolution_out, resolution_out)
                        for _ in 1:length(cell_centered_data)]

  for (cell_data, node_data) in zip(cell_centered_data, node_centered_data)
    # Fill center with original data
    tmp[2:end-1, 2:end-1] .= cell_data

    # Fill sides with opposite data (periodic domain)
    # x-direction
    tmp[1,   2:end-1] .= cell_data[end, :]
    tmp[end, 2:end-1] .= cell_data[1,   :]
    # y-direction
    tmp[2:end-1, 1, ] .= cell_data[:, end]
    tmp[2:end-1, end] .= cell_data[:, 1, ]
    # Corners
    tmp[1,   1, ] = cell_data[end, end]
    tmp[end, 1, ] = cell_data[1,   end]
    tmp[1,   end] = cell_data[end, 1, ]
    tmp[end, end] = cell_data[1,   1, ]

    # Obtain node-centered value by averaging over neighboring cell-centered values
    for j in 1:resolution_out
      for i in 1:resolution_out
        node_data[i, j] = (tmp[i,   j, ] +
                           tmp[i+1, j, ] +
                           tmp[i,   j+1] +
                           tmp[i+1, j+1]) / 4
      end
    end
  end

  # Transpose
  for (index, data) in enumerate(node_centered_data)
    node_centered_data[index] = permutedims(data)
  end

  return node_centered_data
end

# Convert 3d unstructured data to 2d data.
# Additional to the new unstructured data updated coordinates, levels and
# center coordinates are returned.
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function unstructured_3d_to_2d(unstructured_data, coordinates, levels,
                               length_level_0, center_level_0, slice,
                               point)
  if slice === :yz
    slice_dimension = 1
    other_dimensions = [2, 3]
  elseif slice === :xz
    slice_dimension = 2
    other_dimensions = [1, 3]
  elseif slice === :xy
    slice_dimension = 3
    other_dimensions = [1, 2]
  else
    error("illegal dimension '$slice', supported dimensions are :yz, :xz, and :xy")
  end

  # Limits of domain in slice dimension
  lower_limit = center_level_0[slice_dimension] - length_level_0 / 2
  upper_limit = center_level_0[slice_dimension] + length_level_0 / 2

  @assert length(point) >= 3 "Point must be three-dimensional."
  if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit
    error(string("Slice plane is outside of domain.",
        " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit"))
  end

  # Extract data shape information
  n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # New unstructured data has one dimension less.
  # The redundant element ids are removed later.
  @views new_unstructured_data = similar(unstructured_data[1, ..])

  # Declare new empty arrays to fill in new coordinates and levels
  new_coordinates = Array{Float64}(undef, 2, n_elements)
  new_levels = Array{eltype(levels)}(undef, n_elements)

  # Counter for new element ids
  new_id = 0

  # Save vandermonde matrices in a Dict to prevent redundant generation
  vandermonde_to_2d = Dict()

  # Permute dimensions such that the slice dimension is always the
  # third dimension of the array. Below we can always interpolate in the
  # third dimension.
  if slice === :yz
    unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5])
  elseif slice === :xz
    unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5])
  end

  for element_id in 1:n_elements
    # Distance from center to border of this element (half the length)
    element_length = length_level_0 / 2^levels[element_id]
    min_coordinate = coordinates[:, element_id] .- element_length / 2
    max_coordinate = coordinates[:, element_id] .+ element_length / 2

    # Check if slice plane and current element intersect.
    # The first check uses a "greater but not equal" to only match one cell if the
    # slice plane lies between two cells.
    # The second check is needed if the slice plane is at the upper border of
    # the domain due to this.
    if !((min_coordinate[slice_dimension] <= point[slice_dimension] &&
          max_coordinate[slice_dimension] > point[slice_dimension]) ||
        (point[slice_dimension] == upper_limit &&
          max_coordinate[slice_dimension] == upper_limit))
      # Continue for loop if they don't intersect
      continue
    end

    # This element is of interest
    new_id += 1

    # Add element to new coordinates and levels
    new_coordinates[:, new_id] = coordinates[other_dimensions, element_id]
    new_levels[new_id] = levels[element_id]

    # Construct vandermonde matrix (or load from Dict if possible)
    normalized_intercept =
        (point[slice_dimension] - min_coordinate[slice_dimension]) /
        element_length * 2 - 1

    if haskey(vandermonde_to_2d, normalized_intercept)
      vandermonde = vandermonde_to_2d[normalized_intercept]
    else
      # Generate vandermonde matrix to interpolate values at nodes_in to one value
      vandermonde = polynomial_interpolation_matrix(nodes_in, [normalized_intercept])
      vandermonde_to_2d[normalized_intercept] = vandermonde
    end

    # 1D interpolation to specified slice plane
    # We permuted the dimensions above such that now the dimension in which
    # we will interpolate is always the third one.
    for i in 1:n_nodes_in
      for ii in 1:n_nodes_in
        # Interpolate in the third dimension
        data = unstructured_data[i, ii, :, element_id, :]

        value = multiply_dimensionwise(vandermonde, permutedims(data))
        new_unstructured_data[i, ii, new_id, :] = value[:, 1]
      end
    end
  end

  # Remove redundant element ids
  unstructured_data = new_unstructured_data[:, :, 1:new_id, :]
  new_coordinates = new_coordinates[:, 1:new_id]
  new_levels = new_levels[1:new_id]

  center_level_0 = center_level_0[other_dimensions]

  return unstructured_data, new_coordinates, new_levels, center_level_0
end

# Convert 2d unstructured data to 1d slice and interpolate them.
function unstructured_2d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, point)

  if slice === :x
    slice_dimension = 2
    other_dimension = 1
  elseif slice === :y
    slice_dimension = 1
    other_dimension = 2
  else
    error("illegal dimension '$slice', supported dimensions are :x and :y")
  end

  # Set up data structures to stroe new 1D data.
  @views new_unstructured_data = similar(unstructured_data[1, ..])
  @views new_nodes = similar(original_nodes[1, 1, ..])

  n_nodes_in, _, n_elements, n_variables = size(unstructured_data)
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Test if point lies in the domain.
  lower_limit = original_nodes[1, 1, 1, 1]
  upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_elements]

  @assert length(point) >= 2 "Point must be two-dimensional."
  if point[slice_dimension] < lower_limit || point[slice_dimension] > upper_limit
    error(string("Slice axis is outside of domain. ",
        " point[$slice_dimension]=$(point[slice_dimension]) must be between $lower_limit and $upper_limit"))
  end

  # Count the amount of new elements.
  new_id = 0

  # Permute dimensions so that the slice dimension is always in the correct place for later use.
  if slice === :y
    original_nodes = permutedims(original_nodes, [1, 3, 2, 4])
    unstructured_data = permutedims(unstructured_data, [2, 1, 3, 4])
  end

  # Iterate over all elements to find the ones that lie on the slice axis.
  for element_id in 1:n_elements
    min_coordinate = original_nodes[:, 1, 1, element_id]
    max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, element_id]
    element_length = max_coordinate - min_coordinate

    # Test if the element is on the slice axis. If not just continue with the next element.
    if !((min_coordinate[slice_dimension] <= point[slice_dimension] &&
        max_coordinate[slice_dimension] > point[slice_dimension]) ||
        (point[slice_dimension] == upper_limit && max_coordinate[slice_dimension] == upper_limit))

        continue
    end

    new_id += 1

    # Construct vandermonde matrix for interpolation of each 2D element to a 1D element.
    normalized_intercept =
          (point[slice_dimension] - min_coordinate[slice_dimension]) /
          element_length[1] * 2 - 1
    vandermonde = polynomial_interpolation_matrix(nodes_in, normalized_intercept)

    # Interpolate to each node of new 1D element.
    for v in 1:n_variables
      for node in 1:n_nodes_in
        new_unstructured_data[node, new_id, v] = (vandermonde*unstructured_data[node, :, element_id, v])[1]
      end
    end

    new_nodes[:, new_id] = original_nodes[other_dimension, :, 1, element_id]
  end

  return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), new_unstructured_data[:, 1:new_id, :], nvisnodes)
end

# Calculate the arc length of a curve given by ndims x npoints point coordinates (piece-wise linear approximation)
function calc_arc_length(coordinates)
  n_points = size(coordinates)[2]
  arc_length = zeros(n_points)
  for i in 1:n_points-1
    arc_length[i+1] = arc_length[i] + sqrt(sum((coordinates[:,i]-coordinates[:,i+1]).^2))
  end
  return arc_length
end

# Convert 2d unstructured data to 1d data at given curve.
function unstructured_2d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, curve, mesh, solver, cache)

  n_points_curve = size(curve)[2]
  n_nodes, _, n_elements, n_variables = size(unstructured_data)
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)

  # Check if input is correct.
  min = original_nodes[:, 1, 1, 1]
  max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_elements]
  @assert size(curve) == (2, size(curve)[2]) "Coordinates along curve must be 2xn dimensional."
  for element in 1:n_points_curve
    @assert (prod(vcat(curve[:, n_points_curve] .>= min, curve[:, n_points_curve]
            .<= max))) "Some coordinates from `curve` are outside of the domain.."
  end

  # Set nodes acording to the length of the curve.
  arc_length = calc_arc_length(curve)

  # Setup data structures.
  data_on_curve = Array{Float64}(undef, n_points_curve, n_variables)
  temp_data = Array{Float64}(undef, n_nodes, n_points_curve, n_variables)

  # For each coordinate find the corresponding element with its id.
  element_ids = get_elements_by_coordinates(curve, mesh, solver, cache)

  # Iterate over all found elements.
  for element in 1:n_points_curve

    min_coordinate = original_nodes[:, 1, 1, element_ids[element]]
    max_coordinate = original_nodes[:, n_nodes, n_nodes, element_ids[element]]
    element_length = max_coordinate - min_coordinate

    normalized_coordinates = (curve[:, element] - min_coordinate)/element_length[1]*2 .-1

    # Interpolate to a single point in each element.
    vandermonde_x = polynomial_interpolation_matrix(nodes_in, normalized_coordinates[1])
    vandermonde_y = polynomial_interpolation_matrix(nodes_in, normalized_coordinates[2])
    for v in 1:n_variables
      for i in 1:n_nodes
        temp_data[i, element, v] = (vandermonde_y*unstructured_data[i, :, element_ids[element], v])[1]
      end
      data_on_curve[element, v] = (vandermonde_x*temp_data[:, element, v])[]
    end
  end

  return arc_length, data_on_curve, nothing
end

# Convert 3d unstructured data to 1d data at given curve.
function unstructured_3d_to_1d_curve(original_nodes, unstructured_data, nvisnodes, curve, mesh, solver, cache)

  n_points_curve = size(curve)[2]
  n_nodes, _, _, n_elements, n_variables = size(unstructured_data)
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes)

  # Check if input is correct.
  min = original_nodes[:, 1, 1, 1, 1]
  max = max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, n_elements]
  @assert size(curve) == (3, n_points_curve) "Coordinates along curve must be 3xn dimensional."
  for element in 1:n_points_curve
    @assert (prod(vcat(curve[:, n_points_curve] .>= min, curve[:, n_points_curve]
            .<= max))) "Some coordinates from `curve` are outside of the domain.."
  end

  # Set nodes acording to the length of the curve.
  arc_length = calc_arc_length(curve)

  # Setup data structures.
  data_on_curve = Array{Float64}(undef, n_points_curve, n_variables)
  temp_data = Array{Float64}(undef, n_nodes, n_nodes+1, n_points_curve, n_variables)

  # For each coordinate find the corresponding element with its id.
  element_ids = get_elements_by_coordinates(curve, mesh, solver, cache)

  # Iterate over all found elements.
  for element in 1:n_points_curve

    min_coordinate = original_nodes[:, 1, 1, 1, element_ids[element]]
    max_coordinate = original_nodes[:, n_nodes, n_nodes, n_nodes, element_ids[element]]
    element_length = max_coordinate - min_coordinate

    normalized_coordinates = (curve[:, element] - min_coordinate)/element_length[1]*2 .-1

    # Interpolate to a single point in each element.
    vandermonde_x = polynomial_interpolation_matrix(nodes_in, normalized_coordinates[1])
    vandermonde_y = polynomial_interpolation_matrix(nodes_in, normalized_coordinates[2])
    vandermonde_z = polynomial_interpolation_matrix(nodes_in, normalized_coordinates[3])
    for v in 1:n_variables
      for i in 1:n_nodes
        for ii in 1:n_nodes
          temp_data[i, ii, element, v] = (vandermonde_z*unstructured_data[i, ii, :, element_ids[element], v])[1]
        end
        temp_data[i, n_nodes+1, element, v] = (vandermonde_y*temp_data[i, 1:n_nodes, element, v])[1]
      end
      data_on_curve[element, v] = (vandermonde_x*temp_data[:, n_nodes+1, element, v])[1]
    end
  end

  return arc_length, data_on_curve, nothing
end

# Convert 3d unstructured data to 1d slice and interpolate them.
function unstructured_3d_to_1d(original_nodes, unstructured_data, nvisnodes, slice, point)

  if slice === :x
    slice_dimension = 1
    other_dimensions = [2,3]
  elseif slice === :y
    slice_dimension = 2
    other_dimensions = [1,3]
  elseif slice === :z
    slice_dimension = 3
    other_dimensions = [1,2]
  else
    error("illegal dimension '$slice', supported dimensions are :x, :y and :z")
  end

  # Set up data structures to stroe new 1D data.
  @views new_unstructured_data = similar(unstructured_data[1, 1, ..])
  @views temp_unstructured_data = similar(unstructured_data[1, ..])
  @views new_nodes = similar(original_nodes[1, 1, 1,..])

  n_nodes_in, _, _, n_elements, n_variables = size(unstructured_data)
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Test if point lies in the domain.
  lower_limit = original_nodes[1, 1, 1, 1, 1]
  upper_limit = original_nodes[1, n_nodes_in, n_nodes_in, n_nodes_in, n_elements]

  @assert length(point) >= 3 "Point must be three-dimensional."
  if prod(point[other_dimensions] .< lower_limit) || prod(point[other_dimensions] .> upper_limit)
    error(string("Slice axis is outside of domain. ",
        " point[$other_dimensions]=$(point[other_dimensions]) must be between $lower_limit and $upper_limit"))
  end

  # Count the amount of new elements.
  new_id = 0

  # Permute dimensions so that the slice dimensions are always the in correct places for later use.
  if slice === :x
    original_nodes = permutedims(original_nodes, [1, 3, 4, 2, 5])
    unstructured_data = permutedims(unstructured_data, [2, 3, 1, 4, 5])
  elseif slice === :y
    original_nodes = permutedims(original_nodes, [1, 2, 4, 3, 5])
    unstructured_data = permutedims(unstructured_data, [1, 3, 2, 4, 5])
  end

  # Iterate over all elements to find the ones that lie on the slice axis.
  for element_id in 1:n_elements
    min_coordinate = original_nodes[:, 1, 1, 1, element_id]
    max_coordinate = original_nodes[:, n_nodes_in, n_nodes_in, n_nodes_in, element_id]
    element_length = max_coordinate - min_coordinate

    # Test if the element is on the slice axis. If not just continue with the next element.
    if !((prod(min_coordinate[other_dimensions] .<= point[other_dimensions]) &&
        prod(max_coordinate[other_dimensions] .> point[other_dimensions])) ||
        (point[other_dimensions] == upper_limit && prod(max_coordinate[other_dimensions] .== upper_limit)))

        continue
    end

    new_id += 1

    # Construct vandermonde matrix for interpolation of each 2D element to a 1D element.
    normalized_intercept =
          (point[other_dimensions] .- min_coordinate[other_dimensions]) /
          element_length[1] * 2 .- 1
    vandermonde_i = polynomial_interpolation_matrix(nodes_in, normalized_intercept[1])
    vandermonde_ii = polynomial_interpolation_matrix(nodes_in, normalized_intercept[2])

    # Interpolate to each node of new 1D element.
    for v in 1:n_variables
      for i in 1:n_nodes_in
        for ii in 1:n_nodes_in
          temp_unstructured_data[i, ii, new_id, v] = (vandermonde_ii*unstructured_data[ii, :, i, element_id, v])[1]
        end
        new_unstructured_data[i, new_id, v] = (vandermonde_i*temp_unstructured_data[i, :, new_id, v])[1]
      end
    end

    new_nodes[:, new_id] = original_nodes[slice_dimension, 1, 1, :, element_id]
  end

  return get_data_1d(reshape(new_nodes[:, 1:new_id], 1, n_nodes_in, new_id), new_unstructured_data[:, 1:new_id, :], nvisnodes)
end

# Interpolate unstructured DG data to structured data (cell-centered)
#
# This function takes DG data in an unstructured, Cartesian layout and converts it to a uniformely
# distributed 2D layout.
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function unstructured2structured(unstructured_data, normalized_coordinates,
                                 levels, resolution, nvisnodes_per_level)
  # Extract data shape information
  n_nodes_in, _, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Calculate interpolation vandermonde matrices for each level
  max_level = length(nvisnodes_per_level) - 1
  vandermonde_per_level = []
  for l in 0:max_level
    n_nodes_out = nvisnodes_per_level[l + 1]
    dx = 2 / n_nodes_out
    nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out))
    push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out))
  end

  # For each element, calculate index position at which to insert data in global data structure
  lower_left_index = element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level)

  # Create output data structure
  structured = [Matrix{Float64}(undef, resolution, resolution) for _ in 1:n_variables]

  # For each variable, interpolate element data and store to global data structure
  for v in 1:n_variables
    # Reshape data array for use in multiply_dimensionwise function
    reshaped_data = reshape(unstructured_data[:, :, :, v], 1, n_nodes_in, n_nodes_in, n_elements)

    for element_id in 1:n_elements
      # Extract level for convenience
      level = levels[element_id]

      # Determine target indices
      n_nodes_out = nvisnodes_per_level[level + 1]
      first = lower_left_index[:, element_id]
      last = first .+ (n_nodes_out - 1)

      # Interpolate data
      vandermonde = vandermonde_per_level[level + 1]
      structured[v][first[1]:last[1], first[2]:last[2]] .= (
          reshape(multiply_dimensionwise(vandermonde, reshaped_data[:, :, :, element_id]),
                  n_nodes_out, n_nodes_out))
    end
  end

  return structured
end

# For a given normalized element coordinate, return the index of its lower left
# contribution to the global data structure
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function element2index(normalized_coordinates, levels, resolution, nvisnodes_per_level)
  @assert size(normalized_coordinates, 1) == 2 "only works in 2D"

  n_elements = length(levels)

  # First, determine lower left coordinate for all cells
  dx = 2 / resolution
  ndim = 2
  lower_left_coordinate = Array{Float64}(undef, ndim, n_elements)
  for element_id in 1:n_elements
    nvisnodes = nvisnodes_per_level[levels[element_id] + 1]
    lower_left_coordinate[1, element_id] = (
        normalized_coordinates[1, element_id] - (nvisnodes - 1)/2 * dx)
    lower_left_coordinate[2, element_id] = (
        normalized_coordinates[2, element_id] - (nvisnodes - 1)/2 * dx)
  end

  # Then, convert coordinate to global index
  indices = coordinate2index(lower_left_coordinate, resolution)

  return indices
end

# Find 2D array index for a 2-tuple of normalized, cell-centered coordinates (i.e., in [-1,1])
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function coordinate2index(coordinate, resolution::Integer)
  # Calculate 1D normalized coordinates
  dx = 2/resolution
  mesh_coordinates = collect(range(-1 + dx/2, 1 - dx/2, length=resolution))

  # Find index
  id_x = searchsortedfirst.(Ref(mesh_coordinates), coordinate[1, :], lt=(x,y)->x .< y .- dx/2)
  id_y = searchsortedfirst.(Ref(mesh_coordinates), coordinate[2, :], lt=(x,y)->x .< y .- dx/2)
  return transpose(hcat(id_x, id_y))
end

# Calculate the vertices for each mesh cell such that it can be visualized as a closed box
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function calc_vertices(coordinates, levels, length_level_0)
  ndim = size(coordinates, 1)
  @assert ndim == 2 "only works in 2D"

  # Initialize output arrays
  n_elements = length(levels)
  n_points_per_element = 2^ndim+2
  x = Vector{Float64}(undef, n_points_per_element*n_elements)
  y = Vector{Float64}(undef, n_points_per_element*n_elements)

  # Calculate vertices for all coordinates at once
  for element_id in 1:n_elements
    length = length_level_0 / 2^levels[element_id]
    index = n_points_per_element*(element_id-1)
    x[index+1] = coordinates[1, element_id] - 1/2 * length
    x[index+2] = coordinates[1, element_id] + 1/2 * length
    x[index+3] = coordinates[1, element_id] + 1/2 * length
    x[index+4] = coordinates[1, element_id] - 1/2 * length
    x[index+5] = coordinates[1, element_id] - 1/2 * length
    x[index+6] = NaN

    y[index+1] = coordinates[2, element_id] - 1/2 * length
    y[index+2] = coordinates[2, element_id] - 1/2 * length
    y[index+3] = coordinates[2, element_id] + 1/2 * length
    y[index+4] = coordinates[2, element_id] + 1/2 * length
    y[index+5] = coordinates[2, element_id] - 1/2 * length
    y[index+6] = NaN
  end

  return x, y
end

# Calculate the vertices to plot each grid line for StructuredMesh
#
# Note: This is a low-level function that is not considered as part of Trixi's interface and may
#       thus be changed in future releases.
function calc_vertices(node_coordinates, mesh)
  @unpack cells_per_dimension = mesh
  @assert size(node_coordinates, 1) == 2 "only works in 2D"

  linear_indices = LinearIndices(size(mesh))

  # Initialize output arrays
  n_lines = sum(cells_per_dimension) + 2
  max_length = maximum(cells_per_dimension)
  n_nodes = size(node_coordinates, 2)

  # Create output as two matrices `x` and `y`, each holding the node locations for each of the `n_lines` grid lines
  # The # of rows in the matrices must be sufficient to store the longest dimension (`max_length`),
  # and for each the node locations without doubling the corner nodes (`n_nodes-1`), plus the final node (`+1`)
  # Rely on Plots.jl to ignore `NaN`s (i.e., they are not plotted) to handle shorter lines
  x = fill(NaN, max_length*(n_nodes-1)+1, n_lines)
  y = fill(NaN, max_length*(n_nodes-1)+1, n_lines)

  line_index = 1
  # Lines in x-direction
  # Bottom boundary
  i = 1
  for cell_x in axes(mesh, 1)
    for node in 1:(n_nodes-1)
      x[i, line_index] = node_coordinates[1, node, 1, linear_indices[cell_x, 1]]
      y[i, line_index] = node_coordinates[2, node, 1, linear_indices[cell_x, 1]]

      i += 1
    end
  end
  # Last point on bottom boundary
  x[i, line_index] = node_coordinates[1, end, 1, linear_indices[end, 1]]
  y[i, line_index] = node_coordinates[2, end, 1, linear_indices[end, 1]]

  # Other lines in x-direction
  line_index += 1
  for cell_y in axes(mesh, 2)
    i = 1
    for cell_x in axes(mesh, 1)
      for node in 1:(n_nodes-1)
        x[i, line_index] = node_coordinates[1, node, end, linear_indices[cell_x, cell_y]]
        y[i, line_index] = node_coordinates[2, node, end, linear_indices[cell_x, cell_y]]

        i += 1
      end
    end
    # Last point on line
    x[i, line_index] = node_coordinates[1, end, end, linear_indices[end, cell_y]]
    y[i, line_index] = node_coordinates[2, end, end, linear_indices[end, cell_y]]

    line_index += 1
  end

  # Lines in y-direction
  # Left boundary
  i = 1
  for cell_y in axes(mesh, 2)
    for node in 1:(n_nodes-1)
      x[i, line_index] = node_coordinates[1, 1, node, linear_indices[1, cell_y]]
      y[i, line_index] = node_coordinates[2, 1, node, linear_indices[1, cell_y]]

      i += 1
    end
  end
  # Last point on left boundary
  x[i, line_index] = node_coordinates[1, 1, end, linear_indices[1, end]]
  y[i, line_index] = node_coordinates[2, 1, end, linear_indices[1, end]]

  # Other lines in y-direction
  line_index +=1
  for cell_x in axes(mesh, 1)
    i = 1
    for cell_y in axes(mesh, 2)
      for node in 1:(n_nodes-1)
        x[i, line_index] = node_coordinates[1, end, node, linear_indices[cell_x, cell_y]]
        y[i, line_index] = node_coordinates[2, end, node, linear_indices[cell_x, cell_y]]

        i += 1
      end
    end
    # Last point on line
    x[i, line_index] = node_coordinates[1, end, end, linear_indices[cell_x, end]]
    y[i, line_index] = node_coordinates[2, end, end, linear_indices[cell_x, end]]

    line_index += 1
  end

  return x, y
end

# Convert `slice` to orientations (1 -> `x`, 2 -> `y`, 3 -> `z`) for the two axes in a 2D plot
function _get_orientations(mesh, slice)
  if ndims(mesh) == 2 || (ndims(mesh) == 3 && slice === :xy)
    orientation_x = 1
    orientation_y = 2
  elseif ndims(mesh) == 3 && slice === :xz
    orientation_x = 1
    orientation_y = 3
  elseif ndims(mesh) == 3 && slice === :yz
    orientation_x = 2
    orientation_y = 3
  else
    orientation_x = 0
    orientation_y = 0
  end
  return orientation_x, orientation_y
end

# Convert `orientation` into a guide label (see also `_get_orientations`)
function _get_guide(orientation::Integer)
  if orientation == 1
    return "\$x\$"
  elseif orientation == 2
    return "\$y\$"
  elseif orientation == 3
    return "\$z\$"
  else
    return ""
  end
end

#   plotting_interpolation_matrix(dg; kwargs...)
#
# Interpolation matrix which maps discretization nodes to a set of plotting nodes.
# Defaults to the identity matrix of size `length(solver.basis.nodes)`, and interpolates
# to equispaced nodes for DGSEM (set by kwarg `nvisnodes` in the plotting function).
#
# Example:
# ```julia
# A = plotting_interpolation_matrix(dg)
# A * dg.basis.nodes # => vector of nodes at which to plot the solution
# ```
#
# Note: we cannot use UniformScaling to define the interpolation matrix since we use it with `kron`
# to define a multi-dimensional interpolation matrix later.
plotting_interpolation_matrix(dg; kwargs...) = I(length(dg.basis.nodes))

function face_plotting_interpolation_matrix(dg::DGSEM; nvisnodes=2*length(dg.basis.nodes))
  return polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes))
end

function plotting_interpolation_matrix(dg::DGSEM; nvisnodes=2*length(dg.basis.nodes))
  Vp1D = polynomial_interpolation_matrix(dg.basis.nodes, LinRange(-1, 1, nvisnodes))
  # For quadrilateral elements, interpolation to plotting nodes involves applying a 1D interpolation
  # operator to each line of nodes. This is equivalent to multiplying the vector containing all node
  # node coordinates on an element by a Kronecker product of the 1D interpolation operator (e.g., a
  # multi-dimensional interpolation operator).
  return kron(Vp1D, Vp1D)
end

function reference_node_coordinates_2d(dg::DGSEM)
  @unpack nodes = dg.basis
  r = vec([nodes[i] for i in eachnode(dg), j in eachnode(dg)])
  s = vec([nodes[j] for i in eachnode(dg), j in eachnode(dg)])
  return r, s
end

# Given a reference plotting triangulation, this function generates a plotting triangulation for
# the entire global mesh. The output can be plotted using `Makie.mesh`.
function global_plotting_triangulation_makie(pds::PlotDataSeries{<:UnstructuredPlotData2D};
                                             set_z_coordinate_zero = false)
  @unpack variable_id = pds
  pd = pds.plot_data
  @unpack x, y, data, t = pd

  makie_triangles = Makie.to_triangles(t)

  # trimesh[i] holds GeometryBasics.Mesh containing plotting information on the ith element.
  # Note: Float32 is required by GeometryBasics
  num_plotting_nodes, num_elements = size(x)
  trimesh = Vector{GeometryBasics.Mesh{3, Float32}}(undef, num_elements)
  coordinates = zeros(Float32, num_plotting_nodes, 3)
  for element in Base.OneTo(num_elements)
    for i in Base.OneTo(num_plotting_nodes)
      coordinates[i, 1] = x[i, element]
      coordinates[i, 2] = y[i, element]
      if set_z_coordinate_zero == false
        coordinates[i, 3] = data[i, element][variable_id]
      end
    end
    trimesh[element] = GeometryBasics.normal_mesh(Makie.to_vertices(coordinates), makie_triangles)
  end
  plotting_mesh = merge([trimesh...]) # merge meshes on each element into one large mesh
  return plotting_mesh
end

# Returns a list of `Makie.Point`s which can be used to plot the mesh, or a solution "wireframe"
# (e.g., a plot of the mesh lines but with the z-coordinate equal to the value of the solution).
function mesh_plotting_wireframe(pds::PlotDataSeries{<:UnstructuredPlotData2D};
                                 set_z_coordinate_zero = false)
  @unpack variable_id = pds
  pd = pds.plot_data
  @unpack x_face, y_face, face_data = pd

  if set_z_coordinate_zero
    sol_f = zeros(eltype(first(face_data)), size(face_data)) # plot 2d surface by setting z coordinate to zero
  else
    sol_f = StructArrays.component(face_data, variable_id)
  end

  # This line separates solution lines on each edge by NaNs to ensure that they are rendered
  # separately. The coordinates `xf`, `yf` and the solution `sol_f`` are assumed to be a matrix
  # whose columns correspond to different elements. We add NaN separators by appending a row of
  # NaNs to this matrix. We also flatten (e.g., apply `vec` to) the result, as this speeds up
  # plotting.
  xyz_wireframe = GeometryBasics.Point.(map(x->vec(vcat(x, fill(NaN, 1, size(x, 2)))), (x_face, y_face, sol_f))...)

  return xyz_wireframe
end

end # @muladd


text:
import TimeSeries.TimeArray

abstract type AbstractQueryOpt <: AbstractDict{Symbol,Any} end

Base.length(::T) where {T<:AbstractQueryOpt} = fieldcount(T)
Base.length(::Type{T}) where {T<:AbstractQueryOpt} = fieldcount(T)

function Base.iterate(aqo::T, state = 1) where {T<:AbstractQueryOpt}
  (state > length(T)) && return nothing
  (fieldname(T, state) => getfield(aqo, state), state + 1)
end

"""
    struct YahooOpt <: AbstractQueryOpt
      period1  # the start time
      period2  # the end time
      interval # "1d", "1wk" or "1mo"
      events   # currently only `:history` supported
    end

The Yahoo Finance HTTP API query object.

# Examples

```jl-repl
julia> t = Dates.now()
2020-08-09T01:38:04.735

julia> YahooOpt(period1 = t - Year(2), period2 = t)
YahooOpt{DateTime} with 4 entries:
  :period1  => 1533778685
  :period2  => 1596937085
  :interval => "1d"
  :events   => :history
```
"""
struct YahooOpt <: AbstractQueryOpt
  period1::DateTime
  period2::DateTime
  interval::String
  events::Symbol

  YahooOpt(; period1::DateTime = DateTime(1971, 2, 8),
          period2::DateTime = Dates.now(),
          interval::String  = "1d",
          events::Symbol    = :history) =
    new(period1, period2, interval, events)
end

function Base.iterate(opt::YahooOpt, state = 1)
  (state > length(YahooOpt)) && return nothing
  k = fieldname(YahooOpt, state)
  v = getfield(opt, state)
  v′ = (k ∈ (:period1, :period2)) ? round(Int, datetime2unix(v)) : v
  (k => v′, state + 1)
end

"""
    yahoo(symbol::AbstractString, opt::YahooOpt = YahooOpt())::TimeArray
    yahoo(symbol::Symbol, opt::YahooOpt = YahooOpt())::TimeArray

This is a wrapper for downloading historical stock prices from Yahoo Finance.

The yahoo method takes a stock name in the form of a string and returns a
`TimeSeries.TimeArray` corresponding to the Yahoo Finance ticker.
With no argument, the default historical time series is the S&P 500.

# Examples

```julia
AAPL = yahoo(:AAPL)
SPX = yahoo("^GSPC")
NQ = yahoo("^IXIC")
```

```jl-repl
julia> start = DateTime(2018, 1, 1)
2018-01-01T00:00:00

julia> yahoo(:AAPL, YahooOpt(period1 = start))
655×6 TimeArray{Float64,2,Date,Array{Float64,2}} 2018-01-02 to 2020-08-07
...
```

# References

  https://finance.yahoo.com

# See Also

  fred() which accesses the St. Louis Federal Reserve financial and economic data sets.
"""
function yahoo(sym::AbstractString = "^GSPC", opt::YahooOpt = YahooOpt())
    host = rand(["query1", "query2"])
    url  = "https://$host.finance.yahoo.com/v7/finance/download/$sym"
    res  = HTTP.get(url, query = opt)
    @assert res.status == 200
    csv = CSV.File(res.body, missingstrings = ["null"])
    sch = TimeSeries.Tables.schema(csv)
    TimeArray(csv, timestamp = first(sch.names)) |> cleanup_colname!
end

yahoo(s::Symbol, opt::YahooOpt = YahooOpt()) = yahoo(string(s), opt)

"""
Description

    The fred() method is a wrapper to download financial and economic time series data from the St. Louis Federal Reserve (FRED).

Usage

    DGS = fred("DGS10")
    CPI = fred()

Method Signature(s)

    fred(data::String="CPIAUCNS")

Details

    The fred() method takes a string argument that corresponds to a series code from the St. Louis Federal
    Reserve (FRED) database. It returns the data in the TimeSeries.TimeArray data structure.  When no argument
    is provided, the default data set is the Consumer Price Index for All Urban Consumers: All Items (CPIAUCNS).

References

    https://research.stlouisfed.org/fred2

See Also

    yahoo() which is a wrapper from downloading financial time series for stocks from Yahoo Finance.
"""
function fred(data::String="CPIAUCNS")
    url = "http://research.stlouisfed.org/fred2/series/$data/downloaddata/$data.csv"
    res = HTTP.get(url)
    @assert res.status == 200
    csv = CSV.File(res.body)
    sch = TimeSeries.Tables.schema(csv)
    TimeArray(csv, timestamp = first(sch.names)) |> cleanup_colname!
end


text:
#
println("Dep) Test of the AutoIonization and PlasmaShift modules with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-AutoIonization-Plasma.sum")
setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")   ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)
wa   = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(10.),
                          initialConfigs=[Configuration("1s 2s^2 2p^6")],
                          finalConfigs  =[Configuration("1s^2 2s^2 2p^4"), Configuration("1s^2 2s 2p^5"), Configuration("1s^2 2p^6")],
                          processSettings = PlasmaShift.AugerSettings(PlasmaShift.DebyeHueckel(), 0.25, 0., 0, true, LineSelection() ) )

wb = perform(wa)
setDefaults("print summary: close", "")