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

export AbstractTopology,
    BrickTopology,
    StackedBrickTopology,
    CubedShellTopology,
    StackedCubedSphereTopology,
    AnalyticalTopography,
    NoTopography,
    DCMIPMountain,
    EquidistantCubedSphere,
    EquiangularCubedSphere,
    isstacked,
    cubed_sphere_topo_warp,
    compute_lat_long,
    compute_analytical_topography,
    cubed_sphere_warp,
    cubed_sphere_unwarp,
    equiangular_cubed_sphere_warp,
    equiangular_cubed_sphere_unwarp,
    equidistant_cubed_sphere_warp,
    equidistant_cubed_sphere_unwarp

export grid1d, SingleExponentialStretching, InteriorStretching

export basic_topology_info

"""
    AbstractTopology{dim, T, nb}

Represents the connectivity of individual elements, with local dimension `dim`
with `nb` boundary conditions types. The element coordinates are of type `T`.
"""
abstract type AbstractTopology{dim, T, nb} end

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

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

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, nb} <: AbstractTopology{dim, T, nb}

    """
    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 it is a boundary face,
    then it is the boundary element index.
    """
    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}

    """
    Tuple of boundary to element. `bndytoelem[b][i]` is the element which faces
    the `i`th boundary element of boundary `b`.
    """
    bndytoelem::NTuple{nb, Array{Int64, 1}}

    """
    Tuple of boundary to element face. `bndytoface[b][i]` is the face number of
    the element which faces the `i`th boundary element of boundary `b`.
    """
    bndytoface::NTuple{nb, Array{Int64, 1}}

    """
    Element to locally unique vertex number; `elemtouvert[v,e]` is the `v`th vertex
    of element `e`
    """
    elemtouvert::Union{Array{Int64, 2}, Nothing}

    """
    Vertex connectivity information for direct stiffness summation;
    `vtconn[vtconnoff[i]:vtconnoff[i+1]-1] == vtconn[lvt1, lelem1, lvt2, lelem2, ....]`
    for each rank-local unique vertex number,
    where `lvt1` is the element-local vertex number of element `lelem1`, etc.
    Vertices, not shared by multiple elements, are ignored.
    """
    vtconn::Union{AbstractArray{Int64, 1}, Nothing}

    """
    Vertex offset for vertex connectivity information
    """
    vtconnoff::Union{AbstractArray{Int64, 1}, Nothing}

    """
    Face connectivity information for direct stiffness summation;
    `fcconn[ufc,:] = [lfc1, lelem1, lfc2, lelem2, ordr]`
    where `ufc` is the rank-local unique face number, and
    `lfc1` is the element-local face number of element `lelem1`, etc.,
    and ordr is the relative orientation. Faces, not shared by multiple
    elements are ignored.
    """
    fcconn::Union{AbstractArray{Int64, 2}, Nothing}

    """
    Edge connectivity information for direct stiffness summation;
    `edgconn[edgconnoff[i]:edgconnoff[i+1]-1] == edgconn[ledg1, orient, lelem1, ledg2, orient, lelem2, ...]`
    for each rank-local unique edge number,
    where `ledg1` is the element-local edge number, `orient1` is orientation (forward/reverse) of dof along edge.
    Edges, not shared by multiple elements, are ignored.
    """
    edgconn::Union{AbstractArray{Int64, 1}, Nothing}

    """
    Edge offset for edge connectivity information
    """
    edgconnoff::Union{AbstractArray{Int64, 1}, Nothing}

    function BoxElementTopology{dim, T, nb}(
        mpicomm,
        elems,
        realelems,
        ghostelems,
        ghostfaces,
        sendelems,
        sendfaces,
        elemtocoord,
        elemtoelem,
        elemtoface,
        elemtoordr,
        elemtobndy,
        nabrtorank,
        nabrtorecv,
        nabrtosend,
        origsendorder,
        bndytoelem,
        bndytoface;
        device_array = ClimateMachine.array_type(),
        elemtouvert = nothing,
        vtconn = nothing,
        vtconnoff = nothing,
        fcconn = nothing,
        edgconn = nothing,
        edgconnoff = nothing,
    ) where {dim, T, nb}

        exteriorelems = sort(unique(sendelems))
        interiorelems = sort(setdiff(realelems, exteriorelems))
        if vtconn isa Array && vtconnoff isa Array
            vtconn = device_array(vtconn)
            vtconnoff = device_array(vtconnoff)
        end
        if fcconn isa Array
            fcconn = device_array(fcconn)
        end
        if edgconn isa Array && edgconnoff isa Array
            edgconn = device_array(edgconn)
            edgconnoff = device_array(edgconnoff)
        end
        return new{dim, T, nb}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            interiorelems,
            exteriorelems,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            origsendorder,
            bndytoelem,
            bndytoface,
            elemtouvert,
            vtconn,
            vtconnoff,
            fcconn,
            edgconn,
            edgconnoff,
        )
    end
end

"""
    hasboundary(topology::AbstractTopology)

query function to check whether a topology has a boundary (i.e., not fully
periodic)
"""
hasboundary(topology::AbstractTopology{dim, T, nb}) where {dim, T, nb} = nb != 0

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, nb} <: AbstractTopology{dim, T, nb}

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

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

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

abstract type AbstractStackedTopology{dim, T, nb} <:
              AbstractTopology{dim, T, nb} end

"""
    StackedBrickTopology{dim, T, nb} <: AbstractStackedTopology{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, nb} <: AbstractStackedTopology{dim, T, nb}
    topology::BoxElementTopology{dim, T, nb}
    stacksize::Int64
    periodicstack::Bool
end
function Base.getproperty(a::StackedBrickTopology, p::Symbol)
    return p in (:stacksize, :periodicstack) ? getfield(a, p) :
           getproperty(getfield(a, :topology), p)
end

"""
    StackedCubedSphereTopology{T, nb} <: AbstractStackedTopology{3, T, nb}

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, nb} <: AbstractStackedTopology{3, T, nb}
    topology::BoxElementTopology{3, T, nb}
    stacksize::Int64
end
function Base.getproperty(a::StackedCubedSphereTopology, p::Symbol)
    if p == :periodicstack
        return false
    else
        return p == :stacksize ? getfield(a, p) :
               getproperty(getfield(a, :topology), p)
    end
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 ClimateMachine.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 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 =
        connectivity == :face ?
        BrickMesh.connectmesh(mpicomm, topology[1:4]...) :
        BrickMesh.connectmeshfull(mpicomm, topology[1:4]...)
    bndytoelem, bndytoface = BrickMesh.enumerateboundaryfaces!(
        topology.elemtoelem,
        topology.elemtobndy,
        periodicity,
        boundary,
    )

    nb = length(bndytoelem)
    dim = length(elemrange)
    T = eltype(topology.elemtocoord)
    return BrickTopology{dim, T, nb}(BoxElementTopology{dim, T, nb}(
        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,
        bndytoelem,
        bndytoface,
        elemtouvert = topology.elemtovert,
    ))
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 ClimateMachine.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 = :full,
    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)
    ]

    bndytoelem, bndytoface = BrickMesh.enumerateboundaryfaces!(
        elemtoelem,
        elemtobndy,
        periodicity,
        boundary,
    )
    nb = length(bndytoelem)

    T = eltype(basetopo.elemtocoord)
    #----setting up DSS----------
    if basetopo.elemtouvert isa Array
        #---setup vertex DSS
        nvertby2 = div(nvert, 2)
        elemtouvert = similar(basetopo.elemtouvert, nvert, length(elems))
        # base level
        for i in 1:length(basetopo.elems)
            e = stacksize * (i - 1) + 1
            for vt in 1:nvertby2
                elemtouvert[vt, e] =
                    (basetopo.elemtouvert[vt, i] - 1) * (stacksize + 1) + 1
            end
        end
        # interior levels
        for i in 1:length(basetopo.elems), j in 1:(stacksize - 1)
            e = stacksize * (i - 1) + j
            for vt in 1:nvertby2
                elemtouvert[nvertby2 + vt, e] = elemtouvert[vt, e] + 1
                elemtouvert[vt, e + 1] = elemtouvert[nvertby2 + vt, e]
            end
        end
        # top level
        if periodicity[dim]
            for i in 1:length(basetopo.elems)
                e = stacksize * i
                for vt in 1:nvertby2
                    elemtouvert[nvertby2 + vt, e] =
                        (basetopo.elemtouvert[vt, i] - 1) * (stacksize + 1) + 1
                end
            end
        else
            for i in 1:length(basetopo.elems)
                e = stacksize * i
                for vt in 1:nvertby2
                    elemtouvert[nvertby2 + vt, e] = elemtouvert[vt, e] + 1
                end
            end
        end
        vtmax = maximum(unique(elemtouvert))
        vtconn = map(j -> zeros(Int, j), zeros(Int, vtmax))
        for el in elems, lvt in 1:nvert
            vt = elemtouvert[lvt, el]
            push!(vtconn[vt], lvt) # local vertex number
            push!(vtconn[vt], el)  # local elem number
        end
        # building the vertex connectivity device array
        vtconnoff = zeros(Int, vtmax + 1)
        vtconnoff[1] = 1
        temp = Int64[]
        for vt in 1:vtmax
            nconn = length(vtconn[vt])
            if nconn > 2
                vtconnoff[vt + 1] = vtconnoff[vt] + nconn
                append!(temp, vtconn[vt])
            else
                vtconnoff[vt + 1] = vtconnoff[vt]
            end
        end
        vtconn = temp
        #---setup face DSS---------------------
        fcmarker = -ones(Int, nface, length(elems))
        fcno = 0
        for el in realelems
            for fc in 1:nface
                if elemtobndy[fc, el] == 0
                    nabrel = elemtoelem[fc, el]
                    nabrfc = elemtoface[fc, el]
                    if fcmarker[fc, el] == -1 && fcmarker[nabrfc, nabrel] == -1
                        fcno += 1
                        fcmarker[fc, el] = fcno
                        fcmarker[nabrfc, nabrel] = fcno
                    end
                end
            end
        end
        fcconn = -ones(Int, fcno, 5)
        for el in realelems
            for fc in 1:nface
                fcno = fcmarker[fc, el]
                ordr = elemtoordr[fc, el]
                if fcno ≠ -1
                    nabrel = elemtoelem[fc, el]
                    nabrfc = elemtoface[fc, el]
                    fcconn[fcno, 1] = fc
                    fcconn[fcno, 2] = el
                    fcconn[fcno, 3] = nabrfc
                    fcconn[fcno, 4] = nabrel
                    fcconn[fcno, 5] = ordr
                    fcmarker[fc, el] = -1
                    fcmarker[nabrfc, nabrel] = -1
                end
            end
        end
        #---setup edge DSS---------------------
        if dim == 3
            nedge = 12
            edgemask = [
                1 3 5 7 1 2 5 6 1 2 3 4
                2 4 6 8 3 4 7 8 5 6 7 8
            ]
            edges = Array{Int64}(undef, 6, nedge * length(elems))
            ledge = zeros(Int64, 2)
            uedge = Dict{Array{Int64, 1}, Int64}()
            orient = 1
            uedgno = Int64(0)
            ctr = 1
            for el in elems
                for edg in 1:nedge
                    orient = 1
                    for i in 1:2
                        ledge[i] = elemtouvert[edgemask[i, edg], el]
                        edges[i, ctr] = ledge[i]                     # edge vertices [1:2]
                    end
                    if ledge[1] > ledge[2]
                        sort!(ledge)
                        orient = 2
                    end
                    edges[3, ctr] = edg    # local edge number
                    edges[4, ctr] = orient # edge orientation
                    edges[5, ctr] = el     # element number
                    if haskey(uedge, ledge[:])
                        edges[6, ctr] = uedge[ledge[:]]  # unique edge number
                    else
                        uedgno += 1
                        uedge[ledge[:]] = uedgno
                        edges[6, ctr] = uedgno
                    end
                    ctr += 1
                end
            end
            edgconn = map(j -> zeros(Int, j), zeros(Int, uedgno))
            ctr = 1
            for el in elems
                for edg in 1:nedge
                    uedg = edges[6, ctr]
                    push!(edgconn[uedg], edg)           # local edge number
                    push!(edgconn[uedg], edges[4, ctr]) # orientation
                    push!(edgconn[uedg], el)            # element number
                    ctr += 1
                end
            end
            # remove edges belonging to a single element and edges
            # belonging exclusively to ghost elements
            # and build edge connectivity device array
            edgconnoff = Int64[]
            temp = Int64[]
            push!(edgconnoff, 1)
            for i in 1:uedgno
                elen = length(edgconn[i])
                encls = Int(elen / 3)
                edgmrk = true
                if encls > 1
                    for j in 1:encls
                        if edgconn[i][3 * j] ≤ nreal
                            edgmrk = false
                        end
                    end
                end
                if edgmrk
                    edgconn[i] = []
                else
                    shift = edgconnoff[end]
                    push!(edgconnoff, (shift + length(edgconn[i])))
                    append!(temp, edgconn[i][:])
                end
            end
            edgconn = temp
        else
            edgconn = nothing
            edgconnoff = nothing
        end
    else
        elemtouvert = nothing
        vtconn = nothing
        vtconnoff = nothing
        fcconn = nothing
        edgconn = nothing
        edgconnoff = nothing
    end
    #---------------
    StackedBrickTopology{dim, T, nb}(
        BoxElementTopology{dim, T, nb}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            basetopo.origsendorder,
            bndytoelem,
            bndytoface,
            elemtouvert = elemtouvert,
            vtconn = vtconn,
            vtconnoff = vtconnoff,
            fcconn = fcconn,
            edgconn = edgconn,
            edgconnoff = edgconnoff,
        ),
        stacksize,
        periodicity[end],
    )
end

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

Generate a cubed shell mesh with the number of elements along each dimension of
the cubes being `Nelem`. This topology actually creates a cube mesh, and the
warping should be done after the grid is created using the `cubed_sphere_warp`
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*10 elements on each cube face, i.e., the
total number of elements is
`10 * 10 * 6 = 600`, with
```jldoctest brickmesh
using ClimateMachine.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.cubed_sphere_warp(EquiangularCubedSphere(), x1[n], x2[n], x3[n])
end

# in case a unitary equiangular cubed sphere is desired, or

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

# in case an equidistant cubed sphere of radius 10 is desired.
```
"""
function CubedShellTopology(
    mpicomm,
    Neside,
    T;
    connectivity = :full,
    ghostsize = 1,
)

    # We cannot handle anything else right now...
    @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 =
        connectivity == :face ?
        BrickMesh.connectmesh(
            mpicomm,
            topology[1],
            elemtocoord,
            topology[3],
            topology[4],
            dim = 2,
        ) :
        BrickMesh.connectmeshfull(
            mpicomm,
            topology[1],
            elemtocoord,
            topology[3],
            topology[4],
            dim = 2,
        )

    CubedShellTopology{T}(BoxElementTopology{2, T, 0}(
        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,
        (),
        (),
        elemtouvert = topology.elemtovert,
    ))
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

abstract type AbstractCubedSphere end
struct EquiangularCubedSphere <: AbstractCubedSphere end
struct EquidistantCubedSphere <: AbstractCubedSphere end

"""
    cubed_sphere_warp(::EquiangularCubedSphere, 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
[Ronchi1996](@cite)
"""
function cubed_sphere_warp(
    ::EquiangularCubedSphere,
    a,
    b,
    c,
    R = max(abs(a), abs(b), abs(c)),
)

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

    fdim = argmax(abs.((a, b, c)))
    if fdim == 1 && a < 0
        # (-R, *, *) : formulas for Face I from Ronchi, Iacono, Paolucci (1996)
        #              but for us face II of the developed net of the cube
        x1, x2, x3 = f(-R, b / a, c / a)
    elseif fdim == 2 && b < 0
        # ( *,-R, *) : formulas for Face II from Ronchi, Iacono, Paolucci (1996)
        #              but for us face III of the developed net of the cube
        x2, x1, x3 = f(-R, a / b, c / b)
    elseif fdim == 1 && a > 0
        # ( R, *, *) : formulas for Face III from Ronchi, Iacono, Paolucci (1996)
        #              but for us face IV of the developed net of the cube
        x1, x2, x3 = f(R, b / a, c / a)
    elseif fdim == 2 && b > 0
        # ( *, R, *) : formulas for Face IV from Ronchi, Iacono, Paolucci (1996)
        #              but for us face I of the developed net of the cube
        x2, x1, x3 = f(R, a / b, c / b)
    elseif fdim == 3 && c > 0
        # ( *, *, R) : formulas for Face V from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        x3, x2, x1 = f(R, b / c, a / c)
    elseif fdim == 3 && c < 0
        # ( *, *,-R) : formulas for Face VI from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        x3, x2, x1 = f(-R, b / c, a / c)
    else
        error("invalid case for cubed_sphere_warp(::EquiangularCubedSphere): $a, $b, $c")
    end

    return x1, x2, x3
end

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

A wrapper function for the cubed_sphere_warp function, when called with the
EquiangularCubedSphere type
"""
equiangular_cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c))) =
    cubed_sphere_warp(EquiangularCubedSphere(), a, b, c, R)

"""
    cubed_sphere_unwarp(x1, x2, x3)

The inverse of [`cubed_sphere_warp`](@ref). This function projects
a given point `(x_1, x_2, x_3)` from the surface of a sphere onto a cube
"""
function cubed_sphere_unwarp(::EquiangularCubedSphere, 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, *, *) : formulas for Face I from Ronchi, Iacono, Paolucci (1996)
        #              but for us face II of the developed net of the cube
        a, b, c = g(-R, x2 / x1, x3 / x1)
    elseif fdim == 2 && x2 < 0
        # ( *,-R, *) : formulas for Face II from Ronchi, Iacono, Paolucci (1996)
        #              but for us face III of the developed net of the cube
        b, a, c = g(-R, x1 / x2, x3 / x2)
    elseif fdim == 1 && x1 > 0
        # ( R, *, *) : formulas for Face III from Ronchi, Iacono, Paolucci (1996)
        #              but for us face IV of the developed net of the cube
        a, b, c = g(R, x2 / x1, x3 / x1)
    elseif fdim == 2 && x2 > 0
        # ( *, R, *) : formulas for Face IV from Ronchi, Iacono, Paolucci (1996)
        #              but for us face I of the developed net of the cube
        b, a, c = g(R, x1 / x2, x3 / x2)
    elseif fdim == 3 && x3 > 0
        # ( *, *, R) : formulas for Face V from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        c, b, a = g(R, x2 / x3, x1 / x3)
    elseif fdim == 3 && x3 < 0
        # ( *, *,-R) : formulas for Face VI from Ronchi, Iacono, Paolucci (1996)
        #              and the same for us on the developed net of the cube
        c, b, a = g(-R, x2 / x3, x1 / x3)
    else
        error("invalid case for cubed_sphere_unwarp(::EquiangularCubedSphere): $a, $b, $c")
    end

    return a, b, c
end

"""
    equiangular_cubed_sphere_unwarp(x1, x2, x3)

A wrapper function for the cubed_sphere_unwarp function, when called with the
EquiangularCubedSphere type
"""
equiangular_cubed_sphere_unwarp(x1, x2, x3) =
    cubed_sphere_unwarp(EquiangularCubedSphere(), x1, x2, x3)

"""
    cubed_sphere_warp(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 equidistant gnomonic grid outlined in
[Rancic1996](@cite) and [Nair2005](@cite)
"""
function cubed_sphere_warp(
    ::EquidistantCubedSphere,
    a,
    b,
    c,
    R = max(abs(a), abs(b), abs(c)),
)

    r = hypot(a, b, c)

    x1 = R * a / r
    x2 = R * b / r
    x3 = R * c / r

    return x1, x2, x3
end

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

A wrapper function for the cubed_sphere_warp function, when called with the
EquidistantCubedSphere type
"""
equidistant_cubed_sphere_warp(a, b, c, R = max(abs(a), abs(b), abs(c))) =
    cubed_sphere_warp(EquidistantCubedSphere(), a, b, c, R)

"""
    cubed_sphere_unwarp(x1, x2, x3)

The inverse of [`cubed_sphere_warp`](@ref). This function projects
a given point `(x_1, x_2, x_3)` from the surface of a sphere onto a cube
"""
function cubed_sphere_unwarp(::EquidistantCubedSphere, x1, x2, x3)

    r = hypot(1, x2 / x1, x3 / x1)
    R = hypot(x1, x2, x3)

    a = r * x1
    b = r * x2
    c = r * x3

    m = max(abs(a), abs(b), abs(c))

    return a * R / m, b * R / m, c * R / m
end

"""
    equidistant_cubed_sphere_unwarp(x1, x2, x3)

A wrapper function for the cubed_sphere_unwarp function, when called with the
EquidistantCubedSphere type
"""
equidistant_cubed_sphere_unwarp(x1, x2, x3) =
    cubed_sphere_unwarp(EquidistantCubedSphere(), x1, x2, x3)

"""
   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 `cubed_sphere_warp` 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 face,
i.e., the total number of elements is `10 * 10 * 5 * 6 = 3000`, with
```jldoctest brickmesh
using ClimateMachine.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.cubed_sphere_warp(EquiangularCubedSphere(),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 = :full,
    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)
    ]

    bndytoelem, bndytoface = BrickMesh.enumerateboundaryfaces!(
        elemtoelem,
        elemtobndy,
        (false,),
        (boundary,),
    )
    nb = length(bndytoelem)
    #----setting up DSS----------
    if basetopo.elemtouvert isa Array
        #---setup vertex DSS
        nvertby2 = div(nvert, 2)
        elemtouvert = similar(basetopo.elemtouvert, nvert, length(elems))
        # base level
        for i in 1:length(basetopo.elems)
            e = stacksize * (i - 1) + 1
            for vt in 1:nvertby2
                elemtouvert[vt, e] =
                    (basetopo.elemtouvert[vt, i] - 1) * (stacksize + 1) + 1
            end
        end
        # interior levels
        for i in 1:length(basetopo.elems), j in 1:(stacksize - 1)
            e = stacksize * (i - 1) + j
            for vt in 1:nvertby2
                elemtouvert[nvertby2 + vt, e] = elemtouvert[vt, e] + 1
                elemtouvert[vt, e + 1] = elemtouvert[nvertby2 + vt, e]
            end
        end
        # top level
        for i in 1:length(basetopo.elems)
            e = stacksize * i
            for vt in 1:nvertby2
                elemtouvert[nvertby2 + vt, e] = elemtouvert[vt, e] + 1
            end
        end
        vtmax = maximum(unique(elemtouvert))
        vtconn = map(j -> zeros(Int, j), zeros(Int, vtmax))

        for el in elems, lvt in 1:nvert
            vt = elemtouvert[lvt, el]
            push!(vtconn[vt], lvt) # local vertex number
            push!(vtconn[vt], el)  # local elem number
        end
        # building the vertex connectivity device array
        vtconnoff = zeros(Int, vtmax + 1)
        vtconnoff[1] = 1
        temp = Int64[]
        for vt in 1:vtmax
            nconn = length(vtconn[vt])
            if nconn > 2
                vtconnoff[vt + 1] = vtconnoff[vt] + nconn
                append!(temp, vtconn[vt])
            else
                vtconnoff[vt + 1] = vtconnoff[vt]
            end
        end
        vtconn = temp
        #---setup face DSS---------------------
        fcmarker = -ones(Int, nface, length(elems))
        fcno = 0
        for el in realelems
            for fc in 1:nface
                if elemtobndy[fc, el] == 0
                    nabrel = elemtoelem[fc, el]
                    nabrfc = elemtoface[fc, el]
                    if fcmarker[fc, el] == -1 && fcmarker[nabrfc, nabrel] == -1
                        fcno += 1
                        fcmarker[fc, el] = fcno
                        fcmarker[nabrfc, nabrel] = fcno
                    end
                end
            end
        end
        fcconn = -ones(Int, fcno, 5)
        for el in realelems
            for fc in 1:nface
                fcno = fcmarker[fc, el]
                ordr = elemtoordr[fc, el]
                if fcno ≠ -1
                    nabrel = elemtoelem[fc, el]
                    nabrfc = elemtoface[fc, el]
                    fcconn[fcno, 1] = fc
                    fcconn[fcno, 2] = el
                    fcconn[fcno, 3] = nabrfc
                    fcconn[fcno, 4] = nabrel
                    fcconn[fcno, 5] = ordr
                    fcmarker[fc, el] = -1
                    fcmarker[nabrfc, nabrel] = -1
                end
            end
        end
        #---setup edge DSS---------------------
        if dim == 3
            nedge = 12
            edgemask = [
                1 3 5 7 1 2 5 6 1 2 3 4
                2 4 6 8 3 4 7 8 5 6 7 8
            ]
            edges = Array{Int64}(undef, 6, nedge * length(elems))
            ledge = zeros(Int64, 2)
            uedge = Dict{Array{Int64, 1}, Int64}()
            orient = 1
            uedgno = Int64(0)
            ctr = 1
            for el in elems
                for edg in 1:nedge
                    orient = 1
                    for i in 1:2
                        ledge[i] = elemtouvert[edgemask[i, edg], el]
                        edges[i, ctr] = ledge[i]                     # edge vertices [1:2]
                    end
                    if ledge[1] > ledge[2]
                        sort!(ledge)
                        orient = 2
                    end
                    edges[3, ctr] = edg    # local edge number
                    edges[4, ctr] = orient # edge orientation
                    edges[5, ctr] = el     # element number
                    if haskey(uedge, ledge[:])
                        edges[6, ctr] = uedge[ledge[:]]  # unique edge number
                    else
                        uedgno += 1
                        uedge[ledge[:]] = uedgno
                        edges[6, ctr] = uedgno
                    end
                    ctr += 1
                end
            end
            edgconn = map(j -> zeros(Int, j), zeros(Int, uedgno))
            ctr = 1
            for el in elems
                for edg in 1:nedge
                    uedg = edges[6, ctr]
                    push!(edgconn[uedg], edg)           # local edge number
                    push!(edgconn[uedg], edges[4, ctr]) # orientation
                    push!(edgconn[uedg], el)            # element number
                    ctr += 1
                end
            end
            # remove edges belonging to a single element and edges
            # belonging exclusively to ghost elements
            # and build edge connectivity device array
            edgconnoff = Int64[]
            temp = Int64[]
            push!(edgconnoff, 1)
            for i in 1:uedgno
                elen = length(edgconn[i])
                encls = Int(elen / 3)
                edgmrk = true
                if encls > 1
                    for j in 1:encls
                        if edgconn[i][3 * j] ≤ nreal
                            edgmrk = false
                        end
                    end
                end
                if edgmrk
                    edgconn[i] = []
                else
                    shift = edgconnoff[end]
                    push!(edgconnoff, (shift + length(edgconn[i])))
                    append!(temp, edgconn[i][:])
                end
            end
            edgconn = temp
        else
            edgconn = nothing
            edgconnoff = nothing
        end
        #-------------------------------------
    else
        elemtouvert = nothing
        vtconn = nothing
        vtconnoff = nothing
        fcconn = nothing
        edgconn = nothing
        edgconnoff = nothing
    end
    #-----------------------------
    StackedCubedSphereTopology{T, nb}(
        BoxElementTopology{3, T, nb}(
            mpicomm,
            elems,
            realelems,
            ghostelems,
            ghostfaces,
            sendelems,
            sendfaces,
            elemtocoord,
            elemtoelem,
            elemtoface,
            elemtoordr,
            elemtobndy,
            nabrtorank,
            nabrtorecv,
            nabrtosend,
            basetopo.origsendorder,
            bndytoelem,
            bndytoface,
            elemtouvert = elemtouvert,
            vtconn = vtconn,
            vtconnoff = vtconnoff,
            fcconn = fcconn,
            edgconn = edgconn,
            edgconnoff = edgconnoff,
        ),
        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

function basic_topology_info(topology::AbstractStackedTopology)
    nelem = length(topology.elems)
    nvertelem = topology.stacksize
    nhorzelem = div(nelem, nvertelem)
    nrealelem = length(topology.realelems)
    nhorzrealelem = div(nrealelem, nvertelem)

    return (
        nelem = nelem,
        nvertelem = nvertelem,
        nhorzelem = nhorzelem,
        nrealelem = nrealelem,
        nhorzrealelem = nhorzrealelem,
    )
end

function basic_topology_info(topology::AbstractTopology)
    return (
        nelem = length(topology.elems),
        nrealelem = length(topology.realelems),
    )
end

### Helper Functions for Topography Calculations
"""
    compute_lat_long(X,Y,δ,faceid)
Helper function to allow computation of latitute and longitude coordinates
given the cubed sphere coordinates X, Y, δ, faceid
"""
function compute_lat_long(X, Y, δ, faceid)
    if faceid == 1
        λ = atan(X)                     # longitude
        ϕ = atan(cos(λ) * Y)            # latitude
    elseif faceid == 2
        λ = atan(X) + π / 2
        ϕ = atan(Y * cos(atan(X)))
    elseif faceid == 3
        λ = atan(X) + π
        ϕ = atan(Y * cos(atan(X)))
    elseif faceid == 4
        λ = atan(X) + (3 / 2) * π
        ϕ = atan(Y * cos(atan(X)))
    elseif faceid == 5
        λ = atan(X, -Y) + π
        ϕ = atan(1 / sqrt(δ - 1))
    elseif faceid == 6
        λ = atan(X, Y)
        ϕ = -atan(1 / sqrt(δ - 1))
    end
    return λ, ϕ
end

"""
    AnalyticalTopography
Abstract type to allow dispatch over different analytical topography prescriptions
in experiments.
"""
abstract type AnalyticalTopography end

function compute_analytical_topography(
    ::AnalyticalTopography,
    λ,
    ϕ,
    sR,
    r_inner,
    r_outer,
)
    return sR
end

"""
    NoTopography <: AnalyticalTopography
Allows definition of fallback methods in case cubed_sphere_topo_warp is used with
no prescribed topography function.
"""
struct NoTopography <: AnalyticalTopography end

### DCMIP Mountain
"""
    DCMIPMountain <: AnalyticalTopography
Topography description based on standard DCMIP experiments.
"""
struct DCMIPMountain <: AnalyticalTopography end
function compute_analytical_topography(
    ::DCMIPMountain,
    λ,
    ϕ,
    sR,
    r_inner,
    r_outer,
)
    #User specified warp parameters
    R_m = π * 3 / 4
    h0 = 2000
    ζ_m = π / 16
    φ_m = 0
    λ_m = π * 3 / 2
    r_m = acos(sin(φ_m) * sin(ϕ) + cos(φ_m) * cos(ϕ) * cos(λ - λ_m))
    # Define mesh decay profile
    Δ = (r_outer - abs(sR)) / (r_outer - r_inner)
    if r_m < R_m
        zs =
            0.5 *
            h0 *
            (1 + cos(π * r_m / R_m)) *
            cos(π * r_m / ζ_m) *
            cos(π * r_m / ζ_m)
    else
        zs = 0.0
    end
    mR = sign(sR) * (abs(sR) + zs * Δ)
    return mR
end

"""
    cubed_sphere_topo_warp(a, b, c, R = max(abs(a), abs(b), abs(c));
                       r_inner = _planet_radius,
                       r_outer = _planet_radius + domain_height,
                       topography = NoTopography())

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
[Ronchi1996](@cite). Assumes a user specified modified radius using the
compute_analytical_topography function. Defaults to smooth cubed sphere unless otherwise specified
via the AnalyticalTopography type.
"""
function cubed_sphere_topo_warp(
    a,
    b,
    c,
    R = max(abs(a), abs(b), abs(c));
    r_inner = _planet_radius,
    r_outer = _planet_radius + domain_height,
    topography::AnalyticalTopography = NoTopography(),
)

    function f(sR, ξ, η, faceid)
        X, Y = tan(π * ξ / 4), tan(π * η / 4)
        δ = 1 + X^2 + Y^2
        λ, ϕ = compute_lat_long(X, Y, δ, faceid)
        mR = compute_analytical_topography(
            topography,
            λ,
            ϕ,
            sR,
            r_inner,
            r_outer,
        )
        x1 = mR / sqrt(δ)
        x2, x3 = X * x1, Y * x1
        x1, x2, x3
    end
    fdim = argmax(abs.((a, b, c)))

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

    return x1, x2, x3
end

end


text:
# ## Literate example
#
# Whatever is in comment can be Franklin-markdown, so for instance you can
# write maths in much the same way:
#
# $$ \exp(i\tau) = 1 $$
#
# This is Euler's identity following the [$\tau$ day manifesto](https://tauday.com/tau-manifesto).
#
# Code:

a = 5
a^2

# Any valid code basically:

v = ["Thom", "Martha", "Smith", "John"]
f = first.(v)

# So the workflow is to just write a standard Literate document using
# Franklin markdown in the comments and normal code.
# You can also hide lines:

using Random # hide
Random.seed!(1) # hide
randn()

# Or hide everything (note that for the line to be taken as comment for the
# code block there should be no whitespace between the `#` and the `h`).

#hideall
using LinearAlgebra
a = dot(randn(5), randn(5))

# It still shows the result though since `showall = true`! If you want to
# suppress it, just use a `;` on the last line:

#hideall
using LinearAlgebra
a = dot(randn(5), randn(5));

# That's basically it.

println("done!")

ext: jl


text:
# include(joinpath(@__DIR__, "..", "src", "qdldl.jl"))
using Altro.Cqdldl
# using Main.Cqdldl
using SparseArrays
using LinearAlgebra
using Test
using Random

function testallocs(solver::Cqdldl.QDLDLSolver, U, b)
    A = Symmetric(U, :U)
    allocs = @allocated resize!(solver, size(U,1))
    allocs += @allocated Cqdldl.eliminationtree!(solver, U)
    allocs += @allocated Cqdldl.factor!(solver)
    x = copy(b)
    allocs += @allocated Cqdldl.solve!(solver, x)
    @test A*x ≈ b
    @test allocs == 0
end

Random.seed!(1)
n, m = 5, 5
H = Diagonal(rand(n))
D = sprand(m, n, 0.5)
U = [H D'; spzeros(m ,n) Diagonal(fill(-1e-8, m))]
A = Symmetric(U, :U)
b = randn(n + m)

nnzA = nnz(U)
solver = Cqdldl.QDLDLSolver{Float64}(n+m, 2nnzA, 2*nnzA)
testallocs(solver, U, b)

# Change the number nonzero elements
U2 = copy(U)
U2[8,9] = 0.4
@test nnz(U2) > nnz(U)
testallocs(solver, U2, b)

# Change the size of the matrix
a = trues(n + m)
a[8] = false
U3 = U[a,a]
n3 = sum(a)
b3 = randn(size(U3,2))
@test istriu(U3)
testallocs(solver, U3, b3)

F = Cqdldl.QDLDLFactorization(solver)
x3 = F \ b3
@test Symmetric(U3, :U) * x3 ≈ b3


text:
# Matrix: create matrix representation of LinearMap
function Base.Matrix{T}(A::LinearMap) where {T}
    M, N = size(A)
    mat = Matrix{T}(undef, (M, N))
    v = fill(zero(T), N)
    @inbounds for i in 1:N
        v[i] = one(T)
        _unsafe_mul!(view(mat, :, i), A, v)
        v[i] = zero(T)
    end
    return mat
end
Base.Matrix(A::LinearMap{T}) where {T} = Matrix{T}(A)
Base.Array(A::LinearMap) = Matrix(A)
Base.convert(::Type{T}, A::LinearMap) where {T<:Matrix} = T(A)
Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A)
Base.convert(::Type{AbstractMatrix}, A::LinearMap) = convert(Matrix, A)
Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A)

# sparse: create sparse matrix representation of LinearMap
function SparseArrays.sparse(A::LinearMap{T}) where {T}
    M, N = size(A)
    rowind = Int[]
    nzval = T[]
    colptr = Vector{Int}(undef, N+1)
    v = fill(zero(T), N)
    Av = Vector{T}(undef, M)

    @inbounds for i in 1:N
        v[i] = one(T)
        _unsafe_mul!(Av, A, v)
        js = findall(!iszero, Av)
        colptr[i] = length(nzval) + 1
        if length(js) > 0
            append!(rowind, js)
            append!(nzval, Av[js])
        end
        v[i] = zero(T)
    end
    colptr[N+1] = length(nzval) + 1

    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
end
Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A)

# special cases

# ScaledMap
Base.Matrix{T}(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) where {T} =
    convert(Matrix{T}, A.λ * A.lmap.lmap)
SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:MatrixMap}) =
    convert(SparseMatrixCSC, A.λ * A.lmap.lmap)

# UniformScalingMap
Base.Matrix{T}(J::UniformScalingMap) where {T} = Matrix{T}(J.λ*I, size(J))
Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, J.M))

# WrappedMap
Base.Matrix{T}(A::WrappedMap) where {T} = Matrix{T}(A.lmap)
Base.convert(::Type{T}, A::WrappedMap) where {T<:Matrix} = convert(T, A.lmap)
Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)

# TransposeMap & AdjointMap
for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose))
    @eval Base.convert(::Type{AbstractMatrix}, A::$T) = $t(convert(AbstractMatrix, A.lmap))
    @eval SparseArrays.sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap))
end

# LinearCombination
function Base.Matrix{T}(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}}) where {T}
    maps = ΣA.maps
    mats = map(A->getfield(A, :lmap), maps)
    return Matrix{T}(sum(mats))
end
function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{MatrixMap}}})
    maps = ΣA.maps
    mats = map(A->getfield(A, :lmap), maps)
    return convert(SparseMatrixCSC, sum(mats))
end

# CompositeMap
function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, LinearMap}}) where {T}
    B, A = AB.maps
    require_one_based_indexing(B)
    Y = Matrix{eltype(AB)}(undef, size(AB))
    @views for i in 1:size(Y, 2)
        _unsafe_mul!(Y[:, i], A, B.lmap[:, i])
    end
    return Y
end
for ((TA, fieldA), (TB, fieldB)) in (((MatrixMap, :lmap), (MatrixMap, :lmap)),
                                     ((MatrixMap, :lmap), (UniformScalingMap, :λ)),
                                     ((UniformScalingMap, :λ), (MatrixMap, :lmap)))
    @eval function Base.convert(::Type{AbstractMatrix},
                                AB::CompositeMap{<:Any,<:Tuple{$TB,$TA}})
        B, A = AB.maps
        return A.$fieldA*B.$fieldB
    end
end
function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}}) where {T}
    B, A = AB.maps
    return convert(Matrix{T}, A.lmap*B.lmap)
end
function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{MatrixMap, MatrixMap}})
    B, A = AB.maps
    return convert(SparseMatrixCSC, A.lmap*B.lmap)
end
function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}}) where {T}
    A, J = λA.maps
    return convert(Matrix{T}, J.λ*A.lmap)
end
function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{MatrixMap, UniformScalingMap}})
    A, J = λA.maps
    return convert(SparseMatrixCSC, J.λ*A.lmap)
end
function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}}) where {T}
    J, A = Aλ.maps
    return convert(Matrix{T}, A.lmap*J.λ)
end
function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, MatrixMap}})
    J, A = Aλ.maps
    return convert(SparseMatrixCSC, A.lmap*J.λ)
end

# BlockMap & BlockDiagonalMap
Base.Matrix{T}(A::BlockMap) where {T} = hvcat(A.rows, convert.(Matrix{T}, A.maps)...)
Base.convert(::Type{AbstractMatrix}, A::BlockMap) =
    hvcat(A.rows, convert.(AbstractMatrix, A.maps)...)
function SparseArrays.sparse(A::BlockMap)
    return hvcat(
        A.rows,
        convert(SparseMatrixCSC, first(A.maps)),
        convert.(AbstractMatrix, Base.tail(A.maps))...
    )
end
Base.Matrix{T}(A::BlockDiagonalMap) where {T} = Base._cat((1,2), convert.(Matrix{T}, A.maps)...)
Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
function SparseArrays.sparse(A::BlockDiagonalMap)
    return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
end

# KroneckerMap & KroneckerSumMap
Base.Matrix{T}(A::KroneckerMap) where {T} = kron(convert.(Matrix{T}, A.maps)...)
Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) =
    kron(convert.(AbstractMatrix, A.maps)...)
function SparseArrays.sparse(A::KroneckerMap)
    return kron(
        convert(SparseMatrixCSC, first(A.maps)),
        convert.(AbstractMatrix, Base.tail(A.maps))...
    )
end

function Base.Matrix{T}(L::KroneckerSumMap) where {T}
    A, B = L.maps
    IA = Diagonal(ones(Bool, size(A, 1)))
    IB = Diagonal(ones(Bool, size(B, 1)))
    return kron(Matrix{T}(A), IB) + kron(IA, Matrix{T}(B))
end
function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap)
    A, B = L.maps
    IA = Diagonal(ones(Bool, size(A, 1)))
    IB = Diagonal(ones(Bool, size(B, 1)))
    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
end
function SparseArrays.sparse(L::KroneckerSumMap)
    A, B = L.maps
    IA = sparse(Diagonal(ones(Bool, size(A, 1))))
    IB = sparse(Diagonal(ones(Bool, size(B, 1))))
    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
end

# FillMap
Base.Matrix{T}(A::FillMap) where {T} = fill(T(A.λ), size(A))


text:
import Base: start, next, done

# file can be a Cmd, e.g. `xzcat sdfkl32KCsd_enTenTen12.vert.xz`

type SketchEngine; file; dict; unk; eostag; maxlen; column;
    function SketchEngine(file; dict=nothing, eostag="</s>", maxlen=40, column=1, o...)
        if dict == nothing
            unk = 0             # we will construct dict from data
            dict = Dict{Any,Int}(eostag=>eos)
        else
            isa(dict,Dict) || (dict = readvocab(dict))
            unk = length(dict)+1
        end
        new(file, dict, unk, eostag, maxlen, column)
    end
end

unk(s::SketchEngine)=s.unk
maxtoken(s::SketchEngine)=(s.unk > 0 ? s.unk : error("Please specify dict for maxtoken"))

start(s::SketchEngine)=open(s.file)
done(s::SketchEngine,io)=(isa(io,Tuple) && (io=io[1]); eof(io) && (close(io); true))

function next(s::SketchEngine, io)
    sent = Int32[]
    while true
        line = readline(isa(io,Tuple) ? io[1] : io)
        if startswith(line, s.eostag)
            if length(sent) <= s.maxlen
                break
            else
                empty!(sent)
            end
        elseif line[1]=='<' && line[2]!='\t'
            continue
        else
            token = split(line)[s.column]
            if s.unk == 0
                push!(sent, get!(s.dict, token, 1+length(s.dict)))
            else
                push!(sent, get(s.dict, token, s.unk))
            end
        end
    end
    return (sent, io)
end


text:
@testset "alias.jl" begin
    prediction = rand(5)
    lmul!(inv(sum(prediction)), prediction)
    sampler = @inferred(
        ConsistencyResampling.AliasSampler(
            prediction,
            Random.Sampler(Random.GLOBAL_RNG, Float64),
            Random.Sampler(Random.GLOBAL_RNG, 1:5),
        ),
    )
    samples = @inferred(rand(sampler, 100_000))
    @test samples isa Vector{Int}
    @test all(1 ≤ s ≤ 5 for s in samples)
    @test proportions(samples, 1:5) ≈ prediction rtol = 1e-2

    # special cases
    for prediction in ([1.0, 0.0], [0.0, 1.0], [0.0, 1.0, 0.0], [0.5, 0.0, 0.5])
        n = length(prediction)
        sampler = @inferred(
            ConsistencyResampling.AliasSampler(
                prediction,
                Random.Sampler(Random.GLOBAL_RNG, Float64),
                Random.Sampler(Random.GLOBAL_RNG, 1:n),
            ),
        )
        samples = @inferred(rand(sampler, 100_000))
        @test samples isa Vector{Int}
        @test all(1 ≤ s ≤ n for s in samples)
        @test all(!iszero(prediction[s]) for s in samples)
        @test proportions(samples, 1:n) ≈ prediction rtol = 1e-2
    end
end


text:
# const Za = NamedTuple{(:a, :b), Tuple{Branch, Int}}

# struct Za
#   branch::Branch
#   blockid::Int64
# end

struct CFGBlock
  incoming::Array{Tuple{Branch, Int64}}
  outgoing::Array{Branch}
end

struct CFG
  cfg::Array{CFGBlock}
  return_blocks::Array{Tuple{Branch, Int64}}
end

function build_cfg(ir::IR)
  cfg = []
  for b in 1:size(ir.blocks, 1)
    push!(cfg, CFGBlock([], []))
  end
  return_blocks = []
  for b in 1:size(ir.blocks, 1)
    for branch in ir.blocks[b].branches
      if branch.block == 0
        push!(return_blocks, (branch, b))
        continue
      end
      push!(cfg[b].outgoing, branch)
      push!(cfg[branch.block].incoming, (branch, b))
    end
  end
  return CFG(cfg, return_blocks)
end

function build_fwdtypes(ir::IR)
  fwdtypes = Dict{IRTools.Variable, Type}()
  for b in 1:size(ir.blocks, 1)
    block = IRTools.block(ir, b)
    for i in 1:size(IRTools.arguments(block), 1)
      arg = IRTools.arguments(block)[i]
      argtype = IRTools.argtypes(block)[i]
      fwdtypes[arg] = unwrap(argtype)
    end
    for var in keys(block)
      stmt = block[var]
      fwdtypes[var] = unwrap(stmt.type)
    end
  end
  # println(fwdtypes)
  fwdtypes
end


text:
struct Extrapolation{T,N,ITPT,IT,ET} <: AbstractExtrapolation{T,N,ITPT,IT}
    itp::ITPT
    et::ET
end

Base.parent(A::Extrapolation) = A.itp
itpflag(etp::Extrapolation) = itpflag(etp.itp)

# DimSpec{Flag} is not enough for extrapolation dispatch, since we allow nested tuples
# However, no tuples should be nested deeper than this; the first level is for different
# schemes in different dimensions, and the second level is for different schemes in
# different directions.
const ExtrapDimSpec = Union{BoundaryCondition,Tuple{Vararg{Union{BoundaryCondition,NTuple{2,BoundaryCondition}}}}}

etptype(::Extrapolation{T,N,ITPT,IT,ET}) where {T,N,ITPT,IT,ET} = ET
etpflag(etp::Extrapolation{T,N,ITPT,IT,ET}) where {T,N,ITPT,IT,ET} = etp.et

BoundsCheckStyle(etp::AbstractExtrapolation) = CheckWillPass()

"""
`extrapolate(itp, scheme)` adds extrapolation behavior to an interpolation object, according to the provided scheme.

The scheme can take any of these values:

* `Throw` - throws a BoundsError for out-of-bounds indices
* `Flat` - for constant extrapolation, taking the closest in-bounds value
* `Line - linear extrapolation (the wrapped interpolation object must support gradient)
* `Reflect` - reflecting extrapolation (indices must support `mod`)
* `Periodic` - periodic extrapolation (indices must support `mod`)

You can also combine schemes in tuples. For example, the scheme `(Line), Flat())` will use linear extrapolation in the first dimension, and constant in the second.

Finally, you can specify different extrapolation behavior in different direction. `((Line),Flat()), Flat())` will extrapolate linearly in the first dimension if the index is too small, but use constant etrapolation if it is too large, and always use constant extrapolation in the second dimension.
"""
extrapolate(itp::AbstractInterpolation{T,N,IT}, et::ET) where {T,N,IT,ET<:ExtrapDimSpec} =
    Extrapolation{T,N,typeof(itp),IT,ET}(itp, et)

count_interp_dims(::Type{<:Extrapolation{T,N,ITPT}}, n) where {T,N,ITPT} = count_interp_dims(ITPT, n)

@propagate_inbounds function (etp::Extrapolation{T,N})(x::Vararg{Number,N}) where {T,N}
    itp = parent(etp)
    eflag = etpflag(etp)
    xs = inbounds_position(eflag, bounds(itp), x, etp, x)
    g = @inbounds gradient(itp, xs...)
    extrapolate_value(eflag, skip_flagged_nointerp(itp, x), skip_flagged_nointerp(itp, xs), Tuple(g), @inbounds(itp(xs...)))
end
@inline function (etp::Extrapolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
    itp = parent(etp)
    Tret = typeof(lispyprod(zero(T), x...))
    ret = zeros(Tret, shape(x...))
    for (i, y) in zip(eachindex(ret), Iterators.product(x...))
        ret[i] = etp(y...)
    end
    return ret
end

@inline function gradient(etp::AbstractExtrapolation{T,N}, x::Vararg{Number,N}) where {T,N}
    itp = parent(etp)
    if checkbounds(Bool, itp, x...)
        gradient(itp, x...)
    else
        eflag = tcollect(etpflag, etp)
        xs = inbounds_position(eflag, bounds(itp), x, etp, x)
        g = @inbounds gradient(itp, xs...)
        skipni = t->skip_flagged_nointerp(itp, t)
        SVector(extrapolate_gradient.(skipni(eflag), skipni(x), skipni(xs), Tuple(g)))
    end
end

checkbounds(::Bool, ::AbstractExtrapolation, I...) = true

# The last two arguments are just for error-reporting
function inbounds_position(eflag, bounds, x, etp, xN)
    item = inbounds_index(getfirst(eflag), bounds[1], x[1], etp, xN)
    (item, inbounds_position(getrest(eflag), Base.tail(bounds), Base.tail(x), etp, xN)...)
end
inbounds_position(::Any, ::Tuple{}, ::Tuple{}, etp, xN) = ()

# By default, convert all calls to 2-sided calls
inbounds_index(flag::Flag, bounds, x, etp, xN) = inbounds_index((flag, flag), bounds, x, etp, xN)
# But some one-sided calls can be handled more efficiently that way
function inbounds_index(::Throw, (l,u), x, etp, xN)
    @boundscheck(l <= x <= u || Base.throw_boundserror(etp, xN))
    x
end
inbounds_index(::Periodic, (l,u), x, etp, xN) = periodic(x, l, u)
inbounds_index(::Reflect, (l,u), x, etp, xN) = reflect(x, l, u)

# Left-then-right implementations
function inbounds_index((flagl,flagu)::Tuple{Throw,Flag}, (l,u), x, etp, xN)
    @boundscheck(l <= x || Base.throw_boundserror(etp, xN))
    inbounds_index((nothing,flagu), (l,u), x, etp, xN)
end
function inbounds_index((flagl,flagu)::Tuple{Nothing,Throw}, (l,u), x, etp, xN)
    @boundscheck(x <= u || Base.throw_boundserror(etp, xN))
    x
end

function inbounds_index((flagl,flagu)::Tuple{Union{Flat,Line},Flag}, (l,u), x, etp, xN)
    inbounds_index((nothing,flagu), (l,u), maxp(x,l), etp, xN)
end
function inbounds_index((flagl,flagu)::Tuple{Nothing,Union{Flat,Line}}, (l,u), x, etp, xN)
    minp(x,u)
end

function inbounds_index((flagl,flagu)::Tuple{Periodic,Flag}, (l,u), x, etp, xN)
    inbounds_index((nothing,flagu), (l,u), periodic(x, l, u), etp, xN)
end
function inbounds_index((flagl,flagu)::Tuple{Nothing,Periodic}, (l,u), x, etp, xN)
    periodic(x, l, u)
end

function inbounds_index((flagl,flagu)::Tuple{Reflect,Flag}, (l,u), x, etp, xN)
    inbounds_index((nothing,flagu), (l,u), reflect(x, l, u), etp, xN)
end
function inbounds_index((flagl,flagu)::Tuple{Nothing,Reflect}, (l,u), x, etp, xN)
    reflect(x, l, u)
end

minp(a::T, b::T) where T = min(a, b)
minp(a, b) = min(promote(a, b)...)
maxp(a::T, b::T) where T = max(a, b)
maxp(a, b) = max(promote(a, b)...)

function reflect(y, l, u)
    yr = mod(y - l, 2(u-l)) + l
    return ifelse(yr > u, 2u-yr, yr)
end

periodic(y, l, u) = mod(y-l, u-l) + l

@inline function extrapolate_value(eflag, x, xs, g, val)
    val = extrapolate_axis(getfirst(eflag), x[1], xs[1], g[1], val)
    extrapolate_value(getrest(eflag), Base.tail(x), Base.tail(xs), Base.tail(g), val)
end
extrapolate_value(::Any, ::Tuple{}, ::Tuple{}, ::Tuple{}, val) = val

extrapolate_axis(::Flag, x, xs, g, val) = val
extrapolate_axis(::Line, x, xs, g, val) = val + (x-xs)*g

extrapolate_axis((flagl,flagu)::Tuple{Flag,Flag}, x, xs, g, val) =
    extrapolate_axis((nothing,flagu), x, xs, g, val)
extrapolate_axis((flagl,flagu)::Tuple{Nothing,Flag}, x, xs, g, val) = val

extrapolate_axis((flagl,flagu)::Tuple{Line, Flag}, x, xs, g, val) =
    extrapolate_axis((nothing,flagu), x, xs, g, ifelse(x < xs, val + (x-xs)*g, val))
extrapolate_axis((flagl,flagu)::Tuple{Nothing, Line}, x, xs, g, val) =
    ifelse(x > xs, val + (x-xs)*g, val)

extrapolate_gradient(::Flat, x, xs, g) = ifelse(x==xs, g, zero(g))
extrapolate_gradient(::Flag, x, xs, g) = g

include("filled.jl")


text:
"""
    isTimedOut(timeLimit::Union{T,Nothing},startTime::U, additionalTimeElapsed::V)where{T<:Real, U<:Real, V<:Real}

Check if a timer has been exceeded, always true if timeLimit is nothing

# Arguments
- `timeLimit::Union{T,Nothing}`: the time limit in seconds (or nothing)
- `startTime::U`: the time the program started (in seconds since the epoch)
- `additionalTimeElapsed::V`: additional time used (in seconds)
"""
function isTimedOut(timeLimit::Union{T,Nothing},startTime::U, additionalTimeElapsed::V)where{T<:Real, U<:Real, V<:Real}
    if isnothing(timeLimit)
        return false
    end
    if (time() - startTime + additionalTimeElapsed) >= timeLimit
        return true
    else
        return false
    end
end


text:
using Plots

# TO solution
x̄, ū = unpack(z̄, prob)
@show tf_nom = sum([ū[t][end] for t = 1:T-1]) # total time
t_nominal = range(0, stop = tf_nom, length = T)
plot(t_nominal, hcat(x̄...)', width = 2.0)
plot(t_nominal[1:end-1], hcat(ū...)[1:1, :]',
	linetype = :steppost, width = 2.0, color = goldenrod_color)

# DPO solution
x, u = unpack(z, prob_dpo.prob.prob.nom)
@show tf_dpo = sum([u[t][end] for t = 1:T-1])
t_dpo = range(0, stop = tf_dpo, length = T)
plot(t_dpo, hcat(x...)', width = 2.0)
plot(t_dpo[1:end-1], hcat(u...)[1:1, :]',
	linetype = :steppost, width = 2.0, color = red_color)


text:
using Plots
using Random

function c_next(c, delta)
    c_new = deepcopy(c)
    for j in eachindex(c)
        c_new[j] = (c[j] + 1im * delta * c[j - 1] + 1im + conj(delta) * c[j + 1])
    end

    N = normalization_factor(c_new)
    c_new / sqrt(N)
    return c_new
end

function normalization_factor(c)

end

anim = @animate for i=1:100
    plot()
end


text:
module LowRankSketch

import Base: *, show

include("dim_redux.jl")
include("sketch.jl")

end # module


text:
#=
title: IsingSimEUI( manipulatable Ising model simulation with Electron UI )
version: 1.0
date: 2020/08/12
author: okimebarun
url: https://github.com/okimebarun/
url: https://qiita.com/oki_mebarun/
=#

module IsingSimEUI

greet() = print("Hello World!")

###################################################################################################
# load libraries
using Electron

# load local modules
include("HexaLatticeModel03a.jl")
include("SquareLatticeModel03a.jl")
include("IsingMCMC04a.jl")

###################################################################################################
# prepare
# setting for energy function

function genInitialJH(d=d)
    Jd = [Dict{Int,Int8}() for i in 1:d.N]
    Hd = zeros(Float32, d.N)
    (Jd, Hd)
end

function setJH!( Jd=Jd, Hd=Hd, d=d)
    J0 = 1
    H0 = 0
    for i in 1:d.N
        for j in getNeighbors(i,d)
            @inbounds Jd[i][j] = J0
        end
        @inbounds Hd[i] = H0
    end
end

function setH!(H0, Hd=Hd, d=d)
    @simd for i in 1:d.N
        @inbounds Hd[i] = H0
    end
end

###################################################################################################
# JS command

function jscmd(cmd, win=win)
    run(win, cmd)
end

function drawBH(H, S, d=d, N=N)
    x = H
    y = sum(S)/N
    jscmd("drawbox02.setrange(-1.1,1.1,-1.1,1.1);")
    jscmd("trajectory.add($(x), $(y) );")
end

function drawSim02(snap, d=d, N=N)
    jscmd("drawbox01.setrange(0,$(d.ncols),0,$(d.nrows));")

    (xs,ys) = genPoly(1, d)
    len = length(xs)
    ps = zeros(len*2*N)
    i1 = 0
    @inbounds @simd for i in 1:N
        # bits
        (xs,ys) = genPoly(i, d)
        if (snap[i] == 1) # up-spin
            for i in 0:(len-1)
                ps[i*2 + 1 + i1*len] = xs[i + 1]
                ps[i*2 + 2 + i1*len] = ys[i + 1]
            end
            i1 += 2
        else # down-spin
        end
    end
    jscmd("drawbox01.clear();")
    jscmd("drawbox01.drawPoly02($(string(ps)), $(string(len)) );")
end

###################################################################################################
# setting for model
d = SquareLatticeModel.ParamDict(128,128,0)
N = d.nrows * d.ncols
d.N = N

###################################################################################################
# simulation main

main_html_uri = string("file:///", replace(joinpath(@__DIR__, "IsingView01b.html"), '\\' => '/'))

function winmain(mdlno)
    ##########
    # setting for simulation
    trial = 10^5 # MCMC trial

    # gen functors as global
    global getNeighbors, genPoly
    if mdlno == 1
        (getNeighbors, genPoly)= SquareLatticeModel.genFunctor(d)
    elseif mdlno == 2
        (getNeighbors, genPoly)= HexaLatticeModel.genFunctor(d)
    end

    # generate Jd, Hd
    (Jd, Hd) = genInitialJH(d)
    setJH!(Jd, Hd)

    # simulation
    global E, initS, MCMC
    (E, initS, MCMC) = IsingMCMC.genFunctor(Jd, Hd)
    S = ones(Int8, N)
    initS(S, N)

    ##########
    global win
    win = Window( URI(main_html_uri))
    ElectronAPI.setBounds(win, Dict("width"=>1050, "height"=>700))
    ch = msgchannel(win)

    while true
        local msg, jscmd
        try
            msg = take!(ch)
        catch
            println("channel closed.")
            break
        end

        m = match(r"^calc T:(?<T>[\d\.]+) H:(?<H>[-\d\.]+)", msg)
        if !(m === nothing)
            T = parse(Float32, m[:T])
            H = parse(Float32, m[:H])
            print("calc with T=$(T) & H=$(H)")
            setH!(H, Hd, d)
            #
            @time (simE) = MCMC(S, T, N, trial);
            drawSim02(S)
            drawBH(H, S)
        end
    end
end

###################################################################################################
# main functions for PackageCompiler

function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

#
function real_main()
    @show ARGS
    @show Base.PROGRAM_FILE

    # do something
    winmain(1)
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end
#
end # module


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PumpParams <: TechnicalParams
        epnt_data::Union{Nothing, EPANETPumpParams}
        norm_coefs::Union{Nothing, NormPumpParams}
        base_price::Union{Nothing, Float64}
        pattern_name::Union{Nothing, String}
        _forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

Pump characteristics, containing both inp-file data and coefficients for a normalized representation. Currently, only single-speed

# Arguments
- `epnt_data::Union{Nothing, EPANETPumpParams}`: Pump specifications from inp-file
- `norm_coefs::Union{Nothing, NormPumpParams}`: Coefficients that characterize pumps with a normalized representation of head, efficiency, and power. Currently only single-speed
- `base_price::Union{Nothing, Float64}`: 'base' price in dollar/J
- `pattern_name::Union{Nothing, String}`: name of forecast pattern array
- `_forecasts::InfrastructureSystems.Forecasts`
- `internal::InfrastructureSystemsInternal`: internal reference, do not modify
"""
mutable struct PumpParams <: TechnicalParams
    "Pump specifications from inp-file"
    epnt_data::Union{Nothing, EPANETPumpParams}
    "Coefficients that characterize pumps with a normalized representation of head, efficiency, and power. Currently only single-speed"
    norm_coefs::Union{Nothing, NormPumpParams}
    "'base' price in dollar/J"
    base_price::Union{Nothing, Float64}
    "name of forecast pattern array"
    pattern_name::Union{Nothing, String}
    _forecasts::InfrastructureSystems.Forecasts
    "internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PumpParams(epnt_data, norm_coefs, base_price, pattern_name, _forecasts=InfrastructureSystems.Forecasts(), )
    PumpParams(epnt_data, norm_coefs, base_price, pattern_name, _forecasts, InfrastructureSystemsInternal(), )
end

function PumpParams(; epnt_data, norm_coefs, base_price, pattern_name, _forecasts=InfrastructureSystems.Forecasts(), )
    PumpParams(epnt_data, norm_coefs, base_price, pattern_name, _forecasts, )
end

# Constructor for demo purposes; non-functional.
function PumpParams(::Nothing)
    PumpParams(;
        epnt_data=nothing,
        norm_coefs=nothing,
        base_price=nothing,
        pattern_name=nothing,
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get PumpParams epnt_data."""
get_epnt_data(value::PumpParams) = value.epnt_data
"""Get PumpParams norm_coefs."""
get_norm_coefs(value::PumpParams) = value.norm_coefs
"""Get PumpParams base_price."""
get_base_price(value::PumpParams) = value.base_price
"""Get PumpParams pattern_name."""
get_pattern_name(value::PumpParams) = value.pattern_name
"""Get PumpParams _forecasts."""
get__forecasts(value::PumpParams) = value._forecasts
"""Get PumpParams internal."""
get_internal(value::PumpParams) = value.internal


text:
using DrWatson
@quickactivate "StatReth"

# %% 3.1
Pr_Positive_Vampire = 0.95
Pr_Positive_Mortal = 0.01
Pr_Vampire = 0.001
Pr_Positive = Pr_Positive_Vampire * Pr_Vampire + Pr_Positive_Mortal * (1 - Pr_Vampire)
Pr_Vampire_Positive = Pr_Positive_Vampire * Pr_Vampire / Pr_Positive


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function lucompletepivCopy!(A)
    n = size(A, 1)
    rowpiv=zeros(Int, n-1)
    colpiv=zeros(Int, n-1)
    for k = 1:n-1
        As = abs(A[k:n, k:n])
        μ, λ = ind2sub(size(As), indmax(As))
        μ += k-1; λ += k-1
        rowpiv[k] = μ
        A[[k,μ], 1:n] = A[[μ,k], 1:n]
        colpiv[k] = λ
        A[1:n, [k,λ]] = A[1:n, [λ,k]]
        if A[k,k] ≠ 0
            ρ = k+1:n
            A[ρ, k] = A[ρ, k]/A[k, k]
            A[ρ, ρ] = A[ρ, ρ] - A[ρ, k] * A[k, ρ]
        end
    end
    return (A, rowpiv, colpiv)
end

function lucompletepivSub!(A)
    n = size(A, 1)
    rowpiv=zeros(Int, n-1)
    colpiv=zeros(Int, n-1)
    for k = 1:n-1
        As = abs(sub(A, k:n, k:n))
        μ, λ = ind2sub(size(As), indmax(As))
        μ += k-1; λ += k-1
        rowpiv[k] = μ
        A[[k,μ], 1:n] = sub(A, [μ,k], 1:n)
        colpiv[k] = λ
        A[1:n, [k,λ]] = sub(A, 1:n, [λ,k])
        if A[k,k] ≠ 0
            ρ = k+1:n
            A[ρ, k] = sub(A, ρ, k)/A[k, k]
            A[ρ, ρ] = sub(A, ρ, ρ) - sub(A, ρ, k) * sub(A, k, ρ)
        end
    end
    return (A, rowpiv, colpiv)
end


text:
using Plots
using LinearAlgebra

function main( np, nstep )

    ϵ = 1/np

    q = [zeros(2) for i in 1:np]
    q[1] = [0.5,0.5]

    for klm = 2:np
        frein = 0
        overlap = 1
        p = zeros(2)
        while ((overlap == 1) && frein < 7000)
            p .= 2ϵ * ones(2) + (1 - 4ϵ) * rand(2)
            overlap2 = 0
            for subh = 1:(klm-1)
                if norm(p - q[subh],2) < 2ϵ
                    overlap2 = 1
                end
            end
            overlap = overlap2
            frein += 1
            if frein == 10000
                frein
                @error "frein trop grand"
                exit
            end
        end
        q[klm] = p
    end

    offset = [[0, 0], [1,  0], [-1,  0],
              [0, 1], [0, -1], [-1,  1],
              [1, 1], [1, -1], [-1, -1]]

    function compute_dt( qi, qj, vi, vj )

        δz = qi .- qj
        δv = vi .- vj

        c = (δv'δz).^2 - (δv'δv) * (δz'δz - (2ϵ)^2)

        if δv'δz >= 0
            δt = Inf
        elseif c < 0
            δt = Inf
        else
            δt = -(δv'δz + sqrt(c)) / (δv'δv)
        end

        return δt

    end

    v = [randn(2) for j = 1:np]

    Collisions = Inf .* ones(np,np)
    Fantome = zeros(Int, np,np)

    for k=1:np, l=k+1:np

        t = [compute_dt(q[l]+offset[i],q[k],v[l],v[k]) for i in 1:9]

        Collisions[k,l] = minimum(t)
        if isinf(Collisions[k,l])
            Fantome[k,l] = 0
        else
            Fantome[k,l] = argmin(t)
        end
    end

    time = 0

    anim = @animate for step in 1:nstep

        m, n = Tuple(argmin(Collisions))
        dt = Collisions[m,n]
        time += dt

        if Fantome[m,n] > 0
            qa = [q[n] + offset[k] + dt*v[n] for k in 1:9]
            qb = q[m] + dt .* v[m]
            vr = argmin([norm(qx-qb) for qx in qa])

            q[n] = qa[vr]
            q[m] = qb

            J = (dot(v[n]-v[m],q[n]-q[m]))/ 2ϵ
            v[m] = v[m] + J * (q[n]-q[m]) / 2ϵ
            v[n] = v[n] - J * (q[n]-q[m]) / 2ϵ
            q[n] = mod.(q[n],1)
            q[m] = mod.(q[m],1)
        end

        for i in [1:min(m,n)-1;min(m,n)+1:max(m,n)-1;max(m,n)+1:np]
            q[i] = q[i] .+ dt .* v[i]
            q[i] = mod.(q[i],1)
        end

        plt = plot(size  = (200,200),
                   xlims = (0,1),
                   ylims = (0,1),
                   grid  = false,
                   axis  = nothing,
                   framestyle= :none,
                   legend=false,
                   widen = false)

        c = trunc(Int, 200ϵ)
        scatter!( getindex.(q,1), getindex.(q,2), markershape  = :circle,
              markersize   = c ,
              aspect_ratio = :equal)

        Collisions .-= dt
        Collisions[m,:] .= Inf
        Collisions[n,:] .= Inf
        Collisions[:,m] .= Inf
        Collisions[:,n] .= Inf

        for j in [m,n], l in [1:j-1;j+1:np]

            t=[compute_dt(q[l]+offset[k],q[j],v[l],v[j]) for k in 1:9]
            Collisions[j,l]=minimum(t)

            if isinf(Collisions[m,l])
                Fantome[j,l]=0
            else
                Fantome[j,l]=argmin(t)
            end

        end

    end

    gif(anim, joinpath(@__DIR__, "hs_periodic.gif"), fps = 10)

    return time

end

@time main(100, 1000)


text:
using DFControl

try
    Meta.parse("5")
catch
    "5"
end

#maybe use .attached calculations


text:
# This file was generated, do not modify it. # hide
using LinearAlgebra
a = [1, 2, 3, 3, 4, 5, 2, 2]
@show dot(a, a)
println(dot(a, a))


text:
# LU back-substitution in Julia (ala Numerical Recipes)
function banbks(a::Array{Float64,2},n::Int64,m1::Int64,m2::Int64,al::Array{Float64,2},indx::Vector{Int64},b::Vector{Float64})
mm=m1+m2+1
@assert(size(a)  == (n,mm))
@assert(size(al) == (n,m1))
@assert(length(indx) == n)
@assert(length(b) == n)
l=m1
for k=1:n
  i=indx[k]
  if (i != k)
    dum  = b[k]
    b[k] = b[i]
    b[i] = dum
  end
  if (l < n)
    l += 1
  end
  for i=k+1:l
    b[i] -= al[k,i-k]*b[k]
  end
end
l=1
for i in n:-1:1
  dum=b[i]
  for k=2:l
    dum -= a[i,k]*b[k+i-1]
  end
  b[i]=dum/a[i,1]
  if (l < mm)
    l += 1
  end
end
return
end

# (C) Copr. 1986-92 Numerical Recipes Software

function banbks(a::Array{Complex{Float64},2},n::Int64,m1::Int64,m2::Int64,al::Array{Complex{Float64},2},indx::Vector{Int64},b::Vector{Complex{Float64}})
mm=m1+m2+1
@assert(size(a)  == (n,mm))
@assert(size(al) == (n,m1))
@assert(length(indx) == n)
@assert(length(b) == n)
l=m1
for k=1:n
  i=indx[k]
  if (i != k)
    dum  = b[k]
    b[k] = b[i]
    b[i] = dum
  end
  if (l < n)
    l += 1
  end
  for i=k+1:l
    b[i] -= al[k,i-k]*b[k]
  end
end
l=1
for i in n:-1:1
  dum=b[i]
  for k=2:l
    dum -= a[i,k]*b[k+i-1]
  end
  b[i]=dum/a[i,1]
  if (l < mm)
    l += 1
  end
end
return
end

# (C) Copr. 1986-92 Numerical Recipes Software