text:
"""
    jensen_alpha(asset_returns, benchmark_returns; risk_free=0.0)

Jensen's alpha, or simply alpha, is a risk-adjusted excess-performance measure. It indicates the average return of an investment above or below that predicted by the capital asset pricing model (CAPM). The measure adjusts the excess returns such that the risk is identical for the investment and benchmark.

# Arguments
- `asset_returns`:      Vector of asset returns.
- `benchmark_returns`:  Vector of benchmark returns (e.g. market portfolio returns for CAPM beta).
- `risk_free`:          Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).

# Returns
Jensen's alpha measure.

# Sources
- Bacon, Carl (2008). Practical Portfolio Performance Measurement and Attribution, 2nd Edition, John Wiley & Sons Ltd. Page 72.
"""
function jensen_alpha(asset_returns, benchmark_returns; risk_free=0.0)
    (α, β) = capm(asset_returns, benchmark_returns; risk_free)
    α
end

"""
    modified_jensen(asset_returns, benchmark_returns; risk_free=0.0)

Divides Jensen’s alpha by the systematic risk, which measures the systematic risk-adjusted return per unit of systematic risk. See also `jensen_alpha`

# Arguments
- `asset_returns`:      Vector of asset returns.
- `benchmark_returns`:  Vector of benchmark returns (e.g. market portfolio returns for CAPM beta).
- `risk_free`:          Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).

# Returns
Modified Jensen's alpha measure.

# Sources
- Bacon, Carl (2008). Practical Portfolio Performance Measurement and Attribution, 2nd Edition, John Wiley & Sons Ltd. Page 77.
"""
function modified_jensen(asset_returns, benchmark_returns; risk_free=0.0)
    (α, β) = capm(asset_returns, benchmark_returns; risk_free)
    α / β
end


text:
N = 10
filenames = readdir(joinpath("test", "data"))
for filename in filenames
	full_filename = Pkg.dir("DataFrames", "test", "data", filename)
	df = benchmark(() -> readtable(full_filename),
	                                "DataFrame I/O",
	                                filename,
	                                N)
	# TODO: Keep permanent record
	printtable(df, header=false)
end

ext: jl


text:
@inline function add_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:AbstractFloat}
    yhi, ylo = y
    yhi, ylo = add_2(x, yhi, ylo)
    return yhi, ylo
end

@inline function sub_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:AbstractFloat}
    yhi, ylo = y
    yhi, ylo = sub_2(x, yhi, ylo)
    return yhi, ylo
end

@inline function mul_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:AbstractFloat}
    hi, lo = y
    hihi, hilo = mul_2(x, hi)
    lohi, lolo = mul_2(x, lo)
    hi, lo = add_2(hihi, hilo, lohi, lolo)
    return hi, lo
end

@inline function dvi_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:AbstractFloat}
    yinv = inv_dd_dd(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end

@inline function dvi_fpdd_dd_fast(x::T, y::Tuple{T,T}) where {T<:AbstractFloat}
    yinv = inv_dd_dd_fast(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end


text:
using Test
using LatticeTools

@testset "UnitCell" begin
    latticevectors = [0.5 0.0; 0.0 1.0]

    @testset "Constructors" begin
        uc = makeunitcell(latticevectors)
        @test isapprox(uc.latticevectors, latticevectors)
        @test isapprox(uc.reducedreciprocallatticevectors, [2.0 0.0; 0.0 1.0])
        @test isapprox(uc.reciprocallatticevectors, [4*pi 0.0; 0.0 2*pi])

        uc2 = makeunitcell([[0.5,0.0], [0.0,1.0]])

        @test_throws ArgumentError makeunitcell(latticevectors; SiteType=Int)
        @test_throws ArgumentError makeunitcell([1.0 0.0;]; SiteType=String)
        @test_throws ArgumentError makeunitcell([1.0 0.0; 1.0 0.0]; SiteType=String)
        @test_throws ArgumentError makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String, tol=-1E-8)

        uc1d = makeunitcell(2.0)
        uc1d.latticevectors == 2.0 * ones(Float64, 1, 1)
    end

    @testset "Equality" begin
        uc1 = makeunitcell([2.0 0.0; 1.0 2.0]; SiteType=String)
        uc2 = makeunitcell([2.0 0.0; 1.0 2.0]; SiteType=String)
        uc3 = makeunitcell([2.0 0.0; 1.0 2.0]; SiteType=Tuple{String, Int})
        uc4 = makeunitcell([2.0 0.0; 0.0 2.0]; SiteType=String)
        @test uc1 == uc2
        @test uc1 != uc3
        @test uc1 != uc4
    end

    @testset "AngledLatticeVectors" begin
        uc1 = makeunitcell([2.0 0.0; 1.0 2.0])
        uc2 = makeunitcell([[2.0, 1.0], [0.0, 2.0]])
        @test uc1 == uc2
        @test isapprox(uc1.reciprocallatticevectors / π, [1.0 -0.5; 0.0 1.0])
        @test isapprox(uc2.reciprocallatticevectors / π, [1.0 -0.5; 0.0 1.0])
    end

    @testset "Methods" begin
        uc = makeunitcell(latticevectors)
        fc1 = FractCoord([0, 0], [0.5, 0.0])
        fc2 = FractCoord([0, 0], [0.0, 0.5])
        index1 = addsite!(uc, "Ox", fc1)
        index2 = addsite!(uc, "Oy", fc2)
        @test index1 == 1
        @test index2 == 2

        @test dimension(uc) == 2
        @test numsite(uc) == 2
        @test hassite(uc, "Ox")
        @test hassite(uc, "Oy")
        @test !hassite(uc, "Oz")
        @test getsiteindex(uc, "Ox") == 1
        @test getsiteindex(uc, "Oy") == 2
        @test getsite(uc, "Ox") == ("Ox", fc1)
        @test getsite(uc, "Oy") == ("Oy", fc2)
        @test getsitecoord(uc, "Ox") == fc1
        @test getsitecoord(uc, "Oy") == fc2
        @test getsiteindexcoord(uc, "Ox") == (1, fc1)
        @test getsiteindexcoord(uc, "Oy") == (2, fc2)

        @test getsite(uc, 1) == ("Ox", fc1)
        @test getsite(uc, 2) == ("Oy", fc2)
        @test getsitename(uc, 1) == "Ox"
        @test getsitename(uc, 2) == "Oy"
        @test getsitecoord(uc, 1) == fc1
        @test getsitecoord(uc, 2) == fc2

        # orbitals
        @test numorbital(uc) == 0
        @test numorbitals(uc) == 0
        @test orbitalcount(uc) == 0

        addorbital!(uc, "px1", 1)
        addorbital!(uc, "py1", "Ox")
        addorbital!(uc, "px2", 2)
        addorbital!(uc, "py2", "Oy")

        @test numorbital(uc) == 4
        @test numorbitals(uc) == 4
        @test orbitalcount(uc) == 4

        @test hasorbital(uc, "px1")
        @test hasorbital(uc, "py1")
        @test !hasorbital(uc, "pz1")
        @test getorbitalindex(uc, "px1") == 1
        @test getorbitalindex(uc, "py1") == 2
        @test getorbital(uc, 2) == ("py1", 1)
        @test getorbital(uc, "px2") == ("px2", 2)
        @test getorbitalname(uc, 1) == "px1"
        @test getorbitalname(uc, 2) == "py1"
        @test getorbitalsiteindex(uc, 1) == 1
        @test getorbitalsiteindex(uc, 2) == 1
        @test getorbitalsiteindex(uc, 3) == 2
        @test getorbitalsiteindex(uc, 4) == 2
        @test getorbitalcoord(uc, 3) == getsitecoord(uc, 2)
    end

    @testset "Methods Exceptions" begin
        uc = makeunitcell(latticevectors; SiteType=String)
        addsite!(uc, "Ox", FractCoord([0, 0], [0.0, 0.0]))
        @test_throws ArgumentError addsite!(uc, "Oy", FractCoord([0], [0.5]))
        @test_throws ArgumentError addsite!(uc, "Ox", FractCoord([0, 0], [0.0, 0.5]))
        @test_throws ArgumentError addsite!(uc, "Ox2", FractCoord([0, 0], [0.0, 0.0]))
        @test_throws ArgumentError addsite!(uc, "Ox2", FractCoord([1, 0], [0.0, 0.0]))
    end

    @testset "Type" begin
        uc = makeunitcell(latticevectors; SiteType=String)
        fc = Dict("Ox" => FractCoord([0, 0], [0.5, 0.0]),
                  "Oy" => FractCoord([0, 0], [0.0, 0.5]))
        for site in ["Ox", "Oy"]
            addsite!(uc, site, fc[site])
        end
        @test getsitename(uc, 1) == "Ox"
        @test getsitename(uc, 2) == "Oy"
    end

    @testset "fract2carte/carte2fract" begin
        latticevectors = [0.5 0.0; 0.0 1.0]
        uc = makeunitcell(latticevectors)

        rawfractcoord = [-1.2, 1.5]

        correctfractcoord = FractCoord([-2, 1], [0.8, 0.5])
        correctcartecoord = [-0.6, 1.5]

        fractcoord = FractCoord(rawfractcoord)
        cartecoord = fract2carte(uc, fractcoord)
        newfractcoord = carte2fract(uc, cartecoord)
        @test isapprox(fractcoord, correctfractcoord)
        @test isapprox(cartecoord, correctcartecoord)
        @test isapprox(newfractcoord, correctfractcoord)

        @testset "tolerance" begin
            latticevectors = [1.0 0.0; 0.0 1.0]
            uc = makeunitcell(latticevectors)

            fc1 = carte2fract(uc, [1.0 - 1E-9, 1.0])
            fc2 = carte2fract(uc, [1.0 - 1E-9, 1.0]; tol=0.0)
            @test fc1.whole == [1, 1]
            @test isapprox(fc1.fraction, [0.0, 0.0])
            @test fc2.whole == [0, 1]
            @test isapprox(fc2.fraction, [1.0, 0.0])
        end
    end

    @testset "fract2carte/carte2fract exception" begin
        latticevectors = [0.5 0.0; 0.0 1.0]
        uc1 = makeunitcell(latticevectors)

        carte2fract(uc1, [1.0, 2.0])
        @test_throws ArgumentError carte2fract(uc1, [1.0])
        @test_throws ArgumentError carte2fract(uc1, [1.0, 2.0, 3.0])

        fract2carte(uc1, FractCoord([1,2], [0.1, 0.2]))
        @test_throws ArgumentError fract2carte(uc1, FractCoord([1], [0.1]))
        @test_throws ArgumentError fract2carte(uc1, FractCoord([1,2,3], [0.1, 0.2, 0.3]))
    end

    @testset "momentumgrid" begin
        uc = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)

        @test_throws ArgumentError momentumgrid(uc, [2,])
        @test_throws ArgumentError momentumgrid(uc, [2,3,4])
        @test_throws ArgumentError momentumgrid(uc, [-2,3])

        kg = momentumgrid(uc, [2,3])
        for i in 1:2, j in 1:3
            @test isapprox(kg[i,j], [2π * (i-1) / 2, 2π * (j-1) / 3])
        end
    end

    @testset "whichunitcell" begin
        uc = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(uc, "A", FractCoord([0,0], [0.1, 0.1]))
        addsite!(uc, "B", FractCoord([0,0], [0.2, 0.2]))
        @test whichunitcell(uc, "A", [1.1, 2.1]) == [1, 2]
        @test whichunitcell(uc, "A", FractCoord([1,2], [0.1, 0.1])) == [1, 2]
        @test_throws ArgumentError whichunitcell(uc, "B", [1.1, 2.1])
        @test_throws ArgumentError whichunitcell(uc, "B", FractCoord([1,2], [0.1, 0.1]))
    end

    @testset "findsiteindex" begin
        uc = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(uc, "A", FractCoord([0, 0], [0.0, 0.0]))
        @test findsiteindex(uc, FractCoord([0,0], [0.0, 0.0])) == (1, [0,0])
        @test findsiteindex(uc, FractCoord([1,0], [0.0, 0.0])) == (1, [1,0])

        @test findsiteindex(uc, FractCoord([0,0], [0.2, 0.2])) == (-1, Int[])
    end

end


text:
module Interpolation

using CUDA
using DocStringExtensions
using LinearAlgebra
using MPI
using OrderedCollections
using StaticArrays
using KernelAbstractions

using ClimateMachine
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Geometry
import ClimateMachine.Mesh.Elements: baryweights
import ClimateMachine.MPIStateArrays: array_device

export dimensions,
    accumulate_interpolated_data,
    accumulate_interpolated_data!,
    InterpolationBrick,
    InterpolationCubedSphere,
    interpolate_local!,
    project_cubed_sphere!,
    InterpolationTopology

abstract type InterpolationTopology end

"""
    InterpolationBrick{
        FT <: AbstractFloat,CuArrays
        UI8AD <: AbstractArray{UInt8, 2},
        UI16VD <: AbstractVector{UInt16},
        I32V <: AbstractVector{Int32},
    } <: InterpolationTopology

This interpolation data structure and the corresponding functions works for a
brick, where stretching/compression happens only along the x1, x2 & x3 axis.
Here x1 = X1(ξ1), x2 = X2(ξ2) and x3 = X3(ξ3).

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    InterpolationBrick(
        grid::DiscontinuousSpectralElementGrid{FT},
        xbnd::Array{FT,2},
        xres,
    ) where FT <: AbstractFloat

This interpolation structure and the corresponding functions works for a brick,
where stretching/compression happens only along the x1, x2 & x3 axis. Here x1
= X1(ξ1), x2 = X2(ξ2) and x3 = X3(ξ3).

# Arguments for the inner constructor
 - `grid`: DiscontinousSpectralElementGrid
 - `xbnd`: Domain boundaries in x1, x2 and x3 directions
 - `xres`: Resolution of the interpolation grid in x1, x2 and x3 directions
"""
struct InterpolationBrick{
    FT <: AbstractFloat,
    I <: Int,
    FTV <: AbstractVector{FT},
    FTVD <: AbstractVector{FT},
    IVD <: AbstractVector{I},
    FTA2 <: Array{FT, 2},
    UI8AD <: AbstractArray{UInt8, 2},
    UI16VD <: AbstractVector{UInt16},
    I32V <: AbstractVector{Int32},
} <: InterpolationTopology
    "Number of elements"
    Nel::I
    "Total number of interpolation points"
    Np::I
    "Total number of interpolation points on local process"
    Npl::I
    "Domain bounds in x1, x2 and x3 directions"
    xbnd::FTA2
    "Interpolation grid in x1 direction"
    x1g::FTV
    "Interpolation grid in x2 direction"
    x2g::FTV
    "Interpolation grid in x3 direction"
    x3g::FTV
    "Unique ξ1 coordinates of interpolation points within each spectral element"
    ξ1::FTVD
    "Unique ξ2 coordinates of interpolation points within each spectral element"
    ξ2::FTVD
    "Unique ξ3 coordinates of interpolation points within each spectral element"
    ξ3::FTVD
    "Flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point"
    flg::UI8AD
    "Normalization factor"
    fac::FTVD
    "x1 interpolation grid index of interpolation points within each element on the local process"
    x1i::UI16VD
    "x2 interpolation grid index of interpolation points within each element on the local process"
    x2i::UI16VD
    "x3 interpolation grid index of interpolation points within each element on the local process"
    x3i::UI16VD
    "Offsets for each element"
    offset::IVD    # offsets for each element for v
    "GLL points in ξ1 direction"
    m_ξ1::FTVD
    "GLL points in ξ2 direction"
    m_ξ2::FTVD
    "GLL points in ξ3 direction"
    m_ξ3::FTVD
    "Barycentric weights"
    wb1::FTVD
    "Barycentric weights"
    wb2::FTVD
    "Barycentric weights"
    wb3::FTVD
    # MPI setup for gathering interpolated variable on proc # 0
    "Number of interpolation points on each of the processes"
    Np_all::I32V

    "x1 interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    x1i_all::UI16VD
    "x2 interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    x2i_all::UI16VD
    "x3 interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    x3i_all::UI16VD

    function InterpolationBrick(
        grid::DiscontinuousSpectralElementGrid{FT},
        xbnd::Array{FT, 2},
        x1g::AbstractArray{FT, 1},
        x2g::AbstractArray{FT, 1},
        x3g::AbstractArray{FT, 1},
    ) where {FT <: AbstractFloat}
        mpicomm = grid.topology.mpicomm
        pid = MPI.Comm_rank(mpicomm)
        npr = MPI.Comm_size(mpicomm)

        DA = arraytype(grid)                    # device array
        device = arraytype(grid) <: Array ? CPU() : CUDADevice()

        qm = polynomialorders(grid) .+ 1
        ndim = 3
        toler = 4 * eps(FT) # tolerance

        n1g = length(x1g)
        n2g = length(x2g)
        n3g = length(x3g)

        Np = n1g * n2g * n3g
        marker = BitArray{3}(undef, n1g, n2g, n3g)
        fill!(marker, true)

        Nel = length(grid.topology.realelems) # # of elements on local process
        offset = Vector{Int}(undef, Nel + 1) # offsets for the interpolated variable
        n123 = zeros(Int, ndim)        # # of unique ξ1, ξ2, ξ3 points in each cell
        xsten = zeros(Int, 2, ndim)        # x1, x2, x3 start and end for each brick element
        xbndl = zeros(FT, 2, ndim)        # x1,x2,x3 limits (min,max) for each brick element

        ξ1 = map(i -> zeros(FT, i), zeros(Int, Nel))
        ξ2 = map(i -> zeros(FT, i), zeros(Int, Nel))
        ξ3 = map(i -> zeros(FT, i), zeros(Int, Nel))

        x1i = map(i -> zeros(UInt16, i), zeros(UInt16, Nel))
        x2i = map(i -> zeros(UInt16, i), zeros(UInt16, Nel))
        x3i = map(i -> zeros(UInt16, i), zeros(UInt16, Nel))

        x = map(i -> zeros(FT, ndim, i), zeros(Int, Nel)) # interpolation grid points embedded in each cell

        offset[1] = 0

        for el in 1:Nel
            for (xg, dim) in zip((x1g, x2g, x3g), 1:ndim)
                xbndl[1, dim], xbndl[2, dim] =
                    extrema(grid.topology.elemtocoord[dim, :, el])

                st = findfirst(xg .≥ xbndl[1, dim] .- toler)
                if st ≠ nothing
                    if xg[st] > (xbndl[2, dim] + toler)
                        st = nothing
                    end
                end

                if st ≠ nothing
                    xsten[1, dim] = st
                    xsten[2, dim] =
                        findlast(temp -> temp .≤ xbndl[2, dim] .+ toler, xg)
                    n123[dim] = xsten[2, dim] - xsten[1, dim] + 1
                else
                    n123[dim] = 0
                end
            end

            if prod(n123) > 0
                for k in xsten[1, 3]:xsten[2, 3],
                    j in xsten[1, 2]:xsten[2, 2],
                    i in xsten[1, 1]:xsten[2, 1]

                    if marker[i, j, k]
                        push!(
                            ξ1[el],
                            2 * (x1g[i] - xbndl[1, 1]) /
                            (xbndl[2, 1] - xbndl[1, 1]) - 1,
                        )
                        push!(
                            ξ2[el],
                            2 * (x2g[j] - xbndl[1, 2]) /
                            (xbndl[2, 2] - xbndl[1, 2]) - 1,
                        )
                        push!(
                            ξ3[el],
                            2 * (x3g[k] - xbndl[1, 3]) /
                            (xbndl[2, 3] - xbndl[1, 3]) - 1,
                        )

                        push!(x1i[el], UInt16(i))
                        push!(x2i[el], UInt16(j))
                        push!(x3i[el], UInt16(k))
                        marker[i, j, k] = false
                    end
                end
                offset[el + 1] = offset[el] + length(ξ1[el])
            end

        end # el loop

        m_ξ1, m_ξ2, m_ξ3 = referencepoints(grid)
        wb1, wb2, wb3 = baryweights(m_ξ1), baryweights(m_ξ2), baryweights(m_ξ3)

        Npl = offset[end]

        ξ1_d = Array{FT}(undef, Npl)
        ξ2_d = Array{FT}(undef, Npl)
        ξ3_d = Array{FT}(undef, Npl)
        x1i_d = Array{UInt16}(undef, Npl)
        x2i_d = Array{UInt16}(undef, Npl)
        x3i_d = Array{UInt16}(undef, Npl)
        fac_d = zeros(FT, Npl)
        flg_d = zeros(UInt8, 3, Npl)

        for i in 1:Nel
            ctr = 1
            for j in (offset[i] + 1):offset[i + 1]
                ξ1_d[j] = ξ1[i][ctr]
                ξ2_d[j] = ξ2[i][ctr]
                ξ3_d[j] = ξ3[i][ctr]
                x1i_d[j] = x1i[i][ctr]
                x2i_d[j] = x2i[i][ctr]
                x3i_d[j] = x3i[i][ctr]
                # set up interpolation
                fac1 = FT(0)
                fac2 = FT(0)
                fac3 = FT(0)
                for ib in 1:qm[1]
                    if abs(m_ξ1[ib] - ξ1_d[j]) < toler
                        @inbounds flg_d[1, j] = UInt8(ib)
                    else
                        @inbounds fac1 += wb1[ib] / (ξ1_d[j] - m_ξ1[ib])
                    end
                end

                for ib in 1:qm[2]
                    if abs(m_ξ2[ib] - ξ2_d[j]) < toler
                        @inbounds flg_d[2, j] = UInt8(ib)
                    else
                        @inbounds fac2 += wb2[ib] / (ξ2_d[j] - m_ξ2[ib])
                    end
                end

                for ib in 1:qm[3]
                    if abs(m_ξ3[ib] - ξ3_d[j]) < toler
                        @inbounds flg_d[3, j] = UInt8(ib)
                    else
                        @inbounds fac3 += wb3[ib] / (ξ3_d[j] - m_ξ3[ib])
                    end
                end

                flg_d[1, j] ≠ UInt8(0) && (fac1 = FT(1))
                flg_d[2, j] ≠ UInt8(0) && (fac2 = FT(1))
                flg_d[3, j] ≠ UInt8(0) && (fac3 = FT(1))

                fac_d[j] = FT(1) / (fac1 * fac2 * fac3)

                ctr += 1
            end
        end
        # MPI setup for gathering data on proc 0
        root = 0
        Np_all = zeros(Int32, npr)
        Np_all[pid + 1] = Npl

        MPI.Allreduce!(Np_all, +, mpicomm)

        if pid ≠ root
            x1i_all = zeros(UInt16, 0)
            x2i_all = zeros(UInt16, 0)
            x3i_all = zeros(UInt16, 0)
            MPI.Gatherv!(x1i_d, nothing, root, mpicomm)
            MPI.Gatherv!(x2i_d, nothing, root, mpicomm)
            MPI.Gatherv!(x3i_d, nothing, root, mpicomm)
        else
            x1i_all = Array{UInt16}(undef, sum(Np_all))
            x2i_all = Array{UInt16}(undef, sum(Np_all))
            x3i_all = Array{UInt16}(undef, sum(Np_all))
            MPI.Gatherv!(x1i_d, VBuffer(x1i_all, Np_all), root, mpicomm)
            MPI.Gatherv!(x2i_d, VBuffer(x2i_all, Np_all), root, mpicomm)
            MPI.Gatherv!(x3i_d, VBuffer(x3i_all, Np_all), root, mpicomm)
        end

        if device isa CUDADevice
            ξ1_d = DA(ξ1_d)
            ξ2_d = DA(ξ2_d)
            ξ3_d = DA(ξ3_d)
            x1i_d = DA(x1i_d)
            x2i_d = DA(x2i_d)
            x3i_d = DA(x3i_d)
            flg_d = DA(flg_d)
            fac_d = DA(fac_d)
            offset = DA(offset)
            m_ξ1 = DA(m_ξ1)
            m_ξ2 = DA(m_ξ2)
            m_ξ3 = DA(m_ξ3)
            wb1 = DA(wb1)
            wb2 = DA(wb2)
            wb3 = DA(wb3)
            x1i_all = DA(x1i_all)
            x2i_all = DA(x2i_all)
            x3i_all = DA(x3i_all)
        end
        return new{
            FT,
            Int,
            typeof(x1g),
            typeof(ξ1_d),
            typeof(offset),
            typeof(xbnd),
            typeof(flg_d),
            typeof(x1i_d),
            typeof(Np_all),
        }(
            Nel,
            Np,
            Npl,
            xbnd,
            x1g,
            x2g,
            x3g,
            ξ1_d,
            ξ2_d,
            ξ3_d,
            flg_d,
            fac_d,
            x1i_d,
            x2i_d,
            x3i_d,
            offset,
            m_ξ1,
            m_ξ2,
            m_ξ3,
            wb1,
            wb2,
            wb3,
            Np_all,
            x1i_all,
            x2i_all,
            x3i_all,
        )

    end

end # struct InterpolationBrick

"""
    interpolate_local!(
        intrp_brck::InterpolationBrick{FT},
        sv::AbstractArray{FT},
        v::AbstractArray{FT},
    ) where {FT <: AbstractFloat}

This interpolation function works for a brick, where stretching/compression
happens only along the x1, x2 & x3 axis.  Here x1 = X1(ξ1), x2 = X2(ξ2) and x3
= X3(ξ3)

# Arguments
 - `intrp_brck`: Initialized InterpolationBrick structure
 - `sv`: State Array consisting of various variables on the discontinuous
   Galerkin grid
 - `v`:  Interpolated variables
"""
function interpolate_local!(
    intrp_brck::InterpolationBrick{FT},
    sv::AbstractArray{FT},
    v::AbstractArray{FT},
) where {FT <: AbstractFloat}

    offset = intrp_brck.offset
    m_ξ1 = intrp_brck.m_ξ1
    m_ξ2 = intrp_brck.m_ξ2
    m_ξ3 = intrp_brck.m_ξ3
    wb1 = intrp_brck.wb1
    wb2 = intrp_brck.wb2
    wb3 = intrp_brck.wb3
    ξ1 = intrp_brck.ξ1
    ξ2 = intrp_brck.ξ2
    ξ3 = intrp_brck.ξ3
    flg = intrp_brck.flg
    fac = intrp_brck.fac

    qm = (length(m_ξ1), length(m_ξ2), length(m_ξ3))
    Nel = length(offset) - 1
    nvars = size(sv, 2)

    device = array_device(sv)
    comp_stream = Event(device)

    workgroup = (qm[2], qm[3])
    ndrange = (qm[2] * Nel, qm[3] * nvars)

    comp_stream = interpolate_local_kernel!(device, workgroup)(
        offset,
        m_ξ1,
        m_ξ2,
        m_ξ3,
        wb1,
        wb2,
        wb3,
        ξ1,
        ξ2,
        ξ3,
        flg,
        fac,
        sv,
        v,
        Val(qm),
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )
    wait(comp_stream)
    return nothing
end
#---------------------------------------------------------------
@kernel function interpolate_local_kernel!(
    offset::AbstractArray{T, 1},
    m_ξ1::AbstractArray{FT, 1},
    m_ξ2::AbstractArray{FT, 1},
    m_ξ3::AbstractArray{FT, 1},
    wb1::AbstractArray{FT, 1},
    wb2::AbstractArray{FT, 1},
    wb3::AbstractArray{FT, 1},
    ξ1::AbstractArray{FT, 1},
    ξ2::AbstractArray{FT, 1},
    ξ3::AbstractArray{FT, 1},
    flg::AbstractArray{UInt8, 2},
    fac::AbstractArray{FT, 1},
    sv::AbstractArray{FT},
    v::AbstractArray{FT},
    ::Val{qm},
) where {qm, T <: Int, FT <: AbstractFloat}

    el, st_idx = @index(Group, NTuple)
    tj, tk = @index(Local, NTuple)

    vout_jk = @localmem FT (qm[2], qm[3])
    m_ξ1_sh = @localmem FT (qm[1],)
    m_ξ2_sh = @localmem FT (qm[2],)
    m_ξ3_sh = @localmem FT (qm[3],)
    wb1_sh = @localmem FT (qm[1],)
    wb2_sh = @localmem FT (qm[2],)
    wb3_sh = @localmem FT (qm[3],)

    np = @private T (1,)
    off = @private T (1,)

    # load shared memory
    if tk == 1
        m_ξ2_sh[tj] = m_ξ2[tj]
        wb2_sh[tj] = wb2[tj]
    end
    if tj == 1
        m_ξ3_sh[tk] = m_ξ3[tk]
        wb3_sh[tk] = wb3[tk]
    end
    if tj == 1 && tk == 1
        for i in 1:qm[1]
            m_ξ1_sh[i] = m_ξ1[i]
            wb1_sh[i] = wb1[i]
        end
    end
    @synchronize

    np[1] = offset[el + 1] - offset[el]
    off[1] = offset[el]

    for i in 1:np[1] # interpolate point-by-point

        ξ1l = ξ1[off[1] + i]
        ξ2l = ξ2[off[1] + i]

        f1 = flg[1, off[1] + i]

        if f1 == 0 # apply phir
            @inbounds vout_jk[tj, tk] =
                sv[
                    1 + (tj - 1) * qm[1] + (tk - 1) * qm[1] * qm[2],
                    st_idx,
                    el,
                ] * wb1_sh[1] / (ξ1l - m_ξ1_sh[1])
            for ii in 2:qm[1]
                @inbounds vout_jk[tj, tk] +=
                    sv[
                        ii + (tj - 1) * qm[1] + (tk - 1) * qm[1] * qm[2],
                        st_idx,
                        el,
                    ] * wb1_sh[ii] / (ξ1l - m_ξ1_sh[ii])
            end
        else
            @inbounds vout_jk[tj, tk] =
                sv[f1 + (tj - 1) * qm[1] + (tk - 1) * qm[1] * qm[2], st_idx, el]
        end

        if flg[2, off[1] + i] == 0 # apply phis
            @inbounds vout_jk[tj, tk] *= (wb2_sh[tj] / (ξ2l - m_ξ2_sh[tj]))
        end
        @synchronize

        f2 = flg[2, off[1] + i]
        ξ3l = ξ3[off[1] + i]

        if tj == 1 # reduction
            if f2 == 0
                for ij in 2:qm[2]
                    @inbounds vout_jk[1, tk] += vout_jk[ij, tk]
                end
            else
                if f2 ≠ 1
                    @inbounds vout_jk[1, tk] = vout_jk[f2, tk]
                end
            end

            if flg[3, off[1] + i] == 0 # apply phit
                @inbounds vout_jk[1, tk] *= (wb3_sh[tk] / (ξ3l - m_ξ3_sh[tk]))
            end
        end
        @synchronize

        f3 = flg[3, off[1] + i]

        if tj == 1 && tk == 1 # reduction
            if f3 == 0
                for ik in 2:qm[3]
                    @inbounds vout_jk[1, 1] += vout_jk[1, ik]
                end
            else
                if f3 ≠ 1
                    @inbounds vout_jk[1, 1] = vout_jk[1, f3]
                end
            end
            @inbounds v[off[1] + i, st_idx] = vout_jk[1, 1] * fac[off[1] + i]
        end

    end
end
#---------------------------------------------------------------
function dimensions(interpol::InterpolationBrick)
    if Array ∈ typeof(interpol.x1g).parameters
        h_x1g = interpol.x1g
        h_x2g = interpol.x2g
        h_x3g = interpol.x3g
    else
        h_x1g = Array(interpol.x1g)
        h_x2g = Array(interpol.x2g)
        h_x3g = Array(interpol.x3g)
    end
    return OrderedDict(
        "x" => (h_x1g, OrderedDict()),
        "y" => (h_x2g, OrderedDict()),
        "z" => (h_x3g, OrderedDict()),
    )
end

"""
    InterpolationCubedSphere{
    FT <: AbstractFloat,
    T <: Int,
    FTV <: AbstractVector{FT},
    FTVD <: AbstractVector{FT},
    TVD <: AbstractVector{T},
    UI8AD <: AbstractArray{UInt8, 2},
    UI16VD <: AbstractVector{UInt16},
    I32V <: AbstractVector{Int32},
    } <: InterpolationTopology

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

-90⁰  ≤ lat  ≤ 90⁰

-180⁰ ≤ long ≤ 180⁰

Rᵢ ≤ r ≤ Rₒ

# Fields

$(DocStringExtensions.FIELDS)

# Usage

    InterpolationCubedSphere(grid::DiscontinuousSpectralElementGrid, vert_range::AbstractArray{FT}, nhor::Int, lat_res::FT, long_res::FT, rad_res::FT) where {FT <: AbstractFloat}

This interpolation structure and the corresponding functions works for a cubed sphere topology. The data is interpolated along a lat/long/rad grid.

-90⁰  ≤ lat  ≤ 90⁰

-180⁰ ≤ long ≤ 180⁰

Rᵢ ≤ r ≤ Rₒ

# Arguments for the inner constructor
 - `grid`: DiscontinousSpectralElementGrid
 - `vert_range`: Vertex range along the radial coordinate
 - `lat_res`: Resolution of the interpolation grid along the latitude coordinate in radians
 - `long_res`: Resolution of the interpolation grid along the longitude coordinate in radians
 - `rad_res`: Resolution of the interpolation grid along the radial coordinate
"""
struct InterpolationCubedSphere{
    FT <: AbstractFloat,
    T <: Int,
    FTV <: AbstractVector{FT},
    FTVD <: AbstractVector{FT},
    TVD <: AbstractVector{T},
    UI8AD <: AbstractArray{UInt8, 2},
    UI16VD <: AbstractVector{UInt16},
    I32V <: AbstractVector{Int32},
} <: InterpolationTopology
    "Number of elements"
    Nel::T
    "Number of interpolation points"
    Np::T
    "Number of interpolation points on local process"
    Npl::T            # # of interpolation points on the local process
    "Number of interpolation points in radial direction"
    n_rad::T
    "Number of interpolation points in lat direction"
    n_lat::T
    "Number of interpolation points in long direction"
    n_long::T
    "Interpolation grid in radial direction"
    rad_grd::FTV
    "Interpolation grid in lat direction"
    lat_grd::FTV
    "Interpolation grid in long direction"
    long_grd::FTV # rad, lat & long locations of interpolation grid
    "Device array containing ξ1 coordinates of interpolation points within each element"
    ξ1::FTVD
    "Device array containing ξ2 coordinates of interpolation points within each element"
    ξ2::FTVD
    "Device array containing ξ3 coordinates of interpolation points within each element"
    ξ3::FTVD
    "flags when ξ1/ξ2/ξ3 interpolation point matches with a GLL point"
    flg::UI8AD
    "Normalization factor"
    fac::FTVD
    "Radial coordinates of interpolation points withing each element"
    radi::UI16VD
    "Latitude coordinates of interpolation points withing each element"
    lati::UI16VD
    "Longitude coordinates of interpolation points withing each element"
    longi::UI16VD
    "Offsets for each element"
    offset::TVD
    "GLL points in ξ1 direction"
    m_ξ1::FTVD
    "GLL points in ξ2 direction"
    m_ξ2::FTVD
    "GLL points in ξ3 direction"
    m_ξ3::FTVD
    "Barycentric weights in ξ1 direction"
    wb1::FTVD
    "Barycentric weights in ξ2 direction"
    wb2::FTVD
    "Barycentric weights in ξ3 direction"
    wb3::FTVD
    # MPI setup for gathering interpolated variable on proc 0
    "Number of interpolation points on each of the processes"
    Np_all::I32V
    "Radial interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    radi_all::UI16VD
    "Latitude interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    lati_all::UI16VD
    "Longitude interpolation grid index of interpolation points within each element on all processes stored only on proc 0"
    longi_all::UI16VD

    function InterpolationCubedSphere(
        grid::DiscontinuousSpectralElementGrid,
        vert_range::AbstractArray{FT},
        nhor::Int,
        lat_grd::AbstractArray{FT, 1},
        long_grd::AbstractArray{FT, 1},
        rad_grd::AbstractArray{FT},
    ) where {FT <: AbstractFloat}
        mpicomm = MPI.COMM_WORLD
        pid = MPI.Comm_rank(mpicomm)
        npr = MPI.Comm_size(mpicomm)

        DA = arraytype(grid)                    # device array
        device = arraytype(grid) <: Array ? CPU() : CUDADevice()

        qm = polynomialorders(grid) .+ 1
        toler1 = FT(eps(FT) * vert_range[1] * 2.0) # tolerance for unwarp function
        toler2 = FT(eps(FT) * 4.0)                 # tolerance
        toler3 = FT(eps(FT) * vert_range[1] * 10.0) # tolerance for Newton-Raphson

        Nel = length(grid.topology.realelems) # # of local elements on the local process

        nvert_range = length(vert_range)
        nvert = nvert_range - 1              # # of elements in vertical direction
        Nel_glob = nvert * nhor * nhor * 6

        nblck = nhor * nhor * nvert
        Δh = 2 / nhor                               # horizontal grid spacing in unwarped grid

        n_lat, n_long, n_rad =
            Int(length(lat_grd)), Int(length(long_grd)), Int(length(rad_grd))

        Np = n_lat * n_long * n_rad

        uw_grd = zeros(FT, 3)
        diffv = zeros(FT, 3)
        ξ = zeros(FT, 3)

        glob_ord = grid.topology.origsendorder # to account for reordering of elements after the partitioning process

        glob_elem_no = zeros(Int, nvert * length(glob_ord))

        for i in 1:length(glob_ord), j in 1:nvert
            glob_elem_no[j + (i - 1) * nvert] = (glob_ord[i] - 1) * nvert + j
        end
        glob_to_loc = Dict(glob_elem_no[i] => Int(i) for i in 1:Nel) # using dictionary for speedup

        ξ1, ξ2, ξ3 = map(i -> zeros(FT, i), zeros(Int, Nel)),
        map(i -> zeros(FT, i), zeros(Int, Nel)),
        map(i -> zeros(FT, i), zeros(Int, Nel))

        radi, lati, longi = map(i -> zeros(UInt16, i), zeros(UInt16, Nel)),
        map(i -> zeros(UInt16, i), zeros(UInt16, Nel)),
        map(i -> zeros(UInt16, i), zeros(UInt16, Nel))

        offset_d = zeros(Int, Nel + 1)

        for i in 1:n_rad
            rad = rad_grd[i]
            if rad ≤ vert_range[1]       # accounting for minor rounding errors from unwarp function at boundaries
                vert_range[1] - rad < toler1 ? l_nrm = 1 :
                error(
                    "fatal error, rad lower than inner radius: ",
                    vert_range[1] - rad,
                    " $rad_grd /// $lat_grd //// $long_grd",
                )
            elseif rad ≥ vert_range[end] # accounting for minor rounding errors from unwarp function at boundaries
                rad - vert_range[end] < toler1 ? l_nrm = nvert :
                error("fatal error, rad greater than outer radius")
            else                         # normal scenario
                for l in 2:nvert_range
                    if vert_range[l] - rad > FT(0)
                        l_nrm = l - 1
                        break
                    end
                end
            end

            for j in 1:n_lat
                @inbounds x3_grd = rad * sind(lat_grd[j])
                for k in 1:n_long
                    @inbounds x1_grd =
                        rad * cosd(lat_grd[j]) * cosd(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.
                    @inbounds x2_grd =
                        rad * cosd(lat_grd[j]) * sind(long_grd[k]) # inclination -> latitude; azimuthal -> longitude.

                    uw_grd[1], uw_grd[2], uw_grd[3] =
                        Topologies.cubedshellunwarp(x1_grd, x2_grd, x3_grd) # unwarping from sphere to cubed shell

                    x1_uw2_grd = uw_grd[1] / rad # unwrapping cubed shell on to a 2D grid (in 3D space, -1 to 1 cube)
                    x2_uw2_grd = uw_grd[2] / rad
                    x3_uw2_grd = uw_grd[3] / rad

                    if abs(x1_uw2_grd + 1) < toler2 # face 1 (x1 == -1 plane)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (nhor - l2) * nvert +
                            (l3 - 1) * nvert * nhor,
                        )
                    elseif abs(x2_uw2_grd + 1) < toler2 # face 2 (x2 == -1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (l1 - 1) * nvert +
                            (l3 - 1) * nvert * nhor +
                            nblck * 1,
                        )
                    elseif abs(x1_uw2_grd - 1) < toler2 # face 3 (x1 == +1 plane)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (l2 - 1) * nvert +
                            (l3 - 1) * nvert * nhor +
                            nblck * 2,
                        )
                    elseif abs(x3_uw2_grd - 1) < toler2 # face 4 (x3 == +1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (l1 - 1) * nvert +
                            (l2 - 1) * nvert * nhor +
                            nblck * 3,
                        )
                    elseif abs(x2_uw2_grd - 1) < toler2 # face 5 (x2 == +1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l3 = min(div(x3_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (l1 - 1) * nvert +
                            (nhor - l3) * nvert * nhor +
                            nblck * 4,
                        )
                    elseif abs(x3_uw2_grd + 1) < toler2 # face 6 (x3 == -1 plane)
                        l1 = min(div(x1_uw2_grd + 1, Δh) + 1, nhor)
                        l2 = min(div(x2_uw2_grd + 1, Δh) + 1, nhor)
                        el_glob = Int(
                            l_nrm +
                            (l1 - 1) * nvert +
                            (nhor - l2) * nvert * nhor +
                            nblck * 5,
                        )
                    else
                        error("error: unwrapped grid does not lie on any of the 6 faces")
                    end

                    el_loc = get(glob_to_loc, el_glob, nothing)
                    if el_loc ≠ nothing # computing inner coordinates for local elements
                        invert_trilear_mapping_hex!(
                            view(grid.topology.elemtocoord, 1, :, el_loc),
                            view(grid.topology.elemtocoord, 2, :, el_loc),
                            view(grid.topology.elemtocoord, 3, :, el_loc),
                            uw_grd,
                            diffv,
                            toler3,
                            ξ,
                        )
                        push!(ξ1[el_loc], ξ[1])
                        push!(ξ2[el_loc], ξ[2])
                        push!(ξ3[el_loc], ξ[3])
                        push!(radi[el_loc], UInt16(i))
                        push!(lati[el_loc], UInt16(j))
                        push!(longi[el_loc], UInt16(k))
                        offset_d[el_loc + 1] += 1
                    end
                end
            end
        end

        for i in 2:(Nel + 1)
            @inbounds offset_d[i] += offset_d[i - 1]
        end

        Npl = offset_d[Nel + 1]

        v = Vector{FT}(undef, offset_d[Nel + 1]) # Allocating storage for interpolation variable

        ξ1_d = Vector{FT}(undef, Npl)
        ξ2_d = Vector{FT}(undef, Npl)
        ξ3_d = Vector{FT}(undef, Npl)

        flg_d = zeros(UInt8, 3, Npl)
        fac_d = ones(FT, Npl)

        rad_d = Vector{UInt16}(undef, Npl)
        lat_d = Vector{UInt16}(undef, Npl)
        long_d = Vector{UInt16}(undef, Npl)

        m_ξ1, m_ξ2, m_ξ3 = referencepoints(grid)
        wb1, wb2, wb3 = baryweights(m_ξ1), baryweights(m_ξ2), baryweights(m_ξ3)
        for i in 1:Nel
            ctr = 1
            for j in (offset_d[i] + 1):offset_d[i + 1]
                @inbounds ξ1_d[j] = ξ1[i][ctr]
                @inbounds ξ2_d[j] = ξ2[i][ctr]
                @inbounds ξ3_d[j] = ξ3[i][ctr]
                @inbounds rad_d[j] = radi[i][ctr]
                @inbounds lat_d[j] = lati[i][ctr]
                @inbounds long_d[j] = longi[i][ctr]
                # set up interpolation
                fac1 = FT(0)
                fac2 = FT(0)
                fac3 = FT(0)
                for ib in 1:qm[1]
                    if abs(m_ξ1[ib] - ξ1_d[j]) < toler2
                        @inbounds flg_d[1, j] = UInt8(ib)
                    else
                        @inbounds fac1 += wb1[ib] / (ξ1_d[j] - m_ξ1[ib])
                    end
                end

                for ib in 1:qm[2]
                    if abs(m_ξ2[ib] - ξ2_d[j]) < toler2
                        @inbounds flg_d[2, j] = UInt8(ib)
                    else
                        @inbounds fac2 += wb2[ib] / (ξ2_d[j] - m_ξ2[ib])
                    end
                end

                for ib in 1:qm[3]
                    if abs(m_ξ3[ib] - ξ3_d[j]) < toler2
                        @inbounds flg_d[3, j] = UInt8(ib)
                    else
                        @inbounds fac3 += wb3[ib] / (ξ3_d[j] - m_ξ3[ib])
                    end
                end

                flg_d[1, j] ≠ 0 && (fac1 = FT(1))
                flg_d[2, j] ≠ 0 && (fac2 = FT(1))
                flg_d[3, j] ≠ 0 && (fac3 = FT(1))

                fac_d[j] = FT(1) / (fac1 * fac2 * fac3)

                ctr += 1
            end
        end
        # MPI setup for gathering data on proc 0
        root = 0
        Np_all = zeros(Int32, npr)
        Np_all[pid + 1] = Int32(Npl)

        MPI.Allreduce!(Np_all, +, mpicomm)

        if pid ≠ root
            radi_all = zeros(UInt16, 0)
            lati_all = zeros(UInt16, 0)
            longi_all = zeros(UInt16, 0)
            MPI.Gatherv!(rad_d, nothing, root, mpicomm)
            MPI.Gatherv!(lat_d, nothing, root, mpicomm)
            MPI.Gatherv!(long_d, nothing, root, mpicomm)
        else
            radi_all = Array{UInt16}(undef, sum(Np_all))
            lati_all = Array{UInt16}(undef, sum(Np_all))
            longi_all = Array{UInt16}(undef, sum(Np_all))
            MPI.Gatherv!(rad_d, VBuffer(radi_all, Np_all), root, mpicomm)
            MPI.Gatherv!(lat_d, VBuffer(lati_all, Np_all), root, mpicomm)
            MPI.Gatherv!(long_d, VBuffer(longi_all, Np_all), root, mpicomm)
        end

        if device isa CUDADevice
            ξ1_d = DA(ξ1_d)
            ξ2_d = DA(ξ2_d)
            ξ3_d = DA(ξ3_d)

            flg_d = DA(flg_d)
            fac_d = DA(fac_d)

            rad_d = DA(rad_d)
            lat_d = DA(lat_d)
            long_d = DA(long_d)

            m_ξ1 = DA(m_ξ1)
            m_ξ2 = DA(m_ξ2)
            m_ξ3 = DA(m_ξ3)
            wb1 = DA(wb1)
            wb2 = DA(wb2)
            wb3 = DA(wb3)

            offset_d = DA(offset_d)

            rad_grd = DA(rad_grd)
            lat_grd = DA(lat_grd)
            long_grd = DA(long_grd)

            radi_all = DA(radi_all)
            lati_all = DA(lati_all)
            longi_all = DA(longi_all)
        end

        return new{
            FT,
            Int,
            typeof(rad_grd),
            typeof(ξ1_d),
            typeof(offset_d),
            typeof(flg_d),
            typeof(rad_d),
            typeof(Np_all),
        }(
            Nel,
            Np,
            Npl,
            n_rad,
            n_lat,
            n_long,
            rad_grd,
            lat_grd,
            long_grd,
            ξ1_d,
            ξ2_d,
            ξ3_d,
            flg_d,
            fac_d,
            rad_d,
            lat_d,
            long_d,
            offset_d,
            m_ξ1,
            m_ξ2,
            m_ξ3,
            wb1,
            wb2,
            wb3,
            Np_all,
            radi_all,
            lati_all,
            longi_all,
        )

    end # Inner constructor InterpolationCubedSphere

end # struct InterpolationCubedSphere

"""
    invert_trilear_mapping_hex!(X1::AbstractArray{FT,1},
                                X2::AbstractArray{FT,1},
                                X3::AbstractArray{FT,1},
                                 x::AbstractArray{FT,1},
                                 d::AbstractArray{FT,1},
                               tol::FT,
                                 ξ::AbstractArray{FT,1}) where FT <: AbstractFloat

This function computes ξ = (ξ1,ξ2,ξ3) given x = (x1,x2,x3) and the (8) vertex coordinates of a Hexahedron. Newton-Raphson method is used.

# Arguments
 - `X1`: X1 coordinates of the (8) vertices of the hexahedron
 - `X2`: X2 coordinates of the (8) vertices of the hexahedron
 - `X3`: X3 coordinates of the (8) vertices of the hexahedron
 - `x`: (x1,x2,x3) coordinates of the point
 - `d`: (x1,x2,x3) coordinates, temporary storage
 - `tol`: Tolerance for convergence
 - `ξ`: (ξ1,ξ2,ξ3) coordinates of the point
"""
function invert_trilear_mapping_hex!(
    X1::AbstractArray{FT, 1},
    X2::AbstractArray{FT, 1},
    X3::AbstractArray{FT, 1},
    x::AbstractArray{FT, 1},
    d::AbstractArray{FT, 1},
    tol::FT,
    ξ::AbstractArray{FT, 1},
) where {FT <: AbstractFloat}
    max_it = 10     # maximum # of iterations
    ξ .= FT(0) #zeros(FT,3,1) # initial guess => cell centroid
    trilinear_map_minus_x!(ξ, X1, X2, X3, x, d)
    err = sqrt(d[1] * d[1] + d[2] * d[2] + d[3] * d[3])
    ctr = 0
    # Newton-Raphson iterations
    while err > tol
        trilinear_map_IJac_x_vec!(ξ, X1, X2, X3, d)
        ξ .-= d
        trilinear_map_minus_x!(ξ, X1, X2, X3, x, d)
        err = sqrt(d[1] * d[1] + d[2] * d[2] + d[3] * d[3]) #norm(d)
        ctr += 1
        if ctr > max_it
            error(
                "invert_trilinear_mapping_hex: Newton-Raphson not converging to desired tolerance after max_it = ",
                max_it,
                " iterations; err = ",
                err,
                "; toler = ",
                tol,
            )
        end
    end

    clamp!(ξ, FT(-1), FT(1))
    return nothing
end

function trilinear_map_minus_x!(
    ξ::AbstractArray{FT, 1},
    x1v::AbstractArray{FT, 1},
    x2v::AbstractArray{FT, 1},
    x3v::AbstractArray{FT, 1},
    x::AbstractArray{FT, 1},
    d::AbstractArray{FT, 1},
) where {FT <: AbstractFloat}
    p1 = 1 + ξ[1]
    p2 = 1 + ξ[2]
    p3 = 1 + ξ[3]
    m1 = 1 - ξ[1]
    m2 = 1 - ξ[2]
    m3 = 1 - ξ[3]

    d[1] =
        (
            m1 * (
                m2 * (m3 * x1v[1] + p3 * x1v[5]) +
                p2 * (m3 * x1v[3] + p3 * x1v[7])
            ) +
            p1 * (
                m2 * (m3 * x1v[2] + p3 * x1v[6]) +
                p2 * (m3 * x1v[4] + p3 * x1v[8])
            )
        ) / 8.0 - x[1]

    d[2] =
        (
            m1 * (
                m2 * (m3 * x2v[1] + p3 * x2v[5]) +
                p2 * (m3 * x2v[3] + p3 * x2v[7])
            ) +
            p1 * (
                m2 * (m3 * x2v[2] + p3 * x2v[6]) +
                p2 * (m3 * x2v[4] + p3 * x2v[8])
            )
        ) / 8.0 - x[2]

    d[3] =
        (
            m1 * (
                m2 * (m3 * x3v[1] + p3 * x3v[5]) +
                p2 * (m3 * x3v[3] + p3 * x3v[7])
            ) +
            p1 * (
                m2 * (m3 * x3v[2] + p3 * x3v[6]) +
                p2 * (m3 * x3v[4] + p3 * x3v[8])
            )
        ) / 8.0 - x[3]

    return nothing
end

function trilinear_map_IJac_x_vec!(
    ξ::AbstractArray{FT, 1},
    x1v::AbstractArray{FT, 1},
    x2v::AbstractArray{FT, 1},
    x3v::AbstractArray{FT, 1},
    v::AbstractArray{FT, 1},
) where {FT <: AbstractFloat}
    p1 = 1 + ξ[1]
    p2 = 1 + ξ[2]
    p3 = 1 + ξ[3]
    m1 = 1 - ξ[1]
    m2 = 1 - ξ[2]
    m3 = 1 - ξ[3]

    Jac11 =
        (
            m2 * (m3 * (x1v[2] - x1v[1]) + p3 * (x1v[6] - x1v[5])) +
            p2 * (m3 * (x1v[4] - x1v[3]) + p3 * (x1v[8] - x1v[7]))
        ) / 8.0

    Jac12 =
        (
            m1 * (m3 * (x1v[3] - x1v[1]) + p3 * (x1v[7] - x1v[5])) +
            p1 * (m3 * (x1v[4] - x1v[2]) + p3 * (x1v[8] - x1v[6]))
        ) / 8.0

    Jac13 =
        (
            m1 * (m2 * (x1v[5] - x1v[1]) + p2 * (x1v[7] - x1v[3])) +
            p1 * (m2 * (x1v[6] - x1v[2]) + p2 * (x1v[8] - x1v[4]))
        ) / 8.0

    Jac21 =
        (
            m2 * (m3 * (x2v[2] - x2v[1]) + p3 * (x2v[6] - x2v[5])) +
            p2 * (m3 * (x2v[4] - x2v[3]) + p3 * (x2v[8] - x2v[7]))
        ) / 8.0

    Jac22 =
        (
            m1 * (m3 * (x2v[3] - x2v[1]) + p3 * (x2v[7] - x2v[5])) +
            p1 * (m3 * (x2v[4] - x2v[2]) + p3 * (x2v[8] - x2v[6]))
        ) / 8.0

    Jac23 =
        (
            m1 * (m2 * (x2v[5] - x2v[1]) + p2 * (x2v[7] - x2v[3])) +
            p1 * (m2 * (x2v[6] - x2v[2]) + p2 * (x2v[8] - x2v[4]))
        ) / 8.0

    Jac31 =
        (
            m2 * (m3 * (x3v[2] - x3v[1]) + p3 * (x3v[6] - x3v[5])) +
            p2 * (m3 * (x3v[4] - x3v[3]) + p3 * (x3v[8] - x3v[7]))
        ) / 8.0

    Jac32 =
        (
            m1 * (m3 * (x3v[3] - x3v[1]) + p3 * (x3v[7] - x3v[5])) +
            p1 * (m3 * (x3v[4] - x3v[2]) + p3 * (x3v[8] - x3v[6]))
        ) / 8.0

    Jac33 =
        (
            m1 * (m2 * (x3v[5] - x3v[1]) + p2 * (x3v[7] - x3v[3])) +
            p1 * (m2 * (x3v[6] - x3v[2]) + p2 * (x3v[8] - x3v[4]))
        ) / 8.0

    # computing cofactor matrix
    C11 = Jac22 * Jac33 - Jac23 * Jac32
    C12 = -Jac21 * Jac33 + Jac23 * Jac31
    C13 = Jac21 * Jac32 - Jac22 * Jac31
    C21 = -Jac12 * Jac33 + Jac13 * Jac32
    C22 = Jac11 * Jac33 - Jac13 * Jac31
    C23 = -Jac11 * Jac32 + Jac12 * Jac31
    C31 = Jac12 * Jac23 - Jac13 * Jac22
    C32 = -Jac11 * Jac23 + Jac13 * Jac21
    C33 = Jac11 * Jac22 - Jac12 * Jac21

    # computing determinant
    det = Jac11 * C11 + Jac12 * C12 + Jac13 * C13

    Jac11 = (C11 * v[1] + C21 * v[2] + C31 * v[3]) / det
    Jac21 = (C12 * v[1] + C22 * v[2] + C32 * v[3]) / det
    Jac31 = (C13 * v[1] + C23 * v[2] + C33 * v[3]) / det

    v[1] = Jac11
    v[2] = Jac21
    v[3] = Jac31

    return nothing
end

"""
    interpolate_local!(intrp_cs::InterpolationCubedSphere{FT},
                             sv::AbstractArray{FT},
                              v::AbstractArray{FT}) where {FT <: AbstractFloat}

This interpolation function works for cubed spherical shell geometry.

# Arguments
 - `intrp_cs`: Initialized cubed sphere structure
 - `sv`: Array consisting of various variables on the discontinuous Galerkin grid
 - `v`:  Array consisting of variables on the interpolated grid
"""
function interpolate_local!(
    intrp_cs::InterpolationCubedSphere{FT},
    sv::AbstractArray{FT},
    v::AbstractArray{FT},
) where {FT <: AbstractFloat}

    offset = intrp_cs.offset
    m_ξ1 = intrp_cs.m_ξ1
    m_ξ2 = intrp_cs.m_ξ2
    m_ξ3 = intrp_cs.m_ξ3
    wb1 = intrp_cs.wb1
    wb2 = intrp_cs.wb2
    wb3 = intrp_cs.wb3
    ξ1 = intrp_cs.ξ1
    ξ2 = intrp_cs.ξ2
    ξ3 = intrp_cs.ξ3
    flg = intrp_cs.flg
    fac = intrp_cs.fac

    qm = (length(m_ξ1), length(m_ξ2), length(m_ξ3))
    nvars = size(sv, 2)
    Nel = length(offset) - 1
    np_tot = size(v, 1)

    device = array_device(sv)
    comp_stream = Event(device)

    workgroup = (qm[2], qm[3])
    ndrange = (qm[2] * Nel, qm[3] * nvars)

    comp_stream = interpolate_local_kernel!(device, workgroup)(
        offset,
        m_ξ1,
        m_ξ2,
        m_ξ3,
        wb1,
        wb2,
        wb3,
        ξ1,
        ξ2,
        ξ3,
        flg,
        fac,
        sv,
        v,
        Val(qm),
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )
    wait(comp_stream)

    return nothing
end

"""
    project_cubed_sphere!(intrp_cs::InterpolationCubedSphere{FT},
                                 v::AbstractArray{FT},
                              uvwi::Tuple{Int,Int,Int}) where {FT <: AbstractFloat}

This function projects the velocity field along unit vectors in radial, lat and long directions for cubed spherical shell geometry.

# Fields
 - `intrp_cs`: Initialized cubed sphere structure
 - `v`: Array consisting of x1, x2 and x3 components of the vector field
 - `uvwi`:  Tuple providing the column numbers for x1, x2 and x3 components of vector field in the array.
            These columns will be replaced with projected vector fields along unit vectors in long, lat and rad directions.
"""
function project_cubed_sphere!(
    intrp_cs::InterpolationCubedSphere{FT},
    v::AbstractArray{FT},
    uvwi::Tuple{Int, Int, Int},
) where {FT <: AbstractFloat}
    # projecting velocity onto unit vectors in long, lat and rad directions
    # assumes u, v and w are located in columns specified in vector uvwi
    @assert length(uvwi) == 3 "length(uvwi) is not 3"
    lati = intrp_cs.lati
    longi = intrp_cs.longi
    lat_grd = intrp_cs.lat_grd
    long_grd = intrp_cs.long_grd
    _ρu = uvwi[1]
    _ρv = uvwi[2]
    _ρw = uvwi[3]
    np_tot = size(v, 1)

    device = array_device(v)
    comp_stream = Event(device)

    thr_x = min(256, np_tot)

    workgroup = (thr_x,)
    ndrange = (np_tot,)

    comp_stream = project_cubed_sphere_kernel!(device, workgroup)(
        lat_grd,
        long_grd,
        lati,
        longi,
        v,
        _ρu,
        _ρv,
        _ρw,
        ndrange = ndrange,
        dependencies = (comp_stream,),
    )
    wait(comp_stream)
    return nothing
end

@kernel function project_cubed_sphere_kernel!(
    lat_grd::AbstractArray{FT, 1},
    long_grd::AbstractArray{FT, 1},
    lati::AbstractVector{UInt16},
    longi::AbstractVector{UInt16},
    v::AbstractArray{FT},
    _ρu::Int,
    _ρv::Int,
    _ρw::Int,
) where {FT <: AbstractFloat}

    idx = @index(Global, Linear) # global thread ids
    # projecting velocity onto unit vectors in long, lat and rad directions
    # assumed u, v and w are located in columns 2, 3 and 4
    deg2rad = FT(π) / FT(180)
    vrad =
        v[idx, _ρu] *
        cos(lat_grd[lati[idx]] * deg2rad) *
        cos(long_grd[longi[idx]] * deg2rad) +
        v[idx, _ρv] *
        cos(lat_grd[lati[idx]] * deg2rad) *
        sin(long_grd[longi[idx]] * deg2rad) +
        v[idx, _ρw] * sin(lat_grd[lati[idx]] * deg2rad)

    vlat =
        -v[idx, _ρu] *
        sin(lat_grd[lati[idx]] * deg2rad) *
        cos(long_grd[longi[idx]] * deg2rad) -
        v[idx, _ρv] *
        sin(lat_grd[lati[idx]] * deg2rad) *
        sin(long_grd[longi[idx]] * deg2rad) +
        v[idx, _ρw] * cos(lat_grd[lati[idx]] * deg2rad)

    vlon =
        -v[idx, _ρu] * sin(long_grd[longi[idx]] * deg2rad) +
        v[idx, _ρv] * cos(long_grd[longi[idx]] * deg2rad)

    v[idx, _ρu] = vlon
    v[idx, _ρv] = vlat
    v[idx, _ρw] = vrad
    # TODO: cosd / sind having issues on GPU. Unable to isolate the issue at this point. Needs to be revisited.
end

function dimensions(interpol::InterpolationCubedSphere)
    if Array ∈ typeof(interpol.rad_grd).parameters
        h_long_grd = interpol.long_grd
        h_lat_grd = interpol.lat_grd
        h_rad_grd = interpol.rad_grd
    else
        h_long_grd = Array(interpol.long_grd)
        h_lat_grd = Array(interpol.lat_grd)
        h_rad_grd = Array(interpol.rad_grd)
    end
    FT = eltype(h_rad_grd)
    return OrderedDict(
        "long" => (
            h_long_grd,
            OrderedDict("units" => "degrees_east", "long_name" => "longitude"),
        ),
        "lat" => (
            h_lat_grd,
            OrderedDict("units" => "degrees_north", "long_name" => "latitude"),
        ),
        "level" =>
            (h_rad_grd, OrderedDict("units" => "m", "long_name" => "level")),
    )
end

"""
    accumulate_interpolated_data!(intrp::InterpolationTopology,
                                     iv::AbstractArray{FT,2},
                                    fiv::AbstractArray{FT,4}) where {FT <: AbstractFloat}

This interpolation function gathers interpolated data onto process # 0.

# Fields
 - `intrp`: Initialized interpolation topology structure
 - `iv`: Interpolated variables on local process
 - `fiv`: Full interpolated variables accumulated on process # 0
"""
function accumulate_interpolated_data!(
    intrp::InterpolationTopology,
    iv::AbstractArray{FT, 2},
    fiv::AbstractArray{FT, 4},
) where {FT <: AbstractFloat}

    device = array_device(iv)
    mpicomm = MPI.COMM_WORLD
    pid = MPI.Comm_rank(mpicomm)
    npr = MPI.Comm_size(mpicomm)
    root = 0
    nvars = size(iv, 2)

    if intrp isa InterpolationCubedSphere
        nx1 = length(intrp.long_grd)
        nx2 = length(intrp.lat_grd)
        nx3 = length(intrp.rad_grd)
        np_tot = length(intrp.radi_all)
        i1 = intrp.longi_all
        i2 = intrp.lati_all
        i3 = intrp.radi_all
    elseif intrp isa InterpolationBrick
        nx1 = length(intrp.x1g)
        nx2 = length(intrp.x2g)
        nx3 = length(intrp.x3g)
        np_tot = length(intrp.x1i_all)
        i1 = intrp.x1i_all
        i2 = intrp.x2i_all
        i3 = intrp.x3i_all
    else
        error("Unsupported topology; only InterpolationCubedSphere and InterpolationBrick supported")
    end

    if pid == 0 && size(fiv) ≠ (nx1, nx2, nx3, nvars)
        error("size of fiv = $(size(fiv)); which does not match with ($nx1,$nx2,$nx3,$nvars) ")
    end

    if npr > 1
        Np_all = intrp.Np_all
        pid == 0 ? v_all = Array{FT}(undef, np_tot, nvars) :
        v_all = Array{FT}(undef, 0, nvars)
        if device isa CPU

            for vari in 1:nvars
                MPI.Gatherv!(
                    view(iv, :, vari),
                    pid == root ? VBuffer(view(v_all, :, vari), Np_all) :
                    nothing,
                    root,
                    mpicomm,
                )
            end

        elseif device isa CUDADevice

            v = Array(iv)
            for vari in 1:nvars
                MPI.Gatherv!(
                    view(v, :, vari),
                    pid == root ? VBuffer(view(v_all, :, vari), Np_all) :
                    nothing,
                    root,
                    mpicomm,
                )
            end
            v_all = CuArray(v_all)

        else
            error("accumulate_interpolate_data: unsupported device, only CPU() and CUDADevice() supported")
        end
    else
        v_all = iv
    end

    if pid == 0
        comp_stream = Event(device)
        thr_x = min(256, np_tot)
        workgroup = (thr_x,)
        ndrange = (np_tot,)

        comp_stream = accumulate_helper_kernel!(device, workgroup)(
            i1,
            i2,
            i3,
            v_all,
            fiv,
            ndrange = ndrange,
            dependencies = (comp_stream,),
        )
        wait(comp_stream)
    end
    MPI.Barrier(mpicomm)
    return nothing
end

@kernel function accumulate_helper_kernel!(
    i1::AbstractArray{UInt16, 1},
    i2::AbstractArray{UInt16, 1},
    i3::AbstractArray{UInt16, 1},
    v_all::AbstractArray{FT, 2},
    fiv::AbstractArray{FT, 4},
) where {FT <: AbstractFloat}
    idx = @index(Global, Linear)
    nvars = size(v_all, 2)

    for vari in 1:nvars
        @inbounds fiv[i1[idx], i2[idx], i3[idx], vari] = v_all[idx, vari]
    end
end

function accumulate_interpolated_data(
    mpicomm::MPI.Comm,
    intrp::InterpolationTopology,
    iv::AbstractArray{FT, 2},
) where {FT <: AbstractFloat}

    mpirank = MPI.Comm_rank(mpicomm)
    numranks = MPI.Comm_size(mpicomm)
    nvars = size(iv, 2)

    if intrp isa InterpolationCubedSphere
        nx1 = length(intrp.long_grd)
        nx2 = length(intrp.lat_grd)
        nx3 = length(intrp.rad_grd)
        np_tot = length(intrp.radi_all)
        i1 = intrp.longi_all
        i2 = intrp.lati_all
        i3 = intrp.radi_all
    elseif intrp isa InterpolationBrick
        nx1 = length(intrp.x1g)
        nx2 = length(intrp.x2g)
        nx3 = length(intrp.x3g)
        np_tot = length(intrp.x1i_all)
        i1 = intrp.x1i_all
        i2 = intrp.x2i_all
        i3 = intrp.x3i_all
    else
        error("Unsupported topology; only InterpolationCubedSphere and InterpolationBrick supported")
    end

    if array_device(iv) isa CPU
        h_iv = iv
        h_i1 = i1
        h_i2 = i2
        h_i3 = i3
    else
        h_iv = Array(iv)
        h_i1 = Array(i1)
        h_i2 = Array(i2)
        h_i3 = Array(i3)
    end

    if numranks == 1
        v_all = h_iv
    else
        v_all = Array{FT}(undef, mpirank == 0 ? np_tot : 0, nvars)
        for vari in 1:nvars
            MPI.Gatherv!(
                view(h_iv, :, vari),
                mpirank == 0 ? VBuffer(view(v_all, :, vari), intrp.Np_all) :
                nothing,
                0,
                mpicomm,
            )
        end
    end

    if mpirank == 0
        fiv = Array{FT}(undef, nx1, nx2, nx3, nvars)
        for i in 1:np_tot
            for vari in 1:nvars
                @inbounds fiv[h_i1[i], h_i2[i], h_i3[i], vari] = v_all[i, vari]
            end
        end
    else
        fiv = nothing
    end

    return fiv
end

end # module Interpolation

max_forks_repo_name: jm-c/CLIMA


text:
#
# Miscellaneous functions
#

# Redefine unsafe_load, unsafe_store!, read, and write so that they pack the type
function define_packed(ty::DataType)
    @assert isbits(ty)
    packed_offsets = cumsum([sizeof(x) for x in ty.types])
    sz = pop!(packed_offsets)
    unshift!(packed_offsets, 0)

    if sz != sizeof(ty)
        @eval begin
            function Base.unsafe_store!(p::Ptr{$ty}, x::$ty)
                $([:(unsafe_store!(convert(Ptr{$(ty.types[i])}, p+$(packed_offsets[i])), getfield(x, $i)))
                   for i = 1:length(packed_offsets)]...)
            end
            function Base.unsafe_load(p::Ptr{$ty})
                $(Expr(:new, ty, [:(unsafe_load(convert(Ptr{$(ty.types[i])}, p+$(packed_offsets[i]))))
                                   for i = 1:length(packed_offsets)]...))
            end
            Base.sizeof(::Union{$ty,Type{$ty}}) = $sz
        end
    end

    @eval begin
        @inline Base.write(io::MmapIO, x::$ty) = _write(io, x)
        @inline Base.read(io::MmapIO, x::Type{$ty}) = _read(io, x)
        function Base.read(io::IO, ::Type{$ty})
            $(Expr(:new, ty, [:(read(io, $(ty.types[i]))) for i = 1:length(packed_offsets)]...))
        end
        function Base.write(io::IO, x::$ty)
            $([:(write(io, getfield(x, $i))) for i = 1:length(packed_offsets)]...)
            nothing
        end
    end
    nothing
end

# Loads a variable-length size according to flags
# Expects that the first two bits of flags mean:
# 0   The size of the Length of Link Name field is 1 byte.
# 1   The size of the Length of Link Name field is 2 bytes.
# 2   The size of the Length of Link Name field is 4 bytes.
# 3   The size of the Length of Link Name field is 8 bytes.
# Returns the size as a UInt and a new pointer that is offset by the
# size of the size field
function read_size(io::IO, flags::UInt8)
    if (flags & 1) == 0 && (flags & 2) == 0
        Int(read(io, UInt8))
    elseif (flags & 1) == 1 && (flags & 2) == 0
        Int(read(io, UInt16))
    elseif (flags & 1) == 0 && (flags & 2) == 2
        Int(read(io, UInt32))
    else
        Int(read(io, UInt64))
    end
end

# Determine what the size flag should be
# Same rules as above
function size_flag(sz::Integer)
    if sz <= typemax(UInt8)
        UInt8(0)
    elseif sz <= typemax(UInt16)
        UInt8(1)
    elseif sz <= typemax(UInt32)
        UInt8(2)
    else
        UInt8(3)
    end
end

# Store a size
function write_size(io::IO, sz::Integer)
    if sz <= typemax(UInt8)
        write(io, UInt8(sz))
    elseif sz <= typemax(UInt16)
        write(io, UInt16(sz))
    elseif sz <= typemax(UInt32)
        write(io, UInt32(sz))
    else
        write(io, UInt64(sz))
    end
end

# Get the size of the size
function size_size(sz::Integer)
    if sz <= typemax(UInt8)
        1
    elseif sz <= typemax(UInt16)
        2
    elseif sz <= typemax(UInt32)
        4
    else
        8
    end
end


text:
# This file contains utility functions to deal with distributed
# computing

"""
    update_best(channel::RemoteChannel, bestx::SharedArray{Float64, 1})

Listen to a `channel` for results found by lmlovo. If there is an
improvement for the objective function, the shared array `bestx` is
updated.

**Attention**: There might be an unstable state if there is a process
  reading `bestx` while this function is updating it. This should not
  be a problem, since it is used as a starting point.

**Attention 2**: this function is currently out of use.

"""
function update_best(channel::RemoteChannel, bestx::SharedArray{Float64, 1})

    @debug("Running updater.")

    n = length(bestx)

    N::Int = 0

    while isopen(channel)

        θ = try

            take!(channel)

        catch e

            if isa(e, InvalidStateException)
                break
            end

            @warn("Something wrong when reading from channel. Will skip.", e)

            continue

        end

        @debug("Updater has read values from channel.")

        for i = 1:n

            bestx[i] = (N * bestx[i] + θ[i]) / (N + 1)

        end

        N += 1

    end

    @debug("Update channel closed. Exiting thread.")

end

"""

    function consume_tqueue(bqueue::RemoteChannel, tqueue::RemoteChannel,
                            squeue::RemoteChannel, model::Function, gmodel!::Function,
                            data::Array{Float64, 2}, n::Int, pliminf::Int,
                            plimsup::Int, MAXMS::Int, seedMS::MersenneTwister)

This function represents one worker, which runs lmlovo in a multistart
fashion.

It takes a job from the RemoteChannel `tqueue` and runs `lmlovo`
function to it. It might run using a multistart strategy, if
`MAXMS>1`. It sends the best results found for each value obtained in
`tqueue` to channel `squeue`, which will be consumed by the main
process. All the other arguments are the same for [`praff`](@ref)
function.

"""
function consume_tqueue(bqueue::RemoteChannel, tqueue::RemoteChannel,
                        squeue::RemoteChannel,
                        model::Function, gmodel!::Function,
                        data::Array{Float64, 2}, n::Int, pliminf::Int,
                        plimsup::Int, MAXMS::Int,
                        seedMS::MersenneTwister, initguess::Vector{Float64})

    @debug("Started worker $(myid())")

    while isopen(tqueue)

        p = try

            take!(tqueue)

        catch e

            if isa(e, InvalidStateException)
                break
            end

            @warn("Something wrong when reading task. Will skip task.", e)

            continue

        end

        @debug("Received task $(p)")

        if (p.start < pliminf) || (p.stop > plimsup) ||
           (length(p) == 0)

            @warn("Invalid value for task: $(p). Will skip task.")

            continue

        end

        for k in p

            wbest = RAFFOutput(k)

            nf = 0
            nj = 0
            ni = 0

            # Multi-start strategy
            for j = 1:MAXMS

                # New random starting point
                θ = randn(seedMS, n)
                θ .= θ .+ initguess

                # Call function and store results
                rout = lmlovo(model, gmodel!, θ, data, n, k)

                nf += rout.nf
                nj += rout.nj
                ni += rout.iter

                (rout.status == 1) && (rout.f < wbest.f) && (wbest = rout)

                # This block is related to a strategy of smart
                # starting points for the multistart
                # process. Currently, it makes no sense to use it.

                # if rout.f < wbest.f

                #     # Send asynchronously the result to channel if success
                #     if rout.status == 1

                #         @async try

                #             put!(bqueue, rout.solution)

                #             @debug("Added new point to queue.", rout.solution, rout.f)

                #         catch e

                #             @warn(string("Problems when saving best point found in queue. ",
                #                          "Will skip this step"), e)

                #         end

                #     end

                # end

            end

            @debug("Finished. p = $(k) and f = $(wbest.f).")

            try

                wbest = RAFFOutput(wbest.status, wbest.solution, ni, wbest.p, wbest.f, nf, nj,
                                   wbest.outliers)

                put!(squeue, wbest)

            catch e

                if isa(e, InvalidStateException)

                    @warn("Solution queue prematurely closed. Unable to save solution for p=$(k).")

                    return

                end

                @warn("Something wrong when sending the solution to queue for p=$(k).", e)

            end

        end

    end

end

"""

    check_and_close(bqueue::RemoteChannel, tqueue::RemoteChannel,
                    squeue::RemoteChannel, futures::Vector{Future};
                    secs::Float64=0.1)

Check if there is at least one worker process in the vector of
`futures` that has not prematurely finished. If there is no alive
worker, close task, solution and best queues, `tqueue`, `squeue` and
`bqueue`, respectively.

"""
function check_and_close(bqueue::RemoteChannel, tqueue::RemoteChannel,
                         squeue::RemoteChannel, futures::Vector{Future};
                         secs::Float64=0.1)

    n_alive = length(futures)

    @debug("Checking worker status.")

    for (i, f) in enumerate(futures)

        if timedwait(()->isready(f), secs) == :ok

            @warn("Worker $(i) seems to have finished prematurely.",
                  fetch(f))

            n_alive -= 1

        end

    end

    @debug("Workers online: $(n_alive)")

    # Only closes all queues if there are tasks to be completed
    if n_alive == 0 && isopen(tqueue)

        @warn("No live worker found. Will close queues and finish.")

        close(bqueue)

        close(tqueue)

        close(squeue)

    end

end


text:
# mkt.jl
# ======
#
# An implementation of the McDonald Kreitman test.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

include("refs_and_graphs/codon_graphs.jl")
include("mst.jl")

function CodonSet{T}(sequences::Vector{Vector{Codon{T}}},
                     position::Integer) where T <: NucleicAcid
    s = CodonSet{T}()
    @inbounds for seq in sequences
        s |= seq[position]
    end
    return s
end

function simplest_mutation_path(cs::CodonSet{T},
                                cg::CodonGraph{Codon{T}},
                                msts::MSTState{Codon{T}}) where T <: NucleicAcid
    reset!(cg, cs)
    return mst(cg, msts)
end

function poly_count(cs::CodonSet{T},
                    cg::CodonGraph{Codon{T}},
                    msts::MSTState{Codon{T}}) where T <: NucleicAcid

    return simplest_mutation_path(cs, cg, msts)
end

function div_count(csa::CodonSet{T},
                   csb::CodonSet{T},
                   ref::CodonGraphReference{Codon{T}}) where T <: NucleicAcid
    DS = DN = 0
    R = 10
    if (csa & csb) == CodonSet{T}()
        for i in csa
            for j in csb
                edge = lookup_edge(ref, i, j)
                @inbounds DS_i = edge[3]
                @inbounds DN_i = edge[4]
                R_i = rankof(DS_i, DN_i)
                lowrank = R_i < R
                DS = ifelse(lowrank, DS_i, DS)
                DN = ifelse(lowrank, DN_i, DN)
                R = ifelse(lowrank, R_i, R)
            end
        end
    end
    return DS, DN
end

function mkt_PSPN(x::CodonSet{T},
                  y::CodonSet{T},
                  cg::CodonGraph{Codon{T}},
                  msts::MSTState{Codon{T}}) where T <: NucleicAcid

    a = poly_count(x, cg, msts)
    b = poly_count(y, cg, msts)
    return a[1] + b[1], a[2] + b[2]
end

mkt_α(PS, PN, DS, DN) = 1 - (DS * PN) / (DN * PS)

"""
    mkt(x, y, ref)

Compute McDonald Kreitman statistics for two sets of codons.
"""
function mkt(x::Vector{Vector{Codon{T}}},
             y::Vector{Vector{Codon{T}}},
             ref::CodonGraphReference{Codon{T}} = DEFAULT_CODON_GRAPH_REFERENCE) where T <: NucleicAcid

    n = min(minimum(length(xi) for xi in x), minimum(length(yi) for yi in y))
    graph = CodonGraph{Codon{T}}(ref)
    msts = MSTState{Codon{T}}()
    PS = PN = DS = DN = 0
    for i in 1:n
        xset = CodonSet{T}(x, i)
        yset = CodonSet{T}(y, i)

        PS_i, PN_i = mkt_PSPN(xset, yset, graph, msts)
        DS_i, DN_i = div_count(xset, yset, ref)

        PS += PS_i
        PN += PN_i
        DS +- DS_i
        DN += DN_i
    end
    α = mkt_α(PS, PN, DS, DN)
    return PS, PN, DS, DN, α
end