text:
"""
    p1b_p1

The manufactured-solution colliding flow example from Elman et al 2014.
Formulation on linear triangles with a cubic bubble for the velocity.

The formulation is the one derived in Reddy, Introduction to the finite element
method, 1993. Page 486 ff.
"""
module p1b_p1

using LinearAlgebra
using StaticArrays
using MeshCore: nrelations, nentities, ir_identity, attribute
using MeshSteward: T3block
using MeshSteward: Mesh, attach!, baseincrel, boundary
using MeshSteward: vselect, geometry, summary, transform
using MeshSteward: vtkwrite
using Elfel.RefShapes: manifdim, manifdimv
using Elfel.FElements: FEH1_T3_BUBBLE, FEH1_T3, refshape, Jacobian
using Elfel.FESpaces: FESpace, ndofs, setebc!, nunknowns, doftype
using Elfel.FESpaces: numberfreedofs!, numberdatadofs!, numberdofs!
using Elfel.FESpaces: scattersysvec!, makeattribute, gathersysvec!, edofcompnt
using Elfel.FESpaces: highestfreedofnum, highestdatadofnum
using Elfel.FEIterators: FEIterator, ndofsperel, elnodes, eldofs, eldofvals
using Elfel.FEIterators: jacjac, location
using Elfel.QPIterators: QPIterator, bfun, bfungrad, weight
using Elfel.Assemblers: SysmatAssemblerSparse, start!, finish!, assemble!
using Elfel.Assemblers: SysvecAssembler
using Elfel.LocalAssemblers: LocalMatrixAssembler, LocalVectorAssembler, init!
using UnicodePlots

mu = 1.0 # dynamic viscosity
A = 1.0 # half of the length of the side of the square
trueux = (x, y) -> 20 * x * y ^ 3
trueuy = (x, y) -> 5 * x ^ 4 - 5 * y ^ 4
truep = (x, y) -> 60 * x ^ 2 * y - 20 * y ^ 3

function genmesh(N)
    # This mesh will be both for the velocities and for the pressure
    mesh = Mesh()
    attach!(mesh, T3block(2 * A, 2 * A, N, N), "velocity+pressure")
    ir = baseincrel(mesh)
    transform(ir, x -> x .- A)
    eidir = ir_identity(ir)
    attach!(mesh, eidir)
    return mesh
end

function assembleK(Uxh, Uyh, Ph, tndof, mu)
    function integrateK!(ass, elits, qpits, mu)
        uxnedof, uynedof, pnedof = ndofsperel.(elits)
        kuxux = LocalMatrixAssembler(uxnedof, uxnedof, 0.0)
        kuyuy = LocalMatrixAssembler(uynedof, uynedof, 0.0)
        kuxuy = LocalMatrixAssembler(uxnedof, uynedof, 0.0)
        kuxp = LocalMatrixAssembler(uxnedof, pnedof, 0.0)
        kuyp = LocalMatrixAssembler(uynedof, pnedof, 0.0)
        for el in zip(elits...)
            uxel, uyel, pel = el
            init!(kuxux, eldofs(uxel), eldofs(uxel))
            init!(kuyuy, eldofs(uyel), eldofs(uyel))
            init!(kuxuy, eldofs(uxel), eldofs(uyel))
            init!(kuxp, eldofs(uxel), eldofs(pel))
            init!(kuyp, eldofs(uyel), eldofs(pel))
            for qp in zip(qpits...)
                uxqp, uyqp, pqp = qp
                Jac, J = jacjac(pel, pqp)
                JxW = J * weight(pqp)
                gradNp = bfungrad(pqp, Jac)
                gradNux = bfungrad(uxqp, Jac)
                gradNuy = bfungrad(uyqp, Jac)
                Np = bfun(pqp)
                for j in 1:uxnedof, i in 1:uxnedof
                    kuxux[i, j] += (mu * JxW) * (2 * gradNux[i][1] * gradNux[j][1] + gradNux[i][2] * gradNux[j][2])
                end
                for j in 1:uynedof, i in 1:uynedof
                    kuyuy[i, j] += (mu * JxW) * (gradNuy[i][1] * gradNuy[j][1] + 2 * gradNuy[i][2] * gradNuy[j][2])
                end
                for j in 1:uynedof, i in 1:uxnedof
                    kuxuy[i, j] += (mu * JxW) * (gradNux[i][1] * gradNuy[j][2])
                end
                for j in 1:pnedof, i in 1:uxnedof
                    kuxp[i, j] += (-JxW) * (gradNux[i][1] * Np[j])
                end
                for j in 1:pnedof, i in 1:uynedof
                    kuyp[i, j] += (-JxW) * (gradNuy[i][2] * Np[j])
                end
            end
            assemble!(ass, kuxux)
            assemble!(ass, kuxuy)
            assemble!(ass, transpose(kuxuy))
            assemble!(ass, kuyuy)
            assemble!(ass, kuxp)
            assemble!(ass, transpose(kuxp))
            assemble!(ass, kuyp)
            assemble!(ass, transpose(kuyp))
        end
        return ass
    end

    elits = (FEIterator(Uxh), FEIterator(Uyh), FEIterator(Ph))
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uxh, qargs), QPIterator(Uyh, qargs), QPIterator(Ph, qargs))
    ass = SysmatAssemblerSparse(0.0)
    start!(ass, tndof, tndof)
    integrateK!(ass, elits, qpits, mu)
    return finish!(ass)
end

function solve!(U, K, F, nu)
    KT = K * U
    U[1:nu] = K[1:nu, 1:nu] \ (F[1:nu] - KT[1:nu])
end

function evaluate_pressure_error(Ph)
    function integrate!(elit, qpit, truep)
        pnedof = ndofsperel(elit)
        E = 0.0
        for el in elit
            dofvals = eldofvals(el)
            for qp in qpit
                Jac, J = jacjac(el, qp)
                JxW = J * weight(qp)
                Np = bfun(qp)
                pt = truep(location(el, qp)...)
                pa = 0.0
                for j in 1:pnedof
                    pa += (dofvals[j] * Np[j])
                end
                E += (JxW) * (pa - pt)^2
            end
        end
        return sqrt(E)
    end

    elit = FEIterator(Ph)
    qargs = (kind = :default, npts = 3,)
    qpit = QPIterator(Ph, qargs)
    return integrate!(elit, qpit, truep)
end

function evaluate_velocity_error(Uxh, Uyh)
    function integrate!(elits, qpits, trueux, trueuy)
        uxnedof, uynedof = ndofsperel.(elits)
        E = 0.0
        for el in zip(elits...)
            uxel, uyel = el
            uxdofvals = eldofvals(uxel)
            uydofvals = eldofvals(uyel)
            for qp in zip(qpits...)
                uxqp, uyqp = qp
                Jac, J = jacjac(uxel, uxqp)
                JxW = J * weight(uxqp)
                Np = bfun(uxqp)
                uxt = trueux(location(uxel, uxqp)...)
                uyt = trueuy(location(uyel, uyqp)...)
                uxa = 0.0
                uya = 0.0
                for j in 1:uxnedof
                    uxa += (uxdofvals[j] * Np[j])
                    uya += (uydofvals[j] * Np[j])
                end
                E += (JxW) * ((uxa - uxt)^2 + (uya - uyt)^2)
            end
        end
        return sqrt(E)
    end

    elits = (FEIterator(Uxh), FEIterator(Uyh),)
    qargs = (kind = :default, npts = 3,)
    qpits = (QPIterator(Uxh, qargs), QPIterator(Uyh, qargs),)
    return integrate!(elits, qpits, trueux, trueuy)
end

function run(N)
    mesh = genmesh(N)
    # Velocity spaces
    Uxh = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 1)
    Uyh = FESpace(Float64, mesh, FEH1_T3_BUBBLE(), 1)
    locs = geometry(mesh)
    inflate = A / N / 100
    # The entire boundary
    boxes = [[-A A -A -A], [-A -A -A A], [A A -A A], [-A A A A]]
    for box in boxes
        vl = vselect(locs; box = box, inflate = inflate)
        for i in vl
            setebc!(Uxh, 0, i, 1, trueux(locs[i]...))
            setebc!(Uyh, 0, i, 1, trueuy(locs[i]...))
        end
    end
    # Pressure space
    Ph = FESpace(Float64, mesh, FEH1_T3(), 1)
    atcenter = vselect(geometry(Ph.mesh); nearestto = [0.0, 0.0])
    setebc!(Ph, 0, atcenter[1], 1, 0.0)
    # Number the degrees of freedom
    numberdofs!([Uxh, Uyh, Ph])
    tndof = ndofs(Uxh) + ndofs(Uyh) + ndofs(Ph)
    tnunk = nunknowns(Uxh) + nunknowns(Uyh) + nunknowns(Ph)
    # Assemble the coefficient matrix
    K = assembleK(Uxh, Uyh, Ph, tndof, mu)
    # p = spy(K, canvas = DotCanvas)
    # display(p)
    # Solve the system
    U = fill(0.0, tndof)
    gathersysvec!(U, [Uxh, Uyh, Ph])
    F = fill(0.0, tndof)
    solve!(U, K, F, tnunk)
    scattersysvec!([Uxh, Uyh, Ph], U)
    # Postprocessing
    makeattribute(Ph, "p", 1)
    makeattribute(Uxh, "ux", 1)
    makeattribute(Uyh, "uy", 1)
    ep = evaluate_pressure_error(Ph)
        ev = evaluate_velocity_error(Uxh, Uyh)
    vtkwrite("p1b_p1-p", baseincrel(mesh), [(name = "p",), ])
    vtkwrite("p1b_p1-v", baseincrel(mesh), [(name = "ux",), (name = "uy",)])
    return (ep, ev)
end

end

let
    N = 4
    for loop in 1:5
        ep, ev = p1b_p1.run(N)
        @show N, ep, ev
        N = N * 2
    end
end


text:
export files, mdfiles, svgfiles, adocfiles

files(suffix) = filter(endswith(suffix), readdir())
mdfiles() = files(".md")
svgfiles() = files(".svg")
adocfiles() = files(".adoc")

ext: jl

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
abstract type BatchedSamplerCache{T} <: SamplerCache{T} end

# If constructing the cache for a BatchedNet, then build a batched sampler
#cache(s::Sampler, v::State, net::CachedNet{N,C}, par=NotParallel()) where {N,C<:NNBatchedCache} =
#    _batched_sampler_cache(s, v, net, par)

cache(s::Sampler, v, net::CachedNet{N,C}, par_cache=NotParallel()) where {N,C<:NNBatchedCache} =
    _batched_sampler_cache(s, v, net, par_cache)

_batched_sampler_cache(s, v, net, par) =
    _sampler_cache(s, v, net, par)


text:
module APIUtils

using ..CUDA

using LLVM
using LLVM.Interop

using Libdl

# helpers that facilitate working with CUDA APIs
include("call.jl")
include("enum.jl")
include("cache.jl")
include("threading.jl")
include("memoization.jl")

end


text:
module DateSelectors

using Base.Iterators
using Dates
using Intervals
using Random

export DateSelector, NoneSelector, PeriodicSelector, RandomSelector, partition

include("common.jl")
include("NoneSelector.jl")
include("PeriodicSelector.jl")
include("RandomSelector.jl")
end


text:
using GPUAcceleratedTracking, DrWatson, ProgressMeter
@quickactivate "GPUAcceleratedTracking"

allparams = Dict(
    "num_samples" => 2 .^ (11:15),
    "num_ants" => [4],
    "num_correlators" => [3],
    "algorithm" => ["pure", "cplx", "cplx_multi"]
)

dicts = dict_list(allparams)

@showprogress 0.5 "Benchmarking reduction algorithms (pure vs cplx vs cplx_multi)" for (_, d) in enumerate(dicts)
    benchmark_results = run_reduction_benchmark(d)
    @tagsave(
        datadir("benchmarks/reduction", savename("Reduction", d, "jld2")),
        benchmark_results
    )
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
struct ICPController{T, I, Z}
    kxy::T
    zgains::PDGains{T, T}
    m::T
    gz::T
    icptraj::I
    ztraj::Z
    contactmode::PlanarContactMode{T}
end

function ICPController(mechanism::Mechanism{T}, icptraj, z_des::Number) where T
    ztraj = t -> (z_des, zero(z_des), zero(z_des))
    ICPController(
        T(3.0),
        critically_damped_gains(10.0),
        mass(mechanism),
        norm(mechanism.gravitational_acceleration),
        icptraj,
        ztraj,
        PlanarContactMode{T}(root_frame(mechanism))
    )
end

function (controller::ICPController)(t, c::Point3D, cd::FreeVector3D, active_contact_points, state)
    update!(controller.contactmode, active_contact_points, state)
    support_polygon = controller.contactmode.hull

    # Equation (27) in "Design of a momentum-based control framework and application to the humanoid robot atlas"
    # minus the integral term.
    kxy = controller.kxy
    zgains = controller.zgains
    m = controller.m
    gz = controller.gz
    ξ_des_2d, ξd_des_2d = controller.icptraj(t)
    ξ_des = Point3D(c.frame, ξ_des_2d[1], ξ_des_2d[2], zero(eltype(ξ_des_2d))) # TODO: kind of nasty
    ξd_des = FreeVector3D(c.frame, ξd_des_2d[1], ξd_des_2d[2], zero(eltype(ξd_des_2d))) # TODO: kind of nasty
    z_des, zd_des, zdd_des = controller.ztraj(t)
    z = c.v[3]
    zd = cd.v[3]
    ω = sqrt(abs(gz / z))

    ξ = icp(c, cd, ω)
    cmp_des = ξ - ξd_des / ω + kxy * (ξ - ξ_des)

    cmp_des_projected_xy = closest_point(horizontal_projection(cmp_des.v), support_polygon)
    cmp_des_projected = Point3D(cmp_des.frame, cmp_des_projected_xy[1], cmp_des_projected_xy[2], cmp_des.v[3])

    ld_des_xy = (c - cmp_des_projected) * (m * gz) / z
    ld_des_z = m * (pd(zgains, z, z_des, zd, zd_des) + zdd_des)
    FreeVector3D(c.frame, ld_des_xy.v[1], ld_des_xy.v[2], ld_des_z)
end


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solvers
using GeometricIntegrators.Tableaus
using GeometricIntegrators.TestProblems.KuboOscillator
using Test

set_config(:nls_solver, NewtonSolver)

using GeometricIntegrators.TestProblems.KuboOscillator: Δt, nt

include("utils.jl")

sde1   = kubo_oscillator_sde_1()
sde2   = kubo_oscillator_sde_2()
sde3   = kubo_oscillator_sde_3()
psde1  = kubo_oscillator_psde_1()
psde2  = kubo_oscillator_psde_2()
psde3  = kubo_oscillator_psde_3()
spsde1  = kubo_oscillator_spsde_1()
spsde2  = kubo_oscillator_spsde_2()
spsde3  = kubo_oscillator_spsde_3()

@test typeof(Integrator(sde1, getTableauBurrageE1(), Δt)) <: IntegratorSERK
@test typeof(Integrator(sde1, getTableauStochasticGLRK(1), Δt)) <: IntegratorSIRK
@test typeof(Integrator(sde1, getTableauRosslerRS1(), Δt)) <: IntegratorWERK
@test typeof(Integrator(sde1, getTableauSRKw1(), Δt)) <: IntegratorWIRK
@test typeof(Integrator(psde1, getTableauStochasticStormerVerlet(), Δt)) <: IntegratorSIPRK
@test typeof(Integrator(spsde1, getTableauModifiedStochasticStormerVerlet(), Δt)) <: IntegratorSISPRK

### SERK Integrators ###
int = Integrator(sde1, getTableauBurrageE1(), Δt)
sol = Solution(sde1, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 2E-6

int = Integrator(sde2, getTableauBurrageE1(), Δt)
sol = Solution(sde2, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 2E-6

int = Integrator(sde3, getTableauBurrageE1(), Δt)
sol = Solution(sde3, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 2E-6

### SIRK Integrators ###
int = Integrator(sde1, getTableauStochasticGLRK(1), Δt)
sol = Solution(sde1, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

int = Integrator(sde2, getTableauStochasticGLRK(1), Δt)
sol = Solution(sde2, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

int = Integrator(sde3, getTableauStochasticGLRK(1), Δt)
sol = Solution(sde3, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

### WERK Integrators ###
int = Integrator(sde1, getTableauRosslerRS1(), Δt)
sol = Solution(sde1, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-5

int = Integrator(sde2, getTableauRosslerRS1(), Δt)
sol = Solution(sde2, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-5

int = Integrator(sde3, getTableauRosslerRS1(), Δt)
sol = Solution(sde3, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-5

### WIRK Integrators ###
int = Integrator(sde1, getTableauSRKw1(), Δt)
sol = Solution(sde1, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

int = Integrator(sde2, getTableauSRKw1(), Δt)
sol = Solution(sde2, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

int = Integrator(sde3, getTableauSRKw1(), Δt)
sol = Solution(sde3, Δt, nt, conv="weak")
integrate!(int, sol)

@test rel_energy_err_sde(sol) < 1E-14

### SIPRK Integrators ###
int = Integrator(psde1, getTableauStochasticStormerVerlet(), Δt)
sol = Solution(psde1, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 1E-5

int = Integrator(psde2, getTableauStochasticStormerVerlet(), Δt)
sol = Solution(psde2, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 1E-5

int = Integrator(psde3, getTableauStochasticStormerVerlet(), Δt)
sol = Solution(psde3, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 1E-5

### SISPRK Integrators ###
int = Integrator(spsde1, getTableauModifiedStochasticStormerVerlet(), Δt)
sol = Solution(spsde1, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 0.02

int = Integrator(spsde2, getTableauModifiedStochasticStormerVerlet(), Δt)
sol = Solution(spsde2, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 0.02

int = Integrator(spsde3, getTableauModifiedStochasticStormerVerlet(), Δt)
sol = Solution(spsde3, Δt, nt, conv="strong")
integrate!(int, sol)

@test rel_energy_err_psde(sol) < 0.02


text:
"""
    gsdec(sys; job = "finite", smarg, fast = true,  atol = 0,  atol1 = atol, atol2 = atol, rtol = nϵ) -> (sys1, sys2)

Compute for the descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix
`G(λ)`, the additive spectral decomposition `G(λ) = G1(λ) + G2(λ)` such that `G1(λ)`,
the transfer function matrix of the descriptor system `sys1 = (A1-λE1,B1,C1,D1)`,
has only poles in a certain domain of interest `Cg` of the complex plane and `G2(λ)`,
the transfer function matrix of the descriptor system `sys2 = (A2-λE2,B2,C2,0)`, has
only poles outside of `Cg`.

The keyword argument `smarg`, if provided, specifies the stability margin for the
stable eigenvalues of `A-λE`, such that, in the continuous-time case,
the stable eigenvalues have real parts less than or equal to `smarg`, and
in the discrete-time case, the stable eigenvalues
have moduli less than or equal to `smarg`. If `smarg = missing`, the used default values
are: `smarg = -sqrt(ϵ)`, for a continuous-time system, and `smarg = 1-sqrt(ϵ)`,
for a discrete-time system), where `ϵ` is the machine precision of the working accuracy.

The keyword argument `job`, in conjunction with `smarg`, defines the domain of
interest `Cg`, as follows:

for `job = "finite"`, `Cg` is the whole complex plane without the point at infinity, and
   `sys1` has only finite poles and `sys2` has only infinite poles (default);
   the resulting `A2` is nonsingular and upper triangular, while the
   resulting `E2` is nilpotent and upper triangular;

for `job = "infinite"`, `Cg` is the point at infinity, and
   `sys1` has only infinite poles and `sys2` has only finite poles and
   is the strictly proper part of `sys`;
   the resulting `A1` is nonsingular and upper triangular, while the
   resulting `E1` is nilpotent and upper triangular;

for `job = "stable"`, `Cg` is the stability domain of eigenvalues defined by `smarg`, and
    `sys1` has only stable poles and `sys2` has only unstable and infinite poles;
    the resulting pairs `(A1,E1)` and `(A2,E2)` are in generalized Schur form with
    `E1` upper triangular and nonsingular and `E2` upper triangular;

for `job = "unstable"`, `Cg` is the complement of the stability domain of the
    eigenvalues defined by `smarg`, and
    `sys1` has only unstable and infinite poles and `sys2` has only stable poles;
    the resulting pairs `(A1,E1)` and `(A2,E2)` are in generalized Schur form with
    `E1` upper triangular and `E2` upper triangular  and nonsingular.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `A`, the absolute tolerance for the nonzero elements of `E`,  and the relative tolerance
for the nonzero elements of `A` and `E`. The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the order of the system `sys`. The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol`.

The separation of the finite and infinite eigenvalues is performed using
rank decisions based on rank revealing QR-decompositions with column pivoting
if `fast = true` or the more reliable SVD-decompositions if `fast = false`.
"""
function gsdec(SYS::DescriptorStateSpace{T}; job::String = "finite", smarg::Union{Real,Missing} = missing,
               fast::Bool = true,  atol::Real = zero(real(T)),  atol1::Real = atol, atol2::Real = atol,
               rtol::Real = (SYS.nx*eps(real(float(one(T)))))*iszero(min(atol1,atol2))) where T
    disc = !iszero(SYS.Ts)
    if SYS.E == I
       if job == "finite"
          return SYS, dss(zeros(T,SYS.ny,SYS.nu), Ts = SYS.Ts)
       elseif job == "infinite"
          return dss(SYS.D, Ts = SYS.Ts),
                 dss(SYS.A, SYS.B, SYS.C, zeros(T,SYS.ny,SYS.nu), Ts = SYS.Ts)
       elseif job == "stable" || job == "unstable"
          stable_unstable = (job == "stable")
          A, B, C, _, _, blkdims, = ssblkdiag(SYS.A, SYS.B, SYS.C; smarg = smarg, disc = disc, stable_unstable = stable_unstable,
                                              withQ = false, withZ = false)
          n1 = blkdims[1];
          i1 = 1:n1; i2 = n1+1:SYS.nx
          return dss(A[i1,i1], B[i1,:], C[:,i1], SYS.D, Ts = SYS.Ts),
                 dss(A[i2,i2], B[i2,:], C[:,i2], zeros(T,SYS.ny,SYS.nu), Ts = SYS.Ts)
       else
          error("No such job option")
       end
    else
        if job == "finite"
           A, E, B, C, _, _, _, blkdims = fiblkdiag(SYS.A, SYS.E, SYS.B, SYS.C; fast = fast, finite_infinite = true, trinv = false,
                                                    atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
           n1 = blkdims[1];
        elseif job == "infinite"
           A, E, B, C, _, _, _, blkdims = fiblkdiag(SYS.A, SYS.E, SYS.B, SYS.C; fast = fast, finite_infinite = false, trinv = false,
                                                    atol1 = atol1, atol2 = atol2, rtol= rtol, withQ = false, withZ = false)
           n1 = blkdims[1];
        elseif job == "stable"
           A, E, B, C, _, _, _, blkdims, = gsblkdiag(SYS.A, SYS.E, SYS.B, SYS.C; smarg = smarg, disc = disc, fast = fast,
                                                      finite_infinite = true, stable_unstable = true,
                                                      atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
           n1 = blkdims[1];
        elseif job == "unstable"
           A, E, B, C, _, _, _, blkdims, = gsblkdiag(SYS.A, SYS.E, SYS.B, SYS.C; smarg = smarg, disc = disc, fast = fast,
                                                      finite_infinite = false, stable_unstable = false,
                                                      atol1 = atol1, atol2 = atol2, rtol = rtol, withQ = false, withZ = false)
           n1 = blkdims[1]+blkdims[2];
        else
            error("No such job option")
        end
        i1 = 1:n1; i2 = n1+1:SYS.nx
        return dss(A[i1,i1], E[i1,i1], B[i1,:], C[:,i1], SYS.D, Ts = SYS.Ts),
               dss(A[i2,i2], E[i2,i2], B[i2,:], C[:,i2], zeros(T,size(SYS.D)...), Ts = SYS.Ts)
    end
end


text:
module TestArray
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, leveltype

@testset "conversion ordered=$ordered reftype=$R" for ordered in (false, true),
    R in (CategoricalArrays.DefaultRefType, UInt8, UInt, Int8, Int)
    # Vector
    a = ["b", "a", "b"]
    x = CategoricalVector{String, R}(a, ordered=ordered)

    @test x == a
    @test leveltype(typeof(x)) === String
    @test leveltype(x) === String
    @test eltype(x) === CategoricalValue{String, R}
    @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{String}, x) === x
    @test convert(CategoricalArray{String, 1}, x) === x
    @test convert(CategoricalArray{String, 1, R}, x) === x
    @test convert(CategoricalArray{String, 1, DefaultRefType}, x) == x
    @test convert(CategoricalArray{String, 1, UInt8}, x) == x

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

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

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

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

    @test CategoricalVector(a, ordered=ordered) == x
    @test CategoricalVector{String}(a, ordered=ordered) == x
    @test CategoricalVector{String, R}(a, ordered=ordered) == x
    @test CategoricalVector{String, DefaultRefType}(a, ordered=ordered) == x
    @test CategoricalVector{String, UInt8}(a, ordered=ordered) == 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{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 CategoricalArray{Union{String, Missing}}(x, ordered=ordered) == x
    @test CategoricalArray{Union{String, Missing}, 1}(x, ordered=ordered) == x
    @test CategoricalArray{Union{String, Missing}, 1, R}(x, ordered=ordered) == x
    @test CategoricalArray{Union{String, Missing}, 1, DefaultRefType}(x, ordered=ordered) == x
    @test CategoricalArray{Union{String, Missing}, 1, UInt8}(x, ordered=ordered) == x

    @test CategoricalVector{Union{String, Missing}}(x, ordered=ordered) == x
    @test CategoricalVector{Union{String, Missing}, R}(x, ordered=ordered) == x
    @test CategoricalVector{Union{String, Missing}, DefaultRefType}(x, ordered=ordered) == x
    @test CategoricalVector{Union{String, Missing}, UInt8}(x, 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 == x
        @test isa(x2, CategoricalVector{String, R1})
        @test isordered(x2) === ordered
        @test leveltype(x2) === String
        @test eltype(x2) === CategoricalValue{String, R1}

        x2 = categorical(y, compress=comp, ordered=ordered)
        @test x2 == x
        @test isa(x2, CategoricalVector{String, R2})
        @test isordered(x2) === ordered
        @test leveltype(x2) === String
        @test eltype(x2) === CategoricalValue{String, R2}
    end

    x2 = compress(x)
    @test x2 == x
    @test isa(x2, CategoricalVector{String, 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)

    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] === CategoricalValue(x.pool, 2)
    @test x[2] === CategoricalValue(x.pool, 1)
    @test x[3] === CategoricalValue(x.pool, 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[1:2]
    @test typeof(x2) === typeof(x)
    @test x2 == ["b", "a"]
    @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] === CategoricalValue(x.pool, 1)
    @test x[2] === CategoricalValue(x.pool, 1)
    @test x[3] === CategoricalValue(x.pool, 2)

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

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

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

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

    @test_throws ArgumentError levels!(x, ["a"])
    # check that x is restored correctly when dropping levels is not allowed
    @test x == ["a", "b", "b"]
    @test levels(x) == ["b", "a"]

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

    @test_throws ArgumentError levels!(x, ["e", "a", "b", "a"])
    # once again check that x is restored correctly when dropping levels is not allowed
    @test x == ["a", "b", "b"]
    @test levels(x) == ["b", "a"]

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

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

    y = copy(x)

    push!(x, "a")
    @test length(x) == 4
    @test x[end] == "a"
    @test levels(x) == ["e", "a", "b", "c"]

    push!(x, "zz")
    @test length(x) == 5
    @test x[end] == "zz"
    @test levels(x) == ["e", "a", "b", "c", "zz"]

    push!(x, x[1])
    @test length(x) == 6
    @test x[1] == x[end]
    @test levels(x) == ["e", "a", "b", "c", "zz"]

    x2 = copy(x)
    @test_throws MethodError push!(x, 1)
    @test x == x2
    @test x.pool.levels == x2.pool.levels
    @test x.pool.invindex == x2.pool.invindex

    x = y
    insert!(x, 1, "a")
    @test length(x) == 4
    @test x[1] == "a"
    @test levels(x) == ["e", "a", "b", "c"]

    insert!(x, length(x) + 1, "zz")
    @test length(x) == 5
    @test x[end] == "zz"
    @test levels(x) == ["e", "a", "b", "c", "zz"]

    insert!(x, 2, x[1])
    @test length(x) == 6
    @test x[1] == x[2]
    @test levels(x) == ["e", "a", "b", "c", "zz"]

    x2 = copy(x)
    @test_throws MethodError insert!(x, 1, 1)
    @test_throws ArgumentError insert!(x, true, "a")
    @test_throws BoundsError insert!(x, 0, "a")
    @test_throws BoundsError insert!(x, 100, "a")
    @test x == x2
    @test x.pool.levels == x2.pool.levels
    @test x.pool.invindex == x2.pool.invindex

    empty!(x)
    @test length(x) == 0
    @test levels(x) == ["e", "a", "b", "c", "zz"]

    @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{Float64, 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) <: CategoricalValue{Float64}

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

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

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

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

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

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

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

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

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

        @test CategoricalVector{Union{Float64, Missing}}(x, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}}(x, ordered=ordered) == x
        @test CategoricalVector{Union{Float64, Missing}, R}(x, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}, R}(x, ordered=ordered) == x
        @test CategoricalVector{Union{Float64, Missing}, DefaultRefType}(x, ordered=ordered) == x
        @test CategoricalVector{Union{Float32, Missing}, DefaultRefType}(x, 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 == x
            @test isa(x2, CategoricalVector{Float64, R1})
            @test isordered(x2) === ordered
            @test leveltype(x2) === Float64
            @test eltype(x2) === CategoricalValue{Float64, R1}

            x2 = categorical(y, compress=comp, ordered=ordered)
            @test x2 == x
            @test isa(x2, CategoricalVector{Float64, R2})
            @test isordered(x2) === ordered
            @test leveltype(x2) === Float64
            @test eltype(x2) === CategoricalValue{Float64, R2}
        end

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

        @test x[1] === CategoricalValue(x.pool, 1)
        @test x[2] === CategoricalValue(x.pool, 2)
        @test x[3] === CategoricalValue(x.pool, 3)
        @test x[4] === CategoricalValue(x.pool, 4)
        @test_throws BoundsError x[5]

        x2 = x[:]
        @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: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] === CategoricalValue(x.pool, 1)
        @test x[2] === CategoricalValue(x.pool, 3)
        @test x[3] === CategoricalValue(x.pool, 3)
        @test x[4] === CategoricalValue(x.pool, 4)
        @test levels(x) == unique(a)
        @test unique(x) == unique(collect(x))

        x[1:2] .= -1
        @test x[1] === CategoricalValue(x.pool, 5)
        @test x[2] === CategoricalValue(x.pool, 5)
        @test x[3] === CategoricalValue(x.pool, 3)
        @test x[4] === CategoricalValue(x.pool, 4)
        @test levels(x) == vcat(unique(a), -1)
        @test unique(x) == unique(collect(x))

        push!(x, 2.0)
        @test length(x) == 5
        @test x[end] == 2.0
        @test isordered(x) === false
        @test levels(x) == [0.0,  0.5,  1.0,  1.5, -1.0,  2.0]

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

        ordered!(x, ordered)
        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" begin
        a = ["a" "b" "c"; "b" "a" "c"]
        x = CategoricalMatrix{String, 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{String}, x) === x
        @test convert(CategoricalArray{String, 2}, x) === x
        @test convert(CategoricalArray{String, 2, R}, x) === x
        @test convert(CategoricalArray{String, 2, DefaultRefType}, x) == x
        @test convert(CategoricalArray{String, 2, UInt8}, x) == x

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

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

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

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

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

        @test convert(CategoricalArray{Union{String, Missing}}, 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{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 CategoricalArray{Union{String, Missing}}(x, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2}(x, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2}(x, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, R}(x, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, DefaultRefType}(x, ordered=ordered) == x
        @test CategoricalArray{Union{String, Missing}, 2, UInt8}(x, ordered=ordered) == x

        @test CategoricalMatrix{Union{String, Missing}}(x, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, R}(x, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, DefaultRefType}(x, ordered=ordered) == x
        @test CategoricalMatrix{Union{String, Missing}, UInt8}(x, 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 == x
            @test isa(x2, CategoricalMatrix{String, R1})
            @test isordered(x2) === ordered

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

        x2 = compress(x)
        @test x2 == x
        @test isa(x2, CategoricalMatrix{String, 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] === CategoricalValue(x.pool, 1)
        @test x[2] === CategoricalValue(x.pool, 2)
        @test x[3] === CategoricalValue(x.pool, 2)
        @test x[4] === CategoricalValue(x.pool, 1)
        @test x[5] === CategoricalValue(x.pool, 3)
        @test x[6] === CategoricalValue(x.pool, 3)
        @test_throws BoundsError x[7]

        @test x[1,1] === CategoricalValue(x.pool, 1)
        @test x[2,1] === CategoricalValue(x.pool, 2)
        @test x[1,2] === CategoricalValue(x.pool, 2)
        @test x[2,2] === CategoricalValue(x.pool, 1)
        @test x[1,3] === CategoricalValue(x.pool, 3)
        @test x[2,3] === CategoricalValue(x.pool, 3)
        @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" "c"]
        @test levels(x2) == levels(x)
        @test levels(x2) !== levels(x)
        @test isordered(x2) == isordered(x)

        x2 = x[1:1,2]
        @test isa(x2, CategoricalVector{String, R})
        @test x2 == ["b"]
        @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, :]

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

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

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

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

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

    @testset "compress($(typeof(x))) and setindex!()" for x in v
        @test !isassigned(x, 1) && isdefined(x, 1)
        @test !isassigned(x, 2) && isdefined(x, 2)
        @test_throws UndefRefError x[1]
        @test_throws UndefRefError x[2]
        @test isordered(x) === ordered
        @test levels(x) == []

        x2 = compress(x)
        @test isa(x2, CategoricalArray{leveltype(x), ndims(x), UInt8})
        @test !isassigned(x2, 1) && isdefined(x2, 1)
        @test !isassigned(x2, 2) && isdefined(x2, 2)
        @test_throws UndefRefError x2[1]
        @test_throws UndefRefError x2[2]
        @test levels(x2) == []

        x[1] = "c"
        @test x[1] === CategoricalValue(x.pool, 1)
        @test !isassigned(x, 2) && isdefined(x, 2)
        @test_throws UndefRefError x[2]
        @test isordered(x) === false
        @test levels(x) == ["c"]

        x[1] = "a"
        @test x[1] === CategoricalValue(x.pool, 2)
        @test !isassigned(x, 2) && isdefined(x, 2)
        @test_throws UndefRefError x[2]
        @test isordered(x) === false
        @test levels(x) == ["c", "a"]

        x[2] = "c"
        @test x[1] === CategoricalValue(x.pool, 2)
        @test x[2] === CategoricalValue(x.pool, 1)
        @test isordered(x) === false
        @test levels(x) == ["c", "a"]

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

        ordered!(x, ordered)
        v = CategoricalValue(2, CategoricalPool(["xyz", "b"]))
        x[1] = v
        @test x[1] === CategoricalValue(x.pool, 4)
        @test x[2] === CategoricalValue(x.pool, 1)
        @test isordered(x) === false
        @test levels(x) == ["c", "a", "xyz", "b"]
    end
end

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

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

    # To test short-circuiting
    x = CategoricalArray(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
    @test levels(x) == [19:-1:1; 20]
    @test unique(x) == collect(1:10)
end

end


text:
module TestExamplesP4estMesh3D

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "p4est_3d_dgsem")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "P4estMesh3D" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [0.00016263963870641478],
      linf = [0.0014537194925779984])
  end

  @trixi_testset "elixir_advection_unstructured_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_curved.jl"),
      l2   = [0.0004750004258546538],
      linf = [0.026527551737137167])
  end

  @trixi_testset "elixir_advection_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonconforming.jl"),
      l2   = [0.00253595715323843],
      linf = [0.016486952252155795])
  end

  @trixi_testset "elixir_advection_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr.jl"),
      # Expected errors are exactly the same as with TreeMesh!
      l2   = [9.773852895157622e-6],
      linf = [0.0005853874124926162],
      coverage_override = (maxiters=6, initial_refinement_level=1, base_level=1, med_level=2, max_level=3))
  end

  @trixi_testset "elixir_advection_amr_unstructured_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_amr_unstructured_curved.jl"),
      l2   = [1.6236411810065552e-5],
      linf = [0.0010554006923731395],
      tspan = (0.0, 1.0),
      coverage_override = (maxiters=6, initial_refinement_level=0, base_level=0, med_level=1, max_level=2))
  end

  @trixi_testset "elixir_advection_cubed_sphere.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_cubed_sphere.jl"),
      l2   = [0.002006918015656413],
      linf = [0.027655117058380085])
  end

  @trixi_testset "elixir_advection_restart.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
      l2   = [0.002590388934758452],
      linf = [0.01840757696885409])
  end

  @trixi_testset "elixir_euler_source_terms_nonconforming_unstructured_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_unstructured_curved.jl"),
      l2   = [4.070355207909268e-5, 4.4993257426833716e-5, 5.10588457841744e-5, 5.102840924036687e-5, 0.00019986264001630542],
      linf = [0.0016987332417202072, 0.003622956808262634, 0.002029576258317789, 0.0024206977281964193, 0.008526972236273522],
      tspan = (0.0, 0.01))
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [0.0015106060984283647, 0.0014733349038567685, 0.00147333490385685, 0.001473334903856929, 0.0028149479453087093],
      linf = [0.008070806335238156, 0.009007245083113125, 0.009007245083121784, 0.009007245083102688, 0.01562861968368434],
      tspan = (0.0, 1.0))
  end

  @trixi_testset "elixir_euler_free_stream.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream.jl"),
      l2   = [5.162664597942288e-15, 1.941857343642486e-14, 2.0232366394187278e-14, 2.3381518645408552e-14, 7.083114561232324e-14],
      linf = [7.269740365245525e-13, 3.289868377720495e-12, 4.440087186807773e-12, 3.8686831516088205e-12, 9.412914891981927e-12],
      tspan = (0.0, 0.03))
  end

  @trixi_testset "elixir_euler_free_stream_extruded.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
      l2   = [8.444868392439035e-16, 4.889826056731442e-15, 2.2921260987087585e-15, 4.268460455702414e-15, 1.1356712092620279e-14],
      linf = [7.749356711883593e-14, 2.8792246364872653e-13, 1.1121659149182506e-13, 3.3228975127030935e-13, 9.592326932761353e-13],
      tspan=(0.0, 0.1))
  end

  @trixi_testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.010380390326164493, 0.006192950051354618, 0.005970674274073704, 0.005965831290564327, 0.02628875593094754],
      linf = [0.3326911600075694, 0.2824952141320467, 0.41401037398065543, 0.45574161423218573, 0.8099577682187109],
      tspan = (0.0, 0.2),
      coverage_override = (polydeg=3,)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_sedov.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov.jl"),
      l2   = [7.82070951e-02, 4.33260474e-02, 4.33260474e-02, 4.33260474e-02, 3.75260911e-01],
      linf = [7.45329845e-01, 3.21754792e-01, 3.21754792e-01, 3.21754792e-01, 4.76151527e+00],
      tspan = (0.0, 0.3),
      coverage_override = (polydeg=3,)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_source_terms_nonconforming_earth.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonconforming_earth.jl"),
      l2 = [6.040180337738628e-6, 5.4254175153621895e-6, 5.677698851333843e-6, 5.8017136892469794e-6, 1.3637854615117974e-5],
      linf = [0.00013996924184311865, 0.00013681539559939893, 0.00013681539539733834, 0.00013681539541021692, 0.00016833038543762058],
      # Decrease tolerance of adaptive time stepping to get similar results across different systems
      abstol=1.0e-11, reltol=1.0e-11,
      coverage_override = (trees_per_cube_face=(1, 1), polydeg=3)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_circular_wind_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_circular_wind_nonconforming.jl"),
      l2   = [1.573832094977477e-7, 3.863090659429634e-5, 3.867293305754584e-5, 3.686550296950078e-5, 0.05508968493733932],
      linf = [2.2695202613887133e-6, 0.0005314968179916946, 0.0005314969614147458, 0.0005130280733059617, 0.7944959432352334],
      tspan = (0.0, 2e2),
      coverage_override = (trees_per_cube_face=(1, 1), polydeg=3)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_baroclinic_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_baroclinic_instability.jl"),
      l2   = [6.725065410642336e-7, 0.00021710117340245454, 0.000438679759422352, 0.00020836356588024185, 0.07602006689579247],
      linf = [1.9101671995258585e-5, 0.029803626911022396, 0.04847630924006063, 0.022001371349740104, 4.847761006938526],
      tspan = (0.0, 1e2),
      # Decrease tolerance of adaptive time stepping to get similar results across different systems
      abstol=1.0e-9, reltol=1.0e-9,
      coverage_override = (trees_per_cube_face=(1, 1), polydeg=3)) # Prevent long compile time in CI
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic_hohqmesh.jl"),
      l2   = [0.0042023406458005464, 0.004122532789279737, 0.0042448149597303616, 0.0036361316700401765, 0.007389845952982495],
      linf = [0.04530610539892499, 0.02765695110527666, 0.05670295599308606, 0.048396544302230504, 0.1154589758186293])
  end

  @trixi_testset "elixir_mhd_alfven_wave_nonconforming.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_nonconforming.jl"),
      l2   = [0.00019018725889431733, 0.0006523517707148006, 0.0002401595437705759, 0.0007796920661427565,
              0.0007095787460334334, 0.0006558819731628876, 0.0003565026134076906, 0.0007904654548841712,
              9.437300326448332e-7],
      linf = [0.0012482306861187897, 0.006408776208178299, 0.0016845452099629663, 0.0068711236542984555,
              0.004626581522263695, 0.006614624811393632, 0.0030068344747734566, 0.008277825749754025,
              1.3475027166309006e-5],
      tspan = (0.0, 0.25),
      coverage_override = (trees_per_dimension=(1, 1, 1),))
  end

  @trixi_testset "elixir_mhd_shockcapturing_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_shockcapturing_amr.jl"),
      l2   = [0.006298541670176575, 0.0064368506652601265, 0.007108729762852636, 0.006530420607206385,
              0.02061185869237284, 0.005562033787605515, 0.007571716276627825, 0.005571862660453231,
              3.909755063709152e-6],
      linf = [0.20904054009050665, 0.18622917151105936, 0.2347957890323218, 0.19432508025509926,
              0.6858860133405615, 0.15172116633332622, 0.22432820727833747, 0.16805989780225183,
              0.000535219040687628],
      tspan = (0.0, 0.04),
      coverage_override = (maxiters=6, initial_refinement_level=1, base_level=1, max_level=2))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end # module


text:
#!/usr/bin/env julia

@static if VERSION >= v"0.7.0-DEV.3406"
    import Random
    srand(seed) = Random.seed!(seed)
end

import GR

function main()
    srand(0);
    xd = -2 .+ 4 * rand(100)
    yd = -2 .+ 4 * rand(100)
    zd = [xd[i] * exp(-xd[i]^2 - yd[i]^2) for i = 1:100]

    GR.setviewport(0.1, 0.95, 0.1, 0.95)
    GR.setwindow(-2, 2, -2, 2)
    GR.setmarkersize(1)
    GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
    GR.setcharheight(0.024)
    GR.settextalign(2, 0)
    GR.settextfontprec(3, 0)

    x, y, z = GR.gridit(xd, yd, zd, 200, 200)
    h = -0.6:0.05:0.6
    GR.contourf(x, y, h, z, 2)
    GR.polymarker(xd, yd)
    GR.axes(0.25, 0.25, -2, -2, 2, 2, 0.01)

    GR.updatews()
end

main()


text:
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, -5 ≤ x ≤ 10)
@disjunction(
    m,
    0 ≤ x ≤ 3,
    5 ≤ x ≤ 9,
    reformulation=:big_m,
    name=:y
)
@proposition(m, y[1] ∨ y[2]) #this is a redundant proposition

print(m)

# ┌ Warning: disj_y[1] : x in [0.0, 3.0] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
# ┌ Warning: disj_y[2] : x in [5.0, 9.0] uses the `MOI.Interval` set. Each instance of the interval set has been split into two constraints, one for each bound.
# Feasibility
# Subject to
#  XOR(disj_y) : y[1] + y[2] == 1.0         <- XOR constraint
#  y[1] ∨ y[2] : y[1] + y[2] >= 1.0         <- reformulated logical proposition (name is the proposition)
#  disj_y[1,lb] : -x + 5 y[1] <= 5.0       <- left-side of constraint in 1st disjunct (name is assigned to disj_y[1][lb])
#  disj_y[1,ub] : x + 7 y[1] <= 10.0       <- right-side of constraint in 1st disjunct (name is assigned to disj_y[1][ub])
#  disj_y[2,lb] : -x + 10 y[2] <= 5.0      <- left-side of constraint in 2nd disjunct (name is assigned to disj_y[2][lb])
#  disj_y[2,ub] : x + y[2] <= 10.0         <- right-side of constraint in 2nd disjunct (name is assigned to disj_y[2][ub])
#  x >= -5.0                                <- variable lower bound
#  x <= 10.0                                <- variable upper bound
#  y[1] >= 0.0                              <- lower bound on binary
#  y[2] >= 0.0                              <- lower bound on binary
#  y[1] <= 1.0                              <- upper bound on binary
#  y[2] <= 1.0                              <- upper bound on binary
#  y[1] binary                              <- indicator variable (1st disjunct) is binary
#  y[2] binary                              <- indicator variable (2nd disjunct) is binary


text:
function _copy_values(target::Dict{VariableID,Float64},
                      source::Dict{VariableID,Float64}
                      )::Nothing

    for (vid, value) in pairs(source)
        target[vid] = value
    end

    return

end

function _execute_callbacks(phd::PHData, winf::WorkerInf, niter::Int)::Bool

    running = true

    for cb in phd.callbacks
        running &= cb.h(cb.ext, phd, winf, niter)
    end

    return running

end

function _process_reports(phd::PHData,
                          winf::WorkerInf,
                          report_type::Type{R}
                          )::Nothing where R <: Report

    waiting_on = copy(scenarios(phd))

    while !isempty(waiting_on)

        msg = _retrieve_message_type(winf, report_type)
        _verify_report(msg)
        _update_values(phd, msg)

        delete!(waiting_on, msg.scen)
    end

    return

end

function _report(niter::Int,
                 nsqrt::Float64,
                 residual::Float64,
                 xmax::Float64,
                 xhat_sq::Float64,
                 x_sq::Float64
                 )::Nothing

    @printf("Iter: %4d   AbsR: %12.6e   RelR: %12.6e   Xhat: %12.6e   X: %12.6e\n",
            niter, residual, residual / xmax,
            sqrt(xhat_sq)/nsqrt, sqrt(x_sq)/nsqrt
            )
    flush(stdout)

    return

end

function _report_lower_bound(niter::Int, bound::Float64, gap::Float64)::Nothing

    @printf("Iter: %4d   Bound: %12.4e   Abs Gap: %12.4e   Rel Gap: %8.4g\n",
            niter, bound, gap, abs(gap/bound)
            )
    flush(stdout)

    return

end

function _save_iterate(phd::PHData,
                       iter::Int,
                       )::Nothing

    xhd = Dict{XhatID,Float64}()
    xd = Dict{VariableID,Float64}()
    wd = Dict{VariableID,Float64}()

    for (xhid, xhat) in pairs(phd.xhat)
        xhd[xhid] = value(xhat)
        for vid in variables(xhat)
            xd[vid] = branch_value(phd, vid)
            wd[vid] = w_value(phd, vid)
        end
    end

    _save_iterate(phd.history,
                  iter,
                  PHIterate(xhd, xd, wd)
                  )

    return
end

function _save_residual(phd::PHData,
                        iter::Int,
                        xhat_sq::Float64,
                        x_sq::Float64,
                        absr::Float64,
                        relr::Float64,
                        )::Nothing
    _save_residual(phd.history, iter, PHResidual(absr, relr, xhat_sq, x_sq))
    return
end

function _send_solve_commands(phd::PHData,
                              winf::WorkerInf,
                              niter::Int
                              )::Nothing

    @sync for (scen, sinfo) in pairs(phd.scenario_map)

        (w_dict, xhat_dict) = create_ph_dicts(sinfo)

        @async _send_message(winf,
                             sinfo.pid,
                             Solve(scen,
                                   w_dict,
                                   xhat_dict,
                                   niter
                                   )
                             )

    end

    return

end

function _send_lb_solve_commands(phd::PHData,
                                 winf::WorkerInf,
                                 )::Nothing

    @sync for (scen, sinfo) in pairs(phd.scenario_map)
        (w_dict, xhat_dict) = create_ph_dicts(sinfo)
        @async _send_message(winf, sinfo.pid, SolveLowerBound(scen, w_dict))
    end

    return
end

function _update_si_xhat(phd::PHData)::Nothing

    for (xhid, xhat) in pairs(phd.xhat)
        for vid in variables(xhat)
            sinfo = phd.scenario_map[scenario(vid)]
            sinfo.xhat_vars[vid] = value(xhat)
        end
    end

    return
end

function _update_values(phd::PHData,
                        msg::ReportBranch
                        )::Nothing

    sinfo = phd.scenario_map[msg.scen]

    pd = sinfo.problem_data
    pd.obj = msg.obj
    pd.sts = msg.sts
    pd.time = msg.time

    _copy_values(sinfo.branch_vars, msg.vals)

    return
end

function _update_values(phd::PHData,
                        msg::ReportLeaf
                        )::Nothing
    _copy_values(phd.scenario_map[msg.scen].leaf_vars, msg.vals)
    return
end

function _update_values(phd::PHData,
                        msg::ReportLowerBound
                        )::Nothing

    pd = phd.scenario_map[msg.scen].problem_data

    pd.lb_obj = msg.obj
    pd.lb_sts = msg.sts
    pd.lb_time = msg.time

    return

end

function _verify_report(msg::ReportBranch)::Nothing

    if (msg.sts != MOI.OPTIMAL &&
        msg.sts != MOI.LOCALLY_SOLVED &&
        msg.sts != MOI.ALMOST_LOCALLY_SOLVED)
        # TODO: Create user adjustable/definable behavior for when this happens
        @error("Scenario $(msg.scen) subproblem returned $(msg.sts).")
    end

    return

end

function _verify_report(msg::ReportLeaf)::Nothing
    # Intentional no-op
    return
end

function _verify_report(msg::ReportLowerBound)::Nothing

    if (msg.sts != MOI.OPTIMAL &&
        msg.sts != MOI.LOCALLY_SOLVED &&
        msg.sts != MOI.ALMOST_LOCALLY_SOLVED)
        # TODO: Create user adjustable/definable behavior for when this happens
        @error("Scenario $(msg.scen) lower-bound subproblem returned $(msg.sts).")
    end

    return

end

function compute_and_save_xhat(phd::PHData)::Float64

    xhat_res = 0.0

    for (xhid, xhat_var) in pairs(phd.xhat)

        xhat = 0.0
        norm = 0.0

        for vid in variables(xhat_var)

            s = scenario(vid)
            p = phd.scenario_map[s].prob
            x = branch_value(phd, vid)

            xhat += p * x
            norm += p

        end

        xhat_new = xhat / norm
        xhat_old = value(xhat_var)

        xhat_var.value = xhat_new
        xhat_res += (xhat_new - xhat_old)^2
    end

    return xhat_res

end

function compute_and_save_w(phd::PHData)::Float64

    kxsq = 0.0

    for (xhid, xhat_var) in pairs(phd.xhat)

        xhat = value(xhat_var)

        exp = 0.0
        norm = 0.0

        for vid in variables(xhat_var)

            s = scenario(vid)
            p = phd.scenario_map[s].prob
            kx = branch_value(phd, vid) - xhat
            r = get_penalty_value(phd.r, xhid)

            phd.scenario_map[s].w_vars[vid] += r * kx

            kxsq += p * kx^2

            exp += p * w_value(phd, vid)
            norm += p

        end

        if abs(exp) > 1e-6
            vname = name(phd, xhid)
            @warn("Conditional expectation of " *
                  "W[$(vname)] occurring in stage " * stage_id(phd, xhid) *
                  " for scenarios " * stringify(scenario_bundle(phd, xhid)) *
                  " is non-zero: " * string(exp/norm))
        end
    end

    return kxsq

end

function compute_gap(phd::PHData)::NTuple{2,Float64}

    gap = 0.0
    lb = 0.0

    for s in scenarios(phd)

        sinfo = phd.scenario_map[s]
        pd = sinfo.problem_data
        p = sinfo.prob

        gap += p * abs(pd.obj - pd.lb_obj)
        lb += p * pd.lb_obj

    end

    return (lb, gap)

end

function finish(phd::PHData,
                winf::WorkerInf,
                )::Nothing

    # Send shutdown command to workers
    _shutdown(winf)

    # Wait for and process replies
    _process_reports(phd, winf, ReportLeaf)

    # Create hat variables for all leaf variables so
    # that `retrieve_soln` picks up the values
    for (scid, sinfo) in pairs(phd.scenario_map)
        for (vid, vinfo) in pairs(sinfo.leaf_vars)
            xhid = convert_to_xhat_id(phd, vid)
            @assert(!haskey(phd.xhat, xhid))
            # TODO: Correctly assign whether leaf variables are integers. For now,
            # just say they aren't because it doesn't really matter.
            phd.xhat[xhid] = HatVariable(value(phd, vid), vid, false)
        end
    end

    # Wait for workers to exit cleanly
    _wait_for_shutdown(winf)

    return

end

function update_gap(phd::PHData, winf::WorkerInf, niter::Int)::NTuple{2,Float64}

    @timeit(phd.time_info,
            "Issuing lower bound solve commands",
            _send_lb_solve_commands(phd, winf)
            )

    @timeit(phd.time_info,
            "Collecting lower bound results",
            _process_reports(phd, winf, ReportLowerBound)
            )

    (lb, gap) = @timeit(phd.time_info,
                        "Computing gap",
                        compute_gap(phd)
                        )

    _save_lower_bound(phd.history,
                      niter,
                      PHLowerBound(lb, gap, abs(gap/lb))
                      )

    return (lb, gap)

end

function update_ph_variables(phd::PHData)::NTuple{2,Float64}

    xhat_sq = compute_and_save_xhat(phd)
    x_sq = compute_and_save_w(phd)

    return (xhat_sq, x_sq)

end

function solve_subproblems(phd::PHData,
                           winf::WorkerInf,
                           niter::Int,
                           )::Nothing

    # Copy hat values to scenario base structure for easy dispersal
    @timeit(phd.time_info,
            "Other",
            _update_si_xhat(phd))

    # Send solve command to workers
    @timeit(phd.time_info,
            "Issuing solve commands",
            _send_solve_commands(phd, winf, niter))

    # Wait for and process replies
    @timeit(phd.time_info,
            "Collecting results",
            _process_reports(phd, winf, ReportBranch))

    return

end

function hedge(ph_data::PHData,
               worker_inf::WorkerInf,
               max_iter::Int,
               atol::Float64,
               rtol::Float64,
               gap_tol::Float64,
               report::Int,
               save_iter::Int,
               save_res::Int,
               lower_bound::Int,
               )::Tuple{Int,Float64,Float64}

    niter = 0
    report_flag = (report > 0)
    save_iter_flag = (save_iter > 0)
    save_res_flag = (save_res > 0)
    lb_flag = (lower_bound > 0)

    cr = ph_data.history.residuals[-1]
    delete!(ph_data.history.residuals, -1)
    xhat_res_sq = cr.xhat_sq
    x_res_sq = cr.x_sq

    nsqrt = max(sqrt(length(ph_data.xhat)), 1.0)
    xmax = (length(ph_data.xhat) > 0
            ? max(maximum(abs.(value.(values(ph_data.xhat)))), 1e-12)
            : 1.0)
    residual = sqrt(xhat_res_sq + x_res_sq) / nsqrt

    user_continue = @timeit(ph_data.time_info,
                            "User Callbacks",
                            _execute_callbacks(ph_data, worker_inf, niter)
                            )

    if report_flag
        _report(niter, nsqrt, residual, xmax, xhat_res_sq, x_res_sq)
    end

    if lb_flag
        (lb, gap) = @timeit(ph_data.time_info,
                            "Update Gap",
                            update_gap(ph_data, worker_inf, niter)
                            )
        if report_flag
            _report_lower_bound(niter, lb, gap)
        end
    end

    if save_iter_flag
        _save_iterate(ph_data, 0)
    end

    if save_res_flag
        _save_residual(ph_data, 0, xhat_res_sq, x_res_sq, residual, residual/xmax)
    end

    running = (user_continue
               && niter < max_iter
               && residual > atol
               && residual > rtol * xmax
               && (lb_flag ? gap > lb * gap_tol : true)
               )

    while running

        niter += 1

        # Solve subproblems
        @timeit(ph_data.time_info,
                "Solve subproblems",
                solve_subproblems(ph_data, worker_inf, niter)
                )

        # Update xhat and w
        (xhat_res_sq, x_res_sq) = @timeit(ph_data.time_info,
                                          "Update PH Vars",
                                          update_ph_variables(ph_data)
                                          )

        # Update stopping criteria -- xhat_res_sq measures the movement of
        # xhat values from k^th iteration to the (k+1)^th iteration while
        # x_res_sq measures the disagreement between the x variables and
        # its corresponding xhat variable (so lack of consensus amongst the
        # subproblems or violation of the nonanticipativity constraint)
        residual = sqrt(xhat_res_sq + x_res_sq) / nsqrt
        xmax = max(maximum(abs.(value.(values(ph_data.xhat)))), 1e-12)

        user_continue = @timeit(ph_data.time_info,
                                "User Callbacks",
                                _execute_callbacks(ph_data, worker_inf, niter)
                                )

        running = (user_continue
                   && niter < max_iter
                   && residual > atol
                   && residual > rtol * xmax
                   && (lb_flag ? gap > lb * gap_tol : true)
                   )

        if report_flag && (niter % report == 0 || !running)
            _report(niter, nsqrt, residual, xmax, xhat_res_sq, x_res_sq)
        end

        if lb_flag && (niter % lower_bound == 0 || !running)
            (lb, gap) = @timeit(ph_data.time_info,
                                "Update Gap",
                                update_gap(ph_data, worker_inf, niter)
                                )
            if report_flag
                _report_lower_bound(niter, lb, gap)
            end
        end

        if save_iter_flag && niter % save_iter == 0
            _save_iterate(ph_data, niter)
        end

        if save_res_flag && niter % save_res == 0
            _save_residual(ph_data,
                           niter,
                           xhat_res_sq,
                           x_res_sq,
                           residual,
                           residual/xmax
                           )
        end

    end

    @timeit(ph_data.time_info,
            "Finishing",
            finish(ph_data, worker_inf)
            )

    if niter >= max_iter && residual > atol && residual > rtol * xmax
        @warn("Performed $niter iterations without convergence. " *
              "Consider increasing max_iter from $max_iter.")
    end

    return (niter, residual, residual/xmax)
end


text:
using Documenter, XLSXasJSON

makedocs(
   modules = [XLSXasJSON],
 checkdocs = :all,
   authors = "YongHee Kim",
  sitename = "XLSXasJSON.jl",
     pages = Any[
              "Index" => "index.md",
              "User Guide" => "userguide.md",
             ]
)

deploydocs(
    repo   = "github.com/devsisters/XLSXasJSON.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
using Libdl

P4EST_ARTIFACT = false
if VERSION < v"1.3"
    DEFAULT_P4EST_ROOT_DIR = "/usr"
else
    using P4est_jll
    P4EST_ARTIFACT = true
    DEFAULT_P4EST_ROOT_DIR = P4est_jll.artifact_dir
end

P4EST_FOUND        = true
P4EST_ROOT_DIR     = haskey(ENV,"P4EST_ROOT_DIR") ? ENV["P4EST_ROOT_DIR"] : DEFAULT_P4EST_ROOT_DIR
P4EST_DIR          = haskey(ENV,"P4EST_DIR") ? ENV["P4EST_DIR"] : P4EST_ROOT_DIR
P4EST_LIB_DIR      = haskey(ENV,"P4EST_LIB_DIR") ? ENV["P4EST_LIB_DIR"] : joinpath(P4EST_DIR,"lib")
P4EST_INCLUDE_DIR  = haskey(ENV,"P4EST_INCLUDE_DIR") ? ENV["P4EST_INCLUDE_DIR"] : joinpath(P4EST_DIR,"include")
P4EST_LIB_NAME     = haskey(ENV,"P4EST_LIB_NAME") ? ENV["P4EST_LIB_NAME"] : "libp4est"
P4EST_LIB          = haskey(ENV,"P4EST_LIB") ? ENV["P4EST_LIB"] : ""

DEFAULT_P4EST_ENABLE_MPI = !haskey(ENV,"P4EST_ROOT_DIR") && P4EST_ARTIFACT ? false : true
P4EST_ENABLE_MPI   = haskey(ENV,"P4EST_ENABLE_MPI") ? Bool(ENV["P4EST_ENABLE_MPI"]) : DEFAULT_P4EST_ENABLE_MPI

# Check P4EST_DIR exists
if isdir(P4EST_DIR)
    @info "p4est directory found at: $P4EST_DIR"

    # Check P4EST_LIB_DIR (.../lib directory) exists
    if isdir(P4EST_LIB_DIR)
        @info "p4est lib directory found at: $P4EST_LIB_DIR"
    else
        @warn "P4EST lib directory not found: $P4EST_LIB_DIR"
    end

    if (cmp(P4EST_LIB,"") == 0)
    # Check P4EST_LIB (.../libp4est.so file) exists
        P4EST_LIB = string(Libdl.find_library(P4EST_LIB_NAME,[P4EST_LIB_DIR,P4EST_ROOT_DIR]),".",Libdl.dlext)
    end

    if isfile(P4EST_LIB)
        @info "P4EST library found at: $P4EST_LIB"
    else
        P4EST_FOUND = false
        @warn "P4EST lib ($P4EST_LIB) not found at: $P4EST_LIB_DIR"
    end

    # Check P4EST_INCLUDE_DIR (.../include directory) exists
    if isdir(P4EST_INCLUDE_DIR)
        @info "P4EST include directory found at: $P4EST_INCLUDE_DIR"
    else
        @warn "P4EST include directory not found: $P4EST_INCLUDE_DIR"
    end

else
    if (cmp(P4EST_LIB,"") == 0)
        # Check P4EST_LIB (.../libp4est.so file) exists
        P4EST_LIB = string(Libdl.find_library(P4EST_LIB_NAME,[P4EST_LIB_DIR,P4EST_ROOT_DIR]),".",Libdl.dlext)
    end

    if isfile(P4EST_LIB)
        @info "P4EST library found at: $P4EST_LIB"
    else
        P4EST_FOUND = false
        @warn "P4EST lib ($P4EST_LIB) not found at: $P4EST_ROOT_DIR"
    end
end

@info """
P4EST configuration:
==============================================
  - P4EST_FOUND           = $P4EST_FOUND
  - P4EST_DIR             = $P4EST_DIR
  - P4EST_LIB_DIR         = $P4EST_LIB_DIR
  - P4EST_INCLUDE_DIR     = $P4EST_INCLUDE_DIR
  - P4EST_LIB             = $P4EST_LIB
  - P4EST_ENABLE_MPI      = $P4EST_ENABLE_MPI
"""

if !P4EST_FOUND
    error("P4EST library not found")
end

# Write P4EST configuration to deps.jl file
deps_jl = "deps.jl"

if isfile(deps_jl)
  rm(deps_jl)
end

open(deps_jl,"w") do f
  println(f, "# This file is automatically generated")
  println(f, "# Do not edit")
  println(f)
  println(f, :(const P4EST_FOUND           = $P4EST_FOUND))
  println(f, :(const P4EST_DIR             = $P4EST_DIR))
  println(f, :(const P4EST_LIB_DIR         = $P4EST_LIB_DIR))
  println(f, :(const P4EST_INCLUDE_DIR     = $P4EST_INCLUDE_DIR))
  println(f, :(const P4EST_LIB             = $P4EST_LIB))
  println(f, :(const P4EST_ENABLE_MPI      = $P4EST_ENABLE_MPI))
end


text:
using CRlibm
using Compat

using Base.Test

@test cos(0.5, RoundDown) == 0.8775825618903726
@test cos(0.5, RoundUp) == 0.8775825618903728
@test cos(0.5, RoundNearest) == cos(0.5) == 0.8775825618903728
@test cos(1.6, RoundToZero) == -0.029199522301288812
@test cos(1.6, RoundDown) == -0.029199522301288815

function my_eps(prec::Int)
    ldexp(eps(Float64), 53-prec)
end

function my_eps(x::BigFloat)  # only works for precision >= 10?
    my_eps(precision(x) - exponent(x))
end

is_log(f) = string(f)[1:3] == "log"

function do_test(f, val)
    a = f(val, RoundDown)
    b = f(val, RoundUp)
    @test b - a == eps(a) || b - a == eps(b) || b - a == 0
end

function test_CRlibm(function_list)
    @show function_list
    for f in function_list
        #println("Testing CRlibm.$f")

        ff = eval(f)  # the actual Julia function

        for val in (0.51, 103.2, -17.1, -0.00005, 1)

            #@show f, val

            val <= 0.0 && is_log(f) && continue
            abs(val) > 1 && f ∈ (:asin, :acos) && continue

            do_test(ff, val)
        end
    end
end

function test_MPFR()
    for f in CRlibm.MPFR_function_list
        #println("Testing CRlibm.$f")

        ff = eval(f)  # the actual Julia function

        for val in (0.51, 103.2, -17.1, -0.00005)
            #print(val, " ")

            #@show f, val

            val <= 0.0 && is_log(f) && continue
            abs(val) > 1 && f ∈ (:asin, :acos) && continue

            for prec in (20, 100, 1000)

                setprecision(BigFloat, prec) do
                    val = BigFloat(val)
                    do_test(ff, val)
                end

            end
        end
    end
end

println("Testing CRlibm")
test_CRlibm(CRlibm.function_list)
# This will currently fail on the :sinpi etc. functions (that are not defined in MPFR) if MPFR is already enabled because the CRlibm library could not be found

println("Testing shadowing MPFR")
CRlibm.shadow_MPFR()
test_CRlibm(CRlibm.MPFR_function_list)

println("Testing MPFR")
test_MPFR()


text:
using DynamicalSystems
using Plots

# Simulate an discrete dynamical system (Henon map) using DynamicalSystems.jl

@inline @inbounds function hoop(x, p, n)
    dx1 = 1.0 - p[1]*x[1]^2 + x[2]
    dx2 = p[2]*x[1]
    return SVector{2}(dx1, dx2)
end
# Jacobian:
@inline @inbounds function hoop_jac(x, p, n)
    J = @SMatrix [-2*p[1]*x[1]      1.0;
                p[2]            0.0]
    return J
end

ds = DiscreteDynamicalSystem(hoop, [0.0; 0.0], [1.4, 0.3], hoop_jac)

data = trajectory(ds, 10000; dt=1)

plot(data[:, 1], data[:, 2], seriestype=:scatter, markersize=0.5)


text:
include("resnet-util.jl")

# need relu6(x) = clamp.(relu(x), 0, 6)

torch2knet(p::PyObject) = eval(Meta.parse(p._get_name()))(p)

function MobileNetV2(p::PyObject)
    s = torch2knet(p.features)
    push!(s, MobileNetOutput(p.classifier))
end

function MobileNetOutput(p::PyObject)
    Sequential(
        x->pool(x; mode=1, window=size(x)[1:2]),
        x->reshape(x, :, size(x,4)),
        Linear(param(t2a(p[2].weight)); bias=param(t2a(p[2].bias)), dropout=p[1].p),
    )
end

function Sequential(p::PyObject)
    s = Sequential()
    for l in p; push!(s, torch2knet(l)); end
    return s
end

function ConvBNReLU(p::PyObject)
    conv = torch2knet(p[1])
    conv.normalization = torch2knet(p[2])
    conv.activation = torch2knet(p[3])
    return conv
end

function Conv2d(p::PyObject; normalization=nothing, activation=nothing)
    w = param(permutedims(t2a(p.weight), (4,3,2,1)))
    bias = (p.bias === nothing ? nothing :
            param(reshape(t2a(p.bias), (1,1,:,1))))
    Conv(w; bias, normalization, activation,
         p.padding, p.stride, p.dilation, p.groups)
end

function BatchNorm2d(b::PyObject)
    b.track_running_stats = false # so we can do an inference comparison
    bnweight(x) = param(reshape(t2a(x), (1,1,:,1)))
    BatchNorm(
        ; use_estimates = nothing,
        update = b.momentum,
        mean = bnweight(b.running_mean).value,
        var = bnweight(b.running_var).value,
        bias = bnweight(b.bias),
        scale = bnweight(b.weight),
        epsilon = b.eps,
    )
end

function ReLU6(p::PyObject)
    return relu6
end

relu6(x)=clamp(x,0,6)

function InvertedResidual(p::PyObject)
    c = torch2knet(p.conv)
    if c[end] isa BatchNorm && c[end-1] isa Conv && c[end-1].normalization == nothing
        c[end-1].normalization = pop!(c)
    end
    all(c[2].stride .== 1) && size(c[1].w, 3) == size(c[end].w, 4) ? Residual(c, identity) : c
end

mv2 = models.mobilenet_v2(pretrained=true).eval()

a = torch2knet(mv2)

#=
channels: 16,32

julia> a
Sequential

Input:
224x3=>112x32
  Conv(3×3, 3=>32, padding=(1, 1), stride=(2, 2), BatchNorm(), relu6)

[t,c,n,s] = [exp factor, output channels, layers in sequence with same output, first layer stride]

16Channels:
[1,16,1,1]
112x32=>112x16
  Sequential
    Conv(3×3, 1=>32, padding=(1, 1), groups=32, BatchNorm(), relu6)
    Conv(1×1, 32=>16, BatchNorm())

24Channels:
[6, 24, 2, 2]-1
112x16=>56x24
  Sequential
    Conv(1×1, 16=>96, BatchNorm(), relu6)
    Conv(3×3, 1=>96, padding=(1, 1), stride=(2, 2), groups=96, BatchNorm(), relu6)
    Conv(1×1, 96=>24, BatchNorm())

[6, 24, 2, 2]-2
56x24=>56x24
  Residual
    Sequential
      Conv(1×1, 24=>144, BatchNorm(), relu6)
      Conv(3×3, 1=>144, padding=(1, 1), groups=144, BatchNorm(), relu6)
      Conv(1×1, 144=>24, BatchNorm())
    +
    identity

32Channels:
[6, 32, 3, 2]-1
56x24=>28x32
  Sequential
    Conv(1×1, 24=>144, BatchNorm(), relu6)
    Conv(3×3, 1=>144, padding=(1, 1), stride=(2, 2), groups=144, BatchNorm(), relu6)
    Conv(1×1, 144=>32, BatchNorm())

[6, 32, 3, 2]-2
28x32=>28x32-1
  Residual
    Sequential
      Conv(1×1, 32=>192, BatchNorm(), relu6)
      Conv(3×3, 1=>192, padding=(1, 1), groups=192, BatchNorm(), relu6)
      Conv(1×1, 192=>32, BatchNorm())
    +
    identity

[6, 32, 3, 2]-3
28x32=>28x32-2
  Residual
    Sequential
      Conv(1×1, 32=>192, BatchNorm(), relu6)
      Conv(3×3, 1=>192, padding=(1, 1), groups=192, BatchNorm(), relu6)
      Conv(1×1, 192=>32, BatchNorm())
    +
    identity

64Channels:
[6, 64, 4, 2]-1
28x32=>14x64
  Sequential
    Conv(1×1, 32=>192, BatchNorm(), relu6)
    Conv(3×3, 1=>192, padding=(1, 1), stride=(2, 2), groups=192, BatchNorm(), relu6)
    Conv(1×1, 192=>64, BatchNorm())

[6, 64, 4, 2]-2
14x64=>14x64-1
  Residual
    Sequential
      Conv(1×1, 64=>384, BatchNorm(), relu6)
      Conv(3×3, 1=>384, padding=(1, 1), groups=384, BatchNorm(), relu6)
      Conv(1×1, 384=>64, BatchNorm())
    +
    identity

[6, 64, 4, 2]-3
14x64=>14x64-2
  Residual
    Sequential
      Conv(1×1, 64=>384, BatchNorm(), relu6)
      Conv(3×3, 1=>384, padding=(1, 1), groups=384, BatchNorm(), relu6)
      Conv(1×1, 384=>64, BatchNorm())
    +
    identity

[6, 64, 4, 2]-4
14x64=>14x64-3
  Residual
    Sequential
      Conv(1×1, 64=>384, BatchNorm(), relu6)
      Conv(3×3, 1=>384, padding=(1, 1), groups=384, BatchNorm(), relu6)
      Conv(1×1, 384=>64, BatchNorm())
    +
    identity

96Channels:
[6, 96, 3, 1]-1
14x64=>14x96
  Sequential
    Conv(1×1, 64=>384, BatchNorm(), relu6)
    Conv(3×3, 1=>384, padding=(1, 1), groups=384, BatchNorm(), relu6)
    Conv(1×1, 384=>96, BatchNorm())

[6, 96, 3, 1]-2
14x96=>14x96-1
  Residual
    Sequential
      Conv(1×1, 96=>576, BatchNorm(), relu6)
      Conv(3×3, 1=>576, padding=(1, 1), groups=576, BatchNorm(), relu6)
      Conv(1×1, 576=>96, BatchNorm())
    +
    identity

[6, 96, 3, 1]-3
14x96=>14x96-2
  Residual
    Sequential
      Conv(1×1, 96=>576, BatchNorm(), relu6)
      Conv(3×3, 1=>576, padding=(1, 1), groups=576, BatchNorm(), relu6)
      Conv(1×1, 576=>96, BatchNorm())
    +
    identity

160Channels:
[6, 160, 3, 2]-1
14x96=>7x160
  Sequential
    Conv(1×1, 96=>576, BatchNorm(), relu6)
    Conv(3×3, 1=>576, padding=(1, 1), stride=(2, 2), groups=576, BatchNorm(), relu6)
    Conv(1×1, 576=>160, BatchNorm())

[6, 160, 3, 2]-2
7x160=>7x160-1
  Residual
    Sequential
      Conv(1×1, 160=>960, BatchNorm(), relu6)
      Conv(3×3, 1=>960, padding=(1, 1), groups=960, BatchNorm(), relu6)
      Conv(1×1, 960=>160, BatchNorm())
    +
    identity

[6, 160, 3, 2]-3
7x160=>7x160-2
  Residual
    Sequential
      Conv(1×1, 160=>960, BatchNorm(), relu6)
      Conv(3×3, 1=>960, padding=(1, 1), groups=960, BatchNorm(), relu6)
      Conv(1×1, 960=>160, BatchNorm())
    +
    identity

320Channels:
[6, 320, 1, 1]-1
7x160=>7x320
  Sequential
    Conv(1×1, 160=>960, BatchNorm(), relu6)
    Conv(3×3, 1=>960, padding=(1, 1), groups=960, BatchNorm(), relu6)
    Conv(1×1, 960=>320, BatchNorm())

Output:
7x320=>7x1280
  Conv(1×1, 320=>1280, BatchNorm(), relu6)

  Sequential
    x->pool(x; mode=1, window=size(x)[1:2])
    x->reshape(x, :, size(x,4))
    Linear(1280=>1000, bias, dropout=0.2)

=#


text:
using DataFrames
using ExtendedDates
using Test
using TimeDataFrames

include("test_functions.jl")


text:
# # Spherical harmonic addition theorem
# This example confirms numerically that
# ```math
# f(z) = \frac{P_n(z\cdot y) - P_n(x\cdot y)}{z\cdot y - x\cdot y},
# ```
# is actually a degree-$(n-1)$ polynomial on $\mathbb{S}^2$, where $P_n$ is the degree-$n$
# Legendre polynomial, and $x,y,z \in \mathbb{S}^2$.
# To verify, we sample the function on a $N\times M$ equiangular grid
# defined by:
# ```math
# \begin{aligned}
# \theta_n & = (n+\tfrac{1}{2})\pi/N,\quad{\rm for}\quad 0\le n < N,\quad{\rm and}\\
# \varphi_m & = 2\pi m/M,\quad{\rm for}\quad 0\le m < M;
# \end{aligned}
# ```
# we convert the function samples to Fourier coefficients using
# `plan_sph_analysis`; and finally, we transform
# the Fourier coefficients to spherical harmonic coefficients using
# `plan_sph2fourier`.
#
# In the basis of spherical harmonics, it is plain to see the
# addition theorem in action, since $P_n(x\cdot y)$ should only consist of
# exact-degree-$n$ harmonics.
#
# For the storage pattern of the arrays, please consult the
# [documentation](https://MikaelSlevinsky.github.io/FastTransforms).

function threshold!(A::AbstractArray, ϵ)
    for i in eachindex(A)
        if abs(A[i]) < ϵ A[i] = 0 end
    end
    A
end

using FastTransforms, LinearAlgebra, Plots
const GENFIGS = joinpath(pkgdir(FastTransforms), "docs/src/generated")
!isdir(GENFIGS) && mkdir(GENFIGS)
plotlyjs()

# The colatitudinal grid (mod $\pi$):
N = 15
θ = (0.5:N-0.5)/N

# The longitudinal grid (mod $\pi$):
M = 2*N-1
φ = (0:M-1)*2/M

# Arbitrarily, we place $x$ at the North pole:
x = [0,0,1]

# Another vector is completely free:
y = normalize([.123,.456,.789])

# Thus $z \in \mathbb{S}^2$ is our variable vector, parameterized in spherical coordinates:
z = (θ,φ) -> [sinpi(θ)*cospi(φ), sinpi(θ)*sinpi(φ), cospi(θ)]

# On the tensor product grid, the Legendre polynomial $P_n(z\cdot y)$ is:
A = [(2k+1)/(k+1) for k in 0:N-1]
B = zeros(N)
C = [k/(k+1) for k in 0:N]
c = zeros(N); c[N] = 1
pts = vec([z(θ, φ)⋅y for θ in θ, φ in φ])
phi0 = ones(N*M)
F = reshape(FastTransforms.clenshaw!(c, A, B, C, pts, phi0, zeros(N*M)), N, M)

# We superpose a surface plot of $f$ on top of the grid:
X = [sinpi(θ)*cospi(φ) for θ in θ, φ in φ]
Y = [sinpi(θ)*sinpi(φ) for θ in θ, φ in φ]
Z = [cospi(θ) for θ in θ, φ in φ]
scatter3d(vec(X), vec(Y), vec(Z); markersize=1.25, markercolor=:violetred)
surface!(X, Y, Z; surfacecolor=F, legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "sphere1.html"))
###```@raw html
###<object type="text/html" data="../sphere1.html" style="width:100%;height:400px;"></object>
###```

# We show the cut in the surface to help illustrate the definition of the grid.
# In particular, we do not sample the poles.
#
# We precompute a spherical harmonic--Fourier plan:
P = plan_sph2fourier(F)

# And an FFTW Fourier analysis plan on $\mathbb{S}^2$:
PA = plan_sph_analysis(F)

# Its spherical harmonic coefficients demonstrate that it is exact-degree-$n$:
V = PA*F
U = threshold!(P\V, 400*eps())

# The $L^2(\mathbb{S}^2)$ norm of the function is:
nrm1 = norm(U)

# Similarly, on the tensor product grid, our function samples are:
Pnxy = FastTransforms.clenshaw!(c, A, B, C, [x⋅y], [1.0], [0.0])[1]
F = [(F[n, m] - Pnxy)/(z(θ[n], φ[m])⋅y - x⋅y) for n in 1:N, m in 1:M]

# We superpose a surface plot of $f$ on top of the grid:
scatter3d(vec(X), vec(Y), vec(Z); markersize=1.25, markercolor=:violetred)
surface!(X, Y, Z; surfacecolor=F, legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "sphere2.html"))
###```@raw html
###<object type="text/html" data="../sphere2.html" style="width:100%;height:400px;"></object>
###```

# Its spherical harmonic coefficients demonstrate that it is degree-$(n-1)$:
V = PA*F
U = threshold!(P\V, 400*eps())

# Finally, the Legendre polynomial $P_n(z\cdot x)$ is aligned with the grid:
pts = vec([z(θ, φ)⋅x for θ in θ, φ in φ])
F = reshape(FastTransforms.clenshaw!(c, A, B, C, pts, phi0, zeros(N*M)), N, M)

# We superpose a surface plot of $f$ on top of the grid:
scatter3d(vec(X), vec(Y), vec(Z); markersize=1.25, markercolor=:violetred)
surface!(X, Y, Z; surfacecolor=F, legend=false, xlabel="x", ylabel="y", zlabel="f")
savefig(joinpath(GENFIGS, "sphere3.html"))
###```@raw html
###<object type="text/html" data="../sphere3.html" style="width:100%;height:400px;"></object>
###```

# It only has one nonnegligible spherical harmonic coefficient.
# Can you spot it?
V = PA*F
U = threshold!(P\V, 400*eps())

# That nonnegligible coefficient should be
ret = eval("√(2π/($(N-1)+1/2))")

# which is approximately
eval(Meta.parse(ret))

# since the convention in this library is to orthonormalize.
nrm2 = norm(U)

# Note that the integrals of both functions $P_n(z\cdot y)$ and $P_n(z\cdot x)$ and their
# $L^2(\mathbb{S}^2)$ norms are the same because of rotational invariance. The integral of
# either is perhaps not interesting as it is mathematically zero, but the norms
# of either should be approximately the same.
nrm1 ≈ nrm2