text:
"""
International Standard Atmospheric (ISA) model

Only metric units are used.

Source: EUROCONTROL BADA 4 User Manual Chapter 2.2 Atmosphere Model
"""
module Atmosphere

export T, T₀, p, p₀, ρ, ρ₀, a, a₀, g₀, κ,
Vcas2Vtas, Vtas2Vcas, M2Vtas, Vtas2M, Hp_trans, M2Vcas, Vcas2M

"EUROCONTROL BADA 4 User Manual Section 2.2.1/2.2.2"

"Standard atmospheric temperature [K] at Mean Sea Level (MSL)"
const T₀ = 288.15
"Standard atmospheric pressure [Pa] at Mean Sea Level (MSL)"
const p₀ = 101325.0
"Standard atmospheric density [kg/m³] at Mean Sea Level (MSL)"
const ρ₀ = 1.225
"Speed of sound [m/s] at Mean Sea Level (MSL)"
const a₀ = 340.294
"Adiabatic index of air []"
const κ = 1.4
"constant used for Vtas<->Vcas conversion"
const μ = (κ - 1.0)/κ
"Real gas constant for air [M²/(Ks²)]"
const R = 287.05287
"Graviation acceleration [m/s²]"
const g₀ = 9.80665
"ISA temperature gradient [K/m] with altitude below the tropopause"
const βT∇ = -0.0065
"Geopotential pressure altitude [m] of Tropopause"
const Hp_trop = 11000.0

"""
T(Hp::Float64, ΔT::Float64 = 0.0)
Atmospheric temperature [K]
at pressure altitude Hp [m] and with temperature offset ΔT [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-13/16
"""
function T(Hp::Float64, ΔT::Float64 = 0.0)
    if Hp ≤ Hp_trop
        return T₀ + ΔT + βT∇ * Hp
    else
        return T₀ + ΔT + βT∇ * Hp_trop
    end
end

T(Hp::Real, ΔT::Real = zero(Real)) = T(convert(Float64, Hp), convert(Float64, ΔT))

"""
p(Hp::Float64, ΔT::Float64 = 0.0)
Air pressure [Pa]
at pressure altitude Hp [m] and with temperature offset ΔT [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-18/20
"""
function p(Hp::Float64, ΔT::Float64 = 0.0)
    if Hp ≤ Hp_trop
        return p₀*((T(Hp, ΔT) - ΔT) / T₀) ^ (-g₀/(βT∇*R))
    else
        p_trop = p₀*((T(Hp_trop, ΔT) - ΔT) / T₀) ^ (-g₀/(βT∇*R))
        return p_trop * exp(-g₀ / (R*T(Hp_trop)) * (Hp - Hp_trop))
    end
end

p(Hp::Real, ΔT::Real = zero(Real)) = p(convert(Float64, Hp), convert(Float64, ΔT))

"""
ρ(p, T)
Air density [kg/m³] at pressure level [Pa] and Temperature [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-21
"""
function ρ(p::Float64, T::Float64)
    return p / (R * T)
end

ρ(p::Real, T::Real) = ρ(convert(Float64, p), convert(Float64, T))

"""
a(T::Float64)
Speed of sound [m/s]
at pressure level [Pa] and with temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-22
"""
function a(T::Float64)
    return √(κ*R*T)
end

a(T::Real) = a(convert(Float64, T))

"""
Vcas2Vtas(Vcas::Float64, p::Float64, T::Float64)
True airspeed Vtas [m/s] as a function of the calibrated airspeed Vcas [m/s]
at pressure level p [Pa] and with temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-23
"""
function Vcas2Vtas(Vcas::Float64, p::Float64, T::Float64)
    ρᵥ = ρ(p, T)
    return (2.0*p/(μ*ρᵥ)*((1.0+p₀/p*((1.0+μ*ρ₀/(2.0*p₀)*Vcas*Vcas)
            ^(1.0/μ)-1.0))^μ-1.0))^0.5
end

Vcas2Vtas(Vcas::Real, p::Real, T::Real) = Vcas2Vtas(convert(Float64, Vcas),
convert(Float64, p), convert(Float64, T))

"""
Vtas2Vcas(Vtas::Float64, p::Float64, T::Float64)
Calibrated airspeed Vcas [m/s] as a function of the True airspeed Vtas [m/s]
at pressure level p [Pa] and with temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-24
"""
function Vtas2Vcas(Vtas::Float64, p::Float64, T::Float64)
    ρᵥ = ρ(p, T)
    return (2.0*p₀/(μ*ρ₀)*((1.0+p/p₀*((1.0+μ*ρᵥ/(2.0*p)*Vtas*Vtas)
            ^(1.0/μ)-1.0))^μ-1.0))^0.5
end

Vtas2Vcas(Vtas::Real, p::Real, T::Real) = Vtas2Vcas(convert(Float64, Vtas),
convert(Float64, p), convert(Float64, T))

"""
M2Vtas(M::Float64, T::Float64)
True airspeed Vtas [m/s] as a function of the Mach number
and temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-26
"""
function M2Vtas(M::Float64, T::Float64)
    return M*a(T)
end

M2Vtas(M::Real, T::Real) = M2Vtas(convert(Float64, M), convert(Float64, T))

"""
Vtas2M(Vtas::Float64, T::Float64)
Mach as a function of the True airspeed Vtas [m/s]
and temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-26 (reversed)
"""
function Vtas2M(Vtas::Float64, T::Float64)
    return Vtas/a(T)
end

Vtas2M(Vtas::Real, T::Real) = Vtas2M(convert(Float64, Vtas), convert(Float64, T))

"""
M2Vcas(M::Float64, p::Float64, T::Float64)
Calibrated airspeed Vcas [m/s] as a function of the Mach number
at pressure level p [Pa] and with temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-26
"""
function M2Vcas(M::Float64, p::Float64, T::Float64)
    return Vtas2Vcas(M2Vtas(M, T), p, T)
end

M2Vcas(M::Real, p::Real, T::Real) = M2Vcas(convert(Float64, M),
convert(Float64, p), convert(Float64, T))

"""
Vtas2M(Vtas::Float64, T::Float64)
Mach as a function of the True airspeed Vtas [m/s]
at pressure level p [Pa] and with temperature T [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-26 (reversed)
"""
function Vcas2M(Vcas::Float64, p::Float64, T::Float64)
    return Vtas2M(Vcas2Vtas(Vcas, p, T), T)
end

Vcas2M(Vcas::Real, p::Real, T::Real) = Vcas2M(convert(Float64, Vcas),
convert(Float64, p), convert(Float64, T))

"""
Hp_trans(Vcas::Float64, M::Float64, ΔT::Float64 = 0.0)
Transition altitude (also called crossover altitude) [m] between
a given calibrated airspeed [m/s] and a Mach number M
and with temperature offset ΔT [K]
EUROCONTROL BADA 4 User Manual eq. 2.2-27/28/29
"""
function Hp_trans(Vcas::Float64, M::Float64, ΔT::Float64 = 0.0)
    δ_trans = ((1.0+(κ-1.0)/2.0*(Vcas/a₀)^2)
                ^(1.0/μ)-1.0)/((1.0+(κ-1.0)/2.0*M^2)^(1.0/μ)-1.0)
    θ_trans = δ_trans^(-βT∇ * R / g₀)
    return T₀/βT∇*(θ_trans - 1.0)
end

Hp_trans(Vcas::Real, M::Real, ΔT::Real = zero(Real)) = Hp_trans(
convert(Float64, Vcas), convert(Float64, M), convert(Float64, ΔT))

"""
θ(T::Float64)
temperature ratio
EUROCONTROL BADA 4 User Manual eq. 2.2-30
"""
θ(T::Float64) = T / T₀
θ(T::Real) = θ(convert(Float64, T))
"""
δ(p::Float64)
pressure ratio
EUROCONTROL BADA 4 User Manual eq. 2.2-31
"""
δ(p::Float64) = p / p₀
δ(p::Real) = δ(convert(Float64, p))
"""
σ(ρ::Float64)
density ratio
EUROCONTROL BADA 4 User Manual eq. 2.2-32
"""
σ(ρ::Float64) = ρ / ρ₀
σ(ρ::Real) = σ(convert(Float64, ρ))

#TODO We recalculate all items all the time. Is there a faster strategy?
#TODO Interpolate Temperature/Pressure data grid

end # module


text:
@testset "NumField/Elem" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  QasNumberField, _ = NumberField(x - 1)
  Kt, t = PolynomialRing(QasNumberField, "t")
  K1, a1 = NumberField(x^3 - 2)
  K2, (a2, ) = NumberField([x^3 - 2])
  K3, a3 = NumberField(t^3 - 2)
  K4, (a4, ) = NumberField([t^3 - 2])

  @testset for (K, a) in [(K1, a1), (K2, a2), (K3, a3), (K4, a4)]
    b = one(K)
    fl = @inferred isintegral(b)
    @test fl
    fl = @inferred isintegral(a)

    B = @inferred basis(K)
    c = @inferred basis_matrix([one(K), a^4])
    @assert nrows(c) == 2
    @assert ncols(c) == 3
    @assert sum(B[i] * c[1, i] for i in 1:3) == one(K)
    @assert sum(B[i] * c[2, i] for i in 1:3) == a^4

    f = @inferred charpoly(one(K()))
    @test f == (gen(parent(f)) - 1)^3
    f = @inferred minpoly(one(K))
    @test f == (gen(parent(f)) - 1)

    f = @inferred charpoly(a)
    @test f == gen(parent(f))^3 - 2
    f = @inferred minpoly(a)
    @test f == gen(parent(f))^3 - 2

  end

  @testset "NumField/Component" begin
    Qx, x = PolynomialRing(FlintQQ, "x")
    QasNumberField, _ = NumberField(x - 1)
    Kt, t = PolynomialRing(QasNumberField, "t")
    K3, a3 = NumberField(t^3 - 2)
    K4, (a4, ) = NumberField([t^3 - 2])
    K41, mK4 = component(K4, 1)
    @test isisomorphic(K41, K3)[1]
    @test mK4(gen(K41)) == a4

    K1, a1 = NumberField(x^3 - 2)
    K2, (a2, ) = NumberField([x^3 - 2])
    K21, mK2 = component(K2, 1)
    @test isisomorphic(K21, K1)[1]
    @test mK2(gen(K21)) == a2
  end

end

ext: jl


text:
using SymPy

N = 11

b = [symbols("b_$i", real=true) for i in 1:N]
B = diagm(b)

# 2nd order dissipation
D1 = zeros(Int, N, N)
for i in 1:N
    D1[i,i] = 1
end
for i in 1:N-1
    D1[i+1,i] = -1
end
D1[1,2] = 1
display(D1)

B1 = eye(Int, N); B1[1,1] = 0
display(B1)

display(D1' * B * D1)
display(D1' * B1 * D1)

# 4th order dissipation
D2 = zeros(Int, N, N)
for i in 1:N
    D2[i,i] = -2
end
for i in 1:N-1
    D2[i+1,i] = 1
    D2[i,i+1] = 1
end
boundary_width = 1
for i in 1:boundary_width
    D2[i,:] = D2[boundary_width+1,:]
    D2[end-i+1,:] = D2[end-boundary_width,:]
end
display(D2)

B2 = eye(Int, N); B2[1,1] = B2[end,end] = 0
display(B2)

display(D2' * B * D2)
display(D2' * B2 * D2)

# 6th order dissipation
D3 = zeros(Int, N, N)
for i in 1:N
    D3[i,i] = -3
end
for i in 1:N-1
    D3[i+1,i] = 3
    D3[i,i+1] = 1
end
for i in 1:N-2
    D3[i+2,i] = -1
end
boundary_width = 2
for i in 1:boundary_width
    D3[i,:] = D3[boundary_width+1,:]
end
for i in 1:boundary_width-1
    D3[end-i+1,:] = D3[end-boundary_width+1,:]
end
display(D3)

B3 = eye(Int, N); B3[1,1] = B3[2,2]= B3[end,end] = 0
display(B3)

display(D3' * B * D3)
display(D3' * B3 * D3)

# 8th order dissipation
D4 = zeros(Int, N, N)
for i in 1:N
    D4[i,i] = 6
end
for i in 1:N-1
    D4[i+1,i] = -4
    D4[i,i+1] = -4
end
for i in 1:N-2
    D4[i+2,i] = 1
    D4[i,i+2] = 1
end
boundary_width = 2
for i in 1:boundary_width
    D4[i,:] = D4[boundary_width+1,:]
    D4[end-i+1,:] = D4[end-boundary_width,:]
end
display(D4)

B4 = eye(Int, N); B4[1,1] = B4[2,2]= B4[end,end] = B4[end-1,end-1] = 0
display(B4)

display(D4' * B * D4)
display(D4' * B4 * D4)


text:
using QuasinormalModes
using SymEngine

# ------------------------------------------------------------------
# 1. Analytic harmonic oscillator
# ------------------------------------------------------------------

struct HarmonicOscilatorData{N,T} <: QuadraticEigenvalueProblem{N,T}
    nIter::N
    x0::T

    vars::Tuple{Basic,Basic}
    exprs::Tuple{Basic,Basic}
end

function HarmonicOscilatorData(nIter::N, x0::T) where {N,T}

    vars = @vars x ω

    λ0 = 2 * x
    S0 = 1 - ω

    return HarmonicOscilatorData{N,T}(nIter, x0, vars, (λ0, S0))
end

QuasinormalModes.λ0(d::HarmonicOscilatorData{N,T}) where {N,T} = d.exprs[1]
QuasinormalModes.S0(d::HarmonicOscilatorData{N,T}) where {N,T} = d.exprs[2]

QuasinormalModes.get_niter(d::HarmonicOscilatorData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::HarmonicOscilatorData{N,T}) where {N,T} = d.x0

QuasinormalModes.get_ODEvar(d::HarmonicOscilatorData{N,T}) where {N,T} = d.vars[1]
QuasinormalModes.get_ODEeigen(d::HarmonicOscilatorData{N,T}) where {N,T} = d.vars[2]

# ------------------------------------------------------------------
# 2. Numeric harmonic oscillator
# ------------------------------------------------------------------

struct NHarmonicOscilatorData{N,T} <: NumericAIMProblem{N,T}
    nIter::N
    x0::T
end

function NHarmonicOscilatorData(nIter::N, x0::T) where {N,T}
    return NHarmonicOscilatorData{N,T}(nIter, x0)
end

QuasinormalModes.λ0(::NHarmonicOscilatorData{N,T}) where {N,T} = (x, ω) -> 2 * x
QuasinormalModes.S0(::NHarmonicOscilatorData{N,T}) where {N,T} = (x, ω) -> 1 - ω + x - x

QuasinormalModes.get_niter(d::NHarmonicOscilatorData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::NHarmonicOscilatorData{N,T}) where {N,T} = d.x0

# ------------------------------------------------------------------
# 3. Constructing problems and caches
# ------------------------------------------------------------------

p_ana = HarmonicOscilatorData(0x0000A, 0.5);
p_num = NHarmonicOscilatorData(0x0000A, 0.5);

c_ana = AIMCache(p_ana)
c_num = AIMCache(p_num)

# ------------------------------------------------------------------
# 4. Computing quasinormal modes
# ------------------------------------------------------------------

ev_ana = computeEigenvalues(Serial(), p_ana, c_ana)
ev_num = eigenvaluesInGrid(Serial(), p_num, c_num, (0.0, 21.0))

function printEigen(eigenvalues)
    println("--------------------------------------")

    for i in eachindex(eigenvalues)
        println("n = $i, ω = $(eigenvalues[i])")
    end

    println("--------------------------------------")

    return nothing
end

println("Analytic results")
printEigen(reverse!(ev_ana))

println("Numeric results")
printEigen(ev_num)


text:
# lattice.jl
#
# Defines the spatial lattice; supports 1D and 2D square, open and periodic boundaries
# The main data structure is the bond-spin index array, bond_spin[nBond,2]

using Random
Random.seed!(1234)

Dim = 1
nX = 6
PBC = false

if Dim == 2
	nY = nX  #Works for square lattices only right now
end

#one-dimensional lattice
if Dim == 1
    nSpin = nX
    if PBC == true
       nBond = nSpin
    else
       nBond = nSpin-1
    end

    bond_spin = zeros(Int,nBond,2) #assign site indices to bonds
    for i = 1:nBond
        bond_spin[i,1] = i
        bond_spin[i,2] = i+1
    end
    if PBC == true
        bond_spin[nBond,2] = 1
    end

    #println(bond_spin)

#two-dimensional square lattice OBC
elseif Dim == 2 && PBC == false
    nSpin = nX*nY
    nBond = 2*nSpin-nX-nY
    bond_spin = zeros(Int,nBond,2) #assign site indices to bonds

    cnt=1
    for i = 1:div(nBond,2) #horizontal
        #println(i," ",cnt)
        if mod(i,nX-1) != 0
            bond_spin[i,1] = cnt
            bond_spin[i,2] = cnt + 1
            global cnt += 1
        else
            bond_spin[i,1] = cnt
            bond_spin[i,2] = cnt + 1
            global cnt += 2
        end
    end
    global cnt = 1
    for i = (div(nBond,2)+1):nBond #vertical
        bond_spin[i,1] = cnt
        bond_spin[i,2] = cnt + nX
        global cnt += 1
    end

#two-dimensional square lattice PBC
elseif Dim == 2 && PBC == true
    nSpin = nX*nY
    nBond = 2*nSpin
    bond_spin = zeros(Int,nBond,2) #assign site indices to bonds

    for i = 1:nSpin
        bond1 = 2*i - 1 #horizontal
        bond_spin[bond1,1] = i
        bond_spin[bond1,2] = i+1
        if mod(i,nX) == 0
            bond_spin[bond1,2] -= nX
        end

        bond2 = 2*i     #vertical
        bond_spin[bond2,1] = i
        bond_spin[bond2,2] = i+nX
        if i > (nSpin - nX)
            bond_spin[bond2,2] -= nSpin
        end
    end

else
    println("Dimension error")
end


text:
### StochasticPolicy ###
# maintained by @etotheipluspi

mutable struct StochasticPolicy{D, RNG <: AbstractRNG} <: Policy
    distribution::D
    rng::RNG
end
# The constructor below should be used to create the policy so that the action space is initialized correctly
StochasticPolicy(distribution; rng=Base.GLOBAL_RNG) = StochasticPolicy(distribution, rng)

## policy execution ##
function action(policy::StochasticPolicy, s)
    return rand(policy.rng, policy.distribution)
end

## convenience functions ##
updater(policy::StochasticPolicy) = VoidUpdater() # since the stochastic policy does not depend on the belief

# Samples actions uniformly
UniformRandomPolicy(problem, rng=Base.GLOBAL_RNG) = StochasticPolicy(actions(problem), rng)

mutable struct CategoricalTabularPolicy <: Policy
    stochastic::StochasticPolicy
    value::ValuePolicy
end
CategoricalTabularPolicy(mdp::Union{POMDP,MDP}; rng=Base.GLOBAL_RNG) = CategoricalTabularPolicy(StochasticPolicy(Weights(zeros(n_actions(mdp)))), ValuePolicy(mdp))

function action(policy::CategoricalTabularPolicy, s)
    policy.stochastic.distribution = Weights(policy.value.value_table[state_index(policy.value.mdp, s),:])
    return policy.value.act[sample(policy.stochastic.rng, policy.stochastic.distribution)]
end

mutable struct EpsGreedyPolicy <: Policy
    eps::Float64
    val::ValuePolicy
    uni::StochasticPolicy
end

EpsGreedyPolicy(mdp::Union{MDP,POMDP}, eps::Float64;
                rng=Base.GLOBAL_RNG) = EpsGreedyPolicy(eps, ValuePolicy(mdp), UniformRandomPolicy(mdp, rng))

function action(policy::EpsGreedyPolicy, s)
    if rand(policy.uni.rng) > policy.eps
        return action(policy.val, s)
    else
        return action(policy.uni, s)
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "liblsl"
version = v"1.13.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/sccn/liblsl/archive/1.13.0.tar.gz",
    "5b304a5365eba33852da96badfbd9d66556caf4a00c87947a59df2942680a617"),
    DirectorySource("./bundled/patches"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd liblsl-1.13.0

# Add license file
install_license LICENSE

# Link against real time and correct C->C++ library paths on linux
if [[ ${target} == x86_64-linux-* || ${target} == aarch64-linux-* || ${target} == powerpc64le-linux-* ]]; then
    export CXXFLAGS="-lrt"
    export CFLAGS="-lrt -Wl,-rpath-link,/opt/${target}/${target}/lib64"
fi

if [[ ${target} == i686-linux-* || ${target} == arm-linux-* ]]; then
    export CXXFLAGS="-lrt"
    export CFLAGS="-lrt -Wl,-rpath-link,/opt/${target}/${target}/lib"
fi

# Enable C++ 2011 support and patch for MinGW
if [[ ${target} == *-w64-* ]]; then
    atomic_patch -p1 $WORKSPACE/srcdir/liblsl_mingw.diff
    # Line endings are confused on the test internal CMakeLists, fix them before patch
    dos2unix testing/CMakeLists.txt
    atomic_patch -p1 $WORKSPACE/srcdir/lsl_test_internal_mingw.diff
    export CXXFLAGS="-std=c++11"
fi

mkdir build
cd build

cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DLSL_UNIXFOLDERS=1 -DLSL_NO_FANCY_LIBNAME=1 -DLSL_UNITTESTS=1 ../
make -j${nproc}

# We can't run unit-tests as we are cross-compiling
#./lslver
#./testing/lsl_test_internal
#./testing/lsl_test_exported

make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblsl", :liblsl),
    ExecutableProduct("lslver", :lslver)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#md # ---

#md # # Why use Julia language!
#md #
#md # - **You develop in the same language in which you optimize.**
#md # - Packaging system is very efficient (5858 registered packages)
#md # - PyPi (311,500 projects) R (17739 packages)
#md # - It is very easy to create a package (easier than R and Python)
#md # - It is very easy to access to a GPU device.
#md # - Nice interface for Linear Algebra and Differential Equations
#md # - Easy access to BLAS, LAPACK and scientific computing librairies.
#md # - Julia talks to all major Languages - mostly without overhead!

#md # ---

#md # # What's bad
#md #
#md # - It is still hard to build shared library or executable from Julia code.
#md # - Compilation times lattency. Using [Revise.jl](https://github.com/timholy/Revise.jl) helps a lot.
#md # - Plotting takes time (5 seconds for the first plot)
#md # - OpenMP is better than the Julia multithreading library but it is progressing.
#md # - Does not work well with vectorized code, you need to do a lot of inplace computation to avoid memory allocations and use explicit views to avoid copy. There are some packages like [LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl).
#md
#md # [What's Bad About Julia by Jeff Bezanson](https://www.youtube.com/watch?v=TPuJsgyu87U)

#md # ---

#md # ## From zero to Julia!
#md #
#md # https://techytok.com/from-zero-to-julia/
#md #
#md # ## Python-Julia benchmarks by Thierry Dumont
#md #
#md # https://github.com/Thierry-Dumont/BenchmarksPythonJuliaAndCo/wiki
#md #
#md # ## Mailing List
#md # - https://listes.services.cnrs.fr/wws/info/julia
#md #


text:
# sampling functions for simulation
# Copy and modified from pkg StatsBase
# Link: https://github.com/JuliaStats/StatsBase.jl/blob/master/src/sampling.jl
# [MIT License](https://github.com/JuliaStats/StatsBase.jl/blob/master/LICENSE.md)

using Random

"""
    sample([rng], wv::AbstractArray)

Select a single random integer in `wa` with probabilities
proportional to the weights given in `wa`.

Optionally specify a random number generator `rng` as the first argument
(defaults to `Random.GLOBAL_RNG`).

"""
function sample(rng::AbstractRNG, wv::AbstractVector)
    t = rand(rng) * sum(wv)
    n = length(wv)
    i = 1
    cw = wv[1]
    while cw < t && i < n
        i += 1
        @inbounds cw += wv[i]
    end
    i
end

function sample(rng::AbstractRNG, wa::AbstractArray)
    idx = sample(rng, vec(wa))
    ind2sub(axes(wa), idx)
end

sample(wa::AbstractArray) = sample(Random.GLOBAL_RNG, wa)

const ind2sub = Base._ind2sub


text:
#include("convert.jl")
include("rounding.jl")
include("misc.jl")

#include("addsub.jl")
#include("muldiv.jl")
#include("mlfunctions.jl")

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SQLite"
version = v"3.35.5"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.sqlite.org/2021/sqlite-autoconf-3350500.tar.gz", "f52b72a5c319c3e516ed7a92e123139a6e87af08a2dc43d7757724f6132e6db0"),
    FileSource("https://git.archlinux.org/svntogit/packages.git/plain/trunk/license.txt?h=packages/sqlite&id=33cad63ddb1ba86b7c5a47430c98083ce2b4d86b",
               "4e57d9ac979f1c9872e69799c2597eeef4c6ce7224f3ede0bf9dc8d217b1e65d"; filename="LICENSE"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sqlite-autoconf-*/

# Use same flags as
# https://git.archlinux.org/svntogit/packages.git/tree/trunk/PKGBUILD?h=packages/sqlite&id=d8b6ba561152179e943807054388462b7259e6df
export CPPFLAGS="-DSQLITE_ENABLE_COLUMN_METADATA=1 \
                 -DSQLITE_ENABLE_UNLOCK_NOTIFY \
                 -DSQLITE_ENABLE_DBSTAT_VTAB=1 \
                 -DSQLITE_ENABLE_FTS3_TOKENIZER=1 \
                 -DSQLITE_SECURE_DELETE \
                 -DSQLITE_MAX_VARIABLE_NUMBER=250000 \
                 -DSQLITE_MAX_EXPR_DEPTH=10000"

./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=$target \
    --disable-static \
    --disable-amalgamation \
    --enable-fts3 \
    --enable-fts4 \
    --enable-fts5 \
    --enable-rtree \
    --enable-json1
make -j${nproc}
make install
install_license "${WORKSPACE}/srcdir/LICENSE"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libsqlite3", :libsqlite)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: K-JW/Yggdrasil

max_issues_repo_name: K-JW/Yggdrasil

max_forks_repo_name: K-JW/Yggdrasil


text:
module BlockBandedMatrices
using BlockArrays, BandedMatrices, ArrayLayouts, FillArrays, SparseArrays, MatrixFactorizations
using LinearAlgebra

import Base: getindex, setindex!, checkbounds, @propagate_inbounds, convert,
                        +, *, -, /, \, strides, zeros, size,
                        unsafe_convert, fill!, length, first, last,
                        eltype, getindex, to_indices, to_index,
                        reindex, _maybetail, tail, @_propagate_inbounds_meta,
                        ==, axes, copyto!, similar, OneTo

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                        materialize, materialize!

import LinearAlgebra: UniformScaling, isdiag, rmul!, lmul!, ldiv!, rdiv!,
                        AbstractTriangular, AdjOrTrans, HermOrSym, StructuredMatrixStyle,
                        qr, qr!
import LinearAlgebra.BLAS: BlasInt, BlasFloat, @blasfunc, libblas, BlasComplex, BlasReal
import LinearAlgebra.LAPACK: chktrans, chkdiag, liblapack, chklapackerror, checksquare, chkstride1,
                    chkuplo
import MatrixFactorizations: ql, ql!, _ql, QLPackedQ, AdjQRPackedQLayout, AdjQLPackedQLayout, QR, QRPackedQ
import SparseArrays: sparse

import ArrayLayouts: BlasMatLmulVec, MatLmulVec, MatLmulMat,
                    triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
                    triangulardata, sublayout, sub_materialize,
                    AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
                    DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
                    _qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix

import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
                        Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
                        _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                        blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
                        AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks

import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                        inbands_setindex!, inbands_getindex, banded_setindex!,
                        banded_generic_axpy!,
                        BlasFloat, banded_dense_axpy!, MemoryLayout,
                        BandedLayout, BandedColumnMajor, BandedColumns,
                        BandedSubBandedMatrix, bandeddata,
                        _BandedMatrix, colstart, colstop, rowstart, rowstop,
                        BandedStyle, _fill_lmul!, bandshift,
                        _banded_colval, _banded_rowval, _banded_nzval # for sparse

export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockbandwidth, blockbandwidths,
        subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, BlockBidiagonal, isblockbanded

const Block1 = Block{1,Int}
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}

blockcolrange(A...) = blockcolsupport(A...)
blockrowrange(A...) = blockrowsupport(A...)

include("AbstractBlockBandedMatrix.jl")
include("broadcast.jl")
include("BlockSkylineMatrix.jl")
include("BandedBlockBandedMatrix.jl")

include("linalg.jl")
include("blockskylineqr.jl")

include("interfaceimpl.jl")
include("triblockbanded.jl")
include("adjtransblockbanded.jl")

end # module


text:
using ClimateMachine.Mesh.BrickMesh
using ClimateMachine.Mesh.Grids: mappings, commmapping
using Test
using MPI

MPI.Initialized() || MPI.Init()

@testset "Linear Parition" begin
    @test BrickMesh.linearpartition(1, 1, 1) == 1:1
    @test BrickMesh.linearpartition(20, 1, 1) == 1:20
    @test BrickMesh.linearpartition(10, 1, 2) == 1:5
    @test BrickMesh.linearpartition(10, 2, 2) == 6:10
end

@testset "Hilbert Code" begin
    @test BrickMesh.hilbertcode([0, 0], bits = 1) == [0, 0]
    @test BrickMesh.hilbertcode([0, 1], bits = 1) == [0, 1]
    @test BrickMesh.hilbertcode([1, 1], bits = 1) == [1, 0]
    @test BrickMesh.hilbertcode([1, 0], bits = 1) == [1, 1]
    @test BrickMesh.hilbertcode([0, 0], bits = 2) == [0, 0]
    @test BrickMesh.hilbertcode([1, 0], bits = 2) == [0, 1]
    @test BrickMesh.hilbertcode([1, 1], bits = 2) == [0, 2]
    @test BrickMesh.hilbertcode([0, 1], bits = 2) == [0, 3]
    @test BrickMesh.hilbertcode([0, 2], bits = 2) == [1, 0]
    @test BrickMesh.hilbertcode([0, 3], bits = 2) == [1, 1]
    @test BrickMesh.hilbertcode([1, 3], bits = 2) == [1, 2]
    @test BrickMesh.hilbertcode([1, 2], bits = 2) == [1, 3]
    @test BrickMesh.hilbertcode([2, 2], bits = 2) == [2, 0]
    @test BrickMesh.hilbertcode([2, 3], bits = 2) == [2, 1]
    @test BrickMesh.hilbertcode([3, 3], bits = 2) == [2, 2]
    @test BrickMesh.hilbertcode([3, 2], bits = 2) == [2, 3]
    @test BrickMesh.hilbertcode([3, 1], bits = 2) == [3, 0]
    @test BrickMesh.hilbertcode([2, 1], bits = 2) == [3, 1]
    @test BrickMesh.hilbertcode([2, 0], bits = 2) == [3, 2]
    @test BrickMesh.hilbertcode([3, 0], bits = 2) == [3, 3]

    @test BrickMesh.hilbertcode(UInt64.([14, 3, 4])) ==
          UInt64.([0x0, 0x0, 0xe25])
end

@testset "Mesh to Hilbert Code" begin
    let
        etc = Array{Float64}(undef, 2, 4, 6)
        etc[:, :, 1] = [2.0 3.0 2.0 3.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 2] = [3.0 4.0 3.0 4.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 3] = [4.0 5.0 4.0 5.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 4] = [2.0 3.0 2.0 3.0; 5.0 5.0 6.0 6.0]
        etc[:, :, 5] = [3.0 4.0 3.0 4.0; 5.0 5.0 6.0 6.0]
        etc[:, :, 6] = [4.0 5.0 4.0 5.0; 5.0 5.0 6.0 6.0]

        code_exect = UInt64[
            0x0000000000000000 0x1555555555555555 0xffffffffffffffff 0x5555555555555555 0x6aaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa
            0x0000000000000000 0x5555555555555555 0xffffffffffffffff 0x5555555555555555 0xaaaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa
        ]

        code = centroidtocode(MPI.COMM_SELF, etc)

        @test code == code_exect
    end

    let
        nelem = 1
        d = 2

        etc = Array{Float64}(undef, d, d^2, nelem)
        etc[:, :, 1] = [2.0 3.0 2.0 3.0; 4.0 4.0 5.0 5.0]
        code = centroidtocode(MPI.COMM_SELF, etc)

        @test code == zeros(eltype(code), d, nelem)
    end
end

@testset "Vertex Ordering" begin
    @test ((1,), 1) == BrickMesh.vertsortandorder(1)

    @test ((1, 2), 1) == BrickMesh.vertsortandorder(1, 2)
    @test ((1, 2), 2) == BrickMesh.vertsortandorder(2, 1)

    @test ((1, 2, 3), 1) == BrickMesh.vertsortandorder(1, 2, 3)
    @test ((1, 2, 3), 2) == BrickMesh.vertsortandorder(3, 1, 2)
    @test ((1, 2, 3), 3) == BrickMesh.vertsortandorder(2, 3, 1)
    @test ((1, 2, 3), 4) == BrickMesh.vertsortandorder(2, 1, 3)
    @test ((1, 2, 3), 5) == BrickMesh.vertsortandorder(3, 2, 1)
    @test ((1, 2, 3), 6) == BrickMesh.vertsortandorder(1, 3, 2)

    @test_throws ErrorException BrickMesh.vertsortandorder(2, 1, 1)

    @test ((1, 2, 3, 4), 1) == BrickMesh.vertsortandorder(1, 2, 3, 4)
    @test ((1, 2, 3, 4), 2) == BrickMesh.vertsortandorder(1, 3, 2, 4)
    @test ((1, 2, 3, 4), 3) == BrickMesh.vertsortandorder(2, 1, 3, 4)
    @test ((1, 2, 3, 4), 4) == BrickMesh.vertsortandorder(2, 4, 1, 3)
    @test ((1, 2, 3, 4), 5) == BrickMesh.vertsortandorder(3, 1, 4, 2)
    @test ((1, 2, 3, 4), 6) == BrickMesh.vertsortandorder(3, 4, 1, 2)
    @test ((1, 2, 3, 4), 7) == BrickMesh.vertsortandorder(4, 2, 3, 1)
    @test ((1, 2, 3, 4), 8) == BrickMesh.vertsortandorder(4, 3, 2, 1)

    @test_throws ErrorException BrickMesh.vertsortandorder(1, 3, 3, 1)
end

@testset "Mesh" begin
    let
        (etv, etc, etb, fc) = brickmesh((4:7,), (false,))
        etv_expect = [
            1 2 3
            2 3 4
        ]
        etb_expect = [
            1 0 0
            0 0 1
        ]
        fc_expect = Array{Int64, 1}[]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [4 5]
        @test etc[:, :, 2] == [5 6]
        @test etc[:, :, 3] == [6 7]
    end

    let
        (etv, etc, etb, fc) = brickmesh((4:7,), (true,))
        etv_expect = [
            1 2 3
            2 3 4
        ]
        etb_expect = [
            0 0 0
            0 0 0
        ]
        fc_expect = Array{Int64, 1}[[3, 2, 1]]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [4 5]
        @test etc[:, :, 2] == [5 6]
        @test etc[:, :, 3] == [6 7]
    end

    let
        (etv, etc, etb, fc) = brickmesh((2:5, 4:6), (false, true))

        etv_expect = [
            1 2 5 6
            2 3 6 7
            3 4 7 8
            5 6 9 10
            6 7 10 11
            7 8 11 12
        ]'
        etb_expect = [
            1 0 0 1 0 0
            0 0 1 0 0 1
            0 0 0 0 0 0
            0 0 0 0 0 0
        ]
        fc_expect = Array{Int64, 1}[[4, 4, 1, 2], [5, 4, 2, 3], [6, 4, 3, 4]]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [
            2 3 2 3
            4 4 5 5
        ]
        @test etc[:, :, 5] == [
            3 4 3 4
            5 5 6 6
        ]
    end

    let
        (etv, etc, etb, fc) =
            brickmesh((-1:2:1, -1:2:1, -1:1:1), (true, true, true))
        etv_expect = [
            1 5
            2 6
            3 7
            4 8
            5 9
            6 10
            7 11
            8 12
        ]
        etb_expect = zeros(Int64, 6, 2)

        fc_expect = Array{Int64, 1}[
            [1, 2, 1, 3, 5, 7],
            [1, 4, 1, 2, 5, 6],
            [2, 2, 5, 7, 9, 11],
            [2, 4, 5, 6, 9, 10],
            [2, 6, 1, 2, 3, 4],
        ]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect

        @test etc[:, :, 1] == [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            -1 -1 -1 -1 0 0 0 0
        ]

        @test etc[:, :, 2] == [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            0 0 0 0 1 1 1 1
        ]
    end

    let
        (etv, etc, etb, fc) = brickmesh(
            (-1:1, -1:1, -1:1),
            (false, false, false),
            boundary = ((11, 12), (13, 14), (15, 16)),
        )

        @test etb == [
            11 0 11 0 11 0 11 0
            0 12 0 12 0 12 0 12
            13 13 0 0 13 13 0 0
            0 0 14 14 0 0 14 14
            15 15 15 15 0 0 0 0
            0 0 0 0 16 16 16 16
        ]
    end

    let
        x = (1:1000,)
        p = (false,)
        b = ((1, 2),)

        (etv, etc, etb, fc) = brickmesh(x, p, boundary = b)

        n = 50
        (etv_parts, etc_parts, etb_parts, fc_parts) =
            brickmesh(x, p, boundary = b, part = 1, numparts = n)
        for j in 2:n
            (etv_j, etc_j, etb_j, fc_j) =
                brickmesh(x, p, boundary = b, part = j, numparts = n)
            etv_parts = cat(etv_parts, etv_j; dims = 2)
            etc_parts = cat(etc_parts, etc_j; dims = 3)
            etb_parts = cat(etb_parts, etb_j; dims = 2)
        end

        @test etv == etv_parts
        @test etc == etc_parts
        @test etb == etb_parts
    end

    let
        x = (-1:2:10, -1:1:1, -4:1:1)
        p = (true, false, true)
        b = ((1, 2), (3, 4), (5, 6))

        (etv, etc, etb, fc) = brickmesh(x, p, boundary = b)

        n = 50
        (etv_parts, etc_parts, etb_parts, fc_parts) =
            brickmesh(x, p, boundary = b, part = 1, numparts = n)
        for j in 2:n
            (etv_j, etc_j, etb_j, fc_j) =
                brickmesh(x, p, boundary = b, part = j, numparts = n)
            etv_parts = cat(etv_parts, etv_j; dims = 2)
            etc_parts = cat(etc_parts, etc_j; dims = 3)
            etb_parts = cat(etb_parts, etb_j; dims = 2)
        end

        @test etv == etv_parts
        @test etc == etc_parts
        @test etb == etb_parts
    end
end

@testset "Connect" begin
    let
        comm = MPI.COMM_SELF

        mesh = connectmesh(
            comm,
            partition(comm, brickmesh((0:10,), (true,))...)[1:4]...,
        )

        nelem = 10

        @test mesh[:elemtocoord][:, :, 1] == [0 1]
        @test mesh[:elemtocoord][:, :, 2] == [1 2]
        @test mesh[:elemtocoord][:, :, 3] == [2 3]
        @test mesh[:elemtocoord][:, :, 4] == [3 4]
        @test mesh[:elemtocoord][:, :, 5] == [4 5]
        @test mesh[:elemtocoord][:, :, 6] == [5 6]
        @test mesh[:elemtocoord][:, :, 7] == [6 7]
        @test mesh[:elemtocoord][:, :, 8] == [7 8]
        @test mesh[:elemtocoord][:, :, 9] == [8 9]
        @test mesh[:elemtocoord][:, :, 10] == [9 10]

        @test mesh[:elemtoelem] == [
            10 1 2 3 4 5 6 7 8 9
            2 3 4 5 6 7 8 9 10 1
        ]

        @test mesh[:elemtoface] == [
            2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1
        ]

        @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))
        @test mesh[:elemtobndy] == zeros(Int, size(mesh[:elemtoordr]))

        @test mesh[:elems] == 1:nelem
        @test mesh[:realelems] == 1:nelem
        @test mesh[:ghostelems] == nelem .+ (1:0)

        @test length(mesh[:sendelems]) == 0

        @test mesh[:nabrtorank] == Int[]
        @test mesh[:nabrtorecv] == UnitRange{Int}[]
        @test mesh[:nabrtosend] == UnitRange{Int}[]
    end

    let
        comm = MPI.COMM_SELF
        mesh = connectmesh(
            comm,
            partition(comm, brickmesh((0:4, 5:9), (false, true))...)[1:4]...,
        )

        nelem = 16

        @test mesh[:elemtocoord][:, :, 1] == [0 1 0 1; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 2] == [1 2 1 2; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 3] == [1 2 1 2; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 4] == [0 1 0 1; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 5] == [0 1 0 1; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 6] == [0 1 0 1; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 7] == [1 2 1 2; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 8] == [1 2 1 2; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 9] == [2 3 2 3; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 10] == [2 3 2 3; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 11] == [3 4 3 4; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 12] == [3 4 3 4; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 13] == [3 4 3 4; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 14] == [2 3 2 3; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 15] == [2 3 2 3; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 16] == [3 4 3 4; 5 5 6 6]

        @test mesh[:elemtoelem] == [
            1 1 4 4 5 6 6 5 8 7 10 9 14 3 2 15
            2 15 14 3 8 7 10 9 12 11 11 12 13 13 16 16
            6 7 2 1 4 5 8 3 14 9 12 13 16 15 10 11
            4 3 8 5 6 1 2 7 10 15 16 11 12 9 14 13
        ]

        @test mesh[:elemtoface] == [
            1 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
        ]

        @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))

        @test mesh[:elems] == 1:nelem
        @test mesh[:realelems] == 1:nelem
        @test mesh[:ghostelems] == nelem .+ (1:0)

        @test length(mesh[:sendelems]) == 0

        @test mesh[:nabrtorank] == Int[]
        @test mesh[:nabrtorecv] == UnitRange{Int}[]
        @test mesh[:nabrtosend] == UnitRange{Int}[]
    end

    let
        comm = MPI.COMM_SELF
        mesh = connectmeshfull(
            comm,
            partition(comm, brickmesh((0:4, 5:9), (false, true))...)[1:4]...,
        )

        nelem = 16

        @test mesh[:elemtocoord][:, :, 1] == [0 1 0 1; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 2] == [1 2 1 2; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 3] == [1 2 1 2; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 4] == [0 1 0 1; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 5] == [0 1 0 1; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 6] == [0 1 0 1; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 7] == [1 2 1 2; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 8] == [1 2 1 2; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 9] == [2 3 2 3; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 10] == [2 3 2 3; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 11] == [3 4 3 4; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 12] == [3 4 3 4; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 13] == [3 4 3 4; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 14] == [2 3 2 3; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 15] == [2 3 2 3; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 16] == [3 4 3 4; 5 5 6 6]

        @test mesh[:elemtoelem] == [
            1 1 4 4 5 6 6 5 8 7 10 9 14 3 2 15
            2 15 14 3 8 7 10 9 12 11 11 12 13 13 16 16
            6 7 2 1 4 5 8 3 14 9 12 13 16 15 10 11
            4 3 8 5 6 1 2 7 10 15 16 11 12 9 14 13
        ]

        @test mesh[:elemtoface] == [
            1 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
        ]

        @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))

        @test mesh[:elems] == 1:nelem
        @test mesh[:realelems] == 1:nelem
        @test mesh[:ghostelems] == nelem .+ (1:0)

        @test length(mesh[:sendelems]) == 0

        @test mesh[:nabrtorank] == Int[]
        @test mesh[:nabrtorecv] == UnitRange{Int}[]
        @test mesh[:nabrtosend] == UnitRange{Int}[]
    end
end

@testset "Mappings" begin
    let
        comm = MPI.COMM_SELF
        x = (0:4,)
        mesh =
            connectmesh(comm, partition(comm, brickmesh(x, (true,))...)[1:4]...)

        N = 3
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ = mappings(
            ntuple(j -> N, d),
            mesh[:elemtoelem],
            mesh[:elemtoface],
            mesh[:elemtoordr],
        )

        @test vmap⁻ == reshape([1, 4, 5, 8, 9, 12, 13, 16], Nfp, nface, nelem)
        @test vmap⁺ == reshape([16, 5, 4, 9, 8, 13, 12, 1], Nfp, nface, nelem)
    end

    # Single polynomial order
    let
        comm = MPI.COMM_SELF
        x = (-1:1, 0:1)
        p = (false, true)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = 2
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ = mappings(
            ntuple(j -> N, d),
            mesh[:elemtoelem],
            mesh[:elemtoface],
            mesh[:elemtoordr],
        )

        #! format: off
        @test vmap⁻ == reshape(
            [
                 1,  4,  7,  # f=1 e=1
                 3,  6,  9,  # f=2 e=1
                 1,  2,  3,  # f=3 e=1
                 7,  8,  9,  # f=4 e=1
                10, 13, 16,  # f=1 e=2
                12, 15, 18,  # f=2 e=2
                10, 11, 12,  # f=3 e=2
                16, 17, 18,  # f=4 e=2
            ],
            Nfp,
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                 1,  4,  7,  # f=1 e=1
                10, 13, 16,  # f=1 e=2
                 7,  8,  9,  # f=4 e=1
                 1,  2,  3,  # f=3 e=1
                 3,  6,  9,  # f=2 e=1
                12, 15, 18,  # f=2 e=2
                16, 17, 18,  # f=4 e=2
                10, 11, 12,  # f=3 e=2
            ],
            Nfp,
            nface,
            nelem,
        )
        #! format: on
    end

    # Two polynomial orders
    let
        comm = MPI.COMM_SELF
        x = (-1:1, 0:1)
        p = (false, true)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = (2, 3)
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Nfp = div.(prod(N .+ 1), ((N .+ 1) .^ (d - 1)))

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        #! format: off
        @test vmap⁻ == reshape(
            [
                 1,  4,  7, 10,  # f=1 e=1
                 3,  6,  9, 12,  # f=2 e=1
                 1,  2,  3,  0,  # f=3 e=1
                10, 11, 12,  0,  # f=4 e=1
                13, 16, 19, 22,  # f=1 e=2
                15, 18, 21, 24,  # f=2 e=2
                13, 14, 15,  0,  # f=3 e=2
                22, 23, 24,  0,  # f=4 e=2
            ],
            maximum(Nfp),
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                 1,  4,  7, 10,  # f=1 e=1
                13, 16, 19, 22,  # f=2 e=1
                10, 11, 12,  0,  # f=3 e=1
                 1,  2,  3,  0,  # f=4 e=1
                 3,  6,  9, 12,  # f=1 e=2
                15, 18, 21, 24,  # f=2 e=2
                22, 23, 24,  0,  # f=3 e=2
                13, 14, 15,  0,  # f=4 e=2
            ],
            maximum(Nfp),
            nface,
            nelem,
        )
        #! format: on
    end

    # Single polynomial order
    let
        comm = MPI.COMM_SELF
        x = (0:1, 0:1, -1:1)
        p = (false, true, false)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = 2
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Np = (N + 1)^d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ = mappings(
            ntuple(j -> N, d),
            mesh[:elemtoelem],
            mesh[:elemtoface],
            mesh[:elemtoordr],
        )

        #! format: off
        fmask = [
             1  3  1  7 1 19
             4  6  2  8 2 20
             7  9  3  9 3 21
            10 12 10 16 4 22
            13 15 11 17 5 23
            16 18 12 18 6 24
            19 21 19 25 7 25
            22 24 20 26 8 26
            25 27 21 27 9 27
        ]
        #! format: on

        @test vmap⁻ == reshape([fmask[:]; fmask[:] .+ Np], Nfp, nface, nelem)

        @test vmap⁺ == reshape(
            [
                fmask[:, 1]
                fmask[:, 2]
                fmask[:, 4]
                fmask[:, 3]
                fmask[:, 5]
                fmask[:, 5] .+ Np
                fmask[:, 1] .+ Np
                fmask[:, 2] .+ Np
                fmask[:, 4] .+ Np
                fmask[:, 3] .+ Np
                fmask[:, 6]
                fmask[:, 6] .+ Np
            ],
            Nfp,
            nface,
            nelem,
        )
    end

    # Multiple polynomial orders
    let
        comm = MPI.COMM_SELF
        x = (0:1, 0:1, -1:1)
        p = (false, true, false)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = (1, 2, 3)
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Np = prod(N .+ 1)
        Nfp = div.(Np, N .+ 1)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        #! format: off
        fmask =
        (
         [ 1  3  5  7  9 11 13 15 17 19 21 23],
         [ 2  4  6  8 10 12 14 16 18 20 22 24],
         [ 1  2  7  8 13 14 19 20  0  0  0  0],
         [ 5  6 11 12 17 18 23 24  0  0  0  0],
         [ 1  2  3  4  5  6  0  0  0  0  0  0],
         [19 20 21 22 23 24  0  0  0  0  0  0],
        )
        #! format: on

        @test vmap⁻ == reshape(
            [
                fmask[1]
                fmask[2]
                fmask[3]
                fmask[4]
                fmask[5]
                fmask[6]
                fmask[1] .+ Np .* (fmask[1] .> 0)
                fmask[2] .+ Np .* (fmask[2] .> 0)
                fmask[3] .+ Np .* (fmask[3] .> 0)
                fmask[4] .+ Np .* (fmask[4] .> 0)
                fmask[5] .+ Np .* (fmask[5] .> 0)
                fmask[6] .+ Np .* (fmask[6] .> 0)
            ]',
            maximum(Nfp),
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                vmap⁻[:, 1, 1]
                vmap⁻[:, 2, 1]
                vmap⁻[:, 4, 1]
                vmap⁻[:, 3, 1]
                vmap⁻[:, 5, 1]
                vmap⁻[:, 5, 2]
                vmap⁻[:, 1, 2]
                vmap⁻[:, 2, 2]
                vmap⁻[:, 4, 2]
                vmap⁻[:, 3, 2]
                vmap⁻[:, 6, 1]
                vmap⁻[:, 6, 2]
            ],
            maximum(Nfp),
            nface,
            nelem,
        )
    end

    # Test polynomial order 0
    let
        comm = MPI.COMM_SELF
        x = (0:1, -1:1)
        p = (false, true)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = (5, 0)
        Nq = N .+ 1
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        #! format: off
        p = reshape(1:Np, Nq)
        fmask = ntuple(j -> zeros(Int, maximum(Nfp)), 2d)
        fmask[1][1:Nfp[1]] .= p[1, :][:]
        fmask[2][1:Nfp[1]] .= p[end, :][:]
        fmask[3][1:Nfp[2]] .= p[:, 1][:]
        fmask[4][1:Nfp[2]] .= p[:, end][:]
        #! format: on

        @test vmap⁻ == reshape(
            [
                fmask[1]
                fmask[2]
                fmask[3]
                fmask[4]
                fmask[1] .+ Np .* (fmask[1] .> 0)
                fmask[2] .+ Np .* (fmask[2] .> 0)
                fmask[3] .+ Np .* (fmask[3] .> 0)
                fmask[4] .+ Np .* (fmask[4] .> 0)
            ]',
            maximum(Nfp),
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                vmap⁻[:, 1, 1]
                vmap⁻[:, 2, 1]
                vmap⁻[:, 4, 2]
                vmap⁻[:, 3, 2]
                vmap⁻[:, 1, 2]
                vmap⁻[:, 2, 2]
                vmap⁻[:, 4, 1]
                vmap⁻[:, 3, 1]
            ],
            maximum(Nfp),
            nface,
            nelem,
        )
    end

    # Test polynomial order 0
    let
        comm = MPI.COMM_SELF
        x = (0:1, 0:1, -1:1)
        p = (false, true, false)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = (7, 2, 0)
        Nq = N .+ 1
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        #! format: off
        p = reshape(1:Np, Nq)
        fmask = ntuple(j -> zeros(Int, maximum(Nfp)), 2d)
        fmask[1][1:Nfp[1]] .= p[1, :, :][:]
        fmask[2][1:Nfp[1]] .= p[end, :, :][:]
        fmask[3][1:Nfp[2]] .= p[:, 1,  :][:]
        fmask[4][1:Nfp[2]] .= p[:, end,  :][:]
        fmask[5][1:Nfp[3]] .= p[:, :, 1][:]
        fmask[6][1:Nfp[3]] .= p[:, :, end][:]
        #! format: on

        @test vmap⁻ == reshape(
            [
                fmask[1]
                fmask[2]
                fmask[3]
                fmask[4]
                fmask[5]
                fmask[6]
                fmask[1] .+ Np .* (fmask[1] .> 0)
                fmask[2] .+ Np .* (fmask[2] .> 0)
                fmask[3] .+ Np .* (fmask[3] .> 0)
                fmask[4] .+ Np .* (fmask[4] .> 0)
                fmask[5] .+ Np .* (fmask[5] .> 0)
                fmask[6] .+ Np .* (fmask[6] .> 0)
            ]',
            maximum(Nfp),
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                vmap⁻[:, 1, 1]
                vmap⁻[:, 2, 1]
                vmap⁻[:, 4, 1]
                vmap⁻[:, 3, 1]
                vmap⁻[:, 5, 1]
                vmap⁻[:, 5, 2]
                vmap⁻[:, 1, 2]
                vmap⁻[:, 2, 2]
                vmap⁻[:, 4, 2]
                vmap⁻[:, 3, 2]
                vmap⁻[:, 6, 1]
                vmap⁻[:, 6, 2]
            ],
            maximum(Nfp),
            nface,
            nelem,
        )
    end
end

@testset "Get Partition" begin
    let
        Nelem = 150
        (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
        @test so == Nelem:-1:1
        @test ss == [1, Nelem + 1]
        @test rs == [1, Nelem + 1]
    end

    let
        Nelem = 111
        code = [ones(1, Nelem); collect(Nelem:-1:1)']
        (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
        @test so == Nelem:-1:1
        @test ss == [1, Nelem + 1]
        @test rs == [1, Nelem + 1]
    end
end

@testset "Partition" begin
    (etv, etc, etb, fc) =
        brickmesh((-1:2:1, -1:2:1, -2:1:2), (true, true, true))
    (netv, netc, netb, nfc) = partition(MPI.COMM_SELF, etv, etc, etb, fc)[1:4]
    @test etv == netv
    @test etc == netc
    @test etb == netb
    @test fc == nfc
end

@testset "Comm Mappings" begin
    let
        N = 1
        d = 2
        nface = 2d

        commelems = [1, 2, 5]
        commfaces = BitArray(undef, nface, length(commelems))
        commfaces .= false
        nabrtocomm = [1:2, 3:3]

        vmapC, nabrtovmapC =
            commmapping(ntuple(j -> N, d), commelems, commfaces, nabrtocomm)

        @test vmapC == Int[]
        @test nabrtovmapC == UnitRange{Int64}[1:0, 1:0]
    end

    let
        N = 1
        d = 2
        nface = 2d

        commelems = [1, 2, 5]
        commfaces = BitArray([
            false false false
            false true false
            false true false
            false false true
        ])
        nabrtocomm = [1:2, 3:3]

        vmapC, nabrtovmapC =
            commmapping(ntuple(j -> N, d), commelems, commfaces, nabrtocomm)

        @test vmapC == [5, 6, 8, 19, 20]
        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:5]
    end

    # 2D, single polynomial order
    let
        N = 2
        d = 2
        nface = 2d

        commelems = [2, 4, 5]
        commfaces = BitArray([
            true true false
            false false false
            false true true
            false false true
        ])
        nabrtocomm = [1:1, 2:3]

        vmapC, nabrtovmapC =
            commmapping(ntuple(j -> N, d), commelems, commfaces, nabrtocomm)

        @test vmapC == [10, 13, 16, 28, 29, 30, 31, 34, 37, 38, 39, 43, 44, 45]
        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:14]
    end

    # 2D, multiple polynomial orders
    let
        N = (2, 3)
        Np = prod(N .+ 1)
        d = length(N)
        nface = 2d

        commelems = [2, 4, 5]
        #! format: off
        commfaces = BitArray([
             true  true false false # faces of commelems[1] to send
            false false false  true # faces of commelems[2] to send
             true false false  true # faces of commelems[3] to send
        ]')
        #! format: on
        nabrtocomm = [1:1, 2:3]
        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        @test vmapC == [
            (commelems[1] - 1) * Np + 1
            (commelems[1] - 1) * Np + 3
            (commelems[1] - 1) * Np + 4
            (commelems[1] - 1) * Np + 6
            (commelems[1] - 1) * Np + 7
            (commelems[1] - 1) * Np + 9
            (commelems[1] - 1) * Np + 10
            (commelems[1] - 1) * Np + 12
            (commelems[2] - 1) * Np + 10
            (commelems[2] - 1) * Np + 11
            (commelems[2] - 1) * Np + 12
            (commelems[3] - 1) * Np + 1
            (commelems[3] - 1) * Np + 4
            (commelems[3] - 1) * Np + 7
            (commelems[3] - 1) * Np + 10
            (commelems[3] - 1) * Np + 11
            (commelems[3] - 1) * Np + 12
        ]

        @test nabrtovmapC == UnitRange{Int64}[1:8, 9:17]
    end

    # 3D, single polynomial order
    let
        N = 2
        d = 3
        nface = 2d

        commelems = [3, 4, 7, 9]
        commfaces = BitArray([
            true true true false
            false true false false
            false true false false
            false true false false
            false true true true
            false true false true
        ])
        nabrtocomm = [1:1, 2:4]

        vmapC, nabrtovmapC =
            commmapping(ntuple(j -> N, d), commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC == [
            55, 58, 61, 64, 67, 70, 73, 76, 79, 82, 83, 84, 85, 86, 87, 88, 89,
            90, 91, 92, 93, 94, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105,
            106, 107, 108, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172,
            175, 178, 181, 184, 187, 217, 218, 219, 220, 221, 222, 223, 224,
            225, 235, 236, 237, 238, 239, 240, 241, 242, 243,
        ]
        #! format: on
        @test nabrtovmapC == UnitRange{Int64}[1:9, 10:68]
    end

    # 3D, multiple polynomial order
    let
        N = (2, 3, 4)
        Nq = N .+ 1
        Np = prod(Nq)
        d = length(N)
        nface = 2d

        commelems = [3, 4, 7, 9]
        commfaces = BitArray([
            true true true false
            false true false false
            false true false false
            false true false false
            false true true true
            false true false true
        ])

        p = reshape(1:Np, Nq)
        fmask = (
            p[1, :, :][:],     # Face 1
            p[Nq[1], :, :][:], # Face 2
            p[:, 1, :][:],     # Face 3
            p[:, Nq[2], :][:], # Face 4
            p[:, :, 1][:],     # Face 5
            p[:, :, Nq[3]][:], # Face 6
        )

        nabrtocomm = [1:1, 2:4]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC ==
        [
         sort(unique(vcat(fmask[commfaces[:, 1]]...))) .+ (commelems[1] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 2]]...))) .+ (commelems[2] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 3]]...))) .+ (commelems[3] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 4]]...))) .+ (commelems[4] - 1) *Np
        ]
        #! format: on

        @test nabrtovmapC == UnitRange{Int64}[1:20, 21:126]
    end

    # Test mappings with polyorder 0 in one dimension
    let
        N = (0, 1, 2)
        Nq = N .+ 1
        Np = prod(Nq)
        d = length(N)
        nface = 2d

        commelems = [3, 4, 7, 9]
        commfaces = BitArray([
            true true true false
            false true false false
            false true false false
            false true false false
            false true true true
            false true false true
        ])

        p = reshape(1:Np, Nq)
        fmask = (
            p[1, :, :][:],     # Face 1
            p[Nq[1], :, :][:], # Face 2
            p[:, 1, :][:],     # Face 3
            p[:, Nq[2], :][:], # Face 4
            p[:, :, 1][:],     # Face 5
            p[:, :, Nq[3]][:], # Face 6
        )

        nabrtocomm = [1:1, 2:4]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC ==
        [
         sort(unique(vcat(fmask[commfaces[:, 1]]...))) .+ (commelems[1] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 2]]...))) .+ (commelems[2] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 3]]...))) .+ (commelems[3] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 4]]...))) .+ (commelems[4] - 1) *Np
        ]
        #! format: on

        @test nabrtovmapC == UnitRange{Int64}[1:6, 7:22]
    end

    let
        N = (3, 2, 0)
        Nq = N .+ 1
        Np = prod(Nq)
        d = length(N)
        nface = 2d

        commelems = [3, 4, 7, 9, 19]
        commfaces = BitArray([
            true true true false true
            false true false false false
            false true false false true
            false true false false false
            false true true true false
            false true false true false
        ])

        p = reshape(1:Np, Nq)
        fmask = (
            p[1, :, :][:],     # Face 1
            p[Nq[1], :, :][:], # Face 2
            p[:, 1, :][:],     # Face 3
            p[:, Nq[2], :][:], # Face 4
            p[:, :, 1][:],     # Face 5
            p[:, :, Nq[3]][:], # Face 6
        )

        nabrtocomm = [1:1, 2:5]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC ==
        [
         sort(unique(vcat(fmask[commfaces[:, 1]]...))) .+ (commelems[1] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 2]]...))) .+ (commelems[2] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 3]]...))) .+ (commelems[3] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 4]]...))) .+ (commelems[4] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 5]]...))) .+ (commelems[5] - 1) *Np
        ]
        #! format: on

        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:45]
    end

    let
        N = (0, 2)
        Nq = N .+ 1
        Np = prod(Nq)
        d = length(N)
        nface = 2d

        commelems = [3, 4, 7, 9]
        commfaces = BitArray([
            true true true false
            false true false true
            false true false false
            false true false true
        ])

        p = reshape(1:Np, Nq)
        fmask = (
            p[1, :][:],     # Face 1
            p[Nq[1], :][:], # Face 2
            p[:, 1][:],     # Face 3
            p[:, Nq[2]][:], # Face 4
        )

        nabrtocomm = [1:1, 2:4]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC ==
        [
         sort(unique(vcat(fmask[commfaces[:, 1]]...))) .+ (commelems[1] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 2]]...))) .+ (commelems[2] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 3]]...))) .+ (commelems[3] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 4]]...))) .+ (commelems[4] - 1) *Np
        ]
        #! format: on

        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:12]
    end

    let
        N = (3, 0)
        Nq = N .+ 1
        Np = prod(Nq)
        d = length(N)
        nface = 2d

        commelems = [3, 4, 7, 9, 19]
        commfaces = BitArray([
            true true true false true
            false true false false true
            false true false true false
            false true false true false
        ])

        p = reshape(1:Np, Nq)
        fmask = (
            p[1, :][:],     # Face 1
            p[Nq[1], :][:], # Face 2
            p[:, 1][:],     # Face 3
            p[:, Nq[2]][:], # Face 4
        )

        nabrtocomm = [1:1, 2:5]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        #! format: off
        @test vmapC ==
        [
         sort(unique(vcat(fmask[commfaces[:, 1]]...))) .+ (commelems[1] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 2]]...))) .+ (commelems[2] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 3]]...))) .+ (commelems[3] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 4]]...))) .+ (commelems[4] - 1) *Np
         sort(unique(vcat(fmask[commfaces[:, 5]]...))) .+ (commelems[5] - 1) *Np
        ]
        #! format: on

        @test nabrtovmapC == UnitRange{Int64}[1:1, 2:12]
    end

end


text:
################################################################################
#
#   Fields.jl : generic fields
#
################################################################################

include("julia/GF.jl")

include("flint/fmpq.jl")

include("flint/fq.jl")

include("flint/fq_nmod.jl")

include("antic/nf_elem.jl")

include("arb/arb.jl")

include("arb/acb.jl")

//(a::T, b::T) where {T <: FieldElem} = divexact(a, b)

//(x::T, y::Union{Integer, Rational}) where {T <: RingElem} = x//parent(x)(y)

//(x::Union{Integer, Rational}, y::T) where {T <: RingElem} = parent(y)(x)//y

function gcd(x::T, y::T) where {T <: FieldElem}
   check_parent(x, y)
   return iszero(x) && iszero(y) ? zero(parent(y)) : one(parent(y))
end

characteristic(R::Field) = 0


text:
module SimilaritySearch
abstract type Index end

export Index

include("distances/bits.jl")
include("distances/sets.jl")
include("distances/strings.jl")
include("distances/vectors.jl")
include("distances/cos.jl")
include("res/knn.jl")
include("nns/performance.jl")
include("indexes/pivotselection.jl")
include("indexes/seq.jl")
include("indexes/laesa.jl")
include("indexes/pivotselectiontables.jl")
include("knr/SimilarReferences.jl")
include("graph/graph.jl")
end


text:
"""

```
minimize!(sys::System{T},opt::Options)
```

Minimizes the potential energy of `sys.x0`. Modifies `sys.x0`.

### Example
```julia-repl
julia> sys = System(n=100,sides=[100,100])

julia> minimize!(sys,Options(tol=1e-3))
Energy before minimization: 38322.72337856496
Energy after minimization: -74.15646912098042

```

"""
function minimize!(sys::System{T},opt::Options=Options()) where T
  # Simplify code by aliasing common variables
  @unpack n, x0, sides = sys
  x = x0
  u(x) = potential(x,sys,opt)
  f!(f,x) = forces!(f,x,sys,opt)

  # Minimizing the energy of the initial point
  ubest = u(x)
  println("Energy before minimization: ", ubest)

  dx = 1.
  fnorm = 1.
  xtrial = zeros(T,n)
  f = zeros(T,n)
  trial = 0
  while fnorm > 1.e-5 && trial < 20_000
    trial += 1

    # Compute gradient (≡ -forces)
    f!(f,x)
    fnorm = 0.
    for v in f
      fnorm += norm2(v)
    end
    fnorm = fnorm/n

    # Compute trial point ( xtrial = x - (dU/dx)*dx = x + f*dx )
    for i in 1:n
      xtrial[i] = x[i] + f[i]*dx
    end
    utrial = u(xtrial)

    # If energy decreased, accept new point, if not, reject it and decrease dx
    if utrial < ubest
      if opt.minprint
        @printf(
          "ACCEPTED U = %12.5e FNORM = %12.5e DX = %12.5e TRIAL = %10i\n",
          utrial, fnorm, dx, trial
        )
      end
      for i in 1:n
        x[i] = xtrial[i]
      end
      ubest = utrial
      dx = dx * 2
    else
      dx = dx / 2
    end
  end

  for i in 1:n
    x[i] = image(x[i],sides)
  end

  println("Energy after minimization: ", u(x))
  return nothing
end


text:
#
module SEM

using LinearAlgebra,SparseArrays
using Zygote,Plots
using UnPack
using StaticArrays

import FastGaussQuadrature
import Zygote,NNlib
#--------------------------------------#
Base.:*(op::Function,x::Array) = op(x)
#--------------------------------------#
export linspace
linspace(zi::Number,ze::Number,n::Integer) = Array(range(zi,stop=ze,length=n))
#--------------------------------------#
export iscallable
iscallable(op) = !isempty(methods(op))
#--------------------------------------#
sum(A::Array, n::Int) = Base.sum(A, dims=n)
sum(A) = Base.sum(A)
#cumprod(A::Array) = Base.cumprod(A, dims=1)
#cumprod(A::Array, d::Int) = Base.cumprod(A, dims=d)
#flipdim(A, d) = reverse(A, dims=d)
#--------------------------------------#

# SEM building blocks
include("jac.jl")
include("interp.jl")
include("derivMat.jl")
include("semmesh.jl")
include("semq.jl")
include("ndgrid.jl")

include("ABu.jl")
include("geom.jl")
include("mask.jl")

include("time.jl")
include("mesh.jl")

include("pcg.jl")
include("adjoint.jl")

include("gatherScatter.jl")
include("grad.jl")
include("mass.jl")
include("lapl.jl")
include("hlmz.jl")
include("advect.jl")
include("diver.jl")

include("diffusion.jl")
include("convectionDiffusion.jl")
#include("stokes.jl")

include("plt.jl")

end # module
#--------------------------------------#


text:
# tensor.jl
#
# Tensor provides a dense implementation of an AbstractTensor type without any
# symmetry assumptions, i.e. it describes tensors living in the full tensor
# product space of its index spaces.

#++++++++++++++
# Tensor type:
#++++++++++++++
# Type definition and constructors:
#-----------------------------------
struct TensorMap{T<:Number, S<:IndexSpace, N₁, N₂, A} <: AbstractTensorMap{T, S, N₁, N₂}
    data::A
    codom::ProductSpace{S,N₁}
    dom::ProductSpace{S,N₂}
    function TensorMap{T,S,N₁,N₂}(data, spaces::TensorMapSpace{S,N₁,N₂}) where {T<:Number, S<:IndexSpace, N₁, N₂}
        codom = spaces[2]
        dom = spaces[1]
        data2 = validatedata(data, codom, dom, fieldtype(S), sectortype(S))
        new{T,S,N₁,N₂,typeof(data2)}(data2, codom, dom)
    end
    function TensorMap{T,S,N,0}(data, space::TensorSpace{S,N}) where {T<:Number, S<:IndexSpace, N}
        codom = space
        dom = ProductSpace{S,0}(())
        data2 = validatedata(data, codom, dom, fieldtype(S), sectortype(S))
        new{T,S,N,0,typeof(data2)}(data2, codom, dom)
    end
end

const Tensor{T<:Number, S<:IndexSpace, N, A<:AbstractArray{T}} = TensorMap{T, S, N, 0, A}

function validatedata(data::AbstractArray, codom, dom, k::Field, ::Type{Trivial})
    size(data) == (dims(codom)..., dims(dom)...) || throw(DimensionMismatch())
    eltype(data) ⊆ k || warn("eltype(data) = $(eltype(data)) ⊈ $k)")
    return data
end
function validatedata(data::Dict{<:FusionTree{G,N₁,N₂},<:AbstractArray}, codom::ProductSpace{S,N₁}, dom::ProductSpace{S,N₁}, k::Field, ::Type{G}) where {S<:IndexSpace, G<:Sector, N₁,N₂}
    for tree in keys(data)
        sectorcodom = tree.outgoing
        sectordom = tree.incoming
        size(data[tree]) == (dims(codom, sectorcodom)..., dims(dom, sectordom)...) || throw(DimensionMismatch())
        eltype(data[tree]) ⊆ k || warn("eltype(data) = $(eltype(data)) ⊈ $k)")
    end
    return data
end
# TODO: allow to start from full data (a single AbstractArray) and create the dictionary, in the first place for Abelian sectors, or for e.g. SU₂ using Wigner 3j symbols

# Show method:
#-------------
function Base.show(io::IO, t::TensorMap{T,S,N₁,N₂}) where {T<:Number, S<:IndexSpace, N₁, N₂}
    print(io, "TensorMap{", T, "}( ", codom(t), " ← ", dom(t), " )")
    if !get(io, :compact, false)
        println(":")
        Base.showarray(io, t.data, false; header = false)
    end
end
function Base.show(io::IO, t::Tensor{T,S,N}) where {T<:Number, S<:IndexSpace, N}
    print(io," Tensor{", T, "}( ", codom(t), " )")
    if !get(io, :compact, false)
        println(":")
        Base.showarray(io, t.data, false; header = false)
    end
end

# Basic methods for characterising a tensor:
#--------------------------------------------
codom(t::TensorMap) = t.codom
dom(t::TensorMap) = t.dom
space(t::TensorMap) = t.codom ⊗ dual(t.dom)
space(t::TensorMap, n::Int) = space(t::TensorMap)[n]

# General Tensor constructors
#----------------------------
# with data
Tensor{T}(data::AbstractArray{T}, P::TensorSpace{S,N}) where {T<:Number, S<:IndexSpace, N} = Tensor{T,S,N}(data, P)

Tensor(data::AbstractArray{T}, P::TensorSpace) where {T<:Number} = Tensor{T}(data, P)

# without data: uninitialized
Tensor{T}(A::Type{<:AbstractArray{T}}, P::TensorSpace) where {T<:Number} = Tensor{T}(A(dims(P)), P)
Tensor(A::Type{<:AbstractArray{T}}, P::TensorSpace) where {T<:Number} = Tensor{T}(A, P)

Tensor{T}(P::TensorSpace) where {T<:Number} = Tensor(Array{T}, P)
Tensor(P::TensorSpace) = Tensor{Float64}(P)

# generic constructor from callable:
Tensor{T}(f, P::TensorSpace) where {T<:Number} = Tensor{T}(f(T, dims(P)), P)
Tensor(f, P::TensorSpace) = Tensor(f(dims(P)), P)

# for a single `IndexSpace`
Tensor{T}(arg, V::IndexSpace) where {T<:Number} = Tensor{T}(arg, TensorSpace(V))
Tensor(arg, V::IndexSpace) = Tensor(arg, TensorSpace(V))

Tensor(data::AbstractArray{T,N}) where {T<:Real, N} =
    Tensor{T}(data, TensorSpace{CartesianSpace,N}(map(CartesianSpace, size(data))))
Tensor(data::AbstractArray{T,N}) where {T<:Complex, N} =
    Tensor{T}(data, TensorSpace{ComplexSpace,N}(map(ComplexSpace, size(data))))

# General TensorMap constructors
#--------------------------------
# with data
TensorMap{T}(data::AbstractArray{T}, P::TensorMapSpace{S,N₁,N₂}) where {T<:Number, S<:IndexSpace, N₁, N₂} = TensorMap{T,S,N₁,N₂}(data, P)
TensorMap(data::AbstractArray{T}, P::TensorMapSpace) where {T<:Number} = Tensor{T}(data, P)

# without data: uninitialized
TensorMap{T}(A::Type{<:AbstractArray{T}}, P::TensorMapSpace) where {T<:Number} = Tensor{T}(A(dims(P)), P)
TensorMap(A::Type{<:AbstractArray{T}}, P::TensorMapSpace) where {T<:Number} = Tensor{T}(A, P)

TensorMap{T}(P::TensorMapSpace) where {T<:Number} = TensorMap{T}(Array{T}, P)
TensorMap(P::TensorMapSpace) = TensorMap{Float64}(P)

# generic constructor from callable:
TensorMap{T}(f, P::TensorMapSpace) where {T<:Number} = Tensor{T}(f(T, dims(P)), P)
TensorMap(f, P::TensorMapSpace) = Tensor(f(dims(P)), P)

# without spaces
TensorMap(data::AbstractMatrix{T}) where {T<:Real} = Tensor{T}(data, CartesianSpace(size(data,1)) ← CartesianSpace(size(data,2)))
TensorMap(data::AbstractMatrix{T}) where {T<:Complex} = Tensor{T}(data, ComplexSpace(size(data,1)) ← ComplexSpace(size(data,2)))

# Similar
#---------
Base.similar(t::TensorMap{T1, S}, ::Type{T2}, P::TensorMapSpace{S} = (dom(t)=>codom(t))) where {T1,T2,S} = Tensor(similar(t.data, T, dims(P)), P)
Base.similar(t::TensorMap{T1, S}, ::Type{T2}, P::TensorSpace{S}) where {T1,T2,S} = Tensor(similar(t.data, T, dims(P)), P)
Base.similar(t::TensorMap{T1, S}, ::Type{T2}, V::S) where {T1,T2,S} = similar(t, T2, TensorSpace(V))

Base.similar(t::TensorMap{T, S}, P::TensorMapSpace{S} = (dom(t)=>codom(t))) where {T,S} = similar(t, eltype(t), P)
Base.similar(t::TensorMap{T, S}, P::TensorSpace{S}) where {T,S} = simlar(t, eltype(t), P)
Base.similar(t::TensorMap{T, S}, V::S) where {T,S} = similar(t, eltype(t), V)

# Copy and fill tensors:
#------------------------
function Base.copy!(tdest::TensorMap, tsource::TensorMap)
    codom(tdest) == codom(tsource) && dom(tdest) == dom(tsource) || throw(SpaceError())
    copy!(tdest.data, tsource.data)
    return tdest
end
Base.fill!(tdest::TensorMap, value::Number)=fill!(tdest.data, value)

# Conversion and promotion:
#---------------------------
Base.promote_rule(::Type{<:TensorMap{T1,S}},::Type{<:TensorMap{T2,S}}) where {T1, T2, S} = TensorMap{promote_type(T1,T2),S}

Base.convert(::Type{TensorMap{T,S}}, t::Tensor{T,S}) where {T,S} = t
Base.convert(::Type{TensorMap{T1,S}}, t::Tensor{T2,S}) where {T1,T2,S} = copy!(similar(t, T1), t)

# TODO: Check whether we need anything of this old stuff
# Base.promote_rule{S,T1,T2,N1,N2}(::Type{Tensor{S,T1,N1}},::Type{Tensor{S,T2,N2}})=Tensor{S,promote_type(T1,T2)}
# Base.promote_rule{S,T1,T2}(::Type{Tensor{S,T1}},::Type{Tensor{S,T2}})=Tensor{S,promote_type(T1,T2)}
#
# Base.promote_rule{S,T1,T2,N}(::Type{AbstractTensor{S,ProductSpace,T1,N}},::Type{Tensor{S,T2,N}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2),N}
# Base.promote_rule{S,T1,T2,N1,N2}(::Type{AbstractTensor{S,ProductSpace,T1,N1}},::Type{Tensor{S,T2,N2}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2)}
# Base.promote_rule{S,T1,T2}(::Type{AbstractTensor{S,ProductSpace,T1}},::Type{Tensor{S,T2}})=AbstractTensor{S,ProductSpace,promote_type(T1,T2)}

# Base.convert{S,T,N}(::Type{Tensor{S,T,N}},t::Tensor{S,T,N})=t
# Base.convert{S,T1,T2,N}(::Type{Tensor{S,T1,N}},t::Tensor{S,T2,N})=copy!(similar(t,T1),t)
# Base.convert{S,T}(::Type{Tensor{S,T}},t::Tensor{S,T})=t
# Base.convert{S,T1,T2}(::Type{Tensor{S,T1}},t::Tensor{S,T2})=copy!(similar(t,T1),t)
#
# Base.float{S,T<:FloatingPoint}(t::Tensor{S,T})=t
# Base.float(t::Tensor)=tensor(float(t.data),space(t))
#
# Base.real{S,T<:Real}(t::Tensor{S,T})=t
# Base.real(t::Tensor)=tensor(real(t.data),space(t))
# Base.complex{S,T<:Complex}(t::Tensor{S,T})=t
# Base.complex(t::Tensor)=tensor(complex(t.data),space(t))
#
# for (f,T) in ((:float32,    Float32),
#               (:float64,    Float64),
#               (:complex64,  Complex64),
#               (:complex128, Complex128))
#     @eval (Base.$f){S}(t::Tensor{S,$T}) = t
#     @eval (Base.$f)(t::Tensor) = tensor(($f)(t.data),space(t))
# end
#
# Basic vector space methods:
# ---------------------------
function Base.scale!(t1::TensorMap, t2::TensorMap, a::Number)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceError())
    t1.data .= a .* t2.data
    return t1
end

function add!(t1::TensorMap, t2::TensorMap, a::Number)
    (codomain(t1)==codomain(t2) && domain(t1) == domain(t2)) || throw(SpaceError())
    t1.data .+= a .* t2.data
end

# Indexing
#----------
# # linear indexing using ProductBasisVector
Base.getindex(t::TensorMap, indices...)=getindex(t.data, indices...)
Base.setindex!(t::TensorMap, value, indices...)=setindex!(t.data, value, indices...)

# # Index methods
# #---------------
# @eval function insertind{S}(t::Tensor{S},ind::Int,V::S)
#     N=numind(t)
#     0<=ind<=N || throw(IndexError("index out of range"))
#     iscnumber(V) || throw(SpaceError("can only insert index with c-number index space"))
#     spacet=space(t)
#     newspace=spacet[1:ind] ⊗ V ⊗ spacet[ind+1:N]
#     return tensor(t.data,newspace)
# end
# @eval function deleteind(t::Tensor,ind::Int)
#     N=numind(t)
#     1<=ind<=N || throw(IndexError("index out of range"))
#     iscnumber(space(t,ind)) || throw(SpaceError("can only delete index with c-number index space"))
#     spacet=space(t)
#     newspace=spacet[1:ind-1] ⊗ spacet[ind+1:N]
#     return tensor(t.data,newspace)
# end
#
# for (S,TT) in ((CartesianSpace,CartesianTensor),(ComplexSpace,ComplexTensor))
#     @eval function fuseind(t::$TT,ind1::Int,ind2::Int,V::$S)
#         N=numind(t)
#         ind2==ind1+1 || throw(IndexError("only neighbouring indices can be fused"))
#         1<=ind1<=N-1 || throw(IndexError("index out of range"))
#         fuse(space(t,ind1),space(t,ind2),V) || throw(SpaceError("index spaces $(space(t,ind1)) and $(space(t,ind2)) cannot be fused to $V"))
#         spacet=space(t)
#         newspace=spacet[1:ind1-1]*V*spacet[ind2+1:N]
#         return tensor(t.data,newspace)
#     end
#     @eval function splitind(t::$TT,ind::Int,V1::$S,V2::$S)
#         1<=ind<=numind(t) || throw(IndexError("index out of range"))
#         fuse(V1,V2,space(t,ind)) || throw(SpaceError("index space $(space(t,ind)) cannot be split into $V1 and $V2"))
#         spacet=space(t)
#         newspace=spacet[1:ind-1]*V1*V2*spacet[ind+1:numind(t)]
#         return tensor(t.data,newspace)
#     end
# end
#

## NOTE: Do we want this?
# # tensors from concatenation
# function tensorcat{S}(catind, X::Tensor{S}...)
#     catind = collect(catind)
#     isempty(catind) && error("catind should not be empty")
#     # length(unique(catdims)) != length(catdims) && error("every dimension should appear only once")
#
#     nargs = length(X)
#     numindX = map(numind, X)
#
#     all(n->(n == numindX[1]), numindX) || throw(SpaceError("all tensors should have the same number of indices for concatenation"))
#
#     numindC = numindX[1]
#     ncatind = setdiff(1:numindC,catind)
#     spaceCvec = Array(S, numindC)
#     for n = 1:numindC
#         spaceCvec[n] = space(X[1],n)
#     end
#     for i = 2:nargs
#         for n in catind
#             spaceCvec[n] = directsum(spaceCvec[n], space(X[i],n))
#         end
#         for n in ncatind
#             spaceCvec[n] == space(X[i],n) || throw(SpaceError("space mismatch for index $n"))
#         end
#     end
#     spaceC = ⊗(spaceCvec...)
#     typeC = mapreduce(eltype, promote_type, X)
#     dataC = zeros(typeC, map(dim,spaceC))
#
#     offset = zeros(Int,numindC)
#     for i=1:nargs
#         currentdims=ntuple(numindC,n->dim(space(X[i],n)))
#         currentrange=[offset[n]+(1:currentdims[n]) for n=1:numindC]
#         dataC[currentrange...] = X[i].data
#         for n in catind
#             offset[n]+=currentdims[n]
#         end
#     end
#     return tensor(dataC,spaceC)
# end


text:
@safetestset "App Config" begin
  @safetestset "Custom server startup overwrites app config" begin
    using Genie

    @test Genie.config.server_port  == 8000
    @test Genie.config.server_host  == "127.0.0.1"
    @test Genie.config.websockets_port == 8000
    @test Genie.config.run_as_server   == false

    up(9000, "0.0.0.0"; open_browser = false)

    @test Genie.config.server_port  == 9000
    @test Genie.config.server_host  == "0.0.0.0"
    @test Genie.config.websockets_port == 9000
    @test Genie.config.run_as_server   == false

    down()
    sleep(1)

    up(9000, "0.0.0.0", ws_port = 9999; open_browser = false)

    @test Genie.config.server_port  == 9000
    @test Genie.config.server_host  == "0.0.0.0"
    @test Genie.config.websockets_port == 9999
    @test Genie.config.run_as_server   == false

    down()
    sleep(1)

    Genie.config.server_port  = 8000
    Genie.config.server_host  = "127.0.0.1"
    Genie.config.websockets_port = 8000
    Genie.config.run_as_server   = false

    down()
    sleep(1)
  end;
end;


text:
macro loop(name, body)
    return _loop(name, body)
end

macro loop(name, initialize, body, finalize=:())
    return _loop(name, initialize, body, finalize)
end

_loop(name::String, body::Expr) = _loop(name, :(), body)
function _loop(name, initialize, body, finalize=:(), target_speed=50)
    return quote
        function (window::Window, target_speed=$target_speed)
            loop = Loop($name, target_speed)

            task = @task try
                @debug "Running $($name) on thread: $(Threads.threadid())"
                delta_times = CircularBuffer{Float64}(ceil(Int, target_speed/4))
                ticks = 0
                t0 = time()
                elapsed_time = time() - t0
                push!(delta_times, 1/target_speed)
                $initialize
                while length(window.scene_stack) > 0
                    $body
                    expected_time = ticks/target_speed
                    push!(delta_times, time() - (t0 + elapsed_time))
                    elapsed_time += delta_times[end]
                    sleep(max(1/target_speed + expected_time - elapsed_time, 1e-3))
                    ticks += 1
                end
                $finalize
            catch e
                @error "$(loop.name)" exception=(e, catch_backtrace())
                rethrow()
            end

            loop.task = task
            return loop
        end
    end |> esc
end

# Default loops

const DEFAULT_EVENT_LOOP = @loop "event" (event_data = zeros(UInt8, 56)) begin
    while SDL.PollEvent(event_data) != 0
        event = parseevent(window, event_data)
        onevent!(window, event)
    end
end

const DEFAULT_UPDATE_LOOP = @loop "update" (t = elapsed_time; dt = 1/target_speed) begin
    if !loop.paused
        t = elapsed_time
    end

    for event in window.event_queue
        onevent!(window, event, t, dt)
    end
    empty!(window.event_queue)

    sort!(window.timer_queue)
    while length(window.timer_queue) > 0 && window.timer_queue[end].t1 <= t
        pop!(window.timer_queue).fn()
    end

    if !loop.paused
        before_update!(window, t, dt)
        update!(window, t, dt)
        after_update!(window, t, dt)
    end
end

const DEFAULT_RENDER_LOOP = @loop "render" () begin
    fps = length(delta_times)/sum(delta_times)
    render!(window, ticks, fps)
end

"""
    run!(window::Window, [loop])

Execute `window`.

"""
function run!(window::Window)
    loops = map(
        loop -> loop(window),
        (DEFAULT_EVENT_LOOP, DEFAULT_UPDATE_LOOP, DEFAULT_RENDER_LOOP)
    )
    for loop in loops
        window.loops[loop.name] = loop
    end
    Threads.@spawn schedule(loops[1])
    Threads.@spawn schedule(loops[2])
    schedule(loops[3])

    return window
end

function run!(window::Window, loops...)
    for loop in map(loop -> loop(window), loops)
        window.loops[loop.name] = loop
        schedule(loop)
    end
    return window
end

function register_loops!(window::Window, loops...)
    for loop in map(loop -> loop(window), loops)
        window.loops[loop.name] = loop
    end
    return window
end

function reload!(window)
    wait(window)
    empty!(window.tasks)
    return run!(window, DEFAULT_EVENT_LOOP, DEFAULT_UPDATE_LOOP, DEFAULT_RENDER_LOOP)
end

settimer!(window::Window, t::Real, fn::Function) = push!(window.timer_queue, Timeout(t, fn))

struct InterruptQuit end

"""
    quit!(window::Window, e::Event)

Gracefully quit `window`, propagating the event, `e`, which caused
`quit!` to be called to each scene in `window` to allow processing and
optionally interrupting the call by returning `InterruptQuit()`.
"""
function quit!(window::Window, e::Event)
    while length(window.scene_stack) > 0
        scene = pop!(window)
        if onevent!(scene, e) === InterruptQuit()
            push!(window, scene)
            return InterruptQuit()
        end
    end
    destroy!(window)
end

function onevent!(::AbstractObject, ::Event) end
function onevent!(::AbstractObject, ::Event, t, dt) end
for fn in [:before_update!, :update!, :after_update!]
    @eval function $fn(::AbstractObject, ::AbstractLayer, t, dt) end
end
function render!(::Window, ::AbstractObject, frame, fps) end
function render!(::AbstractLayer, ::AbstractObject, frame, fps) end

onevent!(window::Window, ::Event{:notsupported}) = window
function onevent!(window::Window, e::Event)
    if length(window.scene_stack) > 0
        onevent!(activescene(window), e)
    end
    push!(window.event_queue, e)
    return window
end
onevent!(window::Window, e::Event{:quit}) = quit!(window, e)

onevent!(window::Window, e::Event, t, dt) = length(window.scene_stack) > 0 ? onevent!(activescene(window), e, t, dt) : window

for fn in [:before_update!, :update!, :after_update!]
    @eval $fn(window::Window, t, dt) = length(window.scene_stack) > 0 ? $fn(activescene(window), window, t, dt) : nothing
end

render!(window::Window, frame, fps) = length(window.scene_stack) > 0 ? render!(window, activescene(window), frame, fps) : nothing

function onevent!(scene::AbstractScene, e::Event)
    for layer in scene.layers
        onevent!(layer, e)
    end
    return scene
end

function onevent!(scene::AbstractScene, e::Event, t, dt)
    for layer in scene.layers
        onevent!(layer, e, t, dt)
    end
    return scene
end

for fn in [:before_update!, :update!, :after_update!]
    @eval function $fn(scene::AbstractScene, ::Window, t, dt)
        for layer in scene.layers
            $fn(layer, scene, t, dt)
        end
        return scene
    end
end

function render!(window::Window, scene::AbstractScene, frame, fps)
    setcolor!(window, scene.color)
    clear!(window)
    for layer in scene.layers
        render!(window, layer, frame, fps)
    end
    present!(window)
end

function onevent!(layer::AbstractLayer, e::Event)
    for obj in layer.objects
        onevent!(obj, e)
    end
    return layer
end

function onevent!(layer::AbstractLayer, e::Event, t, dt)
    for obj in layer.objects
        onevent!(obj, e, t, dt)
    end
    return layer
end

function before_update!(layer::AbstractLayer, ::AbstractScene, t, dt)
    populate!(layer)
    for obj in layer.objects
        before_update!(obj, layer, t, dt)
    end
    return layer
end

function update!(layer::AbstractLayer, ::AbstractScene, t, dt)
    for obj in layer.objects
        update!(obj, layer, t, dt)
    end
    return layer
end

function after_update!(layer::AbstractLayer, ::AbstractScene, t, dt)
    for obj in layer.objects
        after_update!(obj, layer, t, dt)
    end
    cleanup!(layer)
    return layer
end

function render!(window::Window, layer::AbstractLayer, frame, fps)
    if visible(layer)
        # sort!(layer.objects, by=layer.sortby)
        for obj in layer.objects
            render!(layer, obj, frame, fps)
        end

        for r in layer.render_tasks
            render!(window, layer, r)
        end
        empty!(layer.render_tasks)
    end
    return window
end


text:
function sphere_vol(r)
   return 4/3*pi*r^3
end

sphere_vol(2.5)


text:
module TheoCrystalPedestal

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Theo Crystal Pedestal" => Ahorn.EntityPlacement(
        Maple.TheoCrystalPedestal
    )
)

sprite = "characters/theoCrystal/pedestal.png"

function Ahorn.selection(entity::Maple.TheoCrystalPedestal)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y, jx=0.5, jy=1.0)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TheoCrystalPedestal, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0, jx=0.5, jy=1.0)

end


text:
module DeIdentification

export deidentify, ProjectConfig, DeIdDicts, build_config, build_config_from_csv

import YAML
import Glob
import JSON
import Tables
import CSV
import Dates
import SHA: bytes2hex, sha256
import Random: shuffle, randstring, seed!, make_seed
import Memento
import DataStructures: OrderedDict
import REPL
import Parsers
using REPL.TerminalMenus
using DelimitedFiles

include("config_builder.jl")
include("de_identify.jl")
include("exporting.jl")
include("utils.jl")

"""
    deid_file!(dicts, file_config, project_config, logger)

Reads raw file and deidentifies per file configuration and project configurationg.
Writes the deidentified data to a CSV file and updates the global dictionaries
tracking identifier mappings.
"""
function deid_file!(dicts::DeIdDicts, fc::FileConfig, pc::ProjectConfig, logger)
    # Initiate new file
    infile = try
        CSV.File(fc.filename, dateformat = fc.dateformat)
    catch ArgumentError
        CSV.File(fc.filename)
    end

    dicts = DeIdDicts(dicts, fc.dateformat)

    outfile = joinpath(pc.outdir, "deid_" * fc.name * "_" * getcurrentdate() * ".csv")

    ncol = length(infile.names)
    lastcol = infile.names[end]

    new_names = Vector{Symbol}()
    new_types = Vector{Type}()
    pk = false
    pcol = pc.primary_id

    Memento.info(logger, "$(Dates.now()) Renaming file columns")
    for i = 1:ncol
        n = infile.names[i]
        t = infile.types[i]

        if haskey(fc.rename_cols, n)
            push!(new_names, fc.rename_cols[n])
            push!(new_types, t)

            if fc.rename_cols[n] == pc.primary_id
                pk = true
                pcol = n
            end
        elseif get(fc.colmap, n, Missing) == Drop
            continue
        else
            push!(new_names, n)
            push!(new_types, t)
        end

        if n == pc.primary_id
            pk = true
        end
    end

    Memento.info(logger, "$(Dates.now()) Checking for primary column")
    @assert pk==true "Primary ID must be present in file"

    open(outfile, "w") do io
        # write header to file
        header = [string(n) for n in new_names]
        writedlm(io, reshape(header, 1, length(header)), ',')

        # Process each row
        for (i, row) in Iterators.enumerate(infile)
            try
                val = getoutput(dicts, Hash, getproperty(row, pcol), 0)
                pid = setrid(val, dicts)
                columns = Vector{String}()

                for col in infile.names
                    colname = get(fc.rename_cols, col, col)

                    action = get(fc.colmap, colname, Missing) ::Type

                    if action == Drop
                        continue
                    end

                    VAL = getproperty(row, col)

                    # apply pre-processing transform
                    if haskey(fc.preprocess, colname) && !ismissing(VAL)
                        transform = fc.preprocess[colname]
                        transform = replace(transform, "VAL" => "\"$VAL\"")
                        expr = Meta.parse(transform)
                        VAL = Core.eval(@__MODULE__, expr)
                    end

                    VAL = getoutput(dicts, action, VAL, pid)

                    if col == pcol
                        VAL = pid
                    end

                    # apply post-processing transform
                    if haskey(fc.postprocess, colname) && !ismissing(VAL)
                        transform = fc.postprocess[colname]
                        transform = replace(transform, "VAL" => "\"$VAL\"")
                        expr = Meta.parse(transform)
                        VAL = Core.eval(@__MODULE__, expr)
                    end

                    if eltype(VAL) <: String
                        VAL = replace(VAL, "\"" => "\\\"")
                    end

                    if VAL !== nothing && !ismissing(VAL)
                        push!(columns, string(VAL))
                    else
                        push!(columns, "")
                    end
                end

                writedlm(io, reshape(columns, 1, length(columns)), ',')
            catch e
                Memento.error(logger, "$(Dates.now()) Error occurred while processing row $i")
                rethrow(e)
            end
        end

    end

    return nothing
end

"""
    deidentify(cfg::ProjectConfig)
This is the constructor for the `DeIdentified` struct. We use this type to store
arrays of `DeIdDataFrame` variables, while also keeping a common `salt_dict` and
`dateshift_dict` between `DeIdDataFrame`s. The `salt_dict` allows us to track
what salt was used on what cleartext. This is only necessary in the case of doing
re-identification. The `id_dict` argument is a dictionary containing the hash
digest of the original primary ID to our new research IDs.
"""
function deidentify(cfg::ProjectConfig)
    num_files = length(cfg.file_configs)
    dicts = DeIdDicts(cfg.maxdays, cfg.shiftyears, cfg.dateformat)

    if !isdir(cfg.outdir)
        # mkpath also creates any intermediate paths
        mkpath(cfg.outdir)
    end

    logdir = dirname(cfg.logfile)
    if !isdir(logdir)
        mkpath(logdir)
    end

    # Set up our top-level logger
    logger = Memento.getlogger("deidentify")
    logfile_roller = Memento.FileRoller(cfg.logfile)
    push!(logger, Memento.DefaultHandler(logfile_roller))

    Memento.info(logger, "$(Dates.now()) Logging session for project $(cfg.name)")

    Memento.info(logger, "$(Dates.now()) Setting seed for project $(cfg.name)")
    seed!(cfg.seed)

    for (i, fc) in enumerate(cfg.file_configs)
        Memento.info(logger, "$(Dates.now()) ====================== Processing $(fc.name) ======================")

        Memento.info(logger, "$(Dates.now()) Reading data from $(fc.filename)")
        deid_file!(dicts, fc, cfg, logger)

    end

    write_dicts(dicts, logger, cfg.outdir)

    return dicts
end

"""
    deidentify(config_path)

Run entire pipeline: Processes configuration YAML file, de-identifies the data,
and writes the data to disk.  Returns the dictionaries containing the mappings.
"""
function deidentify(cfg_file::String)
    proj_config = ProjectConfig(cfg_file)
    return deidentify(proj_config)
end

end # module


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions to align strings.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export align_string, align_string_per_line, get_padding_for_string_alignment

"""
    align_string(str::AbstractString, field_width::Int, alignment::Symbol; fill::Bool = false )

Align the string `str` in the field with width `field_width` using `alignment`.

`alignment` can be:

- `:l`: Align the string to the left;
- `:c`: Align the string in the center;
- `:r`: Align the string in the right.

!!! note
    If the printable width of `str` is higher than `field_width`, then nothing
    will be changed.

!!! note
    This function treats `\n` as a normal characters. If one wants to align
    every line, then use the function [`align_string_per_line`](@ref).

# Keyword

- `fill::Bool`: If `true`, the string will be filled with spaces to the right so
    that the resulting string has printable width `field_size` if the initial
    string printable width is lower than it.
- `printable_string_width::Int`: Provide the printable string width to reduce
    the computational burden. If this parameters is lower than 0, the printable
    width is compute internally. (**Default** = -1)
"""
function align_string(
    str::AbstractString,
    field_width::Int,
    alignment::Symbol;
    fill::Bool = false,
    printable_string_width::Int = -1
)
    padding = get_padding_for_string_alignment(
        str,
        field_width,
        alignment;
        fill,
        printable_string_width
    )

    if isnothing(padding)
        return str
    else
        lpad, rpad = padding
        return " "^lpad * str * " "^rpad
    end
end

"""
    align_string_per_line(str::AbstractString, field_width::Int, alignment::Symbol; fill::Bool = false)

Align each line of the string `str` in the field with width `field_width` using
`alignment`.

`alignment` can be:

- `:l`: Align the string to the left;
- `:c`: Align the string in the center;
- `:r`: Align the string in the right.

!!! note
    If the printable width of `str` is higher than `field_width`, then nothing
    will be changed.

# Keyword

- `fill::Bool`: If `true`, the each line will be filled with spaces to the right
    so that the resulting string has printable width `field_size` if the initial
    string printable width is lower than it.
"""
function align_string_per_line(
    str::AbstractString,
    field_width::Int,
    alignment::Symbol;
    fill::Bool = false
)
    if field_width ≤ 0
        return str
    end

    # Split the lines.
    lines = split(str, '\n')
    num_lines = length(lines)

    # Align each one of them.
    buf = IOBuffer(sizehint = floor(Int, sizeof(str) + num_lines * div(field_width, 2)))

    @inbounds for i in 1:num_lines
        write(buf, align_string(lines[i], field_width, alignment; fill))
        i != num_lines && write(buf, '\n')
    end

    return String(take!(buf))
end

"""
    get_padding_for_string_alignment(str::AbstractString, field_width::Int, alignment::Symbol; kwargs...)

Return the left and right padding required to align the string `str` in a field
with width `field_width` using the `alignment`.

`alignment` can be:

- `:l`: Align the string to the left;
- `:c`: Align the string in the center;
- `:r`: Align the string in the right.

!!! note
    If the printable width of `str` is higher than `field_width`, then nothing
    is returned.

!!! note
    This function treats `\n` as a normal characters.

# Keyword

- `fill::Bool`: If `true`, the string will be filled with spaces to the right so
    that the resulting string has printable width `field_size` if the initial
    string printable width is lower than it.
- `printable_string_width::Int`: Provide the printable string width to reduce
    the computational burden. If this parameters is lower than 0, the printable
    width is compute internally. (**Default** = -1)
"""
function get_padding_for_string_alignment(
    str::AbstractString,
    field_width::Int,
    alignment::Symbol;
    fill::Bool = false,
    printable_string_width::Int = -1
)

    str_width = printable_string_width < 0 ?
        printable_textwidth(str) :
        printable_string_width

    if field_width ≤ str_width
        return nothing
    end

    # Compute the padding given the alignment type.
    if alignment == :l
        # In this case, we only need to modify the string if `fill` is `true`.
        if fill
            rpad = field_width - str_width
            return 0, rpad
        end

    elseif alignment == :c
        lpad = div(field_width - str_width, 2)
        rpad = fill ? (field_width - str_width - lpad) : 0
        return lpad, rpad

    elseif alignment == :r
        lpad = field_width - str_width
        return lpad, 0
    end

    return nothing
end