text:
module Topologies
import ..BrickMesh
import MPI
using DocStringExtensions

export AbstractTopology,
    BrickTopology,
    StackedBrickTopology,
    CubedShellTopology,
    StackedCubedSphereTopology,
    isstacked,
    cubedshellwarp,
    cubedshellunwarp

export grid1d, SingleExponentialStretching, InteriorStretching

"""
    AbstractTopology{dim}

Represents the connectivity of individual elements, with local dimension `dim`.
"""
abstract type AbstractTopology{dim} end

"""
    BoxElementTopology{dim, T} <: AbstractTopology{dim}

The local topology of a larger MPI-distributed topology, represented by
`dim`-dimensional box elements.

This contains the necessary information for the connectivity elements of the
elements on the local process, along with "ghost" elements from neighbouring
processes.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct BoxElementTopology{dim, T} <: AbstractTopology{dim}

    """
    MPI communicator for communicating with neighbouring processes.
    """
    mpicomm::MPI.Comm

    """
    Range of element indices
    """
    elems::UnitRange{Int64}

    """
    Range of real (aka nonghost) element indices
    """
    realelems::UnitRange{Int64}

    """
    Range of ghost element indices
    """
    ghostelems::UnitRange{Int64}

    """
    Ghost element to face is received; `ghostfaces[f,ge] == true` if face `f` of
    ghost element `ge` is received.
    """
    ghostfaces::BitArray{2}

    """
    Array of send element indices
    """
    sendelems::Array{Int64, 1}

    """
    Send element to face is sent; `sendfaces[f,se] == true` if face `f` of send
    element `se` is sent.
    """
    sendfaces::BitArray{2}

    """
    Array of real elements that do not have a ghost element as a neighbor.
    """
    interiorelems::Array{Int64, 1}

    """
    Array of real elements that have at least on ghost element as a neighbor.

    Note that this is different from `sendelems` because `sendelems` duplicates
    elements that need to be sent to multiple neighboring processes.
    """
    exteriorelems::Array{Int64, 1}

    """
    Element to vertex coordinates; `elemtocoord[d,i,e]` is the `d`th coordinate
    of corner `i` of element `e`

    !!! note
        currently coordinates always are of size 3 for `(x1, x2, x3)`
    """
    elemtocoord::Array{T, 3}

    """
    Element to neighboring element; `elemtoelem[f,e]` is the number of the
    element neighboring element `e` across face `f`.  If there is no neighboring
    element then `elemtoelem[f,e] == e`.
    """
    elemtoelem::Array{Int64, 2}

    """
    Element to neighboring element face; `elemtoface[f,e]` is the face number of
    the element neighboring element `e` across face `f`.  If there is no
    neighboring element then `elemtoface[f,e] == f`."
    """
    elemtoface::Array{Int64, 2}

    """
    element to neighboring element order; `elemtoordr[f,e]` is the ordering
    number of the element neighboring element `e` across face `f`.  If there is
    no neighboring element then `elemtoordr[f,e] == 1`.
    """
    elemtoordr::Array{Int64, 2}

    """
    Element to boundary number; `elemtobndy[f,e]` is the boundary number of face
    `f` of element `e`.  If there is a neighboring element then `elemtobndy[f,e]
    == 0`.
    """
    elemtobndy::Array{Int64, 2}

    """
    List of the MPI ranks for the neighboring processes
    """
    nabrtorank::Array{Int64, 1}

    """
    Range in ghost elements to receive for each neighbor
    """
    nabrtorecv::Array{UnitRange{Int64}, 1}

    """
    Range in `sendelems` to send for each neighbor
    """
    nabrtosend::Array{UnitRange{Int64}, 1}

    """
    original order in partitioning
    """
    origsendorder::Array{Int64, 1}

    """
    boolean for whether or not this topology has a boundary
    """
    hasboundary::Bool

    function BoxElementTopology{dim, T}(
        mpicomm,
        elems,
        realelems,
        ghostelems,
        ghostfaces,
        sendelems,
        sendfaces,
        elemtocoord,
        elemtoelem,
        elemtoface,
        elemtoordr,
        elemtobndy,
        nabrtorank,
        nabrtorecv,
        nabrtosend,
        origsendorder,
        hasboundary,
    ) where {dim, T}

        exteriorelems = sort(unique(sendelems))
        interiorelems = sort(setdiff(realelems, exteriorelems))

        return new{dim, T}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            interiorelems,
            exteriorelems,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            origsendorder,
            hasboundary,
        )
    end
end

"""
    hasboundary(topology::AbstractTopology)

query function to check whether a topology has a boundary (i.e., not fully
periodic)
"""
hasboundary(topology::AbstractTopology) = topology.hasboundary

if VERSION >= v"1.2-"
    isstacked(::T) where {T <: AbstractTopology} = hasfield(T, :stacksize)
else
    isstacked(::T) where {T <: AbstractTopology} =
        Base.fieldindex(T, :stacksize, false) > 0
end

"""
    BrickTopology{dim, T} <: AbstractTopology{dim}

A simple grid-based topology. This is a convenience wrapper around
[`BoxElementTopology`](@ref).
"""
struct BrickTopology{dim, T} <: AbstractTopology{dim}
    topology::BoxElementTopology{dim, T}
end
Base.getproperty(a::BrickTopology, p::Symbol) =
    getproperty(getfield(a, :topology), p)

"""
    CubedShellTopology{T} <: AbstractTopology{2}

A cube-shell topology. This is a convenience wrapper around
[`BoxElementTopology`](@ref).
"""
struct CubedShellTopology{T} <: AbstractTopology{2}
    topology::BoxElementTopology{2, T}
end
Base.getproperty(a::CubedShellTopology, p::Symbol) =
    getproperty(getfield(a, :topology), p)

"""
    StackedBrickTopology{dim, T} <: AbstractTopology{dim}

A simple grid-based topology, where all elements on the trailing dimension are
stacked to be contiguous. This is a convenience wrapper around
[`BoxElementTopology`](@ref).
"""
struct StackedBrickTopology{dim, T} <: AbstractTopology{dim}
    topology::BoxElementTopology{dim, T}
    stacksize::Int64
end
function Base.getproperty(a::StackedBrickTopology, p::Symbol)
    return p == :stacksize ? getfield(a, p) :
           getproperty(getfield(a, :topology), p)
end

"""
    StackedCubedSphereTopology{3, T} <: AbstractTopology{3}

A cube-sphere topology. All elements on the same "vertical" dimension are
stacked to be contiguous. This is a convenience wrapper around
[`BoxElementTopology`](@ref).
"""
struct StackedCubedSphereTopology{T} <: AbstractTopology{3}
    topology::BoxElementTopology{3, T}
    stacksize::Int64
end
function Base.getproperty(a::StackedCubedSphereTopology, p::Symbol)
    return p == :stacksize ? getfield(a, p) :
           getproperty(getfield(a, :topology), p)
end

""" A wrapper for the BrickTopology """
BrickTopology(mpicomm, Nelems::NTuple{N, Integer}; kw...) where {N} =
    BrickTopology(mpicomm, map(Ne -> 0:Ne, Nelems); kw...)

"""
    BrickTopology{dim, T}(mpicomm, elemrange; boundary, periodicity)

Generate a brick mesh topology with coordinates given by the tuple `elemrange`
and the periodic dimensions given by the `periodicity` tuple.

The elements of the brick are partitioned equally across the MPI ranks based
on a space-filling curve.

By default boundary faces will be marked with a one and other faces with a
zero.  Specific boundary numbers can also be passed for each face of the brick
in `boundary`.  This will mark the nonperiodic brick faces with the given
boundary number.

# Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the
\$x2\$-direction with
```jldoctest brickmesh

using CLIMA.Topologies
using MPI
MPI.Init()
topology = BrickTopology(MPI.COMM_SELF, (2:5,4:6);
                         periodicity=(false,true),
                         boundary=((1,2),(3,4)))
```
This returns the mesh structure for

             x2

              ^
              |
             6-  +-----+-----+-----+
              |  |     |     |     |
              |  |  3  |  4  |  5  |
              |  |     |     |     |
             5-  +-----+-----+-----+
              |  |     |     |     |
              |  |  1  |  2  |  6  |
              |  |     |     |     |
             4-  +-----+-----+-----+
              |
              +--|-----|-----|-----|--> x1
                 2     3     4     5

For example, the (dimension by number of corners by number of elements) array
`elemtocoord` gives the coordinates of the corners of each element.
```jldoctest brickmesh
julia> topology.elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
 2  3  2  3
 4  4  5  5

[:, :, 2] =
 3  4  3  4
 4  4  5  5

[:, :, 3] =
 2  3  2  3
 5  5  6  6

[:, :, 4] =
 3  4  3  4
 5  5  6  6

[:, :, 5] =
 4  5  4  5
 5  5  6  6

[:, :, 6] =
 4  5  4  5
 4  4  5  5
```
Note that the corners are listed in Cartesian order.

The (number of faces by number of elements) array `elemtobndy` gives the
boundary number for each face of each element.  A zero will be given for
connected faces.
```jldoctest brickmesh
julia> topology.elemtobndy
4×6 Array{Int64,2}:
 1  0  1  0  0  0
 0  0  0  0  2  2
 0  0  0  0  0  0
 0  0  0  0  0  0
```
Note that the faces are listed in Cartesian order.

"""
function BrickTopology(
    mpicomm,
    elemrange;
    boundary = ntuple(j -> (1, 1), length(elemrange)),
    periodicity = ntuple(j -> false, length(elemrange)),
    connectivity = :face,
    ghostsize = 1,
)

    if boundary isa Matrix
        boundary = tuple(mapslices(x -> tuple(x...), boundary, dims = 1)...)
    end

    # We cannot handle anything else right now...
    @assert connectivity == :face
    @assert ghostsize == 1

    mpirank = MPI.Comm_rank(mpicomm)
    mpisize = MPI.Comm_size(mpicomm)
    topology = BrickMesh.brickmesh(
        elemrange,
        periodicity,
        part = mpirank + 1,
        numparts = mpisize,
        boundary = boundary,
    )
    topology = BrickMesh.partition(mpicomm, topology...)
    origsendorder = topology[5]
    topology = BrickMesh.connectmesh(mpicomm, topology[1:4]...)

    dim = length(elemrange)
    T = eltype(topology.elemtocoord)
    return BrickTopology{dim, T}(BoxElementTopology{dim, T}(
        mpicomm,
        topology.elems,
        topology.realelems,
        topology.ghostelems,
        topology.ghostfaces,
        topology.sendelems,
        topology.sendfaces,
        topology.elemtocoord,
        topology.elemtoelem,
        topology.elemtoface,
        topology.elemtoordr,
        topology.elemtobndy,
        topology.nabrtorank,
        topology.nabrtorecv,
        topology.nabrtosend,
        origsendorder,
        !minimum(periodicity),
    ))
end

""" A wrapper for the StackedBrickTopology """
StackedBrickTopology(mpicomm, Nelems::NTuple{N, Integer}; kw...) where {N} =
    StackedBrickTopology(mpicomm, map(Ne -> 0:Ne, Nelems); kw...)

"""
    StackedBrickTopology{dim, T}(mpicomm, elemrange; boundary, periodicity)

Generate a stacked brick mesh topology with coordinates given by the tuple
`elemrange` and the periodic dimensions given by the `periodicity` tuple.

The elements are stacked such that the elements associated with range
`elemrange[dim]` are contiguous in the element ordering.

The elements of the brick are partitioned equally across the MPI ranks based
on a space-filling curve.  Further, stacks are not split at MPI boundaries.

By default boundary faces will be marked with a one and other faces with a
zero.  Specific boundary numbers can also be passed for each face of the brick
in `boundary`.  This will mark the nonperiodic brick faces with the given
boundary number.

# Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the
\$x2\$-direction with
```jldoctest brickmesh

using CLIMA.Topologies
using MPI
MPI.Init()
topology = StackedBrickTopology(MPI.COMM_SELF, (2:5,4:6);
                                periodicity=(false,true),
                                boundary=((1,2),(3,4)))
```
This returns the mesh structure stacked in the \$x2\$-direction for

             x2

              ^
              |
             6-  +-----+-----+-----+
              |  |     |     |     |
              |  |  2  |  4  |  6  |
              |  |     |     |     |
             5-  +-----+-----+-----+
              |  |     |     |     |
              |  |  1  |  3  |  5  |
              |  |     |     |     |
             4-  +-----+-----+-----+
              |
              +--|-----|-----|-----|--> x1
                 2     3     4     5

For example, the (dimension by number of corners by number of elements) array
`elemtocoord` gives the coordinates of the corners of each element.
```jldoctest brickmesh
julia> topology.elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
 2  3  2  3
 4  4  5  5

[:, :, 2] =
 2  3  2  3
 5  5  6  6

[:, :, 3] =
 3  4  3  4
 4  4  5  5

[:, :, 4] =
 3  4  3  4
 5  5  6  6

[:, :, 5] =
 4  5  4  5
 4  4  5  5

[:, :, 6] =
 4  5  4  5
 5  5  6  6
```
Note that the corners are listed in Cartesian order.

The (number of faces by number of elements) array `elemtobndy` gives the
boundary number for each face of each element.  A zero will be given for
connected faces.
```jldoctest brickmesh
julia> topology.elemtobndy
4×6 Array{Int64,2}:
 1  0  1  0  0  0
 0  0  0  0  2  2
 0  0  0  0  0  0
 0  0  0  0  0  0
```
Note that the faces are listed in Cartesian order.
"""
function StackedBrickTopology(
    mpicomm,
    elemrange;
    boundary = ntuple(j -> (1, 1), length(elemrange)),
    periodicity = ntuple(j -> false, length(elemrange)),
    connectivity = :face,
    ghostsize = 1,
)

    if boundary isa Matrix
        boundary = tuple(mapslices(x -> tuple(x...), boundary, dims = 1)...)
    end

    dim = length(elemrange)

    dim <= 1 && error("Stacked brick topology works for 2D and 3D")

    # Build the base topology
    basetopo = BrickTopology(
        mpicomm,
        elemrange[1:(dim - 1)];
        boundary = boundary[1:(dim - 1)],
        periodicity = periodicity[1:(dim - 1)],
        connectivity = connectivity,
        ghostsize = ghostsize,
    )

    # Use the base topology to build the stacked topology
    stack = elemrange[dim]
    stacksize = length(stack) - 1

    nvert = 2^dim
    nface = 2dim

    nreal = length(basetopo.realelems) * stacksize
    nghost = length(basetopo.ghostelems) * stacksize

    elems = 1:(nreal + nghost)
    realelems = 1:nreal
    ghostelems = nreal .+ (1:nghost)

    sendelems =
        similar(basetopo.sendelems, length(basetopo.sendelems) * stacksize)
    for i in 1:length(basetopo.sendelems), j in 1:stacksize
        sendelems[stacksize * (i - 1) + j] =
            stacksize * (basetopo.sendelems[i] - 1) + j
    end

    ghostfaces = similar(basetopo.ghostfaces, nface, length(ghostelems))
    ghostfaces .= false

    for i in 1:length(basetopo.ghostelems), j in 1:stacksize
        e = stacksize * (i - 1) + j
        for f in 1:(2 * (dim - 1))
            ghostfaces[f, e] = basetopo.ghostfaces[f, i]
        end
    end

    sendfaces = similar(basetopo.sendfaces, nface, length(sendelems))
    sendfaces .= false

    for i in 1:length(basetopo.sendelems), j in 1:stacksize
        e = stacksize * (i - 1) + j
        for f in 1:(2 * (dim - 1))
            sendfaces[f, e] = basetopo.sendfaces[f, i]
        end
    end

    elemtocoord = similar(basetopo.elemtocoord, dim, nvert, length(elems))

    for i in 1:length(basetopo.elems), j in 1:stacksize
        e = stacksize * (i - 1) + j

        for v in 1:(2^(dim - 1))
            for d in 1:(dim - 1)
                elemtocoord[d, v, e] = basetopo.elemtocoord[d, v, i]
                elemtocoord[d, 2^(dim - 1) + v, e] =
                    basetopo.elemtocoord[d, v, i]
            end

            elemtocoord[dim, v, e] = stack[j]
            elemtocoord[dim, 2^(dim - 1) + v, e] = stack[j + 1]
        end
    end

    elemtoelem = similar(basetopo.elemtoelem, nface, length(elems))
    elemtoface = similar(basetopo.elemtoface, nface, length(elems))
    elemtoordr = similar(basetopo.elemtoordr, nface, length(elems))
    elemtobndy = similar(basetopo.elemtobndy, nface, length(elems))

    for e in 1:(length(basetopo.elems) * stacksize), f in 1:nface
        elemtoelem[f, e] = e
        elemtoface[f, e] = f
        elemtoordr[f, e] = 1
        elemtobndy[f, e] = 0
    end

    for i in 1:length(basetopo.realelems), j in 1:stacksize
        e1 = stacksize * (i - 1) + j

        for f in 1:(2 * (dim - 1))
            e2 = stacksize * (basetopo.elemtoelem[f, i] - 1) + j

            elemtoelem[f, e1] = e2
            elemtoface[f, e1] = basetopo.elemtoface[f, i]

            # We assume a simple orientation right now
            @assert basetopo.elemtoordr[f, i] == 1
            elemtoordr[f, e1] = basetopo.elemtoordr[f, i]
        end

        et = stacksize * (i - 1) + j + 1
        eb = stacksize * (i - 1) + j - 1
        ft = 2 * (dim - 1) + 1
        fb = 2 * (dim - 1) + 2
        ot = 1
        ob = 1

        if j == stacksize
            et = periodicity[dim] ? stacksize * (i - 1) + 1 : e1
            ft = periodicity[dim] ? ft : 2 * (dim - 1) + 2
        end
        if j == 1
            eb = periodicity[dim] ? stacksize * (i - 1) + stacksize : e1
            fb = periodicity[dim] ? fb : 2 * (dim - 1) + 1
        end

        elemtoelem[2 * (dim - 1) + 1, e1] = eb
        elemtoelem[2 * (dim - 1) + 2, e1] = et
        elemtoface[2 * (dim - 1) + 1, e1] = fb
        elemtoface[2 * (dim - 1) + 2, e1] = ft
        elemtoordr[2 * (dim - 1) + 1, e1] = ob
        elemtoordr[2 * (dim - 1) + 2, e1] = ot
    end

    for i in 1:length(basetopo.elems), j in 1:stacksize
        e1 = stacksize * (i - 1) + j

        for f in 1:(2 * (dim - 1))
            elemtobndy[f, e1] = basetopo.elemtobndy[f, i]
        end

        bt = bb = 0

        if j == stacksize
            bt = periodicity[dim] ? bt : boundary[dim][2]
        end
        if j == 1
            bb = periodicity[dim] ? bb : boundary[dim][1]
        end

        elemtobndy[2 * (dim - 1) + 1, e1] = bb
        elemtobndy[2 * (dim - 1) + 2, e1] = bt
    end

    nabrtorank = basetopo.nabrtorank
    nabrtorecv = UnitRange{Int}[
        UnitRange(
            stacksize * (first(basetopo.nabrtorecv[n]) - 1) + 1,
            stacksize * last(basetopo.nabrtorecv[n]),
        ) for n in 1:length(nabrtorank)
    ]
    nabrtosend = UnitRange{Int}[
        UnitRange(
            stacksize * (first(basetopo.nabrtosend[n]) - 1) + 1,
            stacksize * last(basetopo.nabrtosend[n]),
        ) for n in 1:length(nabrtorank)
    ]

    T = eltype(basetopo.elemtocoord)

    StackedBrickTopology{dim, T}(
        BoxElementTopology{dim, T}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            basetopo.origsendorder,
            !minimum(periodicity),
        ),
        stacksize,
    )
end

"""
    CubedShellTopology(mpicomm, Nelem, T) <: AbstractTopology{dim}

Generate a cubed shell mesh with the number of elements along each dimension of
the cubes being `Nelem`. This topology actual creates a cube mesh, and the
warping should be done after the grid is created using the `cubedshellwarp`
function. The coordinates of the points will be of type `T`.

The elements of the shell are partitioned equally across the MPI ranks based
on a space-filling curve.

Note that this topology is logically 2-D but embedded in a 3-D space

# Examples

We can build a cubed shell mesh with 10 elements on each cube, total elements is
`10 * 10 * 6 = 600`, with
```jldoctest brickmesh
using CLIMA.Topologies
using MPI
MPI.Init()
topology = CubedShellTopology(MPI.COMM_SELF, 10, Float64)

# Typically the warping would be done after the grid is created, but the cell
# corners could be warped with...

# Shell radius = 1
x1, x2, x3 = ntuple(j->topology.elemtocoord[j, :, :], 3)
for n = 1:length(x1)
   x1[n], x2[n], x3[n] = Topologies.cubedshellwarp(x1[n], x2[n], x3[n])
end

# Shell radius = 10
x1, x2, x3 = ntuple(j->topology.elemtocoord[j, :, :], 3)
for n = 1:length(x1)
  x1[n], x2[n], x3[n] = Topologies.cubedshellwarp(x1[n], x2[n], x3[n], 10)
end
```
"""
function CubedShellTopology(
    mpicomm,
    Neside,
    T;
    connectivity = :face,
    ghostsize = 1,
)

    # We cannot handle anything else right now...
    @assert connectivity == :face
    @assert ghostsize == 1

    mpirank = MPI.Comm_rank(mpicomm)
    mpisize = MPI.Comm_size(mpicomm)

    topology = cubedshellmesh(Neside, part = mpirank + 1, numparts = mpisize)

    topology = BrickMesh.partition(mpicomm, topology...)
    origsendorder = topology[5]
    dim, nvert = 3, 4
    elemtovert = topology[1]
    nelem = size(elemtovert, 2)
    elemtocoord = Array{T}(undef, dim, nvert, nelem)
    ind2vert = CartesianIndices((Neside + 1, Neside + 1, Neside + 1))
    for e in 1:nelem
        for n in 1:nvert
            v = elemtovert[n, e]
            i, j, k = Tuple(ind2vert[v])
            elemtocoord[:, n, e] =
                (2 * [i - 1, j - 1, k - 1] .- Neside) / Neside
        end
    end

    topology = BrickMesh.connectmesh(
        mpicomm,
        topology[1],
        elemtocoord,
        topology[3],
        topology[4];
        dim = 2,
    )

    CubedShellTopology{T}(BoxElementTopology{2, T}(
        mpicomm,
        topology.elems,
        topology.realelems,
        topology.ghostelems,
        topology.ghostfaces,
        topology.sendelems,
        topology.sendfaces,
        topology.elemtocoord,
        topology.elemtoelem,
        topology.elemtoface,
        topology.elemtoordr,
        topology.elemtobndy,
        topology.nabrtorank,
        topology.nabrtorecv,
        topology.nabrtosend,
        origsendorder,
        false,
    ))
end

"""
    cubedshellmesh(T, Ne; part=1, numparts=1)

Generate a cubed mesh with each of the "cubes" has an `Ne X Ne` grid of
elements.

The mesh can optionally be partitioned into `numparts` and this returns
partition `part`.  This is a simple Cartesian partition and further partitioning
(e.g, based on a space-filling curve) should be done before the mesh is used for
computation.

This mesh returns the cubed spehere in a flatten fashion for the vertex values,
and a remapping is needed to embed the mesh in a 3-D space.

The mesh structures for the cubes is as follows:

```
x2
   ^
   |
4Ne-           +-------+
   |           |       |
   |           |   6   |
   |           |       |
3Ne-           +-------+
   |           |       |
   |           |   5   |
   |           |       |
2Ne-           +-------+
   |           |       |
   |           |   4   |
   |           |       |
 Ne-   +-------+-------+-------+
   |   |       |       |       |
   |   |   1   |   2   |   3   |
   |   |       |       |       |
  0-   +-------+-------+-------+
   |
   +---|-------|-------|------|-> x1
       0      Ne      2Ne    3Ne
```

"""
function cubedshellmesh(Ne; part = 1, numparts = 1)
    dim = 2
    @assert 1 <= part <= numparts

    globalnelems = 6 * Ne^2

    # How many vertices and faces per element
    nvert = 2^dim # 4
    nface = 2dim  # 4

    # linearly partition to figure out which elements we own
    elemlocal = BrickMesh.linearpartition(prod(globalnelems), part, numparts)

    # elemen to vertex maps which we own
    elemtovert = Array{Int}(undef, nvert, length(elemlocal))
    elemtocoord = Array{Int}(undef, dim, nvert, length(elemlocal))

    nelemcube = Ne^dim # Ne^2

    etoijb = CartesianIndices((Ne, Ne, 6))
    bx = [0 Ne 2Ne Ne Ne Ne]
    by = [0 0 0 Ne 2Ne 3Ne]

    vertmap = LinearIndices((Ne + 1, Ne + 1, Ne + 1))
    for (le, e) in enumerate(elemlocal)
        i, j, blck = Tuple(etoijb[e])
        elemtocoord[1, :, le] = bx[blck] .+ [i - 1 i i - 1 i]
        elemtocoord[2, :, le] = by[blck] .+ [j - 1 j - 1 j j]

        for n in 1:4
            ix = i + mod(n - 1, 2)
            jx = j + div(n - 1, 2)
            # set the vertices like they are the face vertices of a cube
            if blck == 1
                elemtovert[n, le] = vertmap[1, Ne + 2 - ix, jx]
            elseif blck == 2
                elemtovert[n, le] = vertmap[ix, 1, jx]
            elseif blck == 3
                elemtovert[n, le] = vertmap[Ne + 1, ix, jx]
            elseif blck == 4
                elemtovert[n, le] = vertmap[ix, jx, Ne + 1]
            elseif blck == 5
                elemtovert[n, le] = vertmap[ix, Ne + 1, Ne + 2 - jx]
            elseif blck == 6
                elemtovert[n, le] = vertmap[ix, Ne + 2 - jx, 1]
            end
        end
    end

    # no boundaries for a shell
    elemtobndy = zeros(Int, nface, length(elemlocal))

    # no faceconnections for a shell
    faceconnections = Array{Array{Int, 1}}(undef, 0)

    (elemtovert, elemtocoord, elemtobndy, faceconnections, collect(elemlocal))
end

"""
    cubedshellwarp(a, b, c, R = max(abs(a), abs(b), abs(c)))

Given points `(a, b, c)` on the surface of a cube, warp the points out to a
spherical shell of radius `R` based on the equiangular gnomonic grid proposed by
Ronchi, Iacono, Paolucci (1996) <https://doi.org/10.1006/jcph.1996.0047>

```
@article{RonchiIaconoPaolucci1996,
  title={The ``cubed sphere'': a new method for the solution of partial
         differential equations in spherical geometry},
  author={Ronchi, C. and Iacono, R. and Paolucci, P. S.},
  journal={Journal of Computational Physics},
  volume={124},
  number={1},
  pages={93--114},
  year={1996},
  doi={10.1006/jcph.1996.0047}
}
```

"""
function cubedshellwarp(a, b, c, R = max(abs(a), abs(b), abs(c)))

    function f(sR, ξ, η)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        x1 = sR / sqrt(X^2 + Y^2 + 1)
        x2, x3 = X * x1, Y * x1
        x1, x2, x3
    end

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : Face I from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : Face II from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : Face III from Ronchi, Iacono, Paolucci (1996)
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : Face IV from Ronchi, Iacono, Paolucci (1996)
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : Face V from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : Face VI from Ronchi, Iacono, Paolucci (1996)
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubedshellwarp: $a, $b, $c")
    end

    return x1, x2, x3
end

"""
    cubedshellunwarp(x1, x2, x3)

The inverse of [`cubedshellwarp`](@ref).
"""
function cubedshellunwarp(x1, x2, x3)

    function g(R, X, Y)
        ξ = atan(X) * 4 / pi
        η = atan(Y) * 4 / pi
        R, R * ξ, R * η
    end

    R = hypot(x1, x2, x3)
    fdim = argmax(abs.((x1, x2, x3)))

    if fdim == 1 && x1 < 0
        # (-R, *, *) : Face I from Ronchi, Iacono, Paolucci (1996)
        a, b, c = g(-R, x2 / x1, x3 / x1)
    elseif fdim == 2 && x2 < 0
        # ( *,-R, *) : Face II from Ronchi, Iacono, Paolucci (1996)
        b, a, c = g(-R, x1 / x2, x3 / x2)
    elseif fdim == 1 && x1 > 0
        # ( R, *, *) : Face III from Ronchi, Iacono, Paolucci (1996)
        a, b, c = g(R, x2 / x1, x3 / x1)
    elseif fdim == 2 && x2 > 0
        # ( *, R, *) : Face IV from Ronchi, Iacono, Paolucci (1996)
        b, a, c = g(R, x1 / x2, x3 / x2)
    elseif fdim == 3 && x3 > 0
        # ( *, *, R) : Face V from Ronchi, Iacono, Paolucci (1996)
        c, b, a = g(R, x2 / x3, x1 / x3)
    elseif fdim == 3 && x3 < 0
        # ( *, *,-R) : Face VI from Ronchi, Iacono, Paolucci (1996)
        c, b, a = g(-R, x2 / x3, x1 / x3)
    else
        error("invalid case for cubedshellunwarp: $a, $b, $c")
    end

    return a, b, c
end

"""
   StackedCubedSphereTopology(mpicomm, Nhorz, Rrange;
                              boundary=(1,1)) <: AbstractTopology{3}

Generate a stacked cubed sphere topology with `Nhorz` by `Nhorz` cells for each
horizontal face and `Rrange` is the radius edges of the stacked elements.  This
topology actual creates a cube mesh, and the warping should be done after the
grid is created using the `cubedshellwarp` function. The coordinates of the
points will be of type `eltype(Rrange)`. The inner boundary condition type is
`boundary[1]` and the outer boundary condition type is `boundary[2]`.

The elements are stacked such that the vertical elements are contiguous in the
element ordering.

The elements of the brick are partitioned equally across the MPI ranks based
on a space-filling curve. Further, stacks are not split at MPI boundaries.

# Examples

We can build a cubed sphere mesh with 10 x 10 x 5 elements on each cube, total
elements is `10 * 10 * 5 * 6 = 3000`, with
```jldoctest brickmesh
using CLIMA.Topologies
using MPI
MPI.Init()
Nhorz = 10
Nstack = 5
Rrange = Float64.(accumulate(+,1:Nstack+1))
topology = StackedCubedSphereTopology(MPI.COMM_SELF, Nhorz, Rrange)

x1, x2, x3 = ntuple(j->reshape(topology.elemtocoord[j, :, :],
                            2, 2, 2, length(topology.elems)), 3)
for n = 1:length(x1)
   x1[n], x2[n], x3[n] = Topologies.cubedshellwarp(x1[n], x2[n], x3[n])
end
```
Note that the faces are listed in Cartesian order.
"""
function StackedCubedSphereTopology(
    mpicomm,
    Nhorz,
    Rrange;
    boundary = (1, 1),
    connectivity = :face,
    ghostsize = 1,
)
    T = eltype(Rrange)

    basetopo = CubedShellTopology(
        mpicomm,
        Nhorz,
        T;
        connectivity = connectivity,
        ghostsize = ghostsize,
    )

    dim = 3
    nvert = 2^dim
    nface = 2dim
    stacksize = length(Rrange) - 1

    nreal = length(basetopo.realelems) * stacksize
    nghost = length(basetopo.ghostelems) * stacksize

    elems = 1:(nreal + nghost)
    realelems = 1:nreal
    ghostelems = nreal .+ (1:nghost)

    sendelems =
        similar(basetopo.sendelems, length(basetopo.sendelems) * stacksize)

    for i in 1:length(basetopo.sendelems), j in 1:stacksize
        sendelems[stacksize * (i - 1) + j] =
            stacksize * (basetopo.sendelems[i] - 1) + j
    end

    ghostfaces = similar(basetopo.ghostfaces, nface, length(ghostelems))
    ghostfaces .= false

    for i in 1:length(basetopo.ghostelems), j in 1:stacksize
        e = stacksize * (i - 1) + j
        for f in 1:(2 * (dim - 1))
            ghostfaces[f, e] = basetopo.ghostfaces[f, i]
        end
    end

    sendfaces = similar(basetopo.sendfaces, nface, length(sendelems))
    sendfaces .= false

    for i in 1:length(basetopo.sendelems), j in 1:stacksize
        e = stacksize * (i - 1) + j
        for f in 1:(2 * (dim - 1))
            sendfaces[f, e] = basetopo.sendfaces[f, i]
        end
    end

    elemtocoord = similar(basetopo.elemtocoord, dim, nvert, length(elems))

    for i in 1:length(basetopo.elems), j in 1:stacksize
        # i is base element
        # e is stacked element
        e = stacksize * (i - 1) + j

        # v is base vertex
        for v in 1:(2^(dim - 1))
            for d in 1:dim # dim here since shell is embedded in 3-D
                # v is lower stacked vertex
                elemtocoord[d, v, e] = basetopo.elemtocoord[d, v, i] * Rrange[j]
                # 2^(dim-1) + v is higher stacked vertex
                elemtocoord[d, 2^(dim - 1) + v, e] =
                    basetopo.elemtocoord[d, v, i] * Rrange[j + 1]
            end
        end
    end

    elemtoelem = similar(basetopo.elemtoelem, nface, length(elems))
    elemtoface = similar(basetopo.elemtoface, nface, length(elems))
    elemtoordr = similar(basetopo.elemtoordr, nface, length(elems))
    elemtobndy = similar(basetopo.elemtobndy, nface, length(elems))

    for e in 1:(length(basetopo.elems) * stacksize), f in 1:nface
        elemtoelem[f, e] = e
        elemtoface[f, e] = f
        elemtoordr[f, e] = 1
        elemtobndy[f, e] = 0
    end

    for i in 1:length(basetopo.realelems), j in 1:stacksize
        e1 = stacksize * (i - 1) + j

        for f in 1:(2 * (dim - 1))
            e2 = stacksize * (basetopo.elemtoelem[f, i] - 1) + j

            elemtoelem[f, e1] = e2
            elemtoface[f, e1] = basetopo.elemtoface[f, i]

            # since the basetopo is 2-D we only need to worry about two
            # orientations
            @assert basetopo.elemtoordr[f, i] ∈ (1, 2)
            #=
            orientation 1:
            2---3     2---3
            |   | --> |   |
            0---1     0---1
            same:
            (a,b) --> (a,b)

            orientation 3:
            2---3     3---2
            |   | --> |   |
            0---1     1---0
            reverse first index:
            (a,b) --> (N+1-a,b)
            =#
            elemtoordr[f, e1] = basetopo.elemtoordr[f, i] == 1 ? 1 : 3
        end

        # If top or bottom of stack set neighbor to self on respective face
        elemtoelem[2 * (dim - 1) + 1, e1] =
            j == 1 ? e1 : stacksize * (i - 1) + j - 1
        elemtoelem[2 * (dim - 1) + 2, e1] =
            j == stacksize ? e1 : stacksize * (i - 1) + j + 1

        elemtoface[2 * (dim - 1) + 1, e1] =
            j == 1 ? 2 * (dim - 1) + 1 : 2 * (dim - 1) + 2
        elemtoface[2 * (dim - 1) + 2, e1] =
            j == stacksize ? 2 * (dim - 1) + 2 : 2 * (dim - 1) + 1

        elemtoordr[2 * (dim - 1) + 1, e1] = 1
        elemtoordr[2 * (dim - 1) + 2, e1] = 1
    end

    # Set the top and bottom boundary condition
    for i in 1:length(basetopo.elems)
        eb = stacksize * (i - 1) + 1
        et = stacksize * (i - 1) + stacksize

        elemtobndy[2 * (dim - 1) + 1, eb] = boundary[1]
        elemtobndy[2 * (dim - 1) + 2, et] = boundary[2]
    end

    nabrtorank = basetopo.nabrtorank
    nabrtorecv = UnitRange{Int}[
        UnitRange(
            stacksize * (first(basetopo.nabrtorecv[n]) - 1) + 1,
            stacksize * last(basetopo.nabrtorecv[n]),
        ) for n in 1:length(nabrtorank)
    ]
    nabrtosend = UnitRange{Int}[
        UnitRange(
            stacksize * (first(basetopo.nabrtosend[n]) - 1) + 1,
            stacksize * last(basetopo.nabrtosend[n]),
        ) for n in 1:length(nabrtorank)
    ]

    StackedCubedSphereTopology{T}(
        BoxElementTopology{3, T}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            basetopo.origsendorder,
            true,
        ),
        stacksize,
    )
end

"""
    grid1d(a, b[, stretch::AbstractGridStretching]; elemsize, nelem)

Discretize the 1D interval [`a`,`b`] into elements.
Exactly one of the following keyword arguments must be provided:
- `elemsize`: the average element size, or
- `nelem`: the number of elements.

The optional `stretch` argument allows stretching, otherwise the element sizes
will be uniform.

Returns either a range object or a vector containing the element boundaries.
"""
function grid1d(a, b, stretch = nothing; elemsize = nothing, nelem = nothing)
    xor(nelem === nothing, elemsize === nothing) ||
    error("Either `elemsize` or `nelem` arguments must be provided")
    if elemsize !== nothing
        nelem = round(Int, abs(b - a) / elemsize)
    end
    grid1d(a, b, stretch, nelem)
end
function grid1d(a, b, ::Nothing, nelem)
    range(a, stop = b, length = nelem + 1)
end

# TODO: document these
abstract type AbstractGridStretching end

"""
    SingleExponentialStretching(A)

Apply single-exponential stretching: `A > 0` will increase the density of points
at the lower boundary, `A < 0` will increase the density at the upper boundary.

# Reference
* "Handbook of Grid Generation" J. F. Thompson, B. K. Soni, N. P. Weatherill
  (Editors) RCR Press 1999, §3.6.1 Single-Exponential Function
"""
struct SingleExponentialStretching{T} <: AbstractGridStretching
    A::T
end
function grid1d(
    a::A,
    b::B,
    stretch::SingleExponentialStretching,
    nelem,
) where {A, B}
    F = float(promote_type(A, B))
    s = range(zero(F), stop = one(F), length = nelem + 1)
    a .+ (b - a) .* expm1.(stretch.A .* s) ./ expm1(stretch.A)
end

struct InteriorStretching{T} <: AbstractGridStretching
    attractor::T
end
function grid1d(a::A, b::B, stretch::InteriorStretching, nelem) where {A, B}
    F = float(promote_type(A, B))
    coe = F(2.5)
    s = range(zero(F), stop = one(F), length = nelem + 1)
    range(a, stop = b, length = nelem + 1) .+
    coe .* (stretch.attractor .- (b - a) .* s) .* (1 .- s) .* s
end

end


text:
@testset "Analytical Nonlin shrinkage (ref⭒) " begin
    ANS = AnalyticalNonlinearShrinkage()
    for X̂ ∈ test_matrices
        size(X̂, 1) < 12 && continue
        ref_result = matlab_ledoitwolf_analytical_shrinkage(X̂)
        c = cov(ANS, X̂); @test c ≈ ref_result
        @test issymmetric(c)
    end

    # weights are not currently exported but it seems to be a good
    # idea to ensure proper errors
    fw1 = FrequencyWeights(rand(1:10, size(test_matrices[1], 1)))
    @test_throws ErrorException cov(ANS, test_matrices[1], fw1)
    @test_throws ErrorException cov(ANS, test_matrices[1], fw1; dims=2)
    @test_throws ErrorException cov(ANS, test_matrices[1], fw1; mean=nothing)
end

ext: jl


text:
using LightXML

struct AnimatedTile
    name::String
    path::String

    origX::Number
    origY::Number

    posX::Number
    posY::Number

    delay::Number
end

struct Mask
    tiles::UInt8
    ignores::UInt8
end

# Relevant characters of the tile mask string
# We never need to check the centre tile, otherwise we wouldn't be in the mask check
const lookupmask = [1, 2, 3, 5, 7, 9, 10, 11]

# Convert mask string into a bitmask.
# ignoremask is 1 for every tile that should NOT be ignored.
function getMaskFromString(s::String)
    data = s[lookupmask]
    tilemask = UInt8(0)
    ignoremask = UInt8(0)
    for (i, c) in enumerate(data)
        tilemask |= UInt8(c == '1') << (8 - i)
        ignoremask |= UInt8(c != 'x') << (8 - i)
    end
    return Mask(tilemask, ignoremask)
end

function sortScore(mask::Mask)
    # Counts the number of 0-bits in the ignore string
    # by counting the number of 1-bits.
    return 8 - sum([(mask.ignores >> i) & 0b1 for i = 0:7])
end

# Get quads for a given mask
function convertTileString(s::String)
    parts = split(s, ";")
    res = Coord[]

    for part in parts
        rawX, rawY = split(part, ",")
        x, y = parse(Int, rawX), parse(Int, rawY)

        push!(res, Coord(x, y))
    end

    return res
end

function loadAnimatedTilesXML(fn::String)
    xdoc = parse_file(fn)
    xroot = root(xdoc)

    animatedTiles = AnimatedTile[]

    for elem in child_elements(xroot)
        push!(animatedTiles, AnimatedTile(
            attribute(elem, "name"),
            attribute(elem, "path"),

            parseNumber(attribute(elem, "origX")),
            parseNumber(attribute(elem, "origY")),

            parseNumber(attribute(elem, "posX")),
            parseNumber(attribute(elem, "posY")),

            parseNumber(attribute(elem, "delay"))
        ))
    end

    return animatedTiles
end

struct Coord
    x::Int
    y::Int
end

struct MaskData
    mask::Mask
    coords::Array{Coord, 1}
    sprites::Array{String, 1}
end

struct Ignore
    wildcard::Bool
    ignores::String

    Ignore(s::String) = new('*' in s, filter(c -> c != '*', s))
end

getTilesetId(id::String) = getTilesetId(id[1])
getTilesetId(id::Char) = codepoint(id) <= 0xff ? id : nothing

function loadTilesetXML(fn::String)
    if !isfile(fn)
        error("File '$fn' not found")
    end

    xdoc = parse_file(fn)
    xroot = root(xdoc)

    paths = Dict{Char, String}()
    masks = Dict{Char, Array{MaskData, 1}}()
    padding = Dict{Char, Array{Coord, 1}}()
    center = Dict{Char, Array{Coord, 1}}()
    ignores = Dict{Char, Ignore}()

    orderedElements = sort(collect(child_elements(xroot)), by=set -> attribute(set, "copy") !== nothing)

    for set in orderedElements
        rawId = attribute(set, "id")
        id = getTilesetId(rawId)
        path = attribute(set, "path")
        copyTileset = attribute(set, "copy")
        ignore = attribute(set, "ignores")

        if id === nothing
            error("Tileset ID must be ASCII character, got '$rawId'")
        end

        if haskey(paths, id)
            error("Tileset ID '$id' is already in use")
        end

        paths[id] = "tilesets/$path"

        if ignore !== nothing
            ignores[id] = Ignore(ignore)
        end

        if copyTileset !== nothing
            tilesetId = copyTileset[1]
            padding[id] = deepcopy(padding[tilesetId])
            center[id] = deepcopy(center[tilesetId])
            masks[id] = deepcopy(masks[tilesetId])
        end

        currMasks = MaskData[]

        for c in child_elements(set)
            mask = attribute(c, "mask")
            tilesString = attribute(c, "tiles")
            quads = convertTileString(tilesString)

            if mask == "padding"
                padding[id] = quads

            elseif mask == "center"
                center[id] = quads

            else
                spritesRaw = attribute(c, "sprites")
                sprites = isa(spritesRaw, String) ? split(spritesRaw, ",") : String[]

                push!(currMasks, MaskData(
                    getMaskFromString(mask),
                    quads,
                    sprites
                ))
            end
        end

        if !isempty(currMasks)
            masks[id] = currMasks
            sort!(masks[id], by=m -> sortScore(m.mask), rev=false)
        end
    end

    return paths, masks, padding, center, ignores
end

struct TilerMeta
    paths::Dict{Char, String}
    masks::Dict{Char, Array{MaskData, 1}}
    paddings::Dict{Char, Array{Coord, 1}}
    centers::Dict{Char, Array{Coord, 1}}
    ignores::Dict{Char, Ignore}
end

TilerMeta(s::String) = TilerMeta(loadTilesetXML(s)...)

fgTilerMeta = nothing
bgTilerMeta = nothing
animatedTilesMeta = nothing

# Filename => (result, load time)
const loadedXMLCache = Dict{String, Tuple{Any, Number}}()

function displayCustomXMLWarning(title::String, message::String)
    topMostInfoDialog("$title\n$message")
end

function displayCustomXMLWarning(title::String, exception::Exception)
    message = sprint(showerror, exception)

    displayCustomXMLWarning(title, message)
end

function getCustomTilesXMLPath(side::Side, filename::String, key::String)
    meta = get(Dict{String, Any}, side.data, "meta")
    metaPath = get(meta, key, "")

    hasRoot, modRoot = getModRoot(filename)
    xmlPath = joinpath(modRoot, metaPath)

    if !isempty(metaPath) && hasRoot
        return xmlPath
    end
end

getCustomTilesXMLPath(side::Nothing, filename::String, key::String) = nothing

function loadCustomXML(loader::Union{Function, Type}, filename::String, errorTitle::String, force::Bool=false, quiet::Bool=false)
    if !force
        result, modified = get(loadedXMLCache, filename, (nothing, 0))

        if modified >= mtime(filename)
            return result
        end
    end

    try
        result = loader(filename)
        loadedXMLCache[filename] = (result, mtime(filename))

        return result

    catch e
        if !quiet
            displayCustomXMLWarning(errorTitle, e)

            println(Base.stderr, e)

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println(Base.stderr, "")
            end
        end

        loadedXMLCache[filename] = (nothing, mtime(filename))
    end

    return nothing
end

# Good enough to cache and retrieve via the custom XML loader
# We don't need any error message, vanilla assets shouldn't fail
# Never need to force reloads, that is mostly just for error message sake
function loadVanillaXML(loader::Union{Function, Type}, filename::String)
    return loadCustomXML(loader, filename, "", false, true)
end

function loadTilesMeta(side::Union{Side, Nothing}, filename::String, fg::Bool=false, force::Bool=false)
    defaultPath = joinpath(storageDirectory, "XML", fg ? "ForegroundTiles.xml" : "BackgroundTiles.xml")
    customPath = getCustomTilesXMLPath(side, filename, fg ? "ForegroundTiles" : "BackgroundTiles")

    if customPath !== nothing
        xmlType = fg ? "Foreground Tiles" : "Background Tiles"
        errorTitle = "Failed to load custom $xmlType XML"
        customMeta = loadCustomXML(TilerMeta, customPath, errorTitle, force)

        if customMeta !== nothing
            return customMeta
        end
    end

    return loadVanillaXML(TilerMeta, defaultPath)
end

function loadAnimatedTilesMeta(side::Union{Side, Nothing}, filename::String, force::Bool=false)
    defaultPath = joinpath(storageDirectory, "XML",  "AnimatedTiles.xml")
    customPath = getCustomTilesXMLPath(side, filename, "AnimatedTiles")

    if customPath !== nothing
        errorTitle = "Failed to load custom Animated Tiles XML"
        customMeta = loadCustomXML(loadAnimatedTilesXML, customPath, errorTitle, force)

        if customMeta !== nothing
            return customMeta
        end
    end

    return loadVanillaXML(loadAnimatedTilesXML, defaultPath)
end

function loadXMLMeta(force::Bool=true)
    side = loadedState.side
    filename = loadedState.filename

    global fgTilerMeta = loadTilesMeta(side, filename, true, force)
    global bgTilerMeta = loadTilesMeta(side, filename, false, force)
    global animatedTilesMeta = loadAnimatedTilesMeta(side, filename, force)
end

function getTile(tiles::Tiles, x::Int, y::Int)
    return get(tiles.data, (y, x), ' ')
end

function checkPadding(data::Array{Char, 2}, x::Int, y::Int, width::Int, height::Int)
    # Checks for '0' even though getIfInboundsKnownSize might return ' '
    # This is due to quirks/special flags in the actuall autotiler
    return getIfInboundsKnownSize(data, y - 2, x, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y + 2, x, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y, x - 2, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y, x + 2, height, width, ' ') == '0'
end

const ignoresDefault = Ignore("")

# Return the "mask" value for a tile
# Returns as UInt8 for performance
# Manual string iteratation for performance
function checkTile(value::Char, target::Char, ignore::Ignore=ignoresDefault)
    if target == '0' || (ignore.wildcard && value != target)
        return 0b0

    else
        for c in ignore.ignores
            if target == c
                return 0b0
            end
        end
    end

    return 0b1
end

function checkMask(mask::Mask, adj::UInt8)
    return ((adj ⊻ mask.tiles) & mask.ignores) == 0
end

function getIfInboundsKnownSize(data::Array{Char, 2}, y::Int, x::Int, height::Int, width::Int, default::Char='0')::Char
    if 1 <= y <= height && 1 <= x <= width
        return @inbounds data[y, x]

    else
        return default
    end
end

# Unrolled for performance
function getAdjacencyMask(data::Array{Char, 2}, y::Int, x::Int, height::Int, width::Int, default::Char, ignores::Ignore)::UInt8
    return checkTile(default, getIfInboundsKnownSize(data, y - 1, x - 1, height, width, default), ignores) << 7 |
        checkTile(default, getIfInboundsKnownSize(data, y - 1, x, height, width, default), ignores) << 6 |
        checkTile(default, getIfInboundsKnownSize(data, y - 1, x + 1, height, width, default), ignores) << 5 |

        checkTile(default, getIfInboundsKnownSize(data, y, x - 1, height, width, default), ignores) << 4 |
        checkTile(default, getIfInboundsKnownSize(data, y, x + 1, height, width, default), ignores) << 3 |

        checkTile(default, getIfInboundsKnownSize(data, y + 1, x - 1, height, width, default), ignores) << 2 |
        checkTile(default, getIfInboundsKnownSize(data, y + 1, x, height, width, default), ignores) << 1 |
        checkTile(default, getIfInboundsKnownSize(data, y + 1, x + 1, height, width, default), ignores)
end

const defaultPaddingCoords = Coord[Coord(0, 0)]
const maskDataDefault = MaskData[]

function getMaskQuads(x::Int, y::Int, tiles::Tiles, meta::TilerMeta)
    data = tiles.data
    height, width = size(data)
    value = @inbounds tiles.data[y, x]

    masks = get(meta.masks, value, maskDataDefault)::Array{MaskData, 1}
    ignores = get(meta.ignores, value, ignoresDefault)::Ignore

    adjacent = getAdjacencyMask(data, y, x, height, width, value, ignores)

    for data in masks
        if checkMask(data.mask, adjacent)
            return data.coords, data.sprites
        end
    end

    if checkPadding(tiles.data, x, y, width, height)
        return get(meta.paddings, value, defaultPaddingCoords), String[]

    else
        return get(meta.centers, value, defaultPaddingCoords), String[]
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
module lmmsim

using DataFrames, DataStructures, MixedModels
export
    dat,
    makedata,
    makeθ,
    mods,
    simulate

include("makedata.jl")
include("simulate.jl")

end # module


text:
"""
    closest(collection, x)

Find the index of the value in `collection` closest to `x`.
"""
closest(x) = Base.Fix2(closest, x)
function closest(collection, x)
    if eltype(collection) <: Number
        return closest_number(collection, x)
    else
        throw(MethodError("clossest does not work for eltype $(eltype(collection))"))
    end
end

function closest_number(collection, x)
    if issorted(collection)
        return _closest_number_forward(collection, x)
    elseif issorted(collection, order=Base.Reverse)
        return _closest_number_reverse(collection, x)
    else
        return _closest_number_search(collection, x)
    end
end

function _closest_number_forward(collection, x)
    i = findfirst(>(x), collection)
    if i == firstindex(collection)
        return i
    else
        iminone = i - oneunit(i)
        if (x - @inbounds(collection[iminone])) < (@inbounds(collection[i]) - x)
            return iminone
        else
            return i
        end
    end
end

function _closest_number_reverse(collection, x)
    i = findlast(>(x), collection)
    if i == lastindex(collection)
        return i
    else
        ipone = i + oneunit(i)
        if (x - @inbounds(collection[ipone])) < (@inbounds(collection[i]) - x)
            return ipone
        else
            return i
        end
    end
end

function _closest_number_search(collection, x)
    p = pairs(collection)
    i, ival = first(p)
    cmpval = abs(ival - x)
    @inbounds for (k,v) in p
        tmp_cmpval = abs(v - x)
        if (tmp_cmpval < cmpval)
            cmpval = tmp_cmpval
            i = k
        end
    end
    return i
end

@propagate_inbounds function ArrayInterface.to_index(::IndexAxis, axis, arg::Base.Fix2{typeof(closest),<:Any})
    idx = arg(keys(axis))
    @boundscheck if idx === nothing
        throw(BoundsError(axis, arg))
    end
    return Int(idx)
end
function Base.checkindex(::Type{Bool}, axis::AbstractAxis, arg::Base.Fix2{typeof(closest),<:Any})
    return arg(keys(axis)) !== nothing
end


text:
#= ForwardDiff as a dependency is removed from this version
import MechanicalSketch
import MechanicalSketch: color_with_lumin, sethue, O, WI, EM, FS, HE, finish, background
import MechanicalSketch: PALETTE, setfont, settext
import MechanicalSketch: dimension_aligned,line, circle, move, do_action
import MechanicalSketch: generate_complex_potential_source, generate_complex_potential_vortex
import MechanicalSketch: @import_expand, x_y_iterators_at_pixels
import MechanicalSketch: place_image, set_scale_sketch, lenient_min_max
import MechanicalSketch: ∙, ∇_rectangle, empty_figure, BinLegendVector, draw_legend
import Interpolations: interpolate, Gridded, Linear, Flat, extrapolate
# using BenchmarkTools

let
# This file is tested line by line using BenchmarkTools. Tests are commented out.

BACKCOLOR = color_with_lumin(PALETTE[8], 10);
include("test_functions_24.jl")
restart(BACKCOLOR)

# Plot the top figure, flow field from test_23.jl visualized with color for direction.
A = flowfield_23();
set_scale_sketch(PHYSWIDTH_23, round(Int, SCREEN_WIDTH_FRAC_23 * WI))
mi, ma = lenient_min_max(A)
toplegend = BinLegendVector(;operand_example = first(A),
        max_magn_legend = ma, noof_magn_bins = 30, noof_ang_bins = 36,
        name = :Velocity)
colmat = toplegend.(A)
upleftpoint, lowrightpoint = place_image(O + (-2EM, -0.25HE + EM), colmat)
legendpos = lowrightpoint + (EM, 0) + (0.0m, PHYSHEIGHT_23)
draw_legend(legendpos, toplegend)

# Centre of the bottom figure
OB = O + (-2EM, + 0.25HE + 0.5EM )

# We are going to use the velocity vector field in a lot of calculations,
# and interpolate between the calculated pixel values.
# We start by defining iterators for coordinates

xs, ys = x_y_iterators_at_pixels(;physwidth = PHYSWIDTH_23, physheight = PHYSHEIGHT_23)
# Note that the image matrix has maximum y in the first index, minimum in the last first index.
ys = reverse(ys)
# Example image, small scale:
# (x, y ) ∈ ({-2, 0.1, 2}, { -1, 1})
# ... the coordinates would be stored as
# M_xy = [(-2,  1)m  ( 0.1,  1)m (2,  1)m;
#         (-2, -1)m  ( 0.1, -1)m (2, -1)m]
#
# ...but 'interpolations' expect coordinates on this form:
#  M_xy_i =  [ (-2,   -1)m   (-2,   1)m ;
#              (0.1,  -1)m   ( 0.1, 1)m ;
#              ( 2,   -1)m   ( 2,   1)m ]
#
# Let's say the values corresponding to M_xy were x * y / m²:
#    M = [ -2        0.1      2;
#           2       -0.1     -2]
#
#  Then 'interpolations' would expect
#     [  2     -2 ;
#       -0.1    0.1 ;
#       -2      2 ]
#  but we would actually want to [row, col] index into M in this way:
#
#     [ [ 2  , 1]   [ 1 , 1];
#       [ 2  , 2]   [ 1 , 2];
#       [ 2  , 3]   [ 1,  3]]
#  if we instead index into transpose(M):
#     [ [ 1  , 2]   [ 1 , 1];
#       [ 2  , 2]   [ 2 , 1];
#       [ 3  , 2]   [ 3,  1]]
#
# and generally
#  M_xy_test[ii, jj] == M_xy[i, j]
#  ii = end + 1 - j
#  jj = i
# TODO fix this, used to work with earlier versions.
fxy_inter = interpolate((xs, ys), transpose(A)[ : , end:-1:1], Gridded(Linear()));
fxy = extrapolate(fxy_inter, Flat());
fxy(0.0m, 0.0m)
fxy(0.0m, 0.0001m)
fxy(100.0m, 0.0001m)

circle(OB, 0.999m, :stroke) # origo
dimension_aligned(OB + (-PHYSWIDTH_23 / 2, PHYSHEIGHT_23 / 2), OB + (PHYSWIDTH_23 / 2, PHYSHEIGHT_23 / 2))
dimension_aligned(OB, OB + (0.0, 1.0)m)
dimension_aligned(OB + (-PHYSWIDTH_23 / 2, - PHYSHEIGHT_23 / 2 ),  OB +  (-PHYSWIDTH_23 / 2, PHYSHEIGHT_23 / 2 ))

setfont("DejaVu Sans", FS)
str = "RK4: 10.4μs vs Euler: 2.7μs"
settext(str, OB + (-WI/3 + EM, + 4EM), markup = true)
setfont("Calibri", FS)

sethue(PALETTE[1])
vx = collect(range(-0.0m, 0.0m, length = 10))
vy = collect(range(0.5m, 0.0m, length = 10))
# @btime
euler_forwardsteps_1!(fxy, vx, vy, 1.5s) # 2.800 μs (0 allocations: 0 bytes)

ex, ey = vx[end], vy[end]
prline_1(OB, vx, vy)
sethue(PALETTE[3])
vx = collect(range(-0.0m, 0.0m, length = 10))
vy = collect(range(0.5m, 0.0m, length = 10))
# @btime
rrk4_forwardsteps_1!(fxy, vx, vy, 1.5s) #10.899 μs (0 allocations: 0 bytes)

exrk, eyrk = vx[end], vy[end]
dimension_aligned(OB + (ex, ey),  OB + (exrk, eyrk))
prline_1(OB, vx, vy)

# println("\nInlining hint, inner functions:")
# @btime
euler_forwardsteps_2!(fxy, vx, vy, 1.5s) # 2.667 μs (0 allocations: 0 bytes)
# @btime
rrk4_forwardsteps_2!(fxy, vx, vy, 1.5s) #  10.300 μs (0 allocations: 0 bytes)

# println("\nNo assertions, inner functions:")
# @btime
euler_forwardsteps_3!(fxy, vx, vy, 1.5s) # 2.644 μs (0 allocations: 0 bytes)
# @btime
rrk4_forwardsteps_3!(fxy, vx, vy, 1.5s) # 10.399 μs (0 allocations: 0 bytes)

# println("\nNo outer function assertions:")

# @btime
euler_forwardsteps_4!(fxy, vx, vy, 1.5s) # 2.745 μs (0 allocations: 0 bytes)
# @btime
rrk4_forwardsteps_4!(fxy, vx, vy, 1.5s) # 11.4 μs (0 allocations: 0 bytes)

# println("Fastmath:")

# @btime
euler_forwardsteps_5!(fxy, vx, vy, 1.5s) # 2.745 μs (0 allocations: 0 bytes)
# @btime
rrk4_forwardsteps_5!(fxy, vx, vy, 1.5s) # 10.4 μs (0 allocations: 0 bytes)

# Neither fastmath nor inlining has any worthwhile effect.

setfont("DejaVu Sans", FS)
str = """
    Flowfield clamped at $CUTOFF_23.\n
    RK4 is four times slower than\nEuler for the same number of\nsteps.\n
    However, the Euler streamlines\n are visually inaccurate here.
    """
settext(str, O + (1.5m, -3.5m), markup = true)

end # Let
finish()
set_scale_sketch()
=#


text:
# Testing of the graph types

using Graphs
using Base.Test

#################################################
#
#  SimpleDirectedGraph
#
#################################################

sgd = simple_graph(3)
add_vertex!(sgd)

# concept test

@test implements_vertex_list(sgd)    == true
@test implements_edge_list(sgd)      == true
@test implements_vertex_map(sgd)     == true
@test implements_edge_map(sgd)       == true

@test implements_adjacency_list(sgd) == true
@test implements_incidence_list(sgd) == true
@test implements_bidirectional_adjacency_list(sgd) == true
@test implements_bidirectional_incidence_list(sgd) == true
@test implements_adjacency_matrix(sgd) == false

# properties

@test is_directed(sgd) == true
@test num_vertices(sgd) == 4
@test num_edges(sgd) == 0

for u = 1 : 4
    @test out_degree(u, sgd) == 0
    @test in_degree(u, sgd) == 0
    @test isempty(out_neighbors(u, sgd))
    @test isempty(out_edges(u, sgd))
    @test isempty(in_neighbors(u, sgd))
    @test isempty(in_edges(u, sgd))
end

@test collect_edges(sgd) == []

es = [  add_edge!(sgd, 1, 2)
        add_edge!(sgd, 1, 3)
        add_edge!(sgd, 2, 4)
        add_edge!(sgd, 3, 4) ]

@test num_edges(sgd) == 4

# outgoing

@test [out_degree(v, sgd) for v = 1:4] == [2, 1, 1, 0]

@test out_edges(1, sgd) == es[1:2]
@test out_edges(2, sgd) == es[[3]]
@test out_edges(3, sgd) == es[[4]]
@test isempty(out_edges(4, sgd))

@test collect(out_neighbors(1, sgd)) == [2, 3]
@test collect(out_neighbors(2, sgd)) == [4]
@test collect(out_neighbors(3, sgd)) == [4]
@test isempty(out_neighbors(4, sgd))

# incoming

@test [in_degree(v, sgd) for v = 1:4] == [0, 1, 1, 2]

@test isempty(in_edges(1, sgd))
@test in_edges(2, sgd) == es[[1]]
@test in_edges(3, sgd) == es[[2]]
@test in_edges(4, sgd) == es[3:4]

@test isempty(in_neighbors(1, sgd))
@test collect(in_neighbors(2, sgd)) == [1]
@test collect(in_neighbors(3, sgd)) == [1]
@test collect(in_neighbors(4, sgd)) == [2, 3]

@test collect_edges(sgd) == [Edge(1,1,2), Edge(2,1,3), Edge(3,2,4), Edge(4,3,4)]

sgd = simple_graph(0)
v = add_vertex!(sgd)
@test v == 1

#################################################
#
#  SimpleUndirectedGraph
#
#################################################

sgu = simple_graph(4, is_directed=false)

@test is_directed(sgu) == false
@test num_vertices(sgu) == 4
@test num_edges(sgu) == 0

for u = 1 : 4
    @test out_degree(u, sgu) == 0
    @test isempty(out_neighbors(u, sgu))
    @test isempty(out_edges(u, sgu))
end

add_edge!(sgu, 1, 2)
add_edge!(sgu, 1, 3)
add_edge!(sgu, 2, 4)
add_edge!(sgu, 3, 4)
add_edge!(sgu, 4, 1)

es = sgu.edges
rs = [revedge(e) for e in es]

@test num_edges(sgu) == 5

# outgoing

@test [out_degree(v, sgu) for v = 1:4] == [3, 2, 2, 3]

@test out_edges(1, sgu) == [es[1], es[2], rs[5]]
@test out_edges(2, sgu) == [rs[1], es[3]]
@test out_edges(3, sgu) == [rs[2], es[4]]
@test out_edges(4, sgu) == [rs[3], rs[4], es[5]]

@test collect(out_neighbors(1, sgu)) == [2, 3, 4]
@test collect(out_neighbors(2, sgu)) == [1, 4]
@test collect(out_neighbors(3, sgu)) == [1, 4]
@test collect(out_neighbors(4, sgu)) == [2, 3, 1]

# incoming

@test [in_degree(v, sgu) for v = 1:4] == [3, 2, 2, 3]

@test in_edges(1, sgu) == [rs[1], rs[2], es[5]]
@test in_edges(2, sgu) == [es[1], rs[3]]
@test in_edges(3, sgu) == [es[2], rs[4]]
@test in_edges(4, sgu) == [es[3], es[4], rs[5]]

@test collect(in_neighbors(1, sgu)) == [2, 3, 4]
@test collect(in_neighbors(2, sgu)) == [1, 4]
@test collect(in_neighbors(3, sgu)) == [1, 4]
@test collect(in_neighbors(4, sgu)) == [2, 3, 1]

for T in [ExVertex, Compat.ASCIIString]

#################################################
#
#  Extended and General directed graph
#
#################################################

    egd = graph(T[], ExEdge{T}[])
    @test is_directed(egd) == true

    names = ["a", "b", "c", "d"]

    for (i, x) in enumerate(names)
        v = add_vertex!(egd, x)
        @test vertex_index(v, egd) == i
    end
    vs = vertices(egd)

    @test num_vertices(egd) == 4
    @test num_edges(egd) == 0

    add_edge!(egd, vs[1], vs[2])
    add_edge!(egd, vs[1], vs[3])
    add_edge!(egd, vs[2], vs[4])
    add_edge!(egd, vs[3], vs[4])
    es = edges(egd)

    @test num_edges(egd) == 4

    # outgoing

    @test [out_degree(v, egd) for v in vs] == [2, 1, 1, 0]

    @test out_edges(vs[1], egd) == es[1:2]
    @test out_edges(vs[2], egd) == es[[3]]
    @test out_edges(vs[3], egd) == es[[4]]
    @test isempty(out_edges(vs[4], egd))

    @test collect(out_neighbors(vs[1], egd)) == [vs[2], vs[3]]
    @test collect(out_neighbors(vs[2], egd)) == vs[[4]]
    @test collect(out_neighbors(vs[3], egd)) == vs[[4]]
    @test isempty(out_neighbors(vs[4], egd))

    # incoming

    @test [in_degree(v, egd) for v in vs] == [0, 1, 1, 2]

    @test isempty(in_edges(vs[1], egd))
    @test in_edges(vs[2], egd) == es[[1]]
    @test in_edges(vs[3], egd) == es[[2]]
    @test in_edges(vs[4], egd) == es[3:4]

    @test isempty(in_neighbors(vs[1], egd))
    @test collect(in_neighbors(vs[2], egd)) == vs[[1]]
    @test collect(in_neighbors(vs[3], egd)) == vs[[1]]
    @test collect(in_neighbors(vs[4], egd)) == vs[2:3]

#################################################
#
#  Extended and General undirected graph
#
#################################################

    egu = graph(T[], ExEdge{T}[]; is_directed=false)
    @test is_directed(egu) == false

    names = ["a", "b", "c", "d"]

    for (i, x) in enumerate(names)
        v = add_vertex!(egu, x)
        @test vertex_index(v, egu) == i
    end
    vs = vertices(egu)

    @test num_vertices(egu) == 4
    @test num_edges(egu) == 0

    add_edge!(egu, vs[1], vs[2])
    add_edge!(egu, vs[1], vs[3])
    add_edge!(egu, vs[2], vs[4])
    add_edge!(egu, vs[3], vs[4])
    es = edges(egu)
    rs = [revedge(e) for e in es]

    @test num_edges(egu) == 4

    # outgoing

    @test [out_degree(v, egu) for v in vs] == [2, 2, 2, 2]

    @test out_edges(vs[1], egu) == [es[1], es[2]]
    @test out_edges(vs[2], egu) == [rs[1], es[3]]
    @test out_edges(vs[3], egu) == [rs[2], es[4]]
    @test out_edges(vs[4], egu) == [rs[3], rs[4]]

    @test collect(out_neighbors(vs[1], egu)) == [vs[2], vs[3]]
    @test collect(out_neighbors(vs[2], egu)) == [vs[1], vs[4]]
    @test collect(out_neighbors(vs[3], egu)) == [vs[1], vs[4]]
    @test collect(out_neighbors(vs[4], egu)) == [vs[2], vs[3]]

    # incoming

    @test [in_degree(v, egu) for v in vs] == [2, 2, 2, 2]

    @test in_edges(vs[1], egu) == [rs[1], rs[2]]
    @test in_edges(vs[2], egu) == [es[1], rs[3]]
    @test in_edges(vs[3], egu) == [es[2], rs[4]]
    @test in_edges(vs[4], egu) == [es[3], es[4]]

    @test collect(in_neighbors(vs[1], egu)) == [vs[2], vs[3]]
    @test collect(in_neighbors(vs[2], egu)) == [vs[1], vs[4]]
    @test collect(in_neighbors(vs[3], egu)) == [vs[1], vs[4]]
    @test collect(in_neighbors(vs[4], egu)) == [vs[2], vs[3]]

end

#################################################
#
#  graph() constructor
#
#################################################

v = [ExVertex(1,""),ExVertex(2,"")]
e = [ExEdge(1,v[1],v[2])]
g = graph(v,e,is_directed=true)
@test num_edges(g) == 1
@test num_vertices(g) == 2

# Integer vertices and edges
n = 100
m = 1000

# We ensure that vertices have random but different values
function  uniqueRands(range,nb::Int)
    result=Dict{Int,Bool}()
    i= 0
    while i < nb
        r = rand(1:10*n)
        if !haskey( result, r)
           i += 1
           result[r] = true
        end
    end
    collect(keys(result))
end

vs = uniqueRands(1:10*n,n)

es = Edge{Int}[]
for i in 1:m
  push!(es,Edge(i,vs[rand(1:n)],vs[rand(1:n)]))
end
g = graph(vs,es)
for i in 1:m
  add_edge!(g,vs[rand(1:n)],vs[rand(1:n)])
end

@test num_vertices(g) == n
@test num_edges(g) == 2 * m
for v in vs
  @test vertex_index(v,g) == g.indexof[v]
end
for i in 1:m
  @test edge_index(es[i],g) == i
end

# same for undirected graph
n = 100
m = 1000

vs = uniqueRands(1:10*n,n)

es = Edge{Int}[]
for i in 1:m
  push!(es,Edge(i,vs[rand(1:n)],vs[rand(1:n)]))
end
g = graph(vs,es,is_directed=false)
for i in 1:m
  add_edge!(g,vs[rand(1:n)],vs[rand(1:n)])
end

@test num_vertices(g) == n
@test num_edges(g) == 2 * m
for v in vs
  @test vertex_index(v,g) == g.indexof[v]
end
for i in 1:m
  @test edge_index(es[i],g) == i
end


text:
using PQBaseCamp
using DataFrames
using CSV
using Distributions

# load up a sample data file -
path_to_data_file = "$(pwd())/test/data/SPY.csv"
df = CSV.read(path_to_data_file, DataFrame)

# build a computation model -
model = LogReturnComputationModel()
model.ticker = "SPY"
model.data = df
model.map = :timestamp => :close
Δ_table = Δ(model)

# estimate a distribution from this data -
D = 𝒟(Laplace, Δ_table; colkey = :Δ)


text:
# ******************************************************************************************
# REQUIRES PACKAGE: PLOTS (incompatible with StatsPlots, I think)
# ******************************************************************************************
using Waterfall
include(joinpath(WATERFALL_DIR,"src","includes.jl"))

DATA_DIR = joinpath(WATERFALL_DIR,"data","pvrd")
df = CSV.read(joinpath(DATA_DIR,"pvrd2-investment-metrics.csv"), DataFrames.DataFrame)
df[!,:Label] .= string.(titlecase.(df[:,:Process])," (",df[:,:Step],")")
df[!,:Units] .= "Efficiency (%)"

# Define global arguments.
distribution=:normal
fuzziness=(0.01,0.3)
kwargs = (label=:Label, distribution=distribution, fuzziness=fuzziness)

p = Plot(df; samples=1000, ylabel=:Units, kwargs...)
v = collect_value(p.cascade)

covariance = LinearAlgebra.LowerTriangular(Statistics.cov(v, v; dims=2))
correlation = LinearAlgebra.LowerTriangular(Statistics.cor(v; dims=2))
# # variance = Statistics.var(v; dims=2) # == LinearAlgebra.diag(covariance)

labels = df[:,:Label]
steps = 15
global plots = fill(Plots.Plot(), (steps,steps))
for ii in 1:steps
    for jj in 1:ii
        begin xx=jj; yy=ii end

        plots[ii,jj] = if ii!==jj
            Plots.scatter(v[xx,:], v[yy,:],
                # xlabel=labels[xx],
                # ylabel=labels[yy],
                label="",
                # label=round(Statistics.cor(v[xx,:], v[yy,:]); sigdigits=4),
                xaxis=[],
                yaxis=[],
                xguidefontsize=8,
                yguidefontsize=8,
                markersize=2,
                markeralpha=0.5,
                markerstrokewidth=0,
            )
        else
            Plots.histogram(v[xx,:],
                xlabel=labels[xx],
                label="",
                xguidefontsize=8,
                yguidefontsize=8,
            )
        end
    end
end

N=5
plots = plots[1:N,1:N]
# allplots = plot(plots[1,1], plots[2,1], plots[1,2], plots[2,2], layout=(2,2))
allplots = plot(plots..., layout=size(plots), size=(800,800))

# printtype(x) = typeof(x)

# p1 = plot(...)
# p2 = plot(...)
# p3 = plot(...)
# plot(p1, p2, p3, layout = l)

# p = Plot{Parallel}(p)

# Luxor.@png begin
#     Luxor.fontsize(14)
#     Luxor.setmatrix([1 0 0 1 LEFT_BORDER TOP_BORDER])

#     draw(p; distribution=distribution, samples=samples)
#     # _draw_highlight(pdata, highlights[hh])

# end WIDTH+LEFT_BORDER+RIGHT_BORDER HEIGHT+TOP_BORDER+BOTTOM_BORDER "correlation.png"


text:
using ImageMorphology
using Test

@testset "ImageMorphology" begin
    include("convexhull.jl")
    include("connected.jl")
    include("dilation_and_erosion.jl")
    include("thinning.jl")
    include("imfill.jl")
end


text:
#
# generate Artificial Data on several manifolds
#
#
#
@doc raw"""
    artificialIn_SAR_image([pts=500])
generate an artificial InSAR image, i.e. phase valued data, of size `pts` x
`pts` points.

This data set was introduced for the numerical examples in

> Bergmann, R., Laus, F., Steidl, G., Weinmann, A.:
> _Second Order Differences of Cyclic Data and Applications in Variational Denoising_
> SIAM J. Imaging Sci., 7(4), 2916–2953, 2014.
> doi: [10.1137/140969993](https://dx.doi.org/10.1137/140969993)
> arxiv: [1405.5349](https://arxiv.org/abs/1405.5349)

"""
function artificialIn_SAR_image(pts::Integer)
    # variables
    # rotation of ellipse
    aEll = 35.0
    cosE = cosd(aEll)
    sinE = sind(aEll)
    aStep = 45.0
    cosA = cosd(aStep)
    sinA = sind(aStep)
    # main and minor axis of the ellipse
    axes_inv = [6, 25]
    # values for the hyperboloid
    mid_point = [0.275; 0.275]
    radius = 0.18
    values = [range(-0.5, 0.5; length=pts)...]
    # Steps
    aSteps = 60.0
    cosS = cosd(aSteps)
    sinS = sind(aSteps)
    l = 0.075
    midP = [-.475, -.0625]#.125, .55]
    img = zeros(Float64, pts, pts)
    for j in eachindex(values), i in eachindex(values)
        # ellipse
        Xr = cosE * values[i] - sinE * values[j]
        Yr = cosE * values[j] + sinE * values[i]
        v = axes_inv[1] * Xr^2 + axes_inv[2] * Yr^2
        k1 = v <= 1.0 ? 10.0 * pi * Yr : 0.0
        # circle
        Xr = cosA * values[i] - sinA * values[j]
        Yr = cosA * values[j] + sinA * values[i]
        v = ((Xr - mid_point[1])^2 + (Yr - mid_point[2])^2) / radius^2
        k2 = v <= 1.0 ? 4.0 * pi * (1.0 - v) : 0.0
        #
        Xr = cosS * values[i] - sinS * values[j]
        Yr = cosS * values[j] + sinS * values[i]
        k3 = 0.0
        for m in 1:8
            in_range = (abs(Xr + midP[1] + m * l) + abs(Yr + midP[2] + m * l)) ≤ l
            k3 += in_range ? 2 * pi * (m / 8) : 0.0
        end
        img[i, j] = mod(k1 + k2 + k3 - pi, 2 * pi) + pi
    end
    return img
end

@doc raw"""
    artificial_S1_slope_signal([pts=500, slope=4.])

Creates a Signal of (phase-valued) data represented on the
[Circle](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold)` `[Manifold](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) with increasing slope.

# Optional
* `pts` – (`500`) number of points to sample the function.
* `slope` – (`4.0`) initial slope that gets increased afterwards

This data set was introduced for the numerical examples in

> Bergmann, R., Laus, F., Steidl, G., Weinmann, A.:
> _Second Order Differences of Cyclic Data and Applications in Variational Denoising_
> SIAM J. Imaging Sci., 7(4), 2916–2953, 2014.
> doi: [10.1137/140969993](https://dx.doi.org/10.1137/140969993)
> arxiv: [1405.5349](https://arxiv.org/abs/1405.5349)

"""
function artificial_S1_slope_signal(pts::Integer=500, slope::Float64=4.0)
    t = range(0.0, 1.0; length=pts)
    f = zero(t)
    f[t .<= 1 / 6] .= -π / 2 .+ slope * π / 8 * t[t .<= 1 / 6]
    # In the following terms, the first max
    f[(1 / 6 .< t) .& (t .<= 1 / 3)] .=
        max(f[f .!= 0]...) .- slope * π / 4 * 1 / 6 .+
        slope * π / 4 .* t[(1 / 6 .< t) .& (t .<= 1 / 3)]
    f[(1 / 3 .< t) .& (t .<= 1 / 2)] .=
        max(f[f .!= 0]...) .- slope * π / 2 * 1 / 3 .+
        slope * π / 2 * t[(1 / 3 .< t) .& (t .<= 1 / 2)]
    f[(1 / 2 .< t) .& (t .<= 2 / 3)] .=
        max(f[f .!= 0]...) .- slope * π * 1 / 2 .+
        slope * π * t[(1 / 2 .< t) .& (t .<= 2 / 3)]
    f[(2 / 3 .< t) .& (t .<= 5 / 6)] .=
        max(f[f .!= 0]...) .- slope * 2 * π * 2 / 3 .+
        slope * 2 * π * t[(2 / 3 .< t) .& (t .<= 5 / 6)]
    f[5 / 6 .< t] .=
        max(f[f .!= 0]...) .- slope * 4 * π * 5 / 6 .+ slope * 4 * π * t[5 / 6 .< t]
    return mod.(f .+ Float64(π), Ref(2 * π)) .- Float64(π)
end

@doc raw"""
    artificial_S1_signal([pts=500])

generate a real-valued signal having piecewise constant, linear and quadratic
intervals with jumps in between. If the resulting manifold the data lives on,
is the [Circle](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.Manifold) the data is also wrapped
to $[-\pi,\pi)$.

# Optional
* `pts` – (`500`) number of points to sample the function

> Bergmann, R., Laus, F., Steidl, G., Weinmann, A.:
> _Second Order Differences of Cyclic Data and Applications in Variational Denoising_
> SIAM J. Imaging Sci., 7(4), 2916–2953, 2014.
> doi: [10.1137/140969993](https://dx.doi.org/10.1137/140969993)
> arxiv: [1405.5349](https://arxiv.org/abs/1405.5349)

"""
function artificial_S1_signal(pts::Integer=500)
    t = range(0.0, 1.0; length=pts)
    f = artificial_S1_signal.(t)
    return mod.(f .+ Float64(π), Ref(2 * π)) .- Float64(π)
end
@doc raw"""
    artificial_S1_signal(x)
evaluate the example signal $f(x), x ∈  [0,1]$,
of phase-valued data introduces in Sec. 5.1 of

> Bergmann, R., Laus, F., Steidl, G., Weinmann, A.:
> _Second Order Differences of Cyclic Data and Applications in Variational Denoising_
> SIAM J. Imaging Sci., 7(4), 2916–2953, 2014.
> doi: [10.1137/140969993](https://dx.doi.org/10.1137/140969993)
> arxiv: [1405.5349](https://arxiv.org/abs/1405.5349)

for values outside that intervall, this Signal is `missing`.
"""
function artificial_S1_signal(x::Real)
    if x < 0
        y = missing
    elseif x <= 1 / 4
        y = -24 * π * (x - 1 / 4)^2 + 3 / 4 * π
    elseif x <= 3 / 8
        y = 4 * π * x - π / 4
    elseif x <= 1 / 2
        y = -π * x - 3 * π / 8
    elseif x <= (3 * 0 + 19) / 32
        y = -(0 + 7) / 8 * π
    elseif x <= (3 * 1 + 19) / 32
        y = -(1 + 7) / 8 * π
    elseif x <= (3 * 2 + 19) / 32
        y = -(2 + 7) / 8 * π
    elseif x <= (3 * 3 + 19) / 32
        y = -(3 + 7) / 8 * π
    elseif x <= 1
        y = 3 / 2 * π * exp(8 - 1 / (1 - x)) - 3 / 4 * π
    else
        y = missing
    end
    return y
end
@doc raw"""
    artificial_S2_whirl_image([pts=64])
generate an artificial image of data on the 2 sphere,

# Arguments
* `pts` – (`64`) size of the image in `pts`$\times$`pts` pixel.

This example dataset was used in the numerical example in Section 5.5 of

> Laus, F., Nikolova, M., Persch, J., Steidl, G.:
> _A Nonlocal Denoising Algorithm for Manifold-Valued Images Using Second Order Statistics_,
> SIAM J. Imaging Sci., 10(1), 416–448, 2017.
> doi: [ 10.1137/16M1087114](https://dx.doi.org/10.1137/16M1087114)
> arxiv: [1607.08481](https://arxiv.org/abs/1607.08481)

It is based on [`artificial_S2_rotation_image`](@ref) extended by small whirl patches.
"""
function artificial_S2_whirl_image(pts::Int=64)
    M = Sphere(2)
    img = artificial_S2_rotation_image(pts, (0.5, 0.5))
    # Set WhirlPatches
    sc = pts / 64
    patchSizes = floor.(sc .* [9, 9, 9, 9, 11, 11, 11, 15, 15, 15, 17, 21])
    patchCenters =
        Integer.(
            floor.(
                sc .*
                [[35, 7] [25, 41] [32, 25] [7, 60] [10, 5] [41, 58] [11, 41] [23, 56] [
                    38, 45
                ] [16, 28] [55, 42] [51, 16]],
            ),
        )
    patchSigns = [1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1]
    for i in 1:length(patchSizes)
        pS = Integer(patchSizes[i])
        pSH = Integer(floor((patchSizes[i] - 1) / 2))
        pC = patchCenters[:, i]
        pC = max.(Ref(1), pC .- pS) .+ pS
        pSgn = patchSigns[i]
        s = pS % 2 == 0 ? 1 : 0
        r = [pC[1] .+ ((-pSH):(pSH + s)), pC[2] .+ ((-pSH):(pSH + s))]
        patch = artificial_S2_whirl_patch(pS)
        if pSgn == -1 # opposite ?
            patch = -patch
        end
        img[r...] = patch
    end
    return img
end
@doc raw"""
    artificial_S2_rotation_image([pts=64, rotations=(.5,.5)])
creates an image with a rotation on each axis as a parametrization.

# Optional Parameters
* `pts` – (`64`) number of pixels along one dimension
* `rotations` – (`(.5,.5)`) number of total rotations performed on the axes.

This dataset was used in the numerical example of Section 5.1 of

> Bačák, M., Bergmann, R., Steidl, G., Weinmann, A.:
> _A Second Order Non-Smooth Variational Model for Restoring Manifold-Valued Images_
> SIAM J. Sci. Comput. 38(1), A567–A597, 2016.
> doi: [10.1137/15M101988X](https://dx.doi.org/10.1137/15M101988X)
> arxiv: [1506.02409](https://arxiv.org/abs/1506.02409)

"""
function artificial_S2_rotation_image(
    pts::Int=64, rotations::Tuple{Float64,Float64}=(0.5, 0.5)
)
    M = Sphere(2)
    img = fill(zeros(3), pts, pts)
    north = [1.0, 0.0, 0.0]
    Rxy(a) = [cos(a) -sin(a) 0.0; sin(a) cos(a) 0.0; 0.0 0.0 1]
    Rxz(a) = [cos(a) 0.0 -sin(a); 0.0 1.0 0.0; sin(a) 0.0 cos(a)]
    for i in 1:pts
        for j in 1:pts
            x = i / pts * 2π * rotations[1]
            y = j / pts * 2π * rotations[2]
            img[i, j] = Rxy(x + y) * Rxz(x - y) * [0, 0, 1]
        end
    end
    return img
end

@doc raw"""
    artificial_S2_whirl_patch([pts=5])

create a whirl within the `pts`$\times$`pts` patch of
[Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html)(@ref)`(2)`-valued image data.

These patches are used within [`artificial_S2_whirl_image`](@ref).

# Optional Parameters
* `pts` – (`5`) size of the patch. If the number is odd, the center is the north
  pole.
"""
function artificial_S2_whirl_patch(pts::Int=5)
    patch = fill([0.0, 0.0, -1.0], pts, pts)
    scaleFactor = sqrt((pts - 1)^2 / 2) * 3 / π
    for i in 1:pts
        for j in 1:pts
            if i != (pts + 1) / 2 || j != (pts + 1) / 2
                α = atan((j - (pts + 1) / 2), (i - (pts + 1) / 2))
                β = sqrt((j - (pts + 1) / 2)^2 + (i - (pts + 1) / 2)^2) / scaleFactor
                patch[i, j] = [
                    sin(α) * sin(π / 2 - β), -cos(α) * sin(π / 2 - β), cos(π / 2 - β)
                ]
            end
        end
    end
    return patch
end
@doc raw"""
    artificial_S2_composite_bezier_curve()

Create the artificial curve in the `Sphere(2)` consisting of 3 segments between the four
points

````math
p_0 = \begin{bmatrix}0&0&1\end{bmatrix}^{\mathrm{T}},
p_1 = \begin{bmatrix}0&-1&0\end{bmatrix}^{\mathrm{T}},
p_2 = \begin{bmatrix}-1&0&0\end{bmatrix}^{\mathrm{T}},
p_3 = \begin{bmatrix}0&0&-1\end{bmatrix}^{\mathrm{T}},
````

where each segment is a cubic Bezér curve, i.e. each point, except $p_3$ has a first point
within the following segment $b_i^+$, $i=0,1,2$ and a last point within the previous
segment, except for $p_0$, which are denoted by $b_i^-$, $i=1,2,3$.
This curve is differentiable by the conditions $b_i^- = \gamma_{b_i^+,p_i}(2)$, $i=1,2$,
where $\gamma_{a,b}$ is the [`shortest_geodesic`](https://juliamanifolds.github.io/Manifolds.jl/stable/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any}) connecting $a$ and $b$.
The remaining points are defined as

````math
\begin{aligned}
    b_0^+ &= \exp_{p_0}\frac{\pi}{8\sqrt{2}}\begin{pmatrix}1&-1&0\end{pmatrix}^{\mathrm{T}},&
    b_1^+ &= \exp_{p_1}-\frac{\pi}{4\sqrt{2}}\begin{pmatrix}-1&0&1\end{pmatrix}^{\mathrm{T}},\\
    b_2^+ &= \exp_{p_2}\frac{\pi}{4\sqrt{2}}\begin{pmatrix}0&1&-1\end{pmatrix}^{\mathrm{T}},&
    b_3^- &= \exp_{p_3}-\frac{\pi}{8\sqrt{2}}\begin{pmatrix}-1&1&0\end{pmatrix}^{\mathrm{T}}.
\end{aligned}
````

This example was used within minimization of acceleration of the paper

> Bergmann, R., Gousenbourger, P.-Y.:
> _A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve_,
> Front. Appl. Math. Stat. 12, 2018.
> doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059)
> arxiv: [1807.10090](https://arxiv.org/abs/1807.10090)

"""
function artificial_S2_composite_bezier_curve()
    M = Sphere(2)
    d0 = [0.0, 0.0, 1.0]
    d1 = [0.0, -1.0, 0.0]
    d2 = [-1.0, 0.0, 0.0]
    d3 = [0.0, 0.0, -1.0]
    #
    # control points - where b1- and b2- are constructed by the C1 condition
    #
    # We define three segments: 1
    b00 = d0 # also known as p0
    ξ0 = π / (8.0 * sqrt(2.0)) .* [1.0, -1.0, 0.0] # staring direction from d0
    b01 = exp(M, d0, ξ0) # b0+
    ξ1 = π / (4.0 * sqrt(2)) .* [1.0, 0.0, 1.0]
    # b02 or b1- and b11 or b1+ are constructed by this vector with opposing sign
    # to achieve a C1 curve
    b02 = exp(M, d1, ξ1)
    b03 = d1
    # 2
    b10 = d1
    b11 = exp(M, d1, -ξ1) # yields c1 condition
    ξ2 = -π / (4 * sqrt(2)) .* [0.0, 1.0, -1.0]
    b12 = exp(M, d2, ξ2)
    b13 = d2
    # 3
    b20 = d2
    b21 = exp(M, d2, -ξ2)
    ξ3 = π / (8.0 * sqrt(2)) .* [-1.0, 1.0, 0.0]
    b22 = exp(M, d3, ξ3)
    b23 = d3
    # hence the matrix of controlpoints for the curve reads
    return [
        BezierSegment([b00, b01, b02, b03]),
        BezierSegment([b10, b11, b12, b13]),
        BezierSegment([b20, b21, b22, b23]),
    ]
end

@doc raw"""
    artificial_SPD_image([pts=64, stepsize=1.5])

create an artificial image of symmetric positive definite matrices of size
`pts`$\times$`pts` pixel with a jump of size `stepsize`.

This dataset was used in the numerical example of Section 5.2 of

> Bačák, M., Bergmann, R., Steidl, G., Weinmann, A.:
> _A Second Order Non-Smooth Variational Model for Restoring Manifold-Valued Images_
> SIAM J. Sci. Comput. 38(1), A567–A597, 2016.
> doi: [10.1137/15M101988X](https://dx.doi.org/10.1137/15M101988X)
> arxiv: [1506.02409](https://arxiv.org/abs/1506.02409)

"""
function artificial_SPD_image(pts::Int=64, stepsize=1.5)
    r = range(0; stop=1 - 1 / pts, length=pts)
    v1 = abs.(2 * pi .* r .- pi)
    v2 = pi .* r
    v3 = range(0; stop=3 * (1 - 1 / pts), length=2 * pts)
    data = fill(Matrix{Float64}(I, 3, 3), pts, pts)
    for row in 1:pts
        for col in 1:pts
            A = [cos(v1[col]) -sin(v1[col]) 0.0; sin(v1[col]) cos(v1[col]) 0.0; 0.0 0.0 1.0]
            B = [1.0 0.0 0.0; 0.0 cos(v2[row]) -sin(v2[row]); 0.0 sin(v2[row]) cos(v2[row])]
            C = [
                cos(v1[mod(col - row, pts) + 1]) 0 -sin(v1[mod(col - row, pts) + 1])
                0.0 1.0 0.0
                sin(v1[mod(col - row, pts) + 1]) 0.0 cos(v1[mod(col - row, pts) + 1])
            ]
            scale = [
                1 + stepsize / 2 * ((row + col) > pts ? 1 : 0)
                1 + v3[row + col] - stepsize * (col > pts / 2 ? 1 : 0)
                4 - v3[row + col] + stepsize * (row > pts / 2 ? 1 : 0)
            ]
            data[row, col] = Matrix(Symmetric(A * B * C * Diagonal(scale) * C' * B' * A'))
        end
    end
    return data
end
@doc raw"""
    artificial_SPD_image2([pts=64, fraction=.66])

create an artificial image of symmetric positive definite matrices of size
`pts`$\times$`pts` pixel with right hand side `fraction` is moved upwards.

This data set was introduced in the numerical examples of Section of

> Bergmann, R., Persch, J., Steidl, G.:
> _A Parallel Douglas Rachford Algorithm for Minimizing ROF-like Functionals on Images with Values in Symmetric Hadamard Manifolds_
> SIAM J. Imaging. Sci. 9(3), pp. 901-937, 2016.
> doi: [10.1137/15M1052858](https://dx.doi.org/10.1137/15M1052858)
> arxiv: [1512.02814](https://arxiv.org/abs/1512.02814)

"""
function artificial_SPD_image2(pts=64, fraction=0.66)
    Zl = 4.0 * Matrix{Float64}(I, 3, 3)
    # create a first matrix
    α = 2.0 * π / 3
    β = π / 3
    B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)]
    A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0]
    Zo = Matrix(Symmetric(A * B * Diagonal([2.0, 4.0, 8.0]) * B' * A'))
    # create a second matrix
    α = -4.0 * π / 3
    β = -π / 3
    B = [1.0 0.0 0.0; 0.0 cos(β) -sin(β); 0.0 sin(β) cos(β)]
    A = [cos(α) -sin(α) 0.0; sin(α) cos(α) 0.0; 0.0 0.0 1.0]
    Zt = A * B * Diagonal([8.0 / sqrt(2.0), 8.0, sqrt(2.0)]) * B' * A'
    data = fill(Matrix{Float64}(I, 3, 3), pts, pts)
    M = SymmetricPositiveDefinite(3)
    for row in 1:pts
        for col in 1:pts
            # (a) from Zo a part to Zt
            C = Zo
            if (row > 1) # in X direction
                C = exp(
                    M,
                    C,
                    log(M, C, Zt),
                    (row - 1) / (2 * (pts - 1)) + ((row > fraction * pts) ? 1 / 2 : 0.0),
                )
            end
            if (col > 1) # and then in Y direction
                C = exp(
                    M,
                    C,
                    vector_transport_to(
                        M, Symmetric(Zo), log(M, Zo, Zl), Symmetric(C), ParallelTransport()
                    ),
                    (col - 1.0) / (pts - 1),
                )
            end
            data[row, col] = C
        end
    end
    return data
end

@doc raw"""
    artificial_S2_lemniscate(p [,pts=128,a=π/2,interval=[0,2π])

generate a Signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by creating the
[Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli)
in the tangent space of `p` sampled at `pts` points and use `exp` to get a
signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html).

# Input
* `p` – the tangent space the Lemniscate is created in
* `pts` – (`128`) number of points to sample the Lemniscate
* `a` – (`π/2`) defines a half axis of the Lemniscate to cover a
   half sphere.
* `interval` – (`[0,2*π]`) range to sample the lemniscate at, the default value
  refers to one closed curve

This dataset was used in the numerical example of Section 5.1 of

> Bačák, M., Bergmann, R., Steidl, G., Weinmann, A.:
> _A Second Order Non-Smooth Variational Model for Restoring Manifold-Valued Images_
> SIAM J. Sci. Comput. 38(1), A567–A597, 2016.
> doi: [10.1137/15M101988X](https://dx.doi.org/10.1137/15M101988X)
> arxiv: [1506.02409](https://arxiv.org/abs/1506.02409)

"""
function artificial_S2_lemniscate(
    p, pts::Integer=128, a::Float64=π / 2.0, interval::Array{Float64,1}=[0.0, 2.0 * π]
)
    return artificial_S2_lemniscate.(Ref(p), range(interval[1], interval[2]; length=pts), a)
end
@doc raw"""
    artificial_S2_lemniscate(p,t; a=π/2)

generate a point from the signal on the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html) $\mathbb S^2$ by
creating the [Lemniscate of Bernoulli](https://en.wikipedia.org/wiki/Lemniscate_of_Bernoulli)
in the tangent space of `p` sampled at `t` and use èxp` to obtain a point on
the [Sphere](https://juliamanifolds.github.io/Manifolds.jl/stable/manifolds/sphere.html).

# Input
* `p` – the tangent space the Lemniscate is created in
* `t` – value to sample the Lemniscate at

# Optional Values
 * `a` – (`π/2`) defines a half axis of the Lemniscate to cover a
   half sphere.

This dataset was used in the numerical example of Section 5.1 of

> Bačák, M., Bergmann, R., Steidl, G., Weinmann, A.:
> _A Second Order Non-Smooth Variational Model for Restoring Manifold-Valued Images_
> SIAM J. Sci. Comput. 38(1), A567–A597, 2016.
> doi: [10.1137/15M101988X](https://dx.doi.org/10.1137/15M101988X)
> arxiv: [1506.02409](https://arxiv.org/abs/1506.02409)

 """
function artificial_S2_lemniscate(p, t::Float64, a::Float64=π / 2.0)
    M = Sphere(2)
    tP = 2.0 * Float64(p[1] >= 0.0) - 1.0 # Take north or south pole
    base = [0.0, 0.0, tP]
    xc = a * (cos(t) / (sin(t)^2 + 1.0))
    yc = a * (cos(t) * sin(t) / (sin(t)^2 + 1.0))
    tV = vector_transport_to(M, base, [xc, yc, 0.0], p, ParallelTransport())
    return exp(M, p, tV)
end


text:
# ---
# title: 1356. Sort Integers by The Number of 1 Bits
# id: problem1356
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Sort, Bit Manipulation
# link: <https://leetcode.com/problems/sort-integers-by-the-number-of-1-bits/description/>
# hidden: true
# ---
#
# Given an integer array `arr`. You have to sort the integers in the array in
# ascending order by the number of **1 's** in their binary representation and
# in case of two or more integers have the same number of **1 's** you have to
# sort them in ascending order.
#
# Return _the sorted array_.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [0,1,2,3,4,5,6,7,8]
#     Output: [0,1,2,4,8,3,5,6,7]
#     **Explantion:** [0] is the only integer with 0 bits.
#     [1,2,4,8] all have 1 bit.
#     [3,5,6] have 2 bits.
#     [7] has 3 bits.
#     The sorted array by bits is [0,1,2,4,8,3,5,6,7]
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1024,512,256,128,64,32,16,8,4,2,1]
#     Output: [1,2,4,8,16,32,64,128,256,512,1024]
#     **Explantion:** All integers have 1 bit in the binary representation, you should just sort them in ascending order.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [10000,10000]
#     Output: [10000,10000]
#
#
# **Example 4:**
#
#
#
#     Input: arr = [2,3,5,7,11,13,17,19]
#     Output: [2,3,5,17,7,11,13,19]
#
#
# **Example 5:**
#
#
#
#     Input: arr = [10,100,1000,10000]
#     Output: [10,100,10000,1000]
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 500`
#   * `0 <= arr[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    abstract type MCMCAlgorithm <: AbstractSamplingAlgorithm end

!!! note

    The details of the `MCMCIterator` and `MCMCAlgorithm` API (see below)
    currently do not form part of the BAT-API, and may change without notice.

The following methods must be defined for subtypes (e.g.
for `SomeAlgorithm<:MCMCAlgorithm`):

```julia
(spec::MCMCSpec{<:SomeAlgorithm})(chainid::Integer)::MCMCIterator
```

In some cases, these custom methods may be necessary (default methods are
defined for `MCMCAlgorithm`):

```julia
BAT.initial_params!(
    params::Union{AbstractVector{<:Real},VectorOfSimilarVectors{<:Real}},
    rng::AbstractRNG,
    posterior::AbstractPosteriorDensity,
    algorithm::SomeAlgorithm
)
```

To implement a new MCMC algorithm, subtypes of both `MCMCAlgorithm` and
[`MCMCIterator`](@ref) are required.
"""
abstract type MCMCAlgorithm <: AbstractSamplingAlgorithm end
export MCMCAlgorithm

"""
    BAT.initial_params!(
        params::Union{AbstractVector{<:Real},VectorOfSimilarVectors{<:Real}},
        rng::AbstractRNG,
        posterior::AbstractPosteriorDensity,
        algorithm::MCMCAlgorithm
    )::typeof(params)

*BAT-internal, not part of stable public API.*

Fill `params` with random initial parameters suitable for `posterior` and
`algorithm`. The default implementation will try to draw the initial
parameters from the prior of the posterior.
"""
function initial_params! end

initial_params!(
    params::Union{AbstractVector{<:Real},VectorOfSimilarVectors{<:Real}},
    rng::AbstractRNG,
    posterior::AbstractPosteriorDensity,
    algorithm::MCMCAlgorithm
) = rand!(rng, sampler(getprior(posterior)), params)

@doc doc"""
    MCMCSpec{
        A<:MCMCAlgorithm,
        M<:AbstractPosteriorDensity,
        R<:AbstractRNGSeed
    }

*BAT-internal, not part of stable public API.*

Specifies a Bayesian MCMC chain.

Constructor:

```julia
MCMCSpec(
    algorithm::MCMCAlgorithm,
    posterior::AbstractPosteriorDensity,
    rngseed::AbstractRNGSeed = AbstractRNGSeed()
)
```

Markov-chain instances, represented by objects of type [`MCMCIterator`](@ref),
are be created via

```julia
(spec::MCMCSpec)(chainid::Integer)
```
"""
struct MCMCSpec{
    A<:MCMCAlgorithm,
    M<:AbstractPosteriorDensity,
    R<:AbstractRNGSeed
}
    algorithm::A
    posterior::M
    rngseed::R
end

# TODO/Decision: Make MCMCSpec a subtype of StatsBase.Sampleable{Multivariate,Continuous}?

MCMCSpec(
    algorithm::MCMCAlgorithm,
    posterior::AbstractPosteriorDensity
) = MCMCSpec(algorithm, posterior, AbstractRNGSeed())

@deprecate MCMCSpec(
    algorithm::MCMCAlgorithm,
    likelihood::Union{AbstractDensity,Distribution{Multivariate,Continuous}},
    prior::Union{AbstractDensity,ParamVolumeBounds},
    rngseed::AbstractRNGSeed = AbstractRNGSeed()
) MCMCSpec(algorithm, PosteriorDensity(likelihood, prior), rngseed)

@with_kw struct MCMCIteratorInfo
    id::Int64
    cycle::Int
    tuned::Bool
    converged::Bool
end

@doc doc"""
    abstract type MCMCIterator end

Represents the current state of a MCMC chain.

!!! note

    The details of the `MCMCIterator` and `MCMCAlgorithm` API (see below)
    currently do not form part of the BAT-API, and may change without notice.

To implement a new MCMC algorithm, subtypes of both [`MCMCAlgorithm`](@ref)
and `MCMCIterator` are required.

The following methods must be defined for subtypes of `MCMCIterator` (e.g.
`SomeMCMCIter<:MCMCIterator`):

```julia
BAT.mcmc_spec(chain::SomeMCMCIter)::MCMCSpec

BAT.getrng(chain::SomeMCMCIter)::AbstractRNG

BAT.mcmc_info(chain::SomeMCMCIter)::MCMCIteratorInfo

BAT.nsteps(chain::SomeMCMCIter)::Int

BAT.nsamples(chain::SomeMCMCIter)::Int

BAT.current_sample(chain::SomeMCMCIter)::DensitySample

BAT.sample_type(chain::SomeMCMCIter)::Type{<:DensitySample}

BAT.samples_available(chain::SomeMCMCIter, nonzero_weights::Bool = false)::Bool

BAT.get_samples!(appendable, chain::SomeMCMCIter, nonzero_weights::Bool)::typeof(appendable)

BAT.next_cycle!(chain::SomeMCMCIter)::SomeMCMCIter

BAT.mcmc_step!(
    callback::AbstractMCMCCallback,
    chain::SomeMCMCIter
)
```

The following methods are implemented by default:

```julia
algorithm(chain::MCMCIterator)
getposterior(chain::MCMCIterator)
rngseed(chain::MCMCIterator)
nparams(chain::MCMCIterator)
DensitySampleVector(chain::MCMCIterator)
mcmc_iterate!(callback, chain::MCMCIterator, ...)
mcmc_iterate!(callbacks, chains::AbstractVector{<:MCMCIterator}, ...)
```
"""
abstract type MCMCIterator end
export MCMCIterator

function mcmc_spec end

function getrng end

function mcmc_info end

function nsteps end

function nsamples end

function current_sample end

function sample_type end

function samples_available end

function get_samples! end

function next_cycle! end

function mcmc_step! end

algorithm(chain::MCMCIterator) = mcmc_spec(chain).algorithm

getposterior(chain::MCMCIterator) = mcmc_spec(chain).posterior

rngseed(chain::MCMCIterator) = mcmc_spec(chain).rngseed

nparams(chain::MCMCIterator) = nparams(getposterior(chain))

DensitySampleVector(chain::MCMCIterator) = DensitySampleVector(sample_type(chain), nparams(chain))

function mcmc_iterate! end

function mcmc_iterate!(
    callback,
    chain::MCMCIterator;
    max_nsamples::Int64 = Int64(1),
    max_nsteps::Int64 = Int64(1000),
    max_time::Float64 = Inf
)
    @debug "Starting iteration over MCMC chain $(chain.info.id), max_nsamples = $max_nsamples, max_nsteps = $max_nsteps, max_time = $max_time"

    cbfunc = Base.convert(AbstractMCMCCallback, callback)

    start_time = time()
    start_nsteps = nsteps(chain)
    start_nsamples = nsamples(chain)

    while (
        (nsamples(chain) - start_nsamples) < max_nsamples &&
        (nsteps(chain) - start_nsteps) < max_nsteps &&
        (time() - start_time) < max_time
    )
        mcmc_step!(cbfunc, chain)
    end

    end_time = time()
    elapsed_time = end_time - start_time

    @debug "Finished iteration over MCMC chain $(chain.info.id), nsamples = $(nsamples(chain)), nsteps = $(nsteps(chain)), time = $(Float32(elapsed_time))"

    chain
end

function mcmc_iterate!(
    callbacks,
    chains::AbstractVector{<:MCMCIterator};
    kwargs...
)
    if isempty(chains)
        @debug "No MCMC chain(s) to iterate over."
        return chains
    else
        @debug "Starting iteration over $(length(chains)) MCMC chain(s)"
    end

    cbv = mcmc_callback_vector(callbacks, eachindex(chains))

    idxs = eachindex(cbv, chains)
    @sync for i in idxs
        @mt_async mcmc_iterate!(cbv[i], chains[i]; kwargs...)
    end

    chains
end