text:
abstract type HadamardInt <: FracIntAlg end

"""
    fracint(f, α, a, b, mu, h, HadamardMat())

Compute Hadamard fractional integral.

### References

```tex
@article{Xu2019SpectralCM,
  title={Spectral Collocation Method for Fractional Differential/Integral Equations with Generalized Fractional Operator},
  author={Qinwu Xu and Zhoushun Zheng},
  journal={International Journal of Differential Equations},
  year={2019}
}
```
"""
struct HadamardMat <: HadamardInt end

function fracint(g, α, a, b, μ, h, ::HadamardMat)
    N = round(Int, (b-a)/h)
    w(x)=x.^μ
    invz(x)=exp.(x)
    (rx, r, D) = GFracMat(log(a), log(b), N, α, w, invz)
    Ga=a^μ.*g.(a)/gamma(1-α)*rx.^(-μ).*(log.(rx./a)).^(-α)
    G=g.(rx)
    H=Ga+D*G
    return rx, H
end

function GFracMat(za, zb, N, α, w, invz)
    r0=((JacobiGL(0, 0, N).+1)./2)
    col_pts=za.+r0*(zb-za)
    col_pt = deepcopy(col_pts)
    if α>0
    FM=Jac_Frac_Diff_Shift(col_pts,0,0,N,α,za,zb,1)
    else
        FM=(2/(zb-za))^α*Jac_Frac_Int(r0*2-1, 0, 0, N, -α, 1)
        FM[1, :] .= 0
    end

    r1 = r0
    interp_pts = r1.*2 .-1
    Vand = Jac_Va(interp_pts, 0, 0, N, 1)

    rx=invz.(col_pt)
    WL = diagm(1 ./w.(rx))
    WL[1, 1]=0
    WR = diagm(w.(rx))
    M=WL*(FM/Vand)*WR
    return rx, col_pt, M
end

function JacobiGL(α,β,N)
    x = zeros(N+1, 1)
    if N==1
        x[1] = -1.0
        x[2] = 1.0
        return x
    end

    (xint, _) = JacobiGQ(α+1,β+1,N-2)
    x = [-1 xint' 1]'
    return x
end

function JacobiGQ(α, β, N)
    if N==0
        x[1] = -(α-β)/(α+β+2)
        w[1] = 2
        return x, w
    end

    J = zeros(N+1)
    h1 = 2*collect(0:N).+α.+β
    J = diagm(-1/2*(α^2-β^2)./(h1.+2)./h1) +
        diagm(1 => 2 ./(h1[1:N].+2).*sqrt.(collect(1:N).*(collect(1:N).+α.+β).*(collect(1:N).+α).*(collect(1:N).+β)./(h1[1:N].+1)./(h1[1:N].+3)))
    if α+β < 10*eps()
        J[1, 1] = 0.0
    end
    J = J + J'

    (D, V) = eigen(J)
    x = diagm(D)

    w = V[1, :].^2*2^(α+β+1)/(α+β+1)*gamma(α+1)*gamma(β+1)/gamma(α+β+1)
    return diag(x), w
end

function Jac_Frac_Diff(x, a, b,N, α, output_type)
    if size(x, 2) !== 1
        x=x'
    end
    L=length(x)

    k = ceil(Int64, α)
    d = gamma.(collect(0:N).+(k+a+b+1))./(2^k*gamma.(collect(0:N).+(a+b+1)))

    if output_type==0
        if N < k
            JFD = zeros(size(x))
            return JFD
        end
        JFD = d[N+1]*Jac_Frac_Int(x,a+k,b+k,N-k,k-α,0)
    else
        JFD = zeros(L, N+1)
        JFD[:, k+1:N+1] = Jac_Frac_Int(x, a+k, b+k, N-k, k-α, 1)
        for j = k+1:N+1
            JFD[:, j] = d[j]*JFD[:, j]
        end
    end
    return JFD
end

function Jac_Frac_Diff_Shift(x,a,b,N,α,min,Max,output_type)
    L=Max-min
    x = 2 .*(x.-min)./L.-1
    FJ = (2/L)^α*Jac_Frac_Diff(x, a, b, N, α, output_type)
    return FJ
end

function Jac_Frac_Int(x, a, b, N::Int64, α, output_type)
    if size(x,2) !== 1
        x=x'
    end
    L=length(x)
    J = zeros(L, Int64(N+1))
    @. J[:, 1] = (x+1)^α/gamma(α+1)
    if N==0
        return J
    end
    @. J[:, 2] = (a+b+2)/2*(x.*(x+1).^α/gamma(α+1)-α*(x+1).^(α+1)/gamma(α+2))+(a-b)/2*J[:, 1]
    if N==1
        return J
    end
    Jac_1 = Jac_Va(-1, a, b, N, 1)
    Aj(a, b, j) = @. (2*j+a+b+1)*(2*j+a+b+2)/(2*(j+1)*(j+a+b+1))
    Bj(a, b, j) = @. (b^2-a^2)*(2*j+a+b+1)/(2*(j+1)*(j+a+b+1)*(2*j+a+b))
    Cj(a, b, j) = @. (j+a)*(j+b)*(2*j+a+b+2)/((j+1)*(j+a+b+1)*(2*j+a+b))
    Ahatj(a, b, j) = @. -2*(j+a)*(j+b)/((j+a+b)*(2*j+a+b)*(2*j+a+b+1))
    Bhatj(a, b, j) = @. 2*(a-b)/((2*j+a+b)*(2*j+a+b+2))
    Chatj(a, b, j) = @. 2*(j+a+b+1)/((2*j+a+b+1)*(2*j+a+b+2))

    for n=2:N
        @. J[:, n+1] = (Aj(a,b,n-1)*x-Bj(a,b,n-1)-α*Aj(a,b,n-1)*Bhatj(a,b,n-1))*J[:, n]-(Cj(a, b, n-1) + α*Aj(a, b, n-1)*Ahatj(a, b, n-1)).*J[:, n-1] + Aj(a, b, n-1)*(Ahatj(a, b, n-1)*Jac_1[n-1]+Bhatj(a,b,n-1)*Jac_1[n]+Chatj(a,b,n-1)*Jac_1[n+1])*(x+1).^α/gamma(α)
        @. J[:, n+1] =J[:, n+1]/(1+α*Aj(a,b,n-1)*Chatj(a,b,n-1))
    end
    if output_type==0
        J=J[:, N+1]
    end
    return J
end

function Jac_Va(x, a, b, N, output_type)
    if size(x,2) !== 1
        x=x'
    end
    L=length(x)
    J = ones(L, N+1)
    if N==0
        J=J[:, 1]
        return J
    end
    @. J[:, 2]=(1+a)+(a+b+2)*(x/2 -1/2)
    if N==1
        if output_type==0
            J=J[:, 2]
            return J
        else
            J=J[:, 1:2]
            return J
        end
    end
    for n=2:N
        @. J[:, n+1] = ((2*n+a+b-1)*((2*n+a+b)*(2*n+a+b-2)*x+a^2-b^2).*J[:, n]-2*(n+a-1)*(n+b-1)*(2*n+a+b)*J[:, n-1])/(2*n*(n+a+b)*(2*n+a+b-2))
    end
    if output_type==0
        J=J[:, N+1]
    end
    return J
end


text:
# ------------------------------------------------------------------------------
# This would be a convenient system to manage cacheing of dual numbers.
# Unfortunately, at least in Julia 1.5, reinterpreting an array of
# Dual numbers leads to a large performance penalty. Until this is
# resolved, a non-reinterpret-based solution is necessary

# Note that ForwardDiff2 is currently under construction, so that
# could be interesting
# ------------------------------------------------------------------------------

using DiscreteElasticRods
using ForwardDiff
abstract type DualType{V,N} end

struct dual_rod{V,N} <: DualType{V,N}
    num::DER.basic_rod{V}
    dual::DER.basic_rod{ForwardDiff.Dual{nothing,V,N}}
end

struct dual_delta{V,N} <: DualType{V,N}
    num::DER.rod_delta{V}
    dual::DER.rod_delta{ForwardDiff.Dual{nothing,V,N}}
end

struct dual_vector{V,N} <: DualType{V,N}
    num::Vector{V}
    dual::Vector{ForwardDiff.Dual{nothing,V,N}}
end

struct dual_cache{V,N} <: DualType{V,N}
    num::Tuple{Vararg{Array{V}}}
    dual::Tuple{Vararg{Array{ForwardDiff.Dual{nothing,V,N}}}}
end

# Constructors
dual_rod(V::Type, Ns::Int, Nd::Int) = dual_rod(DER.allocate_rod(V, Ns), DER.allocate_rod(ForwardDiff.Dual{nothing,V,Nd}, Ns))
dual_delta(V::Type, Ns::Int, Nd::Int) = dual_delta(DER.allocate_delta(V, Ns), DER.allocate_delta(ForwardDiff.Dual{nothing,V,Nd}, Ns))
dual_vector(V::Type, N::Int, Nd::Int) = dual_vector(Vector{V}(undef, N), Vector{ForwardDiff.Dual{nothing,V,Nd}}(undef, N))

# Get the num or dual component
function resolve_type(dual::DualType{V,N}, ::Type{V}) where {V, N}
    return dual.num
end

function resolve_type(dual::DualType{V,N}, ::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
    return reinterpret(ForwardDiff.Dual{T,V,N}, dual.dual)
end

function resolve_type(dual::dual_cache{V,N}, ::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
    return reinterpret_dual_tuple(ForwardDiff.Dual{T,V,N}, dual.dual)
end

# Extend reinterpret (renamed for tuples, to avoid type piracy)
function Base.reinterpret(T::Type, r::DER.basic_rod)
    DER.basic_rod(reinterpret(T, r.x), reinterpret(T, r.d))
end

function Base.reinterpret(T::Type, Δr::DER.rod_delta)
    DER.basic_rod(reinterpret(T, Δr.Δx), reinterpret(T, Δr.Δθ))
end

function reinterpret_dual_tuple(T::Type{ForwardDiff.Dual{T1, V, N}}, a::Tuple{Vararg{Array{ForwardDiff.Dual{T2,V,N}}}}) where {T1, T2, V, N}
    map(x->reinterpret(T, x), a)
end

# Finally, allow any cache structure in my code to be a dual
function allocate_dual_cache(f::Function, V::Type, Ns::Int, Nd::Int)
    dual_cache{V,Nd}(DER.allocate_cache(f, V, Ns), DER.allocate_cache(f, ForwardDiff.Dual{nothing, V, Nd}, Ns))
end

# Minimal functionality
Ns = 50
Nd = 3
NType = Float64
struct MyTag end
DType = ForwardDiff.Dual{MyTag, NType, Nd}

r = dual_rod(NType, Ns, Nd)
resolve_type(r, NType)
resolve_type(r, DType)

r = dual_delta(NType, Ns, Nd)
resolve_type(r, NType)
resolve_type(r, DType)

v = dual_vector(NType, Ns, Nd)
resolve_type(v, NType)
resolve_type(v, DType)

c = allocate_dual_cache(DER.elastic_forces!, Float64, 50, 3)
resolve_type(c, NType)
resolve_type(c, DType)

# How hard does reinterpret hurt timings?
struct MyTag2 end
DType2 = ForwardDiff.Dual{MyTag2, NType, Nd}

Ns = 500
a0 = rand(Ns, 3)
b0 = rand(Ns, 3)
c0 = similar(a0)
@btime DER.cross3!($c0, $a0, $b0)

a1 = DType.(rand(Ns,3))
b1 = DType.(rand(Ns,3))
c1 = similar(a)
@btime DER.cross3!($c1, $a1, $b1)

a2 = reinterpret(DType2, a1)
b2 = reinterpret(DType2, b1)
c2 = reinterpret(DType2, c1)
@btime DER.cross3!($c2, $a2, $b2)

a0 = rand(Ns, 3)
b0 = rand(Ns, 3)
c0 = rand(Ns)
@btime DER.dot3!($c0, $a0, $b0)

a1 = DType.(rand(Ns,3))
b1 = DType.(rand(Ns,3))
c1 = DType.(rand(Ns))
@btime DER.dot3!($c1, $a1, $b1)

a2 = reinterpret(DType2, a1)
b2 = reinterpret(DType2, b1)
c2 = reinterpret(DType2, c1)
@btime DER.dot3!($c2, $a2, $b2)

ext: jl


text:
using CRRao
using Test

@testset "CRRao.jl" begin
    # Write your tests here.
end


text:
module TestMissingArray
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, leveltype
using Missings

const ≅ = isequal

@testset "conversion ordered=$ordered" for ordered in (false, true)
    @testset "conversion reftype=$R" for R in (CategoricalArrays.DefaultRefType, UInt8, UInt, Int8, Int)
        @testset "conversion of $(typeof(a))" for a in (["b", "a", "b"], Union{String, Missing}["b", "a", "b"])
            @testset "Vector with no missing values" begin
                x = CategoricalVector{Union{String, Missing}, R}(a, ordered=ordered)

                @test x == a
                @test leveltype(typeof(x)) === String
                @test leveltype(x) === String
                @test eltype(x) === Union{CategoricalValue{String, R}, Missing}
                @test isordered(x) === ordered
                @test levels(x) == sort(unique(a))
                @test unique(x) == unique(a)
                @test size(x) === (3,)
                @test length(x) === 3

                @test convert(CategoricalArray, x) === x
                @test convert(CategoricalArray{Union{String, Missing}}, x) === x
                @test convert(CategoricalArray{Union{String, Missing}, 1}, x) === x
                @test convert(CategoricalArray{Union{String, Missing}, 1, R}, x) === x
                @test convert(CategoricalArray{Union{String, Missing}, 1, DefaultRefType}, x) == x
                @test convert(CategoricalArray{Union{String, Missing}, 1, UInt8}, x) == x

                @test convert(CategoricalVector, x) === x
                @test convert(CategoricalVector{Union{String, Missing}}, x) === x
                @test convert(CategoricalVector{Union{String, Missing}, R}, x) === x
                @test convert(CategoricalVector{Union{String, Missing}, DefaultRefType}, x) == x
                @test convert(CategoricalVector{Union{String, Missing}, UInt8}, x) == x

                @test convert(CategoricalArray, a) == x
                @test convert(CategoricalArray{Union{String, Missing}}, a) == x
                @test convert(CategoricalArray{Union{String, Missing}, 1}, a) == x
                @test convert(CategoricalArray{Union{String, Missing}, 1, R}, a) == x
                @test convert(CategoricalArray{Union{String, Missing}, 1, DefaultRefType}, a) == x
                @test convert(CategoricalArray{Union{String, Missing}, 1, UInt8}, a) == x

                @test convert(CategoricalVector, a) == x
                @test convert(CategoricalVector{Union{String, Missing}}, a) == x
                @test convert(CategoricalVector{Union{String, Missing}, R}, a) == x
                @test convert(CategoricalVector{Union{String, Missing}, DefaultRefType}, a) == x
                @test convert(CategoricalVector{Union{String, Missing}, UInt8}, a) == x

                @test CategoricalArray{Union{String, Missing}}(a, ordered=ordered) == x
                @test CategoricalArray{Union{String, Missing}, 1}(a, ordered=ordered) == x
                @test CategoricalArray{Union{String, Missing}, 1, R}(a, ordered=ordered) == x
                @test CategoricalArray{Union{String, Missing}, 1, DefaultRefType}(a, ordered=ordered) == x
                @test CategoricalArray{Union{String, Missing}, 1, UInt8}(a, ordered=ordered) == x

                @test CategoricalVector(a, ordered=ordered) == x
                @test CategoricalVector{Union{String, Missing}}(a, ordered=ordered) == x
                @test CategoricalVector{Union{String, Missing}, R}(a, ordered=ordered) == x
                @test CategoricalVector{Union{String, Missing}, DefaultRefType}(a, ordered=ordered) == x
                @test CategoricalVector{Union{String, Missing}, UInt8}(a, ordered=ordered) == x

                @testset "categorical($(typeof(y)), compress=$comp) R1=$R1 R2=$R2" for (y, R1, R2, comp) in
                    ((a, DefaultRefType, UInt8, true),
                     (a, DefaultRefType, DefaultRefType, false),
                     (x, R, UInt8, true),
                     (x, R, R, false))

                    x2 = @inferred categorical(y, ordered=ordered)
                    @test leveltype(x2) === String
                    @test nonmissingtype(eltype(x2)) === CategoricalValue{String, R1}
                    @test x2 == y
                    if eltype(y) >: Missing
                        @test isa(x2, CategoricalVector{Union{String, Missing}, R1})
                    else
                        @test isa(x2, CategoricalVector{String, R1})
                    end
                    @test isordered(x2) === ordered

                    x2 = categorical(y, compress=comp, ordered=ordered)
                    @test x2 == y
                    @test leveltype(x2) === String
                    @test nonmissingtype(eltype(x2)) === CategoricalValue{String, R2}
                    if eltype(y) >: Missing
                        @test isa(x2, CategoricalVector{Union{String, Missing}, R2})
                    else
                        @test isa(x2, CategoricalVector{String, R2})
                    end
                    @test isordered(x2) === ordered
                end

                x2 = compress(x)
                @test x2 == x
                @test isa(x2, CategoricalVector{Union{String, Missing}, UInt8})
                @test isordered(x2) === isordered(x)
                @test levels(x2) == levels(x)

                x2 = copy(x)
                @test x2 == x
                @test isordered(x2) === isordered(x)
                @test typeof(x2) === typeof(x)
                @test levels(x2) == levels(x)

                if !isordered(x)
                    @test ordered!(x, true) === x
                    @test isordered(x) === true
                end
                @test x[1] > x[2]
                @test x[3] > x[2]

                @test ordered!(x, false) === x
                @test isordered(x) === false
                @test_throws Exception x[1] > x[2]
                @test_throws Exception x[3] > x[2]

                @test x[1] === x.pool.valindex[2]
                @test x[2] === x.pool.valindex[1]
                @test x[3] === x.pool.valindex[2]
                @test_throws BoundsError x[4]

                x2 = x[:]
                @test typeof(x2) === typeof(x)
                @test x2 == x
                @test x2 !== x
                @test levels(x2) == levels(x)
                @test levels(x2) !== levels(x)
                @test isordered(x2) == isordered(x)

                x2 = x[2:3]
                @test typeof(x2) === typeof(x)
                @test x2 == ["a", "b"]
                @test levels(x2) == levels(x)
                @test levels(x2) !== levels(x)
                @test isordered(x2) == isordered(x)

                x2 = x[1:1]
                @test typeof(x2) === typeof(x)
                @test x2 == ["b"]
                @test levels(x2) == levels(x)
                @test levels(x2) !== levels(x)
                @test isordered(x2) == isordered(x)

                x2 = x[2:1]
                @test typeof(x2) === typeof(x)
                @test isempty(x2)
                @test levels(x2) == levels(x)
                @test levels(x2) !== levels(x)
                @test isordered(x2) == isordered(x)

                x[1] = x[2]
                @test x[1] === x.pool.valindex[1]
                @test x[2] === x.pool.valindex[1]
                @test x[3] === x.pool.valindex[2]

                x[3] = "c"
                @test x[1] === x.pool.valindex[1]
                @test x[2] === x.pool.valindex[1]
                @test x[3] === x.pool.valindex[3]
                @test levels(x) == ["a", "b", "c"]

                x[2:3] .= "b"
                @test x[1] === x.pool.valindex[1]
                @test x[2] === x.pool.valindex[2]
                @test x[3] === x.pool.valindex[2]
                @test levels(x) == ["a", "b", "c"]

                @test droplevels!(x) === x
                @test levels(x) == ["a", "b"]
                @test x[1] === x.pool.valindex[1]
                @test x[2] === x.pool.valindex[2]
                @test x[3] === x.pool.valindex[2]
                @test levels(x) == ["a", "b"]

                @test levels!(x, ["b", "a"]) === x
                @test levels(x) == ["b", "a"]
                @test x[1] === x.pool.valindex[2]
                @test x[2] === x.pool.valindex[1]
                @test x[3] === x.pool.valindex[1]
                @test levels(x) == ["b", "a"]

                @test_throws ArgumentError levels!(x, ["a"])
                @test_throws ArgumentError levels!(x, ["e", "b"])
                @test_throws ArgumentError levels!(x, ["e", "a", "b", "a"])

                @test levels!(x, ["e", "a", "b"]) === x
                @test levels(x) == ["e", "a", "b"]
                @test x[1] === x.pool.valindex[2]
                @test x[2] === x.pool.valindex[3]
                @test x[3] === x.pool.valindex[3]
                @test levels(x) == ["e", "a", "b"]

                x[1] = "c"
                @test x[1] === x.pool.valindex[4]
                @test x[2] === x.pool.valindex[3]
                @test x[3] === x.pool.valindex[3]
                @test levels(x) == ["e", "a", "b", "c"]

                @test_throws ArgumentError levels!(x, ["e", "c"])
                @test levels!(x, ["e", "c"], allowmissing=true) === x
                @test levels(x) == ["e", "c"]
                @test x[1] === x.pool.valindex[2]
                @test x[2] === missing
                @test x[3] === missing
                @test levels(x) == ["e", "c"]

                push!(x, "e")
                @test length(x) == 4
                @test x ≅ ["c", missing, missing, "e"]
                @test levels(x) == ["e", "c"]

                push!(x, "zz")
                @test length(x) == 5
                @test x ≅ ["c", missing, missing, "e", "zz"]
                @test levels(x) == ["e", "c", "zz"]

                push!(x, x[1])
                @test length(x) == 6
                @test x ≅ ["c", missing, missing, "e", "zz", "c"]
                @test levels(x) == ["e", "c", "zz"]

                push!(x, missing)
                @test length(x) == 7
                @test x ≅ ["c", missing, missing, "e", "zz", "c", missing]
                @test levels(x) == ["e", "c", "zz"]

                empty!(x)
                @test isordered(x) === false
                @test length(x) == 0
                @test levels(x) == ["e", "c", "zz"]
            end
        end

        @testset "Vector with missing values" begin
            a = ["a", "b", missing]
            x = CategoricalVector{Union{String, Missing}, R}(a, ordered=ordered)

            @test x ≅ a
            @test levels(x) == filter(x->!ismissing(x), unique(a))
            @test unique(x) ≅ unique(a)
            @test size(x) === (3,)
            @test length(x) === 3

            @test convert(CategoricalArray, x) === x
            @test convert(CategoricalArray{Union{String, Missing}}, x) === x
            @test convert(CategoricalArray{Union{String, Missing}, 1}, x) === x
            @test convert(CategoricalArray{Union{String, Missing}, 1, R}, x) === x
            @test convert(CategoricalArray{Union{String, Missing}, 1, DefaultRefType}, x) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}, 1, UInt8}, x) ≅ x

            @test convert(CategoricalVector, x) === x
            @test convert(CategoricalVector{Union{String, Missing}}, x) === x
            @test convert(CategoricalVector{Union{String, Missing}, R}, x) === x
            @test convert(CategoricalVector{Union{String, Missing}, DefaultRefType}, x) ≅ x
            @test convert(CategoricalVector{Union{String, Missing}, UInt8}, x) ≅ x

            @test convert(CategoricalArray, a) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}}, a) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}, 1}, a) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}, 1, R}, a) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}, 1, DefaultRefType}, a) ≅ x
            @test convert(CategoricalArray{Union{String, Missing}, 1, UInt8}, a) ≅ x

            @test convert(CategoricalVector, a) ≅ x
            @test convert(CategoricalVector{Union{String, Missing}}, a) ≅ x
            @test convert(CategoricalVector{Union{String, Missing}, R}, a) ≅ x
            @test convert(CategoricalVector{Union{String, Missing}, DefaultRefType}, a) ≅ x
            @test convert(CategoricalVector{Union{String, Missing}, UInt8}, a) ≅ x

            @test CategoricalArray{Union{String, Missing}}(a, ordered=ordered) ≅ x
            @test CategoricalArray{Union{String, Missing}, 1}(a, ordered=ordered) ≅ x
            @test CategoricalArray{Union{String, Missing}, 1, R}(a, ordered=ordered) ≅ x
            @test CategoricalArray{Union{String, Missing}, 1, DefaultRefType}(a, ordered=ordered) ≅ x
            @test CategoricalArray{Union{String, Missing}, 1, UInt8}(a, ordered=ordered) ≅ x

            @test CategoricalVector(a, ordered=ordered) ≅ x
            @test CategoricalVector{Union{String, Missing}}(a, ordered=ordered) ≅ x
            @test CategoricalVector{Union{String, Missing}, R}(a, ordered=ordered) ≅ x
            @test CategoricalVector{Union{String, Missing}, DefaultRefType}(a, ordered=ordered) ≅ x
            @test CategoricalVector{Union{String, Missing}, UInt8}(a, ordered=ordered) ≅ x

            @testset "categorical($(typeof(y)), compress=$comp) R1=$R1 R2=$R2" for (y, R1, R2, comp) in
                ((a, DefaultRefType, UInt8, true),
                 (a, DefaultRefType, DefaultRefType, false),
                 (x, R, UInt8, true),
                 (x, R, R, false))

                x2 = @inferred categorical(y, ordered=ordered)
                @test x2 ≅ y
                if eltype(y) >: Missing
                    @test isa(x2, CategoricalVector{Union{String, Missing}, R1})
                else
                    @test isa(x2, CategoricalVector{String, R1})
                end
                @test isordered(x2) === ordered

                x2 = categorical(y, compress=comp, ordered=ordered)
                @test x2 ≅ y
                if eltype(y) >: Missing
                    @test isa(x2, CategoricalVector{Union{String, Missing}, R2})
                else
                    @test isa(x2, CategoricalVector{String, R2})
                end
                @test isordered(x2) === ordered
            end

            x2 = compress(x)
            @test x2 ≅ x
            @test isa(x2, CategoricalVector{Union{String, Missing}, UInt8})
            @test levels(x2) == levels(x)

            x2 = copy(x)
            @test x2 ≅ x
            @test typeof(x2) === typeof(x)
            @test levels(x2) == levels(x)

            @test x[1] === x.pool.valindex[1]
            @test x[2] === x.pool.valindex[2]
            @test x[3] === missing
            @test_throws BoundsError x[4]

            x2 = x[:]
            @test typeof(x2) === typeof(x)
            @test x2 ≅ x
            @test x2 !== x
            @test levels(x2) == levels(x)
            @test levels(x2) !== levels(x)
            @test isordered(x2) == isordered(x)

            x2 = x[2:3]
            @test typeof(x2) === typeof(x)
            @test x2 ≅ ["b", missing]
            @test levels(x2) == levels(x)
            @test levels(x2) !== levels(x)
            @test isordered(x2) == isordered(x)

            x2 = x[1:1]
            @test typeof(x2) === typeof(x)
            @test x2 == ["a"]
            @test levels(x2) == levels(x)
            @test levels(x2) !== levels(x)
            @test isordered(x2) == isordered(x)

            x2 = x[2:1]
            @test typeof(x2) === typeof(x)
            @test isempty(x2)
            @test levels(x2) == levels(x)
            @test levels(x2) !== levels(x)
            @test isordered(x2) == isordered(x)

            x[1] = "b"
            @test x[1] === x.pool.valindex[2]
            @test x[2] === x.pool.valindex[2]
            @test x[3] === missing

            if ordered
                @test_throws OrderedLevelsException x[3] = "c"
                levels!(x, [levels(x); "c"])
            end
            x[3] = "c"
            @test x[1] === x.pool.valindex[2]
            @test x[2] === x.pool.valindex[2]
            @test x[3] === x.pool.valindex[3]
            @test levels(x) == ["a", "b", "c"]

            x[1] = missing
            @test x[1] === missing
            @test x[2] === x.pool.valindex[2]
            @test x[3] === x.pool.valindex[3]
            @test levels(x) == ["a", "b", "c"]

            x[2:3] .= missing
            @test x[1] === missing
            @test x[2] === missing
            @test x[3] === missing
            @test levels(x) == ["a", "b", "c"]
        end

        @testset "Vector created from range" begin
        # (i.e. non-Array AbstractArray),
        # direct conversion to a vector with different eltype
        a = 0.0:0.5:1.5
        x = CategoricalVector{Union{Float64, Missing}, R}(a, ordered=ordered)

        @test x == collect(a)
        @test isordered(x) === ordered
        @test levels(x) == unique(x) == unique(a)
        @test size(x) === (4,)
        @test length(x) === 4
        @test leveltype(x) === Float64
        @test eltype(x) == Union{CategoricalValue{Float64, R}, Missing}

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{Union{Float64, Missing}}, x) === x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1}, x) === x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, R}, x) === x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, DefaultRefType}, x) == x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, UInt8}, x) == x

        @test convert(CategoricalVector, x) === x
        @test convert(CategoricalVector{Union{Float64, Missing}}, x) === x
        @test convert(CategoricalVector{Union{Float64, Missing}, R}, x) === x
        @test convert(CategoricalVector{Union{Float64, Missing}, DefaultRefType}, x) == x
        @test convert(CategoricalVector{Union{Float64, Missing}, UInt8}, x) == x

        @test convert(CategoricalArray, a) == x
        @test convert(CategoricalArray{Union{Float64, Missing}}, a) == x
        @test convert(CategoricalArray{Union{Float32, Missing}}, a) == x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1}, a) == x
        @test convert(CategoricalArray{Union{Float32, Missing}, 1}, a) == x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, R}, a) == x
        @test convert(CategoricalArray{Union{Float32, Missing}, 1, R}, a) == x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, DefaultRefType}, a) == x
        @test convert(CategoricalArray{Union{Float32, Missing}, 1, DefaultRefType}, a) == x
        @test convert(CategoricalArray{Union{Float64, Missing}, 1, UInt8}, a) == x
        @test convert(CategoricalArray{Union{Float32, Missing}, 1, UInt8}, a) == x

        @test convert(CategoricalVector, a) == x
        @test convert(CategoricalVector{Union{Float64, Missing}}, a) == x
        @test convert(CategoricalVector{Union{Float32, Missing}}, a) == x
        @test convert(CategoricalVector{Union{Float64, Missing}, R}, a) == x
        @test convert(CategoricalVector{Union{Float32, Missing}, R}, a) == x
        @test convert(CategoricalVector{Union{Float64, Missing}, DefaultRefType}, a) == x
        @test convert(CategoricalVector{Union{Float32, Missing}, DefaultRefType}, a) == x
        @test convert(CategoricalVector{Union{Float64, Missing}, UInt8}, a) == x
        @test convert(CategoricalVector{Union{Float32, Missing}, UInt8}, a) == x

        @test CategoricalArray{Union{Float64, Missing}}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float32, Missing}}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float64, Missing}, 1}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float32, Missing}, 1}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float64, Missing}, 1, R}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float32, Missing}, 1, R}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float64, Missing}, 1, DefaultRefType}(a, ordered=ordered) == x
        @test CategoricalArray{Union{Float32, Missing}, 1, DefaultRefType}(a, ordered=ordered) == x

        @test CategoricalVector(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float64, Missing}}(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}}(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float64, Missing}, R}(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}, R}(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float64, Missing}, DefaultRefType}(a, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}, DefaultRefType}(a, ordered=ordered) == x

        @testset "categorical($(typeof(y)), compress=$comp) R1=$R1 R2=$R2" for (y, R1, R2, comp) in
            ((a, DefaultRefType, UInt8, true),
             (a, DefaultRefType, DefaultRefType, false),
             (x, R, UInt8, true),
             (x, R, R, false))

            x2 = @inferred categorical(y, ordered=ordered)
            @test x2 == collect(y)
            if eltype(y) >: Missing
                @test isa(x2, CategoricalVector{Union{Float64, Missing}, R1})
            else
                @test isa(x2, CategoricalVector{Float64, R1})
            end
            @test isordered(x2) === ordered
            @test leveltype(x2) === Float64
            @test nonmissingtype(eltype(x2)) === CategoricalValue{Float64, R1}

            x2 = categorical(y, compress=comp, ordered=ordered)
            @test x2 == collect(y)
            if eltype(y) >: Missing
                @test isa(x2, CategoricalVector{Union{Float64, Missing}, R2})
            else
                @test isa(x2, CategoricalVector{Float64, R2})
            end
            @test isordered(x2) === ordered
            @test leveltype(x2) === Float64
            @test nonmissingtype(eltype(x2)) === CategoricalValue{Float64, R2}
        end

        x2 = compress(x)
        @test x2 == x
        @test isordered(x2) === isordered(x)
        @test isa(x2, CategoricalVector{Union{Float64, Missing}, UInt8})
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 == x
        @test isordered(x2) === isordered(x)
        @test typeof(x2) === typeof(x)
        @test levels(x2) == levels(x)

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test_throws BoundsError x[5]

        x2 = x[:]
        @test typeof(x2) === typeof(x)
        @test x2 == x
        @test x2 !== x
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[1:2]
        @test typeof(x2) === typeof(x)
        @test x2 == [0.0, 0.5]
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[1:1]
        @test typeof(x2) === typeof(x)
        @test x2 == [0.0]
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[2:1]
        @test typeof(x2) === typeof(x)
        @test isempty(x2)
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x[2] = 1
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[3]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test levels(x) == unique(a)
        @test unique(x) == unique(collect(x))

        if ordered
            @test_throws OrderedLevelsException x[1:2] .= -1
            levels!(x, [levels(x); -1])
        end
        x[1:2] .= -1
        @test x[1] === x.pool.valindex[5]
        @test x[2] === x.pool.valindex[5]
        @test x[3] === x.pool.valindex[3]
        @test x[4] === x.pool.valindex[4]
        @test levels(x) == vcat(unique(a), -1)
        @test unique(x) == unique(collect(x))

        if ordered
            @test_throws OrderedLevelsException push!(x, 2.0)
            levels!(x, [levels(x); 2.0])
        end
        push!(x, 2.0)
        @test length(x) == 5
        @test x == [-1.0, -1.0, 1.0, 1.5, 2.0]
        @test isordered(x) === ordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        push!(x, x[1])
        @test length(x) == 6
        @test x == [-1.0, -1.0, 1.0, 1.5, 2.0, -1.0]
        @test isordered(x) === ordered
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

        empty!(x)
        @test length(x) == 0
        @test isordered(x) === ordered
        @test levels(x) == [0.0, 0.5, 1.0, 1.5, -1.0, 2.0]
    end

    @testset "Matrix $(typeof(a)) with no missing values" for a in
        (["a" "b" "c"; "b" "a" "c"],
         Union{String, Missing}["a" "b" "c"; "b" "a" "c"])

        x = CategoricalMatrix{Union{String, Missing}, R}(a, ordered=ordered)

        @test x == a
        @test isordered(x) === ordered
        @test levels(x) == unique(x) == unique(a)
        @test size(x) === (2, 3)
        @test length(x) === 6

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{Union{String, Missing}}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2, R}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2, DefaultRefType}, x) == x
        @test convert(CategoricalArray{Union{String, Missing}, 2, UInt8}, x) == x

        @test convert(CategoricalMatrix, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}}, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}, R}, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}, DefaultRefType}, x) == x
        @test convert(CategoricalMatrix{Union{String, Missing}, UInt8}, x) == x

        @test convert(CategoricalArray, a) == x
        @test convert(CategoricalArray{Union{String, Missing}}, a) == x
        @test convert(CategoricalArray{Union{String, Missing}, 2, R}, a) == x
        @test convert(CategoricalArray{Union{String, Missing}, 2, DefaultRefType}, a) == x
        @test convert(CategoricalArray{Union{String, Missing}, 2, UInt8}, a) == x

        @test convert(CategoricalMatrix, a) == x
        @test convert(CategoricalMatrix{Union{String, Missing}}, a) == x
        @test convert(CategoricalMatrix{Union{String, Missing}, R}, a) == x
        @test convert(CategoricalMatrix{Union{String, Missing}, DefaultRefType}, a) == x
        @test convert(CategoricalMatrix{Union{String, Missing}, UInt8}, a) == x

        @test CategoricalArray{Union{String, Missing}}(a, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2}(a, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2}(a, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, R}(a, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, DefaultRefType}(a, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, UInt8}(a, ordered=ordered) == x

        @test CategoricalMatrix(a, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}}(a, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, R}(a, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, DefaultRefType}(a, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, UInt8}(a, ordered=ordered) == x

        @testset "categorical($(typeof(y)), compress=$comp) R1=$R1 R2=$R2" for (y, R1, R2, comp) in
            ((a, DefaultRefType, UInt8, true),
             (a, DefaultRefType, DefaultRefType, false),
             (x, R, UInt8, true),
             (x, R, R, false))

            x2 = @inferred categorical(y, ordered=ordered)
            @test x2 == y
            if eltype(y) >: Missing
                @test isa(x2, CategoricalMatrix{Union{String, Missing}, R1})
            else
                @test isa(x2, CategoricalMatrix{String, R1})
            end
            @test isordered(x2) === ordered

            x2 = categorical(y, compress=comp, ordered=ordered)
            @test x2 == y
            if eltype(y) >: Missing
                @test isa(x2, CategoricalMatrix{Union{String, Missing}, R2})
            else
                @test isa(x2, CategoricalMatrix{String, R2})
            end
            @test isordered(x2) === ordered
        end

        x2 = compress(x)
        @test x2 == x
        @test isa(x2, CategoricalMatrix{Union{String, Missing}, UInt8})
        @test isordered(x2) === isordered(x)
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 == x
        @test typeof(x2) === typeof(x)
        @test isordered(x2) === isordered(x)
        @test levels(x2) == levels(x)

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === x.pool.valindex[3]
        @test_throws BoundsError x[7]

        @test x[1,1] === x.pool.valindex[1]
        @test x[2,1] === x.pool.valindex[2]
        @test x[1,2] === x.pool.valindex[2]
        @test x[2,2] === x.pool.valindex[1]
        @test x[1,3] === x.pool.valindex[3]
        @test x[2,3] === x.pool.valindex[3]
        @test_throws BoundsError x[1,4]
        @test_throws BoundsError x[4,1]
        @test_throws BoundsError x[4,4]

        @test x[1:2,:] == x
        @test typeof(x[1:2,:]) === typeof(x)
        @test x[1:2,1] == ["a", "b"]
        @test isa(x[1:2,1], CategoricalVector{Union{String, Missing}, R})

        if ordered
            @test_throws OrderedLevelsException x[1] = "z"
            levels!(x, [levels(x); "z"])
        end
        x[1] = "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[2]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,:] .= "a"
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,1:2] .= "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[4]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === x.pool.valindex[3]
        @test levels(x) == ["a", "b", "c", "z"]
        end

        @testset "Matrix with missing values" begin
        a = ["a" missing "c"; "b" "a" missing]
        x = CategoricalMatrix{Union{String, Missing}, R}(a, ordered=ordered)

        @test x ≅ a
        @test isordered(x) === ordered
        @test levels(x) == filter(x->!ismissing(x), unique(a))
        @test unique(x) ≅ unique(a)
        @test size(x) === (2, 3)
        @test length(x) === 6

        @test convert(CategoricalArray, x) === x
        @test convert(CategoricalArray{Union{String, Missing}}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2, R}, x) === x
        @test convert(CategoricalArray{Union{String, Missing}, 2, DefaultRefType}, x) ≅ x
        @test convert(CategoricalArray{Union{String, Missing}, 2, UInt8}, x) ≅ x

        @test convert(CategoricalMatrix, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}}, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}, R}, x) === x
        @test convert(CategoricalMatrix{Union{String, Missing}, DefaultRefType}, x) ≅ x
        @test convert(CategoricalMatrix{Union{String, Missing}, UInt8}, x) ≅ x

        @test convert(CategoricalArray, a) ≅ x
        @test convert(CategoricalArray{Union{String, Missing}}, a) ≅ x
        @test convert(CategoricalArray{Union{String, Missing}, 2, R}, a) ≅ x
        @test convert(CategoricalArray{Union{String, Missing}, 2, DefaultRefType}, a) ≅ x
        @test convert(CategoricalArray{Union{String, Missing}, 2, UInt8}, a) ≅ x

        @test convert(CategoricalMatrix, a) ≅ x
        @test convert(CategoricalMatrix{Union{String, Missing}}, a) ≅ x
        @test convert(CategoricalMatrix{Union{String, Missing}, R}, a) ≅ x
        @test convert(CategoricalMatrix{Union{String, Missing}, DefaultRefType}, a) ≅ x
        @test convert(CategoricalMatrix{Union{String, Missing}, UInt8}, a) ≅ x

        @test CategoricalArray{Union{String, Missing}}(a, ordered=ordered) ≅ x
        @test CategoricalArray{Union{String, Missing}, 2}(a, ordered=ordered) ≅ x
        @test CategoricalArray{Union{String, Missing}, 2}(a, ordered=ordered) ≅ x
        @test CategoricalArray{Union{String, Missing}, 2, R}(a, ordered=ordered) ≅ x
        @test CategoricalArray{Union{String, Missing}, 2, DefaultRefType}(a, ordered=ordered) ≅ x
        @test CategoricalArray{Union{String, Missing}, 2, UInt8}(a, ordered=ordered) ≅ x

        @test CategoricalMatrix(a, ordered=ordered) ≅ x
        @test CategoricalMatrix{Union{String, Missing}}(a, ordered=ordered) ≅ x
        @test CategoricalMatrix{Union{String, Missing}, R}(a, ordered=ordered) ≅ x
        @test CategoricalMatrix{Union{String, Missing}, DefaultRefType}(a, ordered=ordered) ≅ x
        @test CategoricalMatrix{Union{String, Missing}, UInt8}(a, ordered=ordered) ≅ x

        @testset "categorical($(typeof(y)), compress=$comp) R1=$R1 R2=$R2" for (y, R1, R2, comp) in
            ((a, DefaultRefType, UInt8, true),
             (a, DefaultRefType, DefaultRefType, false),
             (x, R, UInt8, true),
             (x, R, R, false))

            x2 = @inferred categorical(y, ordered=ordered)
            @test x2 ≅ y
            @test isa(x2, CategoricalMatrix{Union{String, Missing}, R1})
            @test isordered(x2) === ordered

            x2 = categorical(y, compress=comp, ordered=ordered)
            @test x2 ≅ y
            @test isa(x2, CategoricalMatrix{Union{String, Missing}, R2})
            @test isordered(x2) === ordered
        end

        x2 = compress(x)
        @test x2 ≅ x
        @test isa(x2, CategoricalMatrix{Union{String, Missing}, UInt8})
        @test isordered(x2) === isordered(x)
        @test levels(x2) == levels(x)

        x2 = copy(x)
        @test x2 ≅ x
        @test typeof(x2) === typeof(x)
        @test isordered(x2) === isordered(x)
        @test levels(x2) == levels(x)

        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === missing
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === missing
        @test_throws BoundsError x[7]

        @test x[1,1] === x.pool.valindex[1]
        @test x[2,1] === x.pool.valindex[2]
        @test x[1,2] === missing
        @test x[2,2] === x.pool.valindex[1]
        @test x[1,3] === x.pool.valindex[3]
        @test x[2,3] === missing
        @test_throws BoundsError x[1,4]
        @test_throws BoundsError x[4,1]
        @test_throws BoundsError x[4,4]

        x2 = x[1:2,:]
        @test typeof(x2) === typeof(x)
        @test x2 ≅ x
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[:,[1, 3]]
        @test typeof(x2) === typeof(x)
        @test x2 ≅ ["a" "c"; "b" missing]
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[1:1,2]
        @test isa(x2, CategoricalVector{Union{String, Missing}, R})
        @test x2 ≅ [missing]
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[1:0,:]
        @test typeof(x2) === typeof(x)
        @test size(x2) == (0,3)
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        @test_throws BoundsError x[1:4, :]
        @test_throws BoundsError x[1:1, -1:1]
        @test_throws BoundsError x[4, :]

        if ordered
            @test_throws OrderedLevelsException x[1] = "z"
            levels!(x, [levels(x); "z"])
        end
        x[1] = "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === missing
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[3]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,:] .= "a"
        @test x[1] === x.pool.valindex[1]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,1:2] .= "z"
        @test x[1] === x.pool.valindex[4]
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[4]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[1] = missing
        @test x[1] === missing
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[4]
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,1:2] .= missing
        @test x[1] === missing
        @test x[2] === x.pool.valindex[2]
        @test x[3] === missing
        @test x[4] === x.pool.valindex[1]
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[:,2] .= missing
        @test x[1] === missing
        @test x[2] === x.pool.valindex[2]
        @test x[3] === missing
        @test x[4] === missing
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[1,2] = "a"
        @test x[1] === missing
        @test x[2] === x.pool.valindex[2]
        @test x[3] === x.pool.valindex[1]
        @test x[4] === missing
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]

        x[2,1] = missing
        @test x[1] === missing
        @test x[2] === missing
        @test x[3] === x.pool.valindex[1]
        @test x[4] === missing
        @test x[5] === x.pool.valindex[1]
        @test x[6] === missing
        @test levels(x) == ["a", "b", "c", "z"]
        end

        # Uninitialized array
        v = Any[CategoricalArray{Union{String, Missing}}(undef, 2, ordered=ordered),
                CategoricalArray{Union{String, Missing}, 1}(undef, 2, ordered=ordered),
                CategoricalArray{Union{String, Missing}, 1, R}(undef, 2, ordered=ordered),
                CategoricalVector{Union{String, Missing}}(undef, 2, ordered=ordered),
                CategoricalVector{Union{String, Missing}, R}(undef, 2, ordered=ordered),
                CategoricalArray{Union{String, Missing}}(undef, 2, 3, ordered=ordered),
                CategoricalArray{Union{String, Missing}, 2}(undef, 2, 3, ordered=ordered),
                CategoricalArray{Union{String, Missing}, 2, R}(undef, 2, 3, ordered=ordered),
                CategoricalMatrix{Union{String, Missing}}(undef, 2, 3, ordered=ordered),
                CategoricalMatrix{Union{String, Missing}, R}(undef, 2, 3, ordered=ordered)]

        @testset "Uninitialized $(typeof(x))" for x in v
        @test isordered(x) === ordered
        @test ismissing(x[1])
        @test ismissing(x[2])
        @test levels(x) == []

        x2 = compress(x)
        @test x2 ≅ x
        @test isa(x2, CategoricalArray{Union{leveltype(x), Missing}, ndims(x), UInt8})
        @test isordered(x2) === isordered(x)
        @test levels(x2) == []

        x2 = copy(x)
        @test x2 ≅ x
        @test typeof(x2) === typeof(x)
        @test isordered(x2) === isordered(x)
        @test levels(x2) == []

        if ordered
            @test_throws OrderedLevelsException x[1] = "c"
            levels!(x, [levels(x); "c"])
        end
        x[1] = "c"
        @test x[1] === x.pool.valindex[1]
        @test ismissing(x[2])
        @test levels(x) == ["c"]

        if ordered
            @test_throws OrderedLevelsException x[1] = "a"
            levels!(x, [levels(x); "a"])
        end
        x[1] = "a"
        @test x[1] === x.pool.valindex[2]
        @test ismissing(x[2])
        @test levels(x) == ["c", "a"]

        x[2] = missing
        @test x[1] === x.pool.valindex[2]
        @test x[2] === missing
        @test levels(x) == ["c", "a"]

        if ordered
            @test_throws OrderedLevelsException x[1] = "b"
            levels!(x, [levels(x); "b"])
        end
        x[1] = "b"
        @test x[1] === x.pool.valindex[3]
        @test x[2] === missing
        @test levels(x) == ["c", "a", "b"]
        end
    end
end

@testset "vcat with missings" begin
    ca1 = CategoricalArray(["a", missing])
    ca2 = CategoricalArray([missing, "a"])
    r = vcat(ca1, ca2)
    @test r ≅ CategoricalArray(["a", missing, missing, "a"])
    @test levels(r) == ["a"]
    @test !isordered(r)
    ordered!(ca1,true)
    @test !isordered(vcat(ca1, ca2))
    ordered!(ca2,true)
    @test isordered(vcat(ca1, ca2))
    ordered!(ca1,false)
    @test !isordered(vcat(ca1, ca2))
end

@testset "vcat with all missings" begin
    ca1 = CategoricalArray(["a", missing])
    ca2 = CategoricalArray([missing, missing])
    r = vcat(ca1, ca2)
    @test r ≅ ["a", missing, missing, missing]
    @test levels(r) == ["a"]
    @test !isordered(r)

    @testset "preserves isordered" begin
    # needed for instance when expanding an array with missings
    # such as vcat of DataFrame with missing columns
    ordered!(ca1, true)
    @test isempty(levels(ca2))
    r = vcat(ca1, ca2)
    @test isordered(r)
    end
end

@testset "vcat with all empty array" begin
    ca1 = CategoricalArray(undef, 0)
    ca2 = CategoricalArray([missing, "b"])
    r = vcat(ca1, ca2)
    @test r ≅ [missing, "b"]
    @test levels(r) == ["b"]
    @test !isordered(r)
end

@testset "vcat with all missings and empty" begin
    ca1 = CategoricalArray(undef, 0)
    ca2 = CategoricalArray([missing, missing])
    r = vcat(ca1, ca2)
    @test r ≅ [missing, missing]
    @test levels(r) == String[]
    @test !isordered(r)

    ordered!(ca1, true)
    @test isempty(levels(ca2))
    r = vcat(ca1, ca2)
    @test isordered(r)

    ca1 = CategoricalArray(["a", missing])
    ca2 = CategoricalArray{Union{String, Missing}}(undef, 2)
    ordered!(ca1, true)
    @test isempty(levels(ca2))
    r = vcat(ca1, ca2)
    @test r ≅ ["a", missing, missing, missing]
    @test isordered(r)
end

@testset "unique() and levels()" begin
    x = CategoricalArray(["Old", "Young", "Middle", missing, "Young"])
    @test levels(x) == ["Middle", "Old", "Young"]
    @test unique(x) ≅ ["Old", "Young", "Middle", missing]
    @test levels!(x, ["Young", "Middle", "Old"]) === x
    @test levels(x) == ["Young", "Middle", "Old"]
    @test unique(x) ≅ ["Old", "Young", "Middle", missing]
    @test levels!(x, ["Young", "Middle", "Old", "Unused"]) === x
    @test levels(x) == ["Young", "Middle", "Old", "Unused"]
    @test unique(x) ≅ ["Old", "Young", "Middle", missing]
    @test levels!(x, ["Unused1", "Young", "Middle", "Old", "Unused2"]) === x
    @test levels(x) == ["Unused1", "Young", "Middle", "Old", "Unused2"]
    @test unique(x) ≅ ["Old", "Young", "Middle", missing]

    x = CategoricalArray((Union{String, Missing})[missing])
    @test isa(levels(x), Vector{String}) && isempty(levels(x))
    @test unique(x) ≅ [missing]
    @test levels!(x, ["Young", "Middle", "Old"]) === x
    @test levels(x) == ["Young", "Middle", "Old"]
    @test unique(x) ≅ [missing]

    # To test short-circuiting
    x = CategoricalArray{Union{Int, Missing}}(repeat(1:10, inner=10))
    @test levels(x) == collect(1:10)
    @test unique(x) == collect(1:10)
    @test levels!(x, [19:-1:1; 20]) === x
    x[3] = missing
    @test levels(x) == [19:-1:1; 20]
    @test unique(x) ≅ [1; missing; 2:10]

    # in
    x = CategoricalArray{Int}(repeat(1:1500, inner=10))
    @test !(missing in x)

    x = CategoricalArray{Union{Int, Missing}}(repeat(1:1500, inner=10))
    x[1] = missing
    @test missing in x
end

@testset "Missings.replace should work on CategoricalArrays" begin
    x = categorical(["a", "b", missing, "a"])
    y = ["a", "b", "", "a"]
    r = Missings.replace(x, "")
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
    @test levels(x) == ["a", "b"]
    @test levels(a) == ["a", "b", ""]

    r = Missings.replace(x, "b")
    y = ["a", "b", "b", "a"]
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
    @test levels(x) == ["a", "b"]
    @test levels(a) == ["a", "b"]

    @test_throws MethodError Missings.replace(x, 1)
end

@testset "Missings.replace should work on CategoricalArrays without missing values" begin
    x = categorical(["a", "b", "", "a"])
    y = ["a", "b", "", "a"]
    r = Missings.replace(x, "dummy")
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
    @test levels(x) == ["", "a", "b"]
    @test levels(a) == ["", "a", "b", "dummy"]

    r = Missings.replace(x, "")
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
    @test levels(x) == ["", "a", "b"]
    @test levels(a) == ["", "a", "b"]
end

@testset "Missings.replace should work on CategoricalArrays with empty pools" begin
    x = categorical(Union{String,Missing}[missing])
    y = [""]
    r = Missings.replace(x, "")
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
end

@testset "Missings.replace should work on empty CategoricalArrays" begin
    x = categorical(Union{String,Missing}[])
    y = String[]
    r = Missings.replace(x, "")
    @test isa(r, Missings.EachReplaceMissing)
    a = collect(r)
    @test isa(a, CategoricalVector{String})
    @test y == a
end

end


text:
# Julia wrapper for header: /opt/intel/compilers_and_libraries_2017.1.132/linux/mkl/include/mkl_dnn.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function dnnLayoutCreate_F32(pLayout,dimension,size,strides)
    check_dnnerror(ccall((:dnnLayoutCreate_F32,libmkl),dnnError_t,(Ptr{dnnLayout_t},Csize_t,Ptr{Csize_t},Ptr{Csize_t}),pLayout,dimension,size,strides))
end

function dnnLayoutCreateFromPrimitive_F32(pLayout,primitive,_type)
    check_dnnerror(ccall((:dnnLayoutCreateFromPrimitive_F32,libmkl),dnnError_t,(Ptr{dnnLayout_t},dnnPrimitive_t,dnnResourceType_t),pLayout,primitive,_type))
end

function dnnLayoutGetMemorySize_F32(layout)
    ccall((:dnnLayoutGetMemorySize_F32,libmkl),Csize_t,(dnnLayout_t,),layout)
end

function dnnLayoutCompare_F32(l1,l2)
    ccall((:dnnLayoutCompare_F32,libmkl),Cint,(dnnLayout_t,dnnLayout_t),l1,l2)
end

function dnnAllocateBuffer_F32(pPtr,layout)
    check_dnnerror(ccall((:dnnAllocateBuffer_F32,libmkl),dnnError_t,(Ptr{Ptr{Cvoid}},dnnLayout_t),pPtr,layout))
end

function dnnReleaseBuffer_F32(ptr)
    check_dnnerror(ccall((:dnnReleaseBuffer_F32,libmkl),dnnError_t,(Ptr{Cvoid},),ptr))
end

function dnnLayoutDelete_F32(layout)
    check_dnnerror(ccall((:dnnLayoutDelete_F32,libmkl),dnnError_t,(dnnLayout_t,),layout))
end

function dnnPrimitiveAttributesCreate_F32(attributes)
    check_dnnerror(ccall((:dnnPrimitiveAttributesCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitiveAttributes_t},),attributes))
end

function dnnPrimitiveAttributesDestroy_F32(attributes)
    check_dnnerror(ccall((:dnnPrimitiveAttributesDestroy_F32,libmkl),dnnError_t,(dnnPrimitiveAttributes_t,),attributes))
end

function dnnPrimitiveGetAttributes_F32(primitive,attributes)
    check_dnnerror(ccall((:dnnPrimitiveGetAttributes_F32,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{dnnPrimitiveAttributes_t}),primitive,attributes))
end

function dnnExecute_F32(primitive,resources)
    check_dnnerror(ccall((:dnnExecute_F32,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Ptr{Cvoid}}),primitive,resources))
end

function dnnExecuteAsync_F32(primitive,resources)
    check_dnnerror(ccall((:dnnExecuteAsync_F32,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Ptr{Cvoid}}),primitive,resources))
end

function dnnWaitFor_F32(primitive)
    check_dnnerror(ccall((:dnnWaitFor_F32,libmkl),dnnError_t,(dnnPrimitive_t,),primitive))
end

function dnnDelete_F32(primitive)
    check_dnnerror(ccall((:dnnDelete_F32,libmkl),dnnError_t,(dnnPrimitive_t,),primitive))
end

function dnnConversionCreate_F32(pConversion,from,to)
    check_dnnerror(ccall((:dnnConversionCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnLayout_t,dnnLayout_t),pConversion,from,to))
end

function dnnConversionExecute_F32(conversion,from,to)
    check_dnnerror(ccall((:dnnConversionExecute_F32,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Cvoid},Ptr{Cvoid}),conversion,from,to))
end

function dnnSumCreate_F32(pSum,attributes,nSummands,layout,coefficients)
    check_dnnerror(ccall((:dnnSumCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,dnnLayout_t,Ptr{Cfloat}),pSum,attributes,nSummands,layout,coefficients))
end

function dnnConcatCreate_F32(pConcat,attributes,nSrcTensors,src)
    check_dnnerror(ccall((:dnnConcatCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{dnnLayout_t}),pConcat,attributes,nSrcTensors,src))
end

function dnnSplitCreate_F32(pSplit,attributes,nDstTensors,layout,dstChannelSize)
    check_dnnerror(ccall((:dnnSplitCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,dnnLayout_t,Ptr{Csize_t}),pSplit,attributes,nDstTensors,layout,dstChannelSize))
end

function dnnScaleCreate_F32(pScale,attributes,dataLayout,alpha)
    check_dnnerror(ccall((:dnnScaleCreate_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cfloat),pScale,attributes,dataLayout,alpha))
end

function dnnConvolutionCreateForward_F32(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateForwardBias_F32(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateForwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardData_F32(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardData_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardFilter_F32(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardFilter_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardBias_F32(pConvolution,attributes,algorithm,dimension,dstSize)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t}),pConvolution,attributes,algorithm,dimension,dstSize))
end

function dnnGroupsConvolutionCreateForward_F32(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateForwardBias_F32(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateForwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardData_F32(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardData_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardFilter_F32(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardFilter_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardBias_F32(pConvolution,attributes,algorithm,groups,dimension,dstSize)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t}),pConvolution,attributes,algorithm,groups,dimension,dstSize))
end

function dnnReLUCreateForward_F32(pRelu,attributes,dataLayout,negativeSlope)
    check_dnnerror(ccall((:dnnReLUCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cfloat),pRelu,attributes,dataLayout,negativeSlope))
end

function dnnReLUCreateBackward_F32(pRelu,attributes,diffLayout,dataLayout,negativeSlope)
    check_dnnerror(ccall((:dnnReLUCreateBackward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,dnnLayout_t,Cfloat),pRelu,attributes,diffLayout,dataLayout,negativeSlope))
end

function dnnLRNCreateForward_F32(pLrn,attributes,dataLayout,kernel_size,alpha,beta,k)
    check_dnnerror(ccall((:dnnLRNCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Csize_t,Cfloat,Cfloat,Cfloat),pLrn,attributes,dataLayout,kernel_size,alpha,beta,k))
end

function dnnLRNCreateBackward_F32(pLrn,attributes,diffLayout,dataLayout,kernel_size,alpha,beta,k)
    check_dnnerror(ccall((:dnnLRNCreateBackward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,dnnLayout_t,Csize_t,Cfloat,Cfloat,Cfloat),pLrn,attributes,diffLayout,dataLayout,kernel_size,alpha,beta,k))
end

function dnnBatchNormalizationCreateForward_F32(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cfloat),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnBatchNormalizationCreateBackwardScaleShift_F32(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateBackwardScaleShift_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cfloat),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnBatchNormalizationCreateBackwardData_F32(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateBackwardData_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cfloat),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnPoolingCreateForward_F32(pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType)
    check_dnnerror(ccall((:dnnPoolingCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,dnnLayout_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType))
end

function dnnPoolingCreateBackward_F32(pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType)
    check_dnnerror(ccall((:dnnPoolingCreateBackward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,dnnLayout_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType))
end

function dnnInnerProductCreateForward_F32(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateForward_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateForwardBias_F32(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateForwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardData_F32(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardData_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardFilter_F32(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardFilter_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardBias_F32(pInnerProduct,attributes,dimensions,dstSize)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardBias_F32,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t}),pInnerProduct,attributes,dimensions,dstSize))
end

function dnnLayoutCreate_F64(pLayout,dimension,size,strides)
    check_dnnerror(ccall((:dnnLayoutCreate_F64,libmkl),dnnError_t,(Ptr{dnnLayout_t},Csize_t,Ptr{Csize_t},Ptr{Csize_t}),pLayout,dimension,size,strides))
end

function dnnLayoutCreateFromPrimitive_F64(pLayout,primitive,_type)
    check_dnnerror(ccall((:dnnLayoutCreateFromPrimitive_F64,libmkl),dnnError_t,(Ptr{dnnLayout_t},dnnPrimitive_t,dnnResourceType_t),pLayout,primitive,_type))
end

function dnnLayoutGetMemorySize_F64(layout)
    ccall((:dnnLayoutGetMemorySize_F64,libmkl),Csize_t,(dnnLayout_t,),layout)
end

function dnnLayoutCompare_F64(l1,l2)
    ccall((:dnnLayoutCompare_F64,libmkl),Cint,(dnnLayout_t,dnnLayout_t),l1,l2)
end

function dnnAllocateBuffer_F64(pPtr,layout)
    check_dnnerror(ccall((:dnnAllocateBuffer_F64,libmkl),dnnError_t,(Ptr{Ptr{Cvoid}},dnnLayout_t),pPtr,layout))
end

function dnnReleaseBuffer_F64(ptr)
    check_dnnerror(ccall((:dnnReleaseBuffer_F64,libmkl),dnnError_t,(Ptr{Cvoid},),ptr))
end

function dnnLayoutDelete_F64(layout)
    check_dnnerror(ccall((:dnnLayoutDelete_F64,libmkl),dnnError_t,(dnnLayout_t,),layout))
end

function dnnPrimitiveAttributesCreate_F64(attributes)
    check_dnnerror(ccall((:dnnPrimitiveAttributesCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitiveAttributes_t},),attributes))
end

function dnnPrimitiveAttributesDestroy_F64(attributes)
    check_dnnerror(ccall((:dnnPrimitiveAttributesDestroy_F64,libmkl),dnnError_t,(dnnPrimitiveAttributes_t,),attributes))
end

function dnnPrimitiveGetAttributes_F64(primitive,attributes)
    check_dnnerror(ccall((:dnnPrimitiveGetAttributes_F64,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{dnnPrimitiveAttributes_t}),primitive,attributes))
end

function dnnExecute_F64(primitive,resources)
    check_dnnerror(ccall((:dnnExecute_F64,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Ptr{Cvoid}}),primitive,resources))
end

function dnnExecuteAsync_F64(primitive,resources)
    check_dnnerror(ccall((:dnnExecuteAsync_F64,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Ptr{Cvoid}}),primitive,resources))
end

function dnnWaitFor_F64(primitive)
    check_dnnerror(ccall((:dnnWaitFor_F64,libmkl),dnnError_t,(dnnPrimitive_t,),primitive))
end

function dnnDelete_F64(primitive)
    check_dnnerror(ccall((:dnnDelete_F64,libmkl),dnnError_t,(dnnPrimitive_t,),primitive))
end

function dnnConversionCreate_F64(pConversion,from,to)
    check_dnnerror(ccall((:dnnConversionCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnLayout_t,dnnLayout_t),pConversion,from,to))
end

function dnnConversionExecute_F64(conversion,from,to)
    check_dnnerror(ccall((:dnnConversionExecute_F64,libmkl),dnnError_t,(dnnPrimitive_t,Ptr{Cvoid},Ptr{Cvoid}),conversion,from,to))
end

function dnnSumCreate_F64(pSum,attributes,nSummands,layout,coefficients)
    check_dnnerror(ccall((:dnnSumCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,dnnLayout_t,Ptr{Cdouble}),pSum,attributes,nSummands,layout,coefficients))
end

function dnnConcatCreate_F64(pConcat,attributes,nSrcTensors,src)
    check_dnnerror(ccall((:dnnConcatCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{dnnLayout_t}),pConcat,attributes,nSrcTensors,src))
end

function dnnSplitCreate_F64(pSplit,attributes,nDstTensors,layout,dstChannelSize)
    check_dnnerror(ccall((:dnnSplitCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,dnnLayout_t,Ptr{Csize_t}),pSplit,attributes,nDstTensors,layout,dstChannelSize))
end

function dnnScaleCreate_F64(pScale,attributes,dataLayout,alpha)
    check_dnnerror(ccall((:dnnScaleCreate_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cdouble),pScale,attributes,dataLayout,alpha))
end

function dnnConvolutionCreateForward_F64(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateForwardBias_F64(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateForwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardData_F64(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardData_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardFilter_F64(pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardFilter_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnConvolutionCreateBackwardBias_F64(pConvolution,attributes,algorithm,dimension,dstSize)
    check_dnnerror(ccall((:dnnConvolutionCreateBackwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Ptr{Csize_t}),pConvolution,attributes,algorithm,dimension,dstSize))
end

function dnnGroupsConvolutionCreateForward_F64(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateForwardBias_F64(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateForwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardData_F64(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardData_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardFilter_F64(pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardFilter_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pConvolution,attributes,algorithm,groups,dimension,srcSize,dstSize,filterSize,convolutionStrides,inputOffset,borderType))
end

function dnnGroupsConvolutionCreateBackwardBias_F64(pConvolution,attributes,algorithm,groups,dimension,dstSize)
    check_dnnerror(ccall((:dnnGroupsConvolutionCreateBackwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,Csize_t,Csize_t,Ptr{Csize_t}),pConvolution,attributes,algorithm,groups,dimension,dstSize))
end

function dnnReLUCreateForward_F64(pRelu,attributes,dataLayout,negativeSlope)
    check_dnnerror(ccall((:dnnReLUCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cdouble),pRelu,attributes,dataLayout,negativeSlope))
end

function dnnReLUCreateBackward_F64(pRelu,attributes,diffLayout,dataLayout,negativeSlope)
    check_dnnerror(ccall((:dnnReLUCreateBackward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,dnnLayout_t,Cdouble),pRelu,attributes,diffLayout,dataLayout,negativeSlope))
end

function dnnLRNCreateForward_F64(pLrn,attributes,dataLayout,kernel_size,alpha,beta,k)
    check_dnnerror(ccall((:dnnLRNCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Csize_t,Cdouble,Cdouble,Cdouble),pLrn,attributes,dataLayout,kernel_size,alpha,beta,k))
end

function dnnLRNCreateBackward_F64(pLrn,attributes,diffLayout,dataLayout,kernel_size,alpha,beta,k)
    check_dnnerror(ccall((:dnnLRNCreateBackward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,dnnLayout_t,Csize_t,Cdouble,Cdouble,Cdouble),pLrn,attributes,diffLayout,dataLayout,kernel_size,alpha,beta,k))
end

function dnnBatchNormalizationCreateForward_F64(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cdouble),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnBatchNormalizationCreateBackwardScaleShift_F64(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateBackwardScaleShift_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cdouble),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnBatchNormalizationCreateBackwardData_F64(pBatchNormalization,attributes,dataLayout,eps)
    check_dnnerror(ccall((:dnnBatchNormalizationCreateBackwardData_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnLayout_t,Cdouble),pBatchNormalization,attributes,dataLayout,eps))
end

function dnnPoolingCreateForward_F64(pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType)
    check_dnnerror(ccall((:dnnPoolingCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,dnnLayout_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType))
end

function dnnPoolingCreateBackward_F64(pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType)
    check_dnnerror(ccall((:dnnPoolingCreateBackward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,dnnAlgorithm_t,dnnLayout_t,Ptr{Csize_t},Ptr{Csize_t},Ptr{Cint},dnnBorder_t),pPooling,attributes,op,srcLayout,kernelSize,kernelStride,inputOffset,borderType))
end

function dnnInnerProductCreateForward_F64(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateForward_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateForwardBias_F64(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateForwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardData_F64(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardData_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardFilter_F64(pInnerProduct,attributes,dimensions,srcSize,outputChannels)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardFilter_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t},Csize_t),pInnerProduct,attributes,dimensions,srcSize,outputChannels))
end

function dnnInnerProductCreateBackwardBias_F64(pInnerProduct,attributes,dimensions,dstSize)
    check_dnnerror(ccall((:dnnInnerProductCreateBackwardBias_F64,libmkl),dnnError_t,(Ptr{dnnPrimitive_t},dnnPrimitiveAttributes_t,Csize_t,Ptr{Csize_t}),pInnerProduct,attributes,dimensions,dstSize))
end


text:
#Collection of utility functions

#Load Constants
const SECS = 86400.0
const MJD_J2000 = 51544.5
const JULIAN_CENTURY = 36525
const GMST_CONSTS = (24110.54841, 8640184.812866, 0.093104, 6.2e-6, 1.002737909350795)
const JD_TO_MJD = 2400000.5
const SECONDS_PER_CENTURY = SECS*JULIAN_CENTURY

"""
```julia
  qdyn, M = Get_Dynamic_Pressure(ρ::T, V_VA::S, ss::P) where {T,S,P <: Real}
```
Computes the dynamic pressure and local Mach Number from the air density,
the airspeed (in the V-frame), and the speed of sound.
"""
function Get_Dynamic_Pressure(ρ::T, V_VA::S, ss::P) where {T,S,P <: Real}
  return 0.5*ρ*V_VA^2, V_VA/ss
end

"""
```julia
Get_GMST(MJD::T) where {T <: Real}
```
This function computes the Greenwich Mean Sidereal Time and outputs it as an
angle (0..2pi).
"""
function Get_GMST(MJD::T) where {T <: Real}
  MJD0 = floor(MJD)
  UT1 = SECS*(MJD - MJD0)
  T0 = (MJD0 - MJD_J2000) / JULIAN_CENTURY
  Tval = (MJD - MJD_J2000) / JULIAN_CENTURY
  temp = (GMST_CONSTS[1] + GMST_CONSTS[2]*T0 + GMST_CONSTS[3]*Tval^2 - GMST_CONSTS[4]*Tval^3 + GMST_CONSTS[5]*UT1)/SECS
  return 2pi*(temp - floor(temp))
end

"""
```julia
Get_GMST(dt::DateTime, simtime::T = 0.0) where {T <: Real}
```
This function computes the Greenwich Mean Sidereal Time and outputs it as an
angle (0..2pi).
"""
function Get_GMST(dt::DateTime, simtime::T = 0.0) where {T <: Real}
  return Get_GMST(datetime2modifiedjulian(dt, simtime))
end

"""
```julia
datetime2modifiedjulian(dt::DateTime, simtime::Float64 = 0.0)
```
This function computes the Modified Julian Data from calender date and (optional) simulation time.
"""
function datetime2modifiedjulian(dt::DateTime, simtime::T = 0.0) where {T <: Real}
  return (Dates.datetime2julian(dt) - JD_TO_MJD + simtime/SECS)
end

function centuriespastJ2000(dt::DateTime, simtime::T = 0.0) where {T <: Real}
  return (datetime2modifiedjulian(dt,simtime) - MJD_J2000)/JULIAN_CENTURY
end


text:
"""
$(TYPEDEF)

Interpolant-basis weighted sum-of-squares polynomial cone of dimension `U`, for
real or real-valued complex polynomials, parametrized by vector of matrices
`Ps` derived from interpolant basis and polynomial domain constraints.

    $(FUNCTIONNAME){T, R}(U::Int, Ps::Vector{Matrix{R}}, use_dual::Bool = false)
"""
mutable struct WSOSInterpNonnegative{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    Ps::Vector{Matrix{R}}
    nu::Int

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    hess_fact_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, Matrix{T}}
    inv_hess::Symmetric{T, Matrix{T}}
    hess_fact_mat::Symmetric{T, Matrix{T}}
    hess_fact::Factorization{T}
    use_hess_prod_slow::Bool
    use_hess_prod_slow_updated::Bool

    tempLL::Vector{Matrix{R}}
    tempLL2::Vector{Matrix{R}}
    tempLU::Vector{Matrix{R}}
    ΛFLP::Vector{Matrix{R}}
    tempUU::Matrix{R}
    ΛF::Vector
    Ps_times::Vector{Float64}
    Ps_order::Vector{Int}

    function WSOSInterpNonnegative{T, R}(
        U::Int,
        Ps::Vector{Matrix{R}};
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        for Pk in Ps
            @assert size(Pk, 1) == U
        end
        cone = new{T, R}()
        cone.use_dual_barrier = !use_dual # using dual barrier
        cone.dim = U
        cone.Ps = Ps
        cone.nu = sum(size(Pk, 2) for Pk in Ps)
        return cone
    end
end

reset_data(cone::WSOSInterpNonnegative) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated =
    cone.use_hess_prod_slow = cone.use_hess_prod_slow_updated = false)

function setup_extra_data!(
    cone::WSOSInterpNonnegative{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    dim = cone.dim
    Ls = [size(Pk, 2) for Pk in cone.Ps]
    cone.tempLL = [zeros(R, L, L) for L in Ls]
    cone.tempLL2 = [zeros(R, L, L) for L in Ls]
    cone.tempLU = [zeros(R, L, dim) for L in Ls]
    cone.ΛFLP = [zeros(R, L, dim) for L in Ls]
    cone.tempUU = zeros(R, dim, dim)
    K = length(Ls)
    cone.ΛF = Vector{Any}(undef, K)
    cone.Ps_times = zeros(K)
    cone.Ps_order = collect(1:K)
    return cone
end

set_initial_point!(arr::AbstractVector, cone::WSOSInterpNonnegative) = (arr .= 1)

function update_feas(cone::WSOSInterpNonnegative)
    @assert !cone.feas_updated
    D = Diagonal(cone.point)

    # order the Ps by how long it takes to check feasibility, to improve efficiency
    sortperm!(cone.Ps_order, cone.Ps_times, initialized = true) # stochastic

    cone.is_feas = true
    for k in cone.Ps_order
        cone.Ps_times[k] = @elapsed @inbounds begin
            Pk = cone.Ps[k]
            LUk = cone.tempLU[k]
            LLk = cone.tempLL[k]

            # Λ = Pk' * Diagonal(point) * Pk
            # TODO for complex case, can do sqrt and outer product
            @. LUk = Pk' * cone.point' # currently faster than mul!(LUk, Pk', D)
            mul!(LLk, LUk, Pk)

            ΛFk = cone.ΛF[k] = cholesky!(Hermitian(LLk, :L), check = false)
            if !isposdef(ΛFk)
                cone.is_feas = false
                break
            end
        end
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::WSOSInterpNonnegative{T}) where {T}
    # condition is necessary but not sufficient for dual feasibility
    return all(>(eps(T)), cone.dual_point)
end

function update_grad(cone::WSOSInterpNonnegative)
    @assert cone.is_feas

    cone.grad .= 0
    @inbounds for k in eachindex(cone.Ps)
        ΛFLPk = cone.ΛFLP[k] # computed here
        ldiv!(ΛFLPk, cone.ΛF[k].L, cone.Ps[k]')
        @views for j in 1:cone.dim
            cone.grad[j] -= sum(abs2, ΛFLPk[:, j])
        end
    end

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::WSOSInterpNonnegative)
    @assert cone.grad_updated
    isdefined(cone, :hess) || alloc_hess!(cone)
    UU = cone.tempUU

    cone.hess .= 0
    @inbounds for k in eachindex(cone.Ps)
        outer_prod!(cone.ΛFLP[k], UU, true, false)
        for j in 1:cone.dim, i in 1:j
            cone.hess.data[i, j] += abs2(UU[i, j])
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod_slow!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::WSOSInterpNonnegative,
    )
    cone.use_hess_prod_slow_updated || update_use_hess_prod_slow(cone)
    @assert cone.hess_updated
    cone.use_hess_prod_slow || return hess_prod!(prod, arr, cone)

    @assert cone.grad_updated
    prod .= 0
    @inbounds for k in eachindex(cone.Ps)
        ΛFLPk = cone.ΛFLP[k]
        LUk = cone.tempLU[k]
        LLk = cone.tempLL2[k]
        @views for j in 1:size(arr, 2)
            partial_lambda!(LUk, arr[:, j], LLk, ΛFLPk)
            for i in 1:cone.dim
                prod[i, j] += real(dot(ΛFLPk[:, i], LUk[:, i]))
            end
        end
    end
    return prod
end

function dder3(cone::WSOSInterpNonnegative, dir::AbstractVector)
    @assert cone.grad_updated
    dder3 = cone.dder3
    dder3 .= 0
    @inbounds for k in eachindex(cone.Ps)
        LUk = partial_lambda!(cone.tempLU[k], dir, cone.tempLL2[k], cone.ΛFLP[k])
        @views for j in 1:cone.dim
            dder3[j] += sum(abs2, LUk[:, j])
        end
    end
    return dder3
end

function partial_lambda!(
    LUk::Matrix,
    dir::AbstractVector,
    LLk::Matrix,
    ΛFLPk::Matrix,
    )
    mul!(LUk, ΛFLPk, Diagonal(dir))
    mul!(LLk, LUk, ΛFLPk')
    mul!(LUk, Hermitian(LLk), ΛFLPk)
    return LUk
end


text:
"""
Check if an expression is a call to a function
"""
is_call_to(e, f) = false

is_call_to(e::Expr, f) =
    e.head == :call && e.args[1] == f

is_call_to(e::Expr, f, nargs) =
    is_call_to(e, f) && length(e.args) == (nargs + 1)

"""
    walk_expr_for_syms(expression, symbol_modification_function)

Walk an AST for symbols and transform them with a given function.

Returns a tuple of a new expression and a boolean indicating whether symbols
were replaced in the expression.
"""
walk_expr_for_syms(e, f) =
    return e, e, missing, false

walk_expr_for_syms(q::QuoteNode, f) =
    walk_expr_for_syms(Meta.quot(q.value), f)

walk_expr_for_syms(expr::Expr, symbol_dict::Dict) =
    walk_expr_for_syms(expr, (e, s) -> get(symbol_dict, s, e))

walk_expr_for_syms(expr::Expr, symbol_func::Function, symbols::Array{Symbol,1}) =
    walk_expr_for_syms(expr, (e, s) -> s in symbols ? symbol_func(e, s) : expr)

function walk_expr_for_syms(e::Expr, symbol_func::Function)
    # handle nested symbolic contexts
    if is_call_to(e, :SymbolContext)
        # anonymous contexts
        e, e, missing, false
    elseif e isa Expr && e.head == :call && is_call_to(e.args[1], :SymbolContext)
        # context that is provided with contextual values
        e, e, missing, false
    elseif is_call_to(e, :^, 1)
        e.args[2], e.args[2], missing, false
    elseif is_call_to(e, :sym)
        e, e, missing, false

    # handle macros
    elseif e.head == :macrocall
        e, e, missing, false

    # handle symbol nodes
    elseif e.head == :. && e.args[2] isa QuoteNode
        e, e, missing, false
    elseif e.args[1] == :(:.)
        new_expr = Expr(e.head, symbol_func(e, e.args[1].value), e.args[2:end]...)
        print_expr = Expr(e.head, Highlighted(e.args[1].value), e.args[2:end]...)
        return new_expr, print_expr, [e.args[1].value], true
    elseif e.head == :quote
        symbol_func(e, e.args[1]), Highlighted(e.args[1]), [e.args[1]], true

    # otherwise continue traversal
    else
        walks = walk_expr_for_syms.(e.args, [symbol_func])
        new_expr = Expr(e.head, getindex.(walks, [1])...)
        print_expr = Expr(e.head, getindex.(walks, [2])...)
        syms = unique(vcat(skipmissing(getindex.(walks, [3]))...))
        any_subs = any(last.(walks))
        return new_expr, print_expr, syms, any_subs
    end
end


text:
# TODO: Adapt python model below to new Julia format
# ################################################################################
# ##                                  HfO2                                      ##
# ################################################################################
#
# def n_HfO2_sym(axis='a'):
#     """This function creates a symbolic representation (using SymPy) of the
#     Sellmeier Equation model for the wavelength dependence
#     of crystalline HfO2's ordinary and extraordinary indices of refraction as
#     well as the index of thin-film amorphous HfO2.
#     Sellmeier coefficients for crystalline HfO2 are taken from
#     Cauchy Equation coefficients for amorphous ~100nm sputtered Hafnia (HfO2)
#     films taken from
#     Khoshman and Kordesch. "Optical properties of a-HfO2 thin films."
#     Surface and Coatings Technology 201.6 (2006)
#     https://doi.org/10.1016/j.surfcoat.2006.08.074
#     Fits measured data from sputtered amorphous Hafnia, I haven't found a fit
#     for ALD Hafnia specifically. They also report loss, with a sharp absorption
#     edge near 5.68 ± 0.09 eV (~218 nm)
#     Eqn. form:
#     n = A + B / lam**2 + C / lam**4
#     ng = A + 3 * B / lam**2 + 5 * C / lam**4
#     This model is then exported to other functions that use it and its
#     derivatives to return index, group index and GVD values as a function
#     of temperature and wavelength.
#     Variable units are lm in [um] and T in [deg C]
#     """
#     # if axis is 'e':
#     #     coeffs = [  (0.41344023,0.03684262), # A_i,λ_i Lorentz oscillator strength [1] and resonance wavelength [μm]
#     #                 (0.50497499,0.09076162),
#     #                 (2.4904862,23.771995),
#     #                 ]
#     # elif axis is 'o':
#     #     coeffs = [  (0.48755108,0.04338408), # A_i,λ_i Lorentz oscillator strength [1] and resonance wavelength [μm]
#     #                 (0.39875031,0.09461442),
#     #                 (2.3120353,23.793604),
#     #                 ]
#     if axis is 'a':
#         lm, T = sp.symbols('lm T')
#         # # fit for spectroscopic ellipsometer measurement for a 250nm thick film, good 300-1400nm
#         # A = 1.85
#         # B = 1.17e-2 # note values here are for λ in μm, vs nm in the paper
#         # C = 0.0
#         # fit for spectroscopic ellipsometer measurement for a 112nm thick film, good 300-1400nm
#         A = 1.86
#         B = 7.16e-3 # note values here are for λ in μm, vs nm in the paper
#         C = 0.0
#         n_sym =  A + B / lm**2 + C / lm**4
#         return lm, T, n_sym
#     if axis is 's':
#         lm, T = sp.symbols('lm T')
#         A = 1.59e4 # μm^-2
#         λ0 = 210.16e-3 # μm, note values here are for λ in μm, vs nm in the paper
#         n_sym = sp.sqrt(   1  + ( A * lm**2 ) / ( lm**2/λ0**2 - 1 ) )
#         return lm, T, n_sym
#     else:
#         raise Exception('unrecognized axis! must be "e", "o" or "a"')
#
#
# def n_HfO2(lm_in,T_in=300*u.degK,axis='a'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the index of refraction of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') HfO2.
#     Sellmeier coefficients for crystalline HfO2 are taken from
#     Cauchy Equation coefficients for amorphous ~100nm sputtered Hafnia (HfO2)
#     films taken from
#     Khoshman and Kordesch. "Optical properties of a-HfO2 thin films."
#     Surface and Coatings Technology 201.6 (2006)
#     https://doi.org/10.1016/j.surfcoat.2006.08.074
#     Fits measured data from sputtered amorphous Hafnia, I haven't found a fit
#     for ALD Hafnia specifically. They also report loss, with a sharp absorption
#     edge near 5.68 ± 0.09 eV (~218 nm)
#     Eqn. form:
#     n = A + B / lam**2 + C / lam**4
#     ng = A + 3 * B / lam**2 + 5 * C / lam**4
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = np.array([Q_(lm_in).to(u.um).magnitude]).flatten()
#     T_C = np.array([Q_(T_in).to(u.degC).magnitude]).flatten()
#     lm, T, n_sym = n_HfO2_sym(axis=axis)
#     n = sp.lambdify([lm,T],n_sym,'numpy')
#     output = np.zeros((T_C.size,lm_um.size))
#     for T_idx, TT in enumerate(T_C):
#         output[T_idx,:] = n(lm_um, T_C[T_idx])
#     return output.squeeze()
#
# def n_g_HfO2(lm_in,T_in=300*u.degK,axis='a'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group index of refraction of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') HfO2.
#     Sellmeier coefficients for crystalline HfO2 are taken from
#     Cauchy Equation coefficients for amorphous ~100nm sputtered Hafnia (HfO2)
#     films taken from
#     Khoshman and Kordesch. "Optical properties of a-HfO2 thin films."
#     Surface and Coatings Technology 201.6 (2006)
#     https://doi.org/10.1016/j.surfcoat.2006.08.074
#     Fits measured data from sputtered amorphous Hafnia, I haven't found a fit
#     for ALD Hafnia specifically. They also report loss, with a sharp absorption
#     edge near 5.68 ± 0.09 eV (~218 nm)
#     Eqn. form:
#     n = A + B / lam**2 + C / lam**4
#     ng = A + 3 * B / lam**2 + 5 * C / lam**4
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_HfO2_sym(axis=axis)
#     n_sym_prime = sp.diff(n_sym,lm)
#     n_g_sym = n_sym - lm*n_sym_prime
#     n_g = sp.lambdify([lm,T],n_g_sym,'numpy')
#     return n_g(lm_um, T_C)
#
# def gvd_HfO2(lm_in,T_in=300*u.degK,axis='a'):
#     """Sellmeier Equation model for the temperature and wavelength dependence
#     of the group velocity dispersion of crystalline (axis='o' or 'e') and
#     amorphous (axis='a') HfO2.
#     Sellmeier coefficients for crystalline HfO2 are taken from
#     Cauchy Equation coefficients for amorphous ~100nm sputtered Hafnia (HfO2)
#     films taken from
#     Khoshman and Kordesch. "Optical properties of a-HfO2 thin films."
#     Surface and Coatings Technology 201.6 (2006)
#     https://doi.org/10.1016/j.surfcoat.2006.08.074
#     Fits measured data from sputtered amorphous Hafnia, I haven't found a fit
#     for ALD Hafnia specifically. They also report loss, with a sharp absorption
#     edge near 5.68 ± 0.09 eV (~218 nm)
#     Eqn. form:
#     n = A + B / lam**2 + C / lam**4
#     ng = A + 3 * B / lam**2 + 5 * C / lam**4
#     Variable units are lm in [um] and T in [deg C]
#     """
#     lm_um = Q_(lm_in).to(u.um).magnitude
#     T_C = Q_(T_in).to(u.degC).magnitude
#     lm, T, n_sym = n_HfO2_sym(axis=axis)
#     n_sym_double_prime = sp.diff(n_sym,lm,lm)
#     c = Q_(3e8,'m/s') # unitful definition of speed of light
#     gvd_sym_no_prefactors = (lm**3)*n_sym_double_prime # symbolic gvd without unitful prefactors, to be made unitful below
#     gvd_no_prefactors = sp.lambdify([lm,T],gvd_sym_no_prefactors,'numpy') # numerical gvd without unitful prefactors
#     gvd = (1 / (2 * np.pi * (c**2))) * Q_(gvd_no_prefactors(lm_um,T_C),'um')
#     return gvd.to('fs**2 / mm')