text:
"""
    loglikelihood(v::IHTVariable{T, M})

Calculates the loglikelihood of observing `y` given mean `μ = E(y) = g^{-1}(xβ)`
and some univariate distribution `d`.

Note that loglikelihood is the sum of the logpdfs for each observation.
"""
function loglikelihood(v::IHTVariable{T, M}) where {T <: Float, M}
    d = v.d
    y = v.y
    μ = v.μ
    wts = v.cv_wts
    logl = zero(T)
    ϕ = MendelIHT.deviance(v) / length(y) # variance in the case of normal
    @inbounds for i in eachindex(y)
        logl += loglik_obs(d, y[i], μ[i], wts[i], ϕ)
    end
    return logl
end

"""
This function is taken from GLM.jl from:
https://github.com/JuliaStats/GLM.jl/blob/956a64e7df79e80405867238781f24567bd40c78/src/glmtools.jl#L445

`wt`: in GLM.jl, this is prior frequency (a.k.a. case) weights for observations.
We use this for cross validation weighting: if sample `i` is being fitted,
wt[i] = 1. Otherwise wt[i] = 0.
"""
function loglik_obs end

loglik_obs(::Bernoulli, y, μ, wt, ϕ) = wt*logpdf(Bernoulli(μ), y)
loglik_obs(::Binomial, y, μ, wt, ϕ) = logpdf(Binomial(Int(wt), μ), Int(y*wt))
loglik_obs(::Gamma, y, μ, wt, ϕ) = wt*logpdf(Gamma(inv(ϕ), μ*ϕ), y)
loglik_obs(::InverseGaussian, y, μ, wt, ϕ) = wt*logpdf(InverseGaussian(μ, inv(ϕ)), y)
loglik_obs(::Normal, y, μ, wt, ϕ) = wt*logpdf(Normal(μ, sqrt(ϕ)), y)
loglik_obs(::Poisson, y, μ, wt, ϕ) = wt*logpdf(Poisson(μ), y)
# We use the following parameterization for the Negative Binomial distribution:
#    (Γ(r+y) / (Γ(r) * y!)) * μ^y * r^r / (μ+r)^{r+y}
# The parameterization of NegativeBinomial(r=r, p) in Distributions.jl is
#    Γ(r+y) / (y! * Γ(r)) * p^r(1-p)^y
# Hence, p = r/(μ+r)
loglik_obs(d::NegativeBinomial, y, μ, wt, ϕ) = wt*logpdf(NegativeBinomial(d.r, d.r/(μ+d.r)), y)

"""
    deviance(d, y, μ, wts)

Calculates the sum of the squared deviance residuals (e.g. (y - μ)^2 for Gaussian case)
Each individual sqared deviance residual is evaluated using `devresid`
which is implemented in GLM.jl
"""
function deviance(d::UnivariateDistribution, y::AbstractVector{T},
    μ::AbstractVector{T}, wts::AbstractVector{T}) where {T <: Float}
    dev = zero(T)
    @inbounds for i in eachindex(y)
        dev += wts[i] * devresid(d, y[i], μ[i])
    end
    return dev
end
deviance(v::IHTVariable{T, M}) where {T <: Float, M} =
    MendelIHT.deviance(v.d, v.y, v.μ, v.cv_wts)

"""
    update_μ!(μ, xb, l)

Update the mean (μ) using the linear predictor `xb` with link `l`.
"""
function update_μ!(μ::AbstractVecOrMat{T}, xb::AbstractVecOrMat{T}, l::Link) where {T <: Float}
    @inbounds for i in eachindex(μ)
        μ[i] = linkinv(l, xb[i])
    end
end

function update_μ!(v::IHTVariable{T, M}) where {T <: Float, M}
    μ = v.μ
    xb = v.xb
    zc = v.zc
    l = v.l
    @inbounds for i in eachindex(μ)
        μ[i] = linkinv(l, xb[i] + zc[i]) #genetic + nongenetic contributions
    end
end

"""
    update_xb!(v::IHTVariable{T, M})

Updates the linear predictors `xb` and `zc` with the new proposed `b` and `c`.
`b` is sparse but `c` (beta for non-genetic covariates) is dense.

We clamp the max value of each entry to (-20, 20) because certain distributions
(e.g. Poisson) have exponential link functions, which causes overflow.
"""
function update_xb!(v::IHTVariable{T, M}) where {T <: Float, M}
    # compute X*beta for genetic component
    if v.memory_efficient
        fill!(v.xb, 0)
        fill!(v.xk, 0)
        Threads.@threads for j in eachindex(v.idx)
            if v.idx[j]
                t = Threads.threadid()
                @inbounds @simd for i in 1:length(v.xb)
                    v.xk[i, t] += v.x[i, j] * v.b[j]
                end
            end
        end
        sum!(v.xb, v.xk) # v.xb .= sum(v.xk, dims=2)
    else
        copyto!(v.xk, @view(v.x[:, v.idx]))
        copyto!(v.gk, view(v.b, v.idx)) # use v.gk as storage
        mul!(v.xb, v.xk, v.gk)
    end
    # compute X*beta for nongenetic component
    mul!(v.zc, v.z, v.c)
    if !(typeof(v.d) <: Normal)
        clamp!(v.xb, -20, 20)
        clamp!(v.zc, -20, 20)
    end
end

"""
    score!(v::IHTVariable{T})

Calculates the score (gradient) `X^T * W * (y - g(x^T b))` for different GLMs.
W is a diagonal matrix where `w[i, i] = dμ/dη / var(μ)` (see documentation)
"""
function score!(v::IHTVariable{T, M}) where {T <: Float, M}
    d, l, x, z, y, cv_wts = v.d, v.l, v.x, v.z, v.y, v.cv_wts
    @inbounds for i in eachindex(y)
        η = v.xb[i] + v.zc[i]
        w = mueta(l, η) / glmvar(d, v.μ[i])
        v.r[i] = w * (y[i] - v.μ[i]) * cv_wts[i] # cv_wts handles sample masking for cross validation
    end
    mul!(v.df, Transpose(x), v.r)
    mul!(v.df2, Transpose(z), v.r)
end

"""
Wrapper function to decide whether to use Newton or MM algorithm for estimating
the nuisance paramter of negative binomial regression.
"""
function mle_for_r(v::IHTVariable{T, M}) where {T <: Float, M}
    method = v.est_r
    if method == :MM
        return update_r_MM(v)
    elseif method == :Newton
        return update_r_newton(v)
    else
        throw(ArgumentError("Only support method is Newton or MM, but got $method"))
    end

    return nothing
end

"""
Performs maximum loglikelihood estimation of the nuisance paramter for negative
binomial model using MM's algorithm.
"""
function update_r_MM(v::IHTVariable{T, M}) where {T <: Float, M}
    y = v.y
    μ = v.μ
    r = v.d.r # estimated r in previous iteration
    num = zero(T)
    den = zero(T)
    for i in eachindex(y)
        for j = 0:y[i] - 1
            num = num + (r /(r + j))  # numerator for r
        end
        p = r / (r + μ[i])
        den = den + log(p)  # denominator for r
    end

    return NegativeBinomial(-num / den, T(0.5))
end

"""
Performs maximum loglikelihood estimation of the nuisance paramter for negative
binomial model using Newton's algorithm. Will run a maximum of `maxIter` and
convergence is defaulted to `convTol`.
"""
function update_r_newton(v::IHTVariable{T, M};
    maxIter=100, convTol=T(1.e-6)) where {T <: Float, M}
    y = v.y
    μ = v.μ
    r = v.d.r # estimated r in previous iteration

    function first_derivative(r::T)
        tmp(yi, μi) = -(yi+r)/(μi+r) - log(μi+r) + 1 + log(r) + digamma(r+yi) - digamma(r)
        return sum(tmp(yi, μi) for (yi, μi) in zip(y, μ))
    end

    function second_derivative(r::T)
        tmp(yi, μi) = (yi+r)/(μi+r)^2 - 2/(μi+r) + 1/r + trigamma(r+yi) - trigamma(r)
        return sum(tmp(yi, μi) for (yi, μi) in zip(y, μ))
    end

    function negbin_loglikelihood(r::T)
        v.d = NegativeBinomial(r, T(0.5))
        return MendelIHT.loglikelihood(v)
    end

    function newton_increment(r::T)
        # use gradient descent if hessian not positive definite
        dx  = first_derivative(r)
        dx2 = second_derivative(r)
        if dx2 < 0
            increment = first_derivative(r) / second_derivative(r)
        else
            increment = first_derivative(r)
        end
        return increment
    end

    new_r    = one(T)
    stepsize = one(T)
    for i in 1:maxIter

        # run 1 iteration of Newton's algorithm
        increment = newton_increment(r)
        new_r = r - stepsize * increment

        # linesearch
        old_logl = negbin_loglikelihood(r)
        for j in 1:20
            if new_r <= 0
                stepsize = stepsize / 2
                new_r = r - stepsize * increment
            else
                new_logl = negbin_loglikelihood(new_r)
                if old_logl >= new_logl
                    stepsize = stepsize / 2
                    new_r = r - stepsize * increment
                else
                    break
                end
            end
        end

        #check convergence
        if abs(r - new_r) <= convTol
            return NegativeBinomial(new_r, T(0.5))
        else
            r = new_r
        end
    end

    return NegativeBinomial(r, T(0.5))
end

"""
This function computes the gradient step v.b = P_k(β + η∇f(β)) and updates idx and idc.
"""
function _iht_gradstep!(v::IHTVariable{T, M}, η::T) where {T <: Float, M}
    lg = length(v.group)
    J = v.J
    k = length(v.ks) > 0 ? v.ks : v.k

    # take gradient step: b = b + ηv, v = score
    BLAS.axpy!(η, v.df, v.b)
    BLAS.axpy!(η, v.df2, v.c)

    # project to sparsity, scaling model by weight vector, if supplied
    if lg == 0
        vectorize!(v.full_b, v.b, v.c, v.weight, v.zkeep)
        project_k!(v.full_b, k + v.zkeepn)
        unvectorize!(v.full_b, v.b, v.c, v.weight, v.zkeep)
    else
        # TODO: enable model selection for non-genetic covariates on group projection
        project_group_sparse!(v.b, v.group, J, k)
    end

    #recombute support
    v.idx .= v.b .!= 0
    v.idc .= v.c .!= 0

    # if more than J*k entries are selected, randomly choose J*k of them
    typeof(k) == Int && _choose!(v)

    # make necessary resizing since grad step might include/exclude non-genetic covariates
    check_covariate_supp!(v)
end

"""
    vectorize!(a::AbstractVector, b::AbstractVector, c::AbstractVector, ckeep::BitVector)

Without allocations, copies `b` into `a` and the copies `c` into
remaining parts of `a`.

`ckeep` tracks entries of `c` that will not be projected. That is, entries of `a`
corresponding to entries in `c[i]` will be filled with `Inf` if `ckeep[i] = true`.
"""
function vectorize!(a::AbstractVector, b::AbstractVector, c::AbstractVector,
    weight::AbstractVector, ckeep::BitVector)
    lb = length(b)
    lw = length(weight)
    la = length(a)

    # scale model by weight vector, if supplied
    if lw == 0
        copyto!(@view(a[1:lb]), b)
        copyto!(@view(a[lb+1:la]), c)
    else
        @inbounds for i in 1:lb
            a[i] = b[i] * weight[i]
        end
        ii = 1
        @inbounds for i in lb+1:la
            a[i] = c[ii] * weight[i]
            ii += 1
        end
    end

    # don't project certain non-genetic covariates
    a_view = @view(a[lb+1:la])
    @view(a_view[ckeep]) .= typemax(eltype(a))
end

"""
    unvectorize!(a::AbstractVector, B::AbstractMatrix, C::AbstractMatrix, Ckeep::BitVector)

Without allocations, copies the first `length(b)` part of `a` into `b`
and the remaining parts of `a` into `c`.

`ckeep` tracks entries of `c` that will not be projected. That is, entries in `c[i]`
will not be touched if `ckeep[i] = true`.
"""
function unvectorize!(a::AbstractVector, b::AbstractVector, c::AbstractVector,
    weight::AbstractVector, ckeep::BitVector)
    lb = length(b)
    lw = length(weight)
    la = length(a)

    # scale model by weight vector, if supplied
    if lw == 0
        copyto!(b, @view(a[1:lb]))
        ii = 1
        @inbounds for i in lb+1:la
            if !ckeep[ii]
                c[ii] = a[i]
            end
            ii += 1
        end
    else
        @inbounds for i in 1:lb
            b[i] = a[i] / weight[i]
        end
        ii = 1
        @inbounds for i in lb+1:la
            if !ckeep[ii]
                c[ii] = a[i] / weight[i]
            end
            ii += 1
        end
    end
end

"""
When initializing the IHT algorithm, take largest elements in magnitude of each
group of the score as nonzero components of b. This function set v.idx = 1 for
those indices. If `init_beta=true`, then beta values will be initialized to
their univariate values (see [`initialize_beta`](@ref)), in which case we will simply
choose top `k` entries

`J` is the maximum number of active groups, and `k` is the maximum number of
predictors per group.
"""
function init_iht_indices!(v::IHTVariable, init_beta::Bool, cv_idx::BitVector, verbose::Bool=false)
    fill!(v.b, 0)
    fill!(v.b0, 0)
    fill!(v.best_b, 0)
    fill!(v.xb, 0)
    fill!(v.xk, 0)
    fill!(v.gk, 0)
    fill!(v.xgk, 0)
    fill!(v.idx, false)
    fill!(v.idx0, false)
    copyto!(v.idc, v.zkeep)
    copyto!(v.idc0, v.zkeep)
    fill!(v.r, 0)
    fill!(v.df, 0)
    fill!(v.df2, 0)
    fill!(v.c, 0)
    fill!(v.best_c, 0)
    fill!(v.c0, 0)
    fill!(v.zc, 0)
    fill!(v.zdf2, 0)
    fill!(v.μ, 0)
    fill!(v.cv_wts, 0)
    fill!(v.full_b, 0)
    v.cv_wts[cv_idx] .= 1

    init_beta && !(typeof(v.d) <: Normal) &&
        throw(ArgumentError("Intializing beta values only work for Gaussian phenotypes! Sorry!"))

    # find the intercept by Newton's method
    ybar = zero(eltype(v.y))
    @inbounds @simd for i in eachindex(v.y)
        ybar += v.y[i] * v.cv_wts[i]
    end
    ybar /= count(!iszero, v.cv_wts)
    for iteration = 1:20
        g1 = linkinv(v.l, v.c[1])
        g2 = mueta(v.l, v.c[1])
        v.c[1] = v.c[1] - clamp((g1 - ybar) / g2, -1.0, 1.0)
        abs(g1 - ybar) < 1e-10 && break
    end
    mul!(v.zc, v.z, v.c)

    # update mean vector and use them to compute score (gradient)
    update_μ!(v)
    score!(v)

    if init_beta
        initialize_beta!(v, cv_idx, verbose)
        project_k!(v)
    else
        # first `k` non-zero entries are chosen based on largest gradient
        vectorize!(v.full_b, v.df, v.df2, v.weight, v.zkeep)
        if length(v.ks) == 0 # no group projection
            project_k!(v.full_b, v.k + v.zkeepn) # project k + number of nongentic covariates to keep
            unvectorize!(v.full_b, v.df, v.df2, v.weight, v.zkeep)
            v.idx .= v.df .!= 0
            v.idc .= v.zkeep

            # Choose randomly if more are selected
            _choose!(v)
        else
            project_group_sparse!(v.df, v.group, v.J, v.ks)
            v.idx .= v.b .!= 0
            fill!(v.idc, true)
        end
    end

    # make necessary resizing when necessary
    check_covariate_supp!(v)

    # store relevant components of x for first iteration
    v.memory_efficient || copyto!(v.xk, @view(v.x[:, v.idx]))
end

"""
if more than J*k entries are selected after projection, randomly select top J*k entries.
This can happen if entries of b are equal to each other.
"""
function _choose!(v::IHTVariable{T}) where {T <: Float}
    sparsity = v.k + v.zkeepn
    groups = (v.J == 0 ? 1 : v.J)

    nonzero = sum(v.idx) + sum(v.idc) - v.zkeepn
    if nonzero > groups * sparsity
        z = zero(eltype(v.b))
        non_zero_idx = findall(!iszero, v.idx)
        excess = nonzero - groups * sparsity
        for pos in sample(non_zero_idx, excess, replace=false)
            v.b[pos]   = z
            v.idx[pos] = false
        end
    end
end

"""
In `_init_iht_indices` and `_iht_gradstep!`, if non-genetic cov got
included/excluded, we must resize `xk` and `gk`.

TODO: Use ElasticArrays.jl
"""
function check_covariate_supp!(v::IHTVariable{T}) where {T <: Float}
    nzidx = sum(v.idx)
    if nzidx != size(v.xk, 2) && !v.memory_efficient
        v.xk = zeros(T, size(v.xk, 1), nzidx)
        v.gk = zeros(T, nzidx)
    end
    # @inbounds for i in eachindex(v.zkeep)
    #     v.zkeep[i] && !v.idc[i] && error("A non-genetic covariate was accidentally set to 0! Shouldn't happen!")
    # end
end

"""
    _iht_backtrack_(logl::T, prev_logl::T, η_step::Int64, nstep::Int64)

Returns true if one of the following conditions is met:
1. New loglikelihood is smaller than the old one
2. Current backtrack (`η_step`) exceeds maximum allowed backtracking (`nstep`, default = 3)
"""
function _iht_backtrack_(logl::T, prev_logl::T, η_step::Int64, nstep::Int64) where {T <: Float}
    (prev_logl > logl) && (η_step < nstep)
end

"""
    standardize!(z::AbstractVecOrMat)

Standardizes each column of `z` to mean 0 and variance 1. Make sure you
do not standardize the intercept.
"""
@inline function standardize!(z::AbstractVecOrMat)
    n, q = size(z)
    μ = _mean(z)
    σ = _std(z, μ)

    @inbounds for j in 1:q
        @simd for i in 1:n
            z[i, j] = (z[i, j] - μ[j]) * σ[j]
        end
    end
end

@inline function _mean(z)
    n, q = size(z)
    μ = zeros(q)
    @inbounds for j in 1:q
        tmp = 0.0
        @simd for i in 1:n
            tmp += z[i, j]
        end
        μ[j] = tmp / n
    end
    return μ
end

function _std(z, μ)
    n, q = size(z)
    σ = zeros(q)

    @inbounds for j in 1:q
        @simd for i in 1:n
            σ[j] += (z[i, j] - μ[j])^2
        end
        σ[j] = 1.0 / sqrt(σ[j] / (n - 1))
    end
    return σ
end

"""
    project_k!(x::AbstractVector, k::Integer)

Sets all but the largest `k` entries of `x` to 0.

# Examples:
```julia-repl
using MendelIHT
x = [1.0; 2.0; 3.0]
project_k!(x, 2) # keep 2 largest entry
julia> x
3-element Array{Float64,1}:
 0.0
 2.0
 3.0
```

# Arguments:
- `x`: the vector to project.
- `k`: the number of components of `x` to preserve.
"""
function project_k!(x::AbstractVector{T}, k::Int64) where {T <: Float}
    k < 0 && throw(DomainError("Attempted to project to sparsity level $k"))
    a = abs(partialsort(x, k, by=abs, rev=true))
    @inbounds for i in eachindex(x)
        abs(x[i]) < a && (x[i] = zero(T))
    end
end

function project_k!(v::IHTVariable)
    v.k < 0 && throw(DomainError("Attempted to project to sparsity level $(v.k)"))

    # copy genetic and non-genetic effects to full_grad, project, and copy back
    vectorize!(v.full_b, v.b, v.c, v.weight, v.zkeep)
    project_k!(v.full_b, v.k + v.zkeepn)
    unvectorize!(v.full_b, v.b, v.c, v.weight, v.zkeep)

    # update support
    v.idx .= v.b .!= 0
    v.idc .= v.c .!= 0
    check_covariate_supp!(v)
end

"""
    project_group_sparse!(y::AbstractVector, group::AbstractVector, J::Integer, k<:Real)

When `k` is an integer, projects the vector `y` onto the set with at most `J` active groups
and at most `k` active predictors per group. To have variable group sparsity level, input `k`
as a vector of integers. We will preserve `k[1]` elements for group 1, `k[2]` predictors for
group 2...etc. This function assumes there are no unknown or overlaping group membership.

Note: In the `group` vector, the first group must be 1, and the second group must be 2...etc.

# Examples
```julia-repl
using MendelIHT
J, k, n = 2, 3, 20
y = collect(1.0:20.0)
y_copy = copy(y)
group = rand(1:5, n)
project_group_sparse!(y, group, J, k)
for i = 1:length(y)
    println(i,"  ",group[i],"  ",y[i],"  ",y_copy[i])
end

J, k, n = 2, 0.9, 20
y = collect(1.0:20.0)
y_copy = copy(y)
group = rand(1:5, n)
project_group_sparse!(y, group, J, k)
for i = 1:length(y)
    println(i,"  ",group[i],"  ",y[i],"  ",y_copy[i])
end
```

# Arguments
- `y`: The vector to project
- `group`: Vector encoding group membership
- `J`: Max number of non-zero group
- `k`: Maximum predictors per group. Can be a positive integer or a vector of integers.
"""
function project_group_sparse!(y::AbstractVector{T}, group::AbstractVector{Int64},
    J::Int64, k::Int64) where {T <: Float}
    groups = maximum(group)          # number of groups
    group_count = zeros(Int, groups) # counts number of predictors in each group
    group_norm = zeros(groups)       # l2 norm of each group
    perm = zeros(Int64, length(y))   # vector holding the permuation vector after sorting
    sortperm!(perm, y, by = abs, rev = true)

    #calculate the magnitude of each group, where only top predictors contribute
    for i in eachindex(y)
        j = perm[i]
        n = group[j]
        if group_count[n] < k
            group_norm[n] = group_norm[n] + y[j]^2
            group_count[n] = group_count[n] + 1
        end
    end

    #go through the top predictors in order. Set predictor to 0 if criteria not met
    group_rank = zeros(Int64, length(group_norm))
    sortperm!(group_rank, group_norm, rev = true)
    group_rank = invperm(group_rank)
    fill!(group_count, 1)
    for i in eachindex(y)
        j = perm[i]
        n = group[j]
        if (group_rank[n] > J) || (group_count[n] > k)
            y[j] = 0.0
        else
            group_count[n] = group_count[n] + 1
        end
    end
end

function project_group_sparse!(y::AbstractVector{T}, group::AbstractVector{Int64},
    J::Int64, k::Vector{Int}) where {T <: Float}
    groups = maximum(group)          # number of groups
    group_count = zeros(Int, groups) # counts number of predictors in each group
    group_norm = zeros(groups)       # l2 norm of each group
    perm = zeros(Int64, length(y))   # vector holding the permuation vector after sorting
    sortperm!(perm, y, by = abs, rev = true)

    #calculate the magnitude of each group, where only top predictors contribute
    for i in eachindex(y)
        j = perm[i]
        n = group[j]
        if group_count[n] < k[n]
            group_norm[n] = group_norm[n] + y[j]^2
            group_count[n] = group_count[n] + 1
        end
    end

    #go through the top predictors in order. Set predictor to 0 if criteria not met
    group_rank = zeros(Int64, length(group_norm))
    sortperm!(group_rank, group_norm, rev = true)
    group_rank = invperm(group_rank)
    fill!(group_count, 1)
    for i in eachindex(y)
        j = perm[i]
        n = group[j]
        if (group_rank[n] > J) || (group_count[n] > k[n])
            y[j] = 0.0
        else
            group_count[n] = group_count[n] + 1
        end
    end
end

"""
    maf_weights(x::SnpArray; max_weight::T = Inf)

Calculates the prior weight based on minor allele frequencies.

Returns an array of weights where `w[i] = 1 / (2 * sqrt(p[i] (1 - p[i]))) ∈ (1, ∞).`
Here `p` is the minor allele frequency computed by `maf()` in SnpArrays.

- `x`: A SnpArray
- `max_weight`: Maximum weight for any predictor. Defaults to `Inf`.
"""
function maf_weights(x::SnpArray; max_weight::T = Inf) where {T <: Float}
    p = SnpArrays.maf(x)
    p .= 1 ./ (2 .* sqrt.(p .* (1 .- p)))
    clamp!(p, 1.0, max_weight)
    return p
end

"""
Function that saves `b`, `xb`, `idx`, `idc`, `c`, and `zc` after each iteration.
"""
function save_prev!(v::IHTVariable{T}, cur_logl::T, best_logl::T) where {T <: Float}
    copyto!(v.b0, v.b)     # b0 = b
    copyto!(v.idx0, v.idx) # idx0 = idx
    copyto!(v.idc0, v.idc) # idc0 = idc
    copyto!(v.c0, v.c)     # c0 = c
    if cur_logl > best_logl
        copyto!(v.best_b, v.b) # best_b = b
        copyto!(v.best_c, v.c) # best_c = c
    end
    return max(cur_logl, best_logl)
end

"""
Computes the best step size η = v'v / v'Jv = v'v / v'X'WXv

Here v is the score and J is the expected information matrix, which is
computed by J = g'(xb) / var(μ), assuming dispersion is 1. Note

Note: cross validation weights are needed in W.
"""
function iht_stepsize!(v::IHTVariable{T, M}) where {T <: Float, M}
    z = v.z # non genetic covariates
    d = v.d # distribution
    l = v.l # link function

    # first compute xgk = X*v using relevant components of gradient
    if v.memory_efficient
        fill!(v.xgk, 0)
        fill!(v.xk, 0)
        Threads.@threads for j in eachindex(v.idx)
            if v.idx[j]
                t = Threads.threadid()
                @inbounds @simd for i in 1:length(v.xgk)
                    v.xk[i, t] += v.x[i, j] * v.df[j]
                end
            end
        end
        sum!(v.xgk, v.xk) # v.xgk .= sum(v.xk, dims=2)
    else
        copyto!(v.gk, view(v.df, v.idx))
        mul!(v.xgk, v.xk, v.gk)
    end
    mul!(v.zdf2, view(z, :, v.idc), view(v.df2, v.idc))
    v.xgk .+= v.zdf2 # xgk = Xv

    # Compute sqrt(W); use zdf2 as storage
    @inbounds @simd for i in eachindex(v.zdf2)
        v.zdf2[i] = sqrt(abs2(mueta(l, v.xb[i] + v.zc[i])) / glmvar(d, v.μ[i])) * v.cv_wts[i]
    end
    v.xgk .*= v.zdf2 # xgk = sqrt(W)Xv

    # now compute and return step size. Note non-genetic covariates are separated from x
    numer = v.memory_efficient ? sum(abs2, @view(v.df[v.idx])) : sum(abs2, v.gk)
    numer += sum(abs2, @view(v.df2[v.idc]))
    denom = dot(v.xgk, v.xgk)
    η = numer / denom

    # for bad boundary cases (sometimes, k = 1 in cross validation generates weird η)
    isinf(η) && (η = T(1e-8))
    isnan(η) && (η = T(1e-8))

    return η :: T
end

"""
    initialize_beta!(v, cv_wts)

Initialze beta to univariate regression values. That is, `β[i]` is set to the estimated
beta with `y` as response, and `x[:, i]` with an intercept term as covariate.

Progress will be displayed if `verbose=true`

TODO: this function assumes quantitative (Gaussian) phenotypes. Make it work for other distributions
"""
function initialize_beta!(
    v::IHTVariable,
    cv_wts::BitVector, # cross validation weights; 1 = sample is present, 0 = not present
    verbose::Bool=true
    )
    y, x, z, β, c, T = v.y, v.x, v.z, v.b, v.c, eltype(v.b)
    xtx_store = [zeros(T, 2, 2) for _ in 1:Threads.nthreads()]
    xty_store = [zeros(T, 2) for _ in 1:Threads.nthreads()]
    xstore = [zeros(T, sum(cv_wts)) for _ in 1:Threads.nthreads()]
    ystore = y[cv_wts]
    pmeter = verbose ? Progress(nsnps(v) + ncovariates(v), 5,
        "Initializing β to univariate regression values...") : nothing
    # genetic covariates
    c0 = zero(eltype(z))
    Threads.@threads for i in 1:nsnps(v)
        id = Threads.threadid()
        copyto!(xstore[id], @view(x[cv_wts, i]))
        linreg!(xstore[id], ystore, xtx_store[id], xty_store[id])
        c0 += xty_store[id][1]
        β[i] = xty_store[id][2]
        verbose && next!(pmeter) # update progress
    end
    # non-genetic covariates
    Threads.@threads for i in 2:ncovariates(v)
        id = Threads.threadid()
        copyto!(xstore[id], @view(z[cv_wts, i]))
        linreg!(xstore[id], ystore, xtx_store[id], xty_store[id])
        c0 += xty_store[id][1]
        c[i] = xty_store[id][2]
        verbose && next!(pmeter) # update progress
    end
    c[1] = c0 / (nsnps(v) + ncovariates(v) - 1)
    clamp!(v.b, -2, 2)
    clamp!(v.c, -2, 2)
    copyto!(v.b0, β)
    copyto!(v.c0, c)
end

"""
    linreg!(x::Vector, y::Vector)

Performs linear regression with `y` as response, `x` and a vector of 1 as
covariate. `β̂` will be stored in `xty_store`.

Code inspired from Doug Bates on Discourse:
https://discourse.julialang.org/t/efficient-way-of-doing-linear-regression/31232/28
"""
function linreg!(
    x::AbstractVector{T},
    y::AbstractVector{T},
    xtx_store::AbstractMatrix{T} = zeros(T, 2, 2),
    xty_store::AbstractVector{T} = zeros(T, 2)
    ) where {T<:AbstractFloat}
    N = length(x)
    N == length(y) || throw(DimensionMismatch())
    xtx_store[1, 1] = N
    xtx_store[1, 2] = sum(x)
    xtx_store[2, 2] = sum(abs2, x)
    xty_store[1] = sum(y)
    xty_store[2] = dot(x, y)
    try # rare SNPs may have 0s everywhere, causing cholesky to fail
        ldiv!(cholesky!(Symmetric(xtx_store, :U)), xty_store)
    catch
        return xty_store
    end
    return xty_store
end

"""
This function initializes 1 instance of a GeneralizedLinearModel(G<:GlmResp, L<:LinPred, Bool).
"""
function initialize_glm_object()
    d = Bernoulli
    l = canonicallink(d())
    x = rand(100, 2)
    y = rand(0:1, 100)
    return fit(GeneralizedLinearModel, x, y, d(), l)
end

"""
    naive_impute(x, destination)

Imputes missing entries of a SnpArray using the mode of each SNP, and
saves the result in a new file called destination in current directory.
Non-missing entries are the same.
"""
function naive_impute(x::SnpArray, destination::String)
    n, p = size(x)
    y = SnpArray(destination, n, p)

    @inbounds for j in 1:p

        #identify mode
        entry0, entry1, entry2 = 0, 0, 0
        for i in 1:n
            if x[i, j] == 0x00
                y[i, j] = 0x00
                entry0 += 1
            elseif x[i, j] == 0x02
                y[i, j] = 0x02
                entry1 += 1
            elseif x[i, j] == 0x03
                y[i, j] = 0x03
                entry2 += 1
            end
        end
        most_often = max(entry0, entry1, entry2)
        missing_entry = 0x00
        if most_often == entry1
            missing_entry = 0x02
        elseif most_often == entry2
            missing_entry = 0x03
        end

        # impute
        for i in 1:n
            if x[i, j] == 0x01
                y[i, j] = missing_entry
            end
        end
    end

    return nothing
end

# small function to check sparsity parameter `k` is reasonable.
function check_group(k, group)
    if typeof(k) <: Vector
        @assert length(group) > 1 "Doubly sparse projection specified (since k" *
            " is a vector) but there are no group information."
        for i in 1:length(k)
            group_member = count(x -> x == i, group)
            group_member > k[i] || throw(DomainError("Maximum predictors for group " *
                "$i was $(k[i]) but there are only $group_member predictors is this " *
                "group. Please choose a smaller number."))
        end
    else
        @assert k >= 0 "Value of k (max predictors per group) must be nonnegative!\n"
    end
end

# helper function from https://discourse.julialang.org/t/how-to-find-out-the-version-of-a-package-from-its-module/37755
pkgversion(m::Module) = Pkg.TOML.parsefile(joinpath(dirname(string(first(methods(m.eval)).file)), "..", "Project.toml"))["version"]

function print_iht_signature(io::IO)
    v = pkgversion(MendelIHT)
    println(io, "****                   MendelIHT Version $v                  ****")
    println(io, "****     Benjamin Chu, Kevin Keys, Chris German, Hua Zhou       ****")
    println(io, "****   Jin Zhou, Eric Sobel, Janet Sinsheimer, Kenneth Lange    ****")
    println(io, "****                                                            ****")
    println(io, "****                 Please cite our paper!                     ****")
    println(io, "****         https://doi.org/10.1093/gigascience/giaa044        ****")
    println(io, "")
    io != stdout && print_iht_signature(stdout)
end
print_iht_signature() = print_iht_signature(stdout)

function print_parameters(io::IO, k, d, l, use_maf, group, debias, tol, max_iter, min_iter)
    regression = typeof(d) <: Normal ? "linear" : typeof(d) <: Bernoulli ?
        "logistic" : typeof(d) <: Poisson ? "Poisson" :
        typeof(d) <: NegativeBinomial ? "NegativeBinomial" :
        typeof(d) <: MvNormal ? "Multivariate Gaussian" : "unknown"
    println(io, "Running sparse $regression regression")
    println(io, "Number of threads = ", Threads.nthreads())
    println(io, "Link functin = $l")
    typeof(k) <: Int && println(io, "Sparsity parameter (k) = $k")
    typeof(k) <: Vector{Int} && println(io, "Sparsity parameter (k) = using group membership specified in k")
    println(io, "Prior weight scaling = ", use_maf ? "on" : "off")
    println(io, "Doubly sparse projection = ", length(group) > 0 ? "on" : "off")
    println(io, "Debias = ", debias ? "on" : "off")
    println(io, "Max IHT iterations = $max_iter")
    println(io, "Converging when tol < $tol and iteration ≥ $min_iter:\n")
    io != stdout && print_parameters(stdout, k, d, l, use_maf, group, debias, tol, max_iter, min_iter)
end
print_parameters(k, d, l, use_maf, group, debias, tol, max_iter, min_iter) =
    print_parameters(stdout, k, d, l, use_maf, group, debias, tol, max_iter, min_iter)

function check_convergence(v::IHTVariable)
    the_norm = max(chebyshev(v.b, v.b0), chebyshev(v.c, v.c0)) #max(abs(x - y))
    scaled_norm = the_norm / (max(norm(v.b0, Inf), norm(v.c0, Inf)) + 1.0)
    return scaled_norm
end

function backtrack!(v::IHTVariable, η::Float)
    # recompute gradient step
    copyto!(v.b, v.b0)
    copyto!(v.c, v.c0)
    _iht_gradstep!(v, η)

    # recompute η = xb, μ = g(η), and loglikelihood to see if we're now increasing
    update_xb!(v)
    update_μ!(v)
    if v.est_r != :None
        v.d = mle_for_r(v)
    end

    return loglikelihood(v)
end

function check_data_dim(y::AbstractVecOrMat, x::AbstractMatrix, z::AbstractVecOrMat)
    if is_multivariate(y)
        r, n1 = size(y)
        p, n2 = size(x)
        q, n3 = size(z)
        n1 == n2 == n3 || throw(DimensionMismatch("Detected multivariate analysis" *
            " but size(y, 2) = $n1, size(x, 2) = $n2, size(z, 2) = $n3 which don't " *
            "match. Recall each column of `y`, `x`, `z` should be sample " *
            "phenotypes/genotypes/covariates."))
    else
        n1 = length(y)
        n2, p = size(x)
        n3, = size(z)
        n1 == n2 == n3 || throw(DimensionMismatch("Detected univariate analysis " *
        "but length(y) = $n1, size(x, 1) = $n2, size(z, 1) = $n3 which don't match. " *
        "Recall each `y` should be a vector of phenotypes, and each row of `x`" *
        " and `z` should be sample genotypes/covariates."))
    end
end

function save_best_model!(v::IHTVariable)
    # compute η = xb with the best estimated model
    copyto!(v.b, v.best_b)
    copyto!(v.c, v.best_c)
    v.idx .= v.b .!= 0
    v.idc .= v.c .!= 0
    check_covariate_supp!(v)
    update_xb!(v)

    # update estimated mean μ with genotype predictors
    update_μ!(v.μ, v.xb, v.l)
end

"""
    debias!(v::IHTVariable)

After each IHT iteration, `β` is sparse. This function solves for the exact
solution `β̂` on the non-zero indices of `β`, a process known as debiasing.
"""
function debias!(v::IHTVariable)
    v.memory_efficient && error("Currently debiasing only works with memory_efficient=false")
    if sum(v.idx) == size(v.xk, 2)
        temp_glm = fit(GeneralizedLinearModel, v.xk, v.y, v.d, v.l)
        view(v.b, v.idx) .= temp_glm.pp.beta0
    end
end


text:
using ClimateMachine.Mesh.Elements
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Metrics
using LinearAlgebra: I
using Test
using Random: MersenneTwister

const _ξ1x1, _ξ2x1, _ξ3x1 = Grids._ξ1x1, Grids._ξ2x1, Grids._ξ3x1
const _ξ1x2, _ξ2x2, _ξ3x2 = Grids._ξ1x2, Grids._ξ2x2, Grids._ξ3x2
const _ξ1x3, _ξ2x3, _ξ3x3 = Grids._ξ1x3, Grids._ξ2x3, Grids._ξ3x3
const _M, _MI = Grids._M, Grids._MI
const _x1, _x2, _x3 = Grids._x1, Grids._x2, Grids._x3
const _JcV = Grids._JcV
const _nvgeo = Grids._nvgeo

const _n1, _n2, _n3 = Grids._n1, Grids._n2, Grids._n3
const _sM, _vMI = Grids._sM, Grids._vMI
const _nsgeo = Grids._nsgeo

@testset "1-D Metric terms" begin
    for FT in (Float32, Float64)
        #{{{
        let
            N = (4,)
            Nq = N .+ 1
            Np = prod(Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(j -> Elements.lglpoints(FT, N[j]), dim)
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)

            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            dim = 1
            e2c = Array{FT, 3}(undef, 1, 2, 2)
            e2c[:, :, 1] = [-1 0]
            e2c[:, :, 2] = [0 10]
            nelem = size(e2c, 3)

            (vgeo, sgeo, _) =
                Grids.computegeometry(e2c, D, ξ, ω, (x...) -> identity(x))
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            @test vgeo[:, _x1, 1] ≈ (ξ[1] .- 1) / 2
            @test vgeo[:, _x1, 2] ≈ 5 * (ξ[1] .+ 1)

            @test vgeo[:, _M, 1] ≈ ω[1] .* ones(FT, Nq) / 2
            @test vgeo[:, _M, 2] ≈ 5 * ω[1] .* ones(FT, Nq)
            @test vgeo[:, _ξ1x1, 1] ≈ 2 * ones(FT, Nq)
            @test vgeo[:, _ξ1x1, 2] ≈ ones(FT, Nq) / 5
            @test sgeo[_n1, 1, 1, :] ≈ -ones(FT, nelem)
            @test sgeo[_n1, 1, 2, :] ≈ ones(FT, nelem)
            @test sgeo[_sM, 1, 1, :] ≈ ones(FT, nelem)
            @test sgeo[_sM, 1, 2, :] ≈ ones(FT, nelem)
        end
        #}}}
    end

    # N = 0 test
    for FT in (Float32, Float64)
        #{{{
        let
            N = (0,)
            Nq = N .+ 1
            Np = prod(Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(j -> Elements.glpoints(FT, N[j]), dim)
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)

            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            dim = 1
            e2c = Array{FT, 3}(undef, 1, 2, 2)
            e2c[:, :, 1] = [-1 0]
            e2c[:, :, 2] = [0 10]
            nelem = size(e2c, 3)

            (vgeo, sgeo, _) =
                Grids.computegeometry(e2c, D, ξ, ω, (x...) -> identity(x))
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            @test vgeo[1, _x1, 1] ≈ sum(e2c[:, :, 1]) / 2
            @test vgeo[1, _x1, 2] ≈ sum(e2c[:, :, 2]) / 2

            @test vgeo[:, _M, 1] ≈ ω[1] .* ones(FT, Nq) / 2
            @test vgeo[:, _M, 2] ≈ 5 * ω[1] .* ones(FT, Nq)
            @test vgeo[:, _ξ1x1, 1] ≈ 2 * ones(FT, Nq)
            @test vgeo[:, _ξ1x1, 2] ≈ ones(FT, Nq) / 5

            @test sgeo[_n1, 1, 1, :] ≈ -ones(FT, nelem)
            @test sgeo[_n1, 1, 2, :] ≈ ones(FT, nelem)
            @test sgeo[_sM, 1, 1, :] ≈ ones(FT, nelem)
            @test sgeo[_sM, 1, 2, :] ≈ ones(FT, nelem)
        end
        #}}}
    end
end

@testset "2-D Metric terms" begin
    # for FT in (Float32, Float64)
    for FT in (Float32, Float64), N in ((4, 4), (4, 6), (6, 4))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        dim = length(N)
        nface = 2dim

        # Create element operators for each polynomial order
        ξω = ntuple(j -> Elements.lglpoints(FT, N[j]), dim)
        ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
        D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

        # linear and rotation test
        #{{{
        let
            e2c = Array{FT, 3}(undef, 2, 4, 4)
            e2c[:, :, 1] = [
                0 2 0 2
                0 0 2 2
            ]
            e2c[:, :, 2] = [
                2 2 0 0
                0 2 0 2
            ]
            e2c[:, :, 3] = [
                2 0 2 0
                2 2 0 0
            ]
            e2c[:, :, 4] = [
                0 0 2 2
                2 0 2 0
            ]
            nelem = size(e2c, 3)

            x_exact = Array{FT, 3}(undef, Nq..., 4)
            x_exact[:, :, 1] .= 1 .+ ξ[1]
            x_exact[:, :, 2] .= 1 .- ξ[2]'
            x_exact[:, :, 3] .= 1 .- ξ[1]
            x_exact[:, :, 4] .= 1 .+ ξ[2]'

            y_exact = Array{FT, 3}(undef, Nq..., 4)
            y_exact[:, :, 1] .= 1 .+ ξ[2]'
            y_exact[:, :, 2] .= 1 .+ ξ[1]
            y_exact[:, :, 3] .= 1 .- ξ[2]'
            y_exact[:, :, 4] .= 1 .- ξ[1]

            M_exact =
                ones(FT, Nq..., 4) .* reshape(kron(reverse(ω)...), Nq..., 1)

            ξ1x1_exact = zeros(FT, Nq..., 4)
            ξ1x1_exact[:, :, 1] .= 1
            ξ1x1_exact[:, :, 3] .= -1

            ξ1x2_exact = zeros(FT, Nq..., 4)
            ξ1x2_exact[:, :, 2] .= 1
            ξ1x2_exact[:, :, 4] .= -1

            ξ2x1_exact = zeros(FT, Nq..., 4)
            ξ2x1_exact[:, :, 2] .= -1
            ξ2x1_exact[:, :, 4] .= 1

            ξ2x2_exact = zeros(FT, Nq..., 4)
            ξ2x2_exact[:, :, 1] .= 1
            ξ2x2_exact[:, :, 3] .= -1

            sM_exact = fill(FT(NaN), maximum(Nfp), nface, nelem)
            sM_exact[1:Nfp[1], 1, :] .= 1 .* ω[2]
            sM_exact[1:Nfp[1], 2, :] .= 1 .* ω[2]
            sM_exact[1:Nfp[2], 3, :] .= 1 .* ω[1]
            sM_exact[1:Nfp[2], 4, :] .= 1 .* ω[1]

            nx_exact = fill(FT(NaN), maximum(Nfp), nface, nelem)
            nx_exact[1:Nfp[1], 1:2, :] .= 0
            nx_exact[1:Nfp[2], 3:4, :] .= 0

            nx_exact[1:Nfp[1], 1, 1] .= -1
            nx_exact[1:Nfp[1], 2, 1] .= 1
            nx_exact[1:Nfp[2], 3, 2] .= 1
            nx_exact[1:Nfp[2], 4, 2] .= -1
            nx_exact[1:Nfp[1], 1, 3] .= 1
            nx_exact[1:Nfp[1], 2, 3] .= -1
            nx_exact[1:Nfp[2], 3, 4] .= -1
            nx_exact[1:Nfp[2], 4, 4] .= 1

            ny_exact = fill(FT(NaN), maximum(Nfp), nface, nelem)
            ny_exact[1:Nfp[1], 1:2, :] .= 0
            ny_exact[1:Nfp[2], 3:4, :] .= 0

            ny_exact[1:Nfp[2], 3, 1] .= -1
            ny_exact[1:Nfp[2], 4, 1] .= 1
            ny_exact[1:Nfp[1], 1, 2] .= -1
            ny_exact[1:Nfp[1], 2, 2] .= 1
            ny_exact[1:Nfp[2], 3, 3] .= 1
            ny_exact[1:Nfp[2], 4, 3] .= -1
            ny_exact[1:Nfp[1], 1, 4] .= 1
            ny_exact[1:Nfp[1], 2, 4] .= -1

            (vgeo, sgeo, _) =
                Grids.computegeometry(e2c, D, ξ, ω, (x...) -> identity(x))
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)

            @test (@view vgeo[:, :, _x1, :]) ≈ x_exact
            @test (@view vgeo[:, :, _x2, :]) ≈ y_exact
            @test (@view vgeo[:, :, _M, :]) ≈ M_exact
            @test (@view vgeo[:, :, _ξ1x1, :]) ≈ ξ1x1_exact
            @test (@view vgeo[:, :, _ξ1x2, :]) ≈ ξ1x2_exact
            @test (@view vgeo[:, :, _ξ2x1, :]) ≈ ξ2x1_exact
            @test (@view vgeo[:, :, _ξ2x2, :]) ≈ ξ2x2_exact
            msk = isfinite.(sM_exact)
            @test sgeo[_sM, :, :, :][msk] ≈ sM_exact[msk]
            @test sgeo[_n1, :, :, :][msk] ≈ nx_exact[msk]
            @test sgeo[_n2, :, :, :][msk] ≈ ny_exact[msk]

            nothing
        end
        #}}}

        # Polynomial 2-D test
        #{{{
        let
            f(ξ1, ξ2) = (
                9 .* ξ1 - (1 .+ ξ1) .* ξ2 .^ 2 +
                (ξ1 .- 1) .^ 2 .* (1 .- ξ2 .^ 2 .+ ξ2 .^ 3),
                10 .* ξ2 .+ ξ1 .^ 4 .* (1 .- ξ2) .+ ξ1 .^ 2 .* ξ2 .* (1 .+ ξ2),
            )
            fx1ξ1(ξ1, ξ2) =
                7 .+ ξ2 .^ 2 .- 2 .* ξ2 .^ 3 .+
                2 .* ξ1 .* (1 .- ξ2 .^ 2 .+ ξ2 .^ 3)
            fx1ξ2(ξ1, ξ2) =
                -2 .* (1 .+ ξ1) .* ξ2 .+
                (-1 .+ ξ1) .^ 2 .* ξ2 .* (-2 .+ 3 .* ξ2)
            fx2ξ1(ξ1, ξ2) =
                -4 .* ξ1 .^ 3 .* (-1 .+ ξ2) .+ 2 .* ξ1 .* ξ2 .* (1 .+ ξ2)
            fx2ξ2(ξ1, ξ2) = 10 .- ξ1 .^ 4 .+ ξ1 .^ 2 .* (1 .+ 2 .* ξ2)

            e2c = Array{FT, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            # Create the metrics
            (x1ξ1, x1ξ2, x2ξ1, x2ξ2) = let
                (vgeo, _) = Grids.computegeometry(
                    e2c,
                    D,
                    ξ,
                    ω,
                    (x...) -> identity(x),
                )
                vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
                ξ1, ξ2 = vgeo[:, :, _x1, :], vgeo[:, :, _x2, :]
                (fx1ξ1(ξ1, ξ2), fx1ξ2(ξ1, ξ2), fx2ξ1(ξ1, ξ2), fx2ξ2(ξ1, ξ2))
            end
            J = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1)
            M = J .* reshape(kron(reverse(ω)...), Nq..., 1)

            meshwarp(ξ1, ξ2, _) = (f(ξ1, ξ2)..., 0)
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            x1 = @view vgeo[:, :, _x1, :]
            x2 = @view vgeo[:, :, _x2, :]

            @test M ≈ (@view vgeo[:, :, _M, :])
            @test (@view vgeo[:, :, _ξ1x1, :]) ≈ x2ξ2 ./ J
            @test (@view vgeo[:, :, _ξ2x1, :]) ≈ -x2ξ1 ./ J
            @test (@view vgeo[:, :, _ξ1x2, :]) ≈ -x1ξ2 ./ J
            @test (@view vgeo[:, :, _ξ2x2, :]) ≈ x1ξ1 ./ J

            # check the normals?
            sM = @view sgeo[_sM, :, :, :]
            n1 = @view sgeo[_n1, :, :, :]
            n2 = @view sgeo[_n2, :, :, :]
            @test all(hypot.(n1[1:Nfp[1], 1:2, :], n2[1:Nfp[1], 1:2, :]) .≈ 1)
            @test all(hypot.(n1[1:Nfp[2], 3:4, :], n2[1:Nfp[2], 3:4, :]) .≈ 1)
            @test sM[1:Nfp[1], 1, :] .* n1[1:Nfp[1], 1, :] ≈
                  -x2ξ2[1, :, :] .* ω[2]
            @test sM[1:Nfp[1], 1, :] .* n2[1:Nfp[1], 1, :] ≈
                  x1ξ2[1, :, :] .* ω[2]
            @test sM[1:Nfp[1], 2, :] .* n1[1:Nfp[1], 2, :] ≈
                  x2ξ2[Nq[1], :, :] .* ω[2]
            @test sM[1:Nfp[1], 2, :] .* n2[1:Nfp[1], 2, :] ≈
                  -x1ξ2[Nq[1], :, :] .* ω[2]
            @test sM[1:Nfp[2], 3, :] .* n1[1:Nfp[2], 3, :] ≈
                  x2ξ1[:, 1, :] .* ω[1]
            @test sM[1:Nfp[2], 3, :] .* n2[1:Nfp[2], 3, :] ≈
                  -x1ξ1[:, 1, :] .* ω[1]
            @test sM[1:Nfp[2], 4, :] .* n1[1:Nfp[2], 4, :] ≈
                  -x2ξ1[:, Nq[2], :] .* ω[1]
            @test sM[1:Nfp[2], 4, :] .* n2[1:Nfp[2], 4, :] ≈
                  x1ξ1[:, Nq[2], :] .* ω[1]
        end
        #}}}

        # Constant preserving test
        #{{{
        let
            rng = MersenneTwister(777)
            f(ξ1, ξ2) = (
                ξ1 + (ξ1 * ξ2 * rand(rng) + rand(rng)) / 10,
                ξ2 + (ξ1 * ξ2 * rand(rng) + rand(rng)) / 10,
            )

            e2c = Array{FT, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            meshwarp(ξ1, ξ2, _) = (f(ξ1, ξ2)..., 0)
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            x1 = @view vgeo[:, :, _x1, :]
            x2 = @view vgeo[:, :, _x2, :]

            (Cx1, Cx2) = (zeros(FT, Nq...), zeros(FT, Nq...))

            J = vgeo[:, :, _M, :] ./ reshape(kron(reverse(ω)...), Nq..., 1)
            ξ1x1 = @view vgeo[:, :, _ξ1x1, :]
            ξ1x2 = @view vgeo[:, :, _ξ1x2, :]
            ξ2x1 = @view vgeo[:, :, _ξ2x1, :]
            ξ2x2 = @view vgeo[:, :, _ξ2x2, :]

            e = 1
            for n in 1:Nq[2]
                Cx1[:, n] += D[1] * (J[:, n, e] .* ξ1x1[:, n, e])
                Cx2[:, n] += D[1] * (J[:, n, e] .* ξ1x2[:, n, e])
            end
            for n in 1:Nq[1]
                Cx1[n, :] += D[2] * (J[n, :, e] .* ξ2x1[n, :, e])
                Cx2[n, :] += D[2] * (J[n, :, e] .* ξ2x2[n, :, e])
            end
            @test maximum(abs.(Cx1)) ≤ 100 * eps(FT)
            @test maximum(abs.(Cx2)) ≤ 100 * eps(FT)
        end
        #}}}
    end

    #N = 0 test
    #{{{
    let
        for FT in (Float32, Float64)
            N = (2, 0)
            Nq = N .+ 1
            Np = prod(Nq)
            Nfp = div.(Np, Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(
                j ->
                    Nq[j] == 1 ? Elements.glpoints(FT, N[j]) :
                    Elements.lglpoints(FT, N[j]),
                dim,
            )
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            fx1(ξ1, ξ2) = ξ1 + (1 + ξ1)^2 * ξ2 / 10
            fx1ξ1(ξ1, ξ2) = 1 + 2 * (1 + ξ1) * ξ2 / 10
            fx1ξ2(ξ1, ξ2) = (1 + ξ1)^2 / 10
            fx2(ξ1, ξ2) = ξ2 - (1 + ξ1)^2
            fx2ξ1(ξ1, ξ2) = -2 * (1 + ξ1)
            fx2ξ2(ξ1, ξ2) = 1

            e2c = Array{FT, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            # Create the metrics
            (x1, x2, x1ξ1, x1ξ2, x2ξ1, x2ξ2) = let
                ξ1 = zeros(FT, Np, nelem)
                ξ2 = zeros(FT, Np, nelem)
                Metrics.creategrid!(ξ1, ξ2, e2c, ξ...)
                (
                    fx1.(ξ1, ξ2),
                    fx2.(ξ1, ξ2),
                    fx1ξ1.(ξ1, ξ2),
                    fx1ξ2.(ξ1, ξ2),
                    fx2ξ1.(ξ1, ξ2),
                    fx2ξ2.(ξ1, ξ2),
                )
            end
            J = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1)

            M = J .* reshape(kron(reverse(ω)...), Nq..., 1)

            meshwarp(ξ1, ξ2, _) = (fx1(ξ1, ξ2), fx2(ξ1, ξ2), 0)
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            @test x1 ≈ vgeo[:, :, _x1, :]
            @test x2 ≈ vgeo[:, :, _x2, :]

            @test M ≈ vgeo[:, :, _M, :]
            @test (@view vgeo[:, :, _ξ2x1, :]) .* J ≈ -x2ξ1
            @test (@view vgeo[:, :, _ξ2x2, :]) .* J ≈ x1ξ1
            @test (@view vgeo[:, :, _ξ1x1, :]) .* J ≈ x2ξ2
            @test (@view vgeo[:, :, _ξ1x2, :]) .* J ≈ -x1ξ2

            # check the normals?
            sM = @view sgeo[_sM, :, :, :]
            n1 = @view sgeo[_n1, :, :, :]
            n2 = @view sgeo[_n2, :, :, :]
            @test all(hypot.(n1[1:Nfp[1], 1:2, :], n2[1:Nfp[1], 1:2, :]) .≈ 1)
            @test all(hypot.(n1[1:Nfp[2], 3:4, :], n2[1:Nfp[2], 3:4, :]) .≈ 1)
            @test sM[1:Nfp[1], 1, :] .* n1[1:Nfp[1], 1, :] ≈
                  -x2ξ2[1, :, :] .* ω[2]
            @test sM[1:Nfp[1], 1, :] .* n2[1:Nfp[1], 1, :] ≈
                  x1ξ2[1, :, :] .* ω[2]
            @test sM[1:Nfp[1], 2, :] .* n1[1:Nfp[1], 2, :] ≈
                  x2ξ2[Nq[1], :, :] .* ω[2]
            @test sM[1:Nfp[1], 2, :] .* n2[1:Nfp[1], 2, :] ≈
                  -x1ξ2[Nq[1], :, :] .* ω[2]

            # for these faces we need the N = 1 metrics
            (x1ξ1, x2ξ1) = let
                @assert Nq[2] == 1 && Nq[1] != 1
                Nq_N1 = max.(2, Nq)
                ξ1 = zeros(FT, Nq_N1..., nelem)
                ξ2 = zeros(FT, Nq_N1..., nelem)
                Metrics.creategrid!(
                    ξ1,
                    ξ2,
                    e2c,
                    ξ[1],
                    Elements.lglpoints(FT, 1)[1],
                )
                (fx1ξ1.(ξ1, ξ2), fx2ξ1.(ξ1, ξ2))
            end
            @test sM[1:Nfp[2], 3, :] .* n1[1:Nfp[2], 3, :] ≈
                  x2ξ1[:, 1, :] .* ω[1]
            @test sM[1:Nfp[2], 3, :] .* n2[1:Nfp[2], 3, :] ≈
                  -x1ξ1[:, 1, :] .* ω[1]
            @test sM[1:Nfp[2], 4, :] .* n1[1:Nfp[2], 4, :] ≈
                  -x2ξ1[:, 2, :] .* ω[1]
            @test sM[1:Nfp[2], 4, :] .* n2[1:Nfp[2], 4, :] ≈
                  x1ξ1[:, 2, :] .* ω[1]
        end
    end
    #}}}

    # Constant preserving test for N = 0
    #{{{
    let
        for FT in (Float32, Float64), N in ((4, 0), (0, 4))
            Nq = N .+ 1
            Np = prod(Nq)
            Nfp = div.(Np, Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(
                j ->
                    Nq[j] == 1 ? Elements.glpoints(FT, N[j]) :
                    Elements.lglpoints(FT, N[j]),
                dim,
            )
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            rng = MersenneTwister(777)
            fx1(ξ1, ξ2) = ξ1 + (ξ1 * ξ2 * rand(rng) + rand(rng)) / 10
            fx2(ξ1, ξ2) = ξ2 + (ξ1 * ξ2 * rand(rng) + rand(rng)) / 10

            e2c = Array{FT, 3}(undef, 2, 4, 1)
            e2c[:, :, 1] = [-1 1 -1 1; -1 -1 1 1]
            nelem = size(e2c, 3)

            meshwarp(ξ1, ξ2, _) = (fx1(ξ1, ξ2), fx2(ξ1, ξ2), 0)
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)

            M = vgeo[:, _M, :]
            ξ1x1 = vgeo[:, _ξ1x1, :]
            ξ2x1 = vgeo[:, _ξ2x1, :]
            ξ1x2 = vgeo[:, _ξ1x2, :]
            ξ2x2 = vgeo[:, _ξ2x2, :]

            I1 = Matrix(I, Nq[1], Nq[1])
            I2 = Matrix(I, Nq[2], Nq[2])
            D1 = kron(I2, D[1])
            D2 = kron(D[2], I1)

            # Face interpolation operators
            L = (
                kron(I2, I1[1, :]'),
                kron(I2, I1[Nq[1], :]'),
                kron(I2[1, :]', I1),
                kron(I2[Nq[2], :]', I1),
            )
            sM = ntuple(f -> sgeo[_sM, 1:Nfp[cld(f, 2)], f, :], nface)
            n1 = ntuple(f -> sgeo[_n1, 1:Nfp[cld(f, 2)], f, :], nface)
            n2 = ntuple(f -> sgeo[_n2, 1:Nfp[cld(f, 2)], f, :], nface)

            # If constant preserving then:
            #   \sum_{j} = D' * M * ξjxk = \sum_{f} L_f' * sM_f * n1_f
            @test D1' * (M .* ξ1x1) + D2' * (M .* ξ2x1) ≈
                  mapreduce((L, sM, n1) -> L' * (sM .* n1), +, L, sM, n1)
            @test D1' * (M .* ξ1x2) + D2' * (M .* ξ2x2) ≈
                  mapreduce((L, sM, n2) -> L' * (sM .* n2), +, L, sM, n2)
        end
    end
    #}}}
end

@testset "3-D Metric terms" begin
    # linear test
    #{{{
    for FT in (Float32, Float64), N in ((2, 2, 2), (2, 3, 4), (4, 3, 2))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        dim = length(N)
        nface = 2dim

        # Create element operators for each polynomial order
        ξω = ntuple(j -> Elements.lglpoints(FT, N[j]), dim)
        ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
        D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

        e2c = Array{FT, 3}(undef, dim, 8, 2)
        e2c[:, :, 1] = [
            0 2 0 2 0 2 0 2
            0 0 2 2 0 0 2 2
            0 0 0 0 2 2 2 2
        ]
        e2c[:, :, 2] = [
            2 2 0 0 2 2 0 0
            0 2 0 2 0 2 0 2
            0 0 0 0 2 2 2 2
        ]

        nelem = size(e2c, 3)

        x_exact = Array{FT, 4}(undef, Nq..., nelem)
        x_exact[:, :, :, 1] .= 1 .+ ξ[1]
        x_exact[:, :, :, 2] .= 1 .- ξ[2]'

        ξ1x1_exact = zeros(Int, Nq..., nelem)
        ξ1x1_exact[:, :, :, 1] .= 1

        ξ1x2_exact = zeros(Int, Nq..., nelem)
        ξ1x2_exact[:, :, :, 2] .= 1

        ξ2x1_exact = zeros(Int, Nq..., nelem)
        ξ2x1_exact[:, :, :, 2] .= -1

        ξ2x2_exact = zeros(Int, Nq..., nelem)
        ξ2x2_exact[:, :, :, 1] .= 1

        ξ3x3_exact = ones(Int, Nq..., nelem)

        y_exact = Array{FT, 4}(undef, Nq..., nelem)
        y_exact[:, :, :, 1] .= 1 .+ ξ[2]'
        y_exact[:, :, :, 2] .= 1 .+ ξ[1]

        z_exact = Array{FT, 4}(undef, Nq..., nelem)
        z_exact[:, :, :, :] .= reshape(1 .+ ξ[3], 1, 1, Nq[3])

        M_exact =
            ones(Int, Nq..., nelem) .* reshape(kron(reverse(ω)...), Nq..., 1)

        sJ_exact = ones(Int, maximum(Nfp), nface, nelem)

        nx_exact = zeros(Int, maximum(Nfp), nface, nelem)
        nx_exact[:, 1, 1] .= -1
        nx_exact[:, 2, 1] .= 1
        nx_exact[:, 3, 2] .= 1
        nx_exact[:, 4, 2] .= -1

        ny_exact = zeros(Int, maximum(Nfp), nface, nelem)
        ny_exact[:, 3, 1] .= -1
        ny_exact[:, 4, 1] .= 1
        ny_exact[:, 1, 2] .= -1
        ny_exact[:, 2, 2] .= 1

        nz_exact = zeros(Int, maximum(Nfp), nface, nelem)
        nz_exact[:, 5, 1:2] .= -1
        nz_exact[:, 6, 1:2] .= 1

        (vgeo, sgeo, _) =
            Grids.computegeometry(e2c, D, ξ, ω, (x...) -> identity(x))
        vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)

        @test (@view vgeo[:, :, :, _x1, :]) ≈ x_exact
        @test (@view vgeo[:, :, :, _x2, :]) ≈ y_exact
        @test (@view vgeo[:, :, :, _x3, :]) ≈ z_exact
        @test (@view vgeo[:, :, :, _M, :]) ≈ M_exact
        @test (@view vgeo[:, :, :, _ξ1x1, :]) ≈ ξ1x1_exact
        @test (@view vgeo[:, :, :, _ξ1x2, :]) ≈ ξ1x2_exact
        @test maximum(abs.(@view vgeo[:, :, :, _ξ1x3, :])) ≤ 100 * eps(FT)
        @test (@view vgeo[:, :, :, _ξ2x1, :]) ≈ ξ2x1_exact
        @test (@view vgeo[:, :, :, _ξ2x2, :]) ≈ ξ2x2_exact
        @test maximum(abs.(@view vgeo[:, :, :, _ξ2x3, :])) ≤ 100 * eps(FT)
        @test maximum(abs.(@view vgeo[:, :, :, _ξ3x1, :])) ≤ 100 * eps(FT)
        @test maximum(abs.(@view vgeo[:, :, :, _ξ3x2, :])) ≤ 100 * eps(FT)
        @test (@view vgeo[:, :, :, _ξ3x3, :]) ≈ ξ3x3_exact
        for d in 1:dim
            for f in (2d - 1):(2d)
                ωf = ntuple(j -> ω[mod1(d + j, dim)], dim - 1)
                if !(dim == 3 && d == 2)
                    ωf = reverse(ωf)
                end
                Mf = kron(1, ωf...)
                @test isapprox(
                    (@view sgeo[_sM, 1:Nfp[d], f, :]),
                    sJ_exact[1:Nfp[d], f, :] .* Mf,
                    atol = √eps(FT),
                    rtol = √eps(FT),
                )
                @test isapprox(
                    (@view sgeo[_n1, 1:Nfp[d], f, :]),
                    nx_exact[1:Nfp[d], f, :];
                    atol = √eps(FT),
                    rtol = √eps(FT),
                )
                @test isapprox(
                    (@view sgeo[_n2, 1:Nfp[d], f, :]),
                    ny_exact[1:Nfp[d], f, :];
                    atol = √eps(FT),
                    rtol = √eps(FT),
                )
                @test isapprox(
                    (@view sgeo[_n3, 1:Nfp[d], f, :]),
                    nz_exact[1:Nfp[d], f, :];
                    atol = √eps(FT),
                    rtol = √eps(FT),
                )
            end
        end
    end
    #}}}

    # Polynomial 3-D test
    #{{{
    for FT in (Float32, Float64), N in ((9, 9, 9), (9, 9, 10), (10, 9, 11))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        dim = length(N)
        nface = 2dim

        # Create element operators for each polynomial order
        ξω = ntuple(j -> Elements.lglpoints(FT, N[j]), dim)
        ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
        D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

        f(ξ1, ξ2, ξ3) = @.( (
            ξ2 + ξ1 * ξ3 - (ξ1^2 * ξ2^2 * ξ3^2) / 4,
            ξ3 - ((ξ1 * ξ2 * ξ3 + 1) / 2)^3 + 1,
            ξ1 + 2 * ((ξ1 + 1) / 2)^6 * ((ξ2 + 1) / 2)^6 * ((ξ3 + 1) / 2)^6,
        ))

        fx1ξ1(ξ1, ξ2, ξ3) = @.(ξ3 - (ξ1 * ξ2^2 * ξ3^2) / 2)
        fx1ξ2(ξ1, ξ2, ξ3) = @.(1 - (ξ1^2 * ξ2 * ξ3^2) / 2)
        fx1ξ3(ξ1, ξ2, ξ3) = @.(ξ1 - (ξ1^2 * ξ2^2 * ξ3) / 2)
        fx2ξ1(ξ1, ξ2, ξ3) = @.(-(3 * ξ2 * ξ3 * ((ξ1 * ξ2 * ξ3 + 1) / 2)^2) / 2)
        fx2ξ2(ξ1, ξ2, ξ3) = @.(-(3 * ξ1 * ξ3 * ((ξ1 * ξ2 * ξ3 + 1) / 2)^2) / 2)
        fx2ξ3(ξ1, ξ2, ξ3) =
            @.(1 - (3 * ξ1 * ξ2 * ((ξ1 * ξ2 * ξ3 + 1) / 2)^2) / 2)
        fx3ξ1(ξ1, ξ2, ξ3) =
            @.(6 * ((ξ1 + 1) / 2)^5 * ((ξ2 + 1) / 2)^6 * ((ξ3 + 1) / 2)^6 + 1)
        fx3ξ2(ξ1, ξ2, ξ3) =
            @.(6 * ((ξ1 + 1) / 2)^6 * ((ξ2 + 1) / 2)^5 * ((ξ3 + 1) / 2)^6)
        fx3ξ3(ξ1, ξ2, ξ3) =
            @.(6 * ((ξ1 + 1) / 2)^6 * ((ξ2 + 1) / 2)^6 * ((ξ3 + 1) / 2)^5)

        e2c = Array{FT, 3}(undef, dim, 8, 1)
        e2c[:, :, 1] = [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            -1 -1 -1 -1 1 1 1 1
        ]

        nelem = size(e2c, 3)

        # Compute exact metrics
        (x1ξ1, x1ξ2, x1ξ3, x2ξ1, x2ξ2, x2ξ3, x3ξ1, x3ξ2, x3ξ3) = let
            (vgeo, _) =
                Grids.computegeometry(e2c, D, ξ, ω, (x...) -> identity(x))
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)
            ξ1 = vgeo[:, :, :, _x1, :]
            ξ2 = vgeo[:, :, :, _x2, :]
            ξ3 = vgeo[:, :, :, _x3, :]
            (
                fx1ξ1(ξ1, ξ2, ξ3),
                fx1ξ2(ξ1, ξ2, ξ3),
                fx1ξ3(ξ1, ξ2, ξ3),
                fx2ξ1(ξ1, ξ2, ξ3),
                fx2ξ2(ξ1, ξ2, ξ3),
                fx2ξ3(ξ1, ξ2, ξ3),
                fx3ξ1(ξ1, ξ2, ξ3),
                fx3ξ2(ξ1, ξ2, ξ3),
                fx3ξ3(ξ1, ξ2, ξ3),
            )
        end
        J = (
            x1ξ1 .* (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) +
            x2ξ1 .* (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) +
            x3ξ1 .* (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2)
        )

        ξ1x1 = (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) ./ J
        ξ1x2 = (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) ./ J
        ξ1x3 = (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2) ./ J
        ξ2x1 = (x2ξ3 .* x3ξ1 - x2ξ1 .* x3ξ3) ./ J
        ξ2x2 = (x3ξ3 .* x1ξ1 - x3ξ1 .* x1ξ3) ./ J
        ξ2x3 = (x1ξ3 .* x2ξ1 - x1ξ1 .* x2ξ3) ./ J
        ξ3x1 = (x2ξ1 .* x3ξ2 - x2ξ2 .* x3ξ1) ./ J
        ξ3x2 = (x3ξ1 .* x1ξ2 - x3ξ2 .* x1ξ1) ./ J
        ξ3x3 = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1) ./ J

        (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, f)
        vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)

        @test (@view vgeo[:, :, :, _M, :]) ≈
              J .* reshape(kron(reverse(ω)...), Nq..., 1)
        @test (@view vgeo[:, :, :, _ξ1x1, :]) ≈ ξ1x1
        @test (@view vgeo[:, :, :, _ξ1x2, :]) ≈ ξ1x2
        @test (@view vgeo[:, :, :, _ξ1x3, :]) ≈ ξ1x3
        @test (@view vgeo[:, :, :, _ξ2x1, :]) ≈ ξ2x1
        @test (@view vgeo[:, :, :, _ξ2x2, :]) ≈ ξ2x2
        @test (@view vgeo[:, :, :, _ξ2x3, :]) ≈ ξ2x3
        @test (@view vgeo[:, :, :, _ξ3x1, :]) ≈ ξ3x1
        @test (@view vgeo[:, :, :, _ξ3x2, :]) ≈ ξ3x2
        @test (@view vgeo[:, :, :, _ξ3x3, :]) ≈ ξ3x3
        n1 = @view sgeo[_n1, :, :, :]
        n2 = @view sgeo[_n2, :, :, :]
        n3 = @view sgeo[_n3, :, :, :]
        sM = @view sgeo[_sM, :, :, :]
        for d in 1:dim
            for f in (2d - 1):(2d)
                @test all(
                    hypot.(
                        n1[1:Nfp[d], f, :],
                        n2[1:Nfp[d], f, :],
                        n3[1:Nfp[d], f, :],
                    ) .≈ 1,
                )
            end
        end
        d, f = 1, 1
        Mf = kron(1, ω[3], ω[2])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (-J[1, :, :, :] .* ξ1x1[1, :, :, :])[:] .* Mf,
            (-J[1, :, :, :] .* ξ1x2[1, :, :, :])[:] .* Mf,
            (-J[1, :, :, :] .* ξ1x3[1, :, :, :])[:] .* Mf,
        ]
        d, f = 1, 2
        Mf = kron(1, ω[3], ω[2])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (J[Nq[d], :, :, :] .* ξ1x1[Nq[d], :, :, :])[:] .* Mf,
            (J[Nq[d], :, :, :] .* ξ1x2[Nq[d], :, :, :])[:] .* Mf,
            (J[Nq[d], :, :, :] .* ξ1x3[Nq[d], :, :, :])[:] .* Mf,
        ]
        d, f = 2, 3
        Mf = kron(1, ω[3], ω[1])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (-J[:, 1, :, :] .* ξ2x1[:, 1, :, :])[:] .* Mf,
            (-J[:, 1, :, :] .* ξ2x2[:, 1, :, :])[:] .* Mf,
            (-J[:, 1, :, :] .* ξ2x3[:, 1, :, :])[:] .* Mf,
        ]
        d, f = 2, 4
        Mf = kron(1, ω[3], ω[1])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (J[:, Nq[d], :, :] .* ξ2x1[:, Nq[d], :, :])[:] .* Mf,
            (J[:, Nq[d], :, :] .* ξ2x2[:, Nq[d], :, :])[:] .* Mf,
            (J[:, Nq[d], :, :] .* ξ2x3[:, Nq[d], :, :])[:] .* Mf,
        ]
        d, f = 3, 5
        Mf = kron(1, ω[2], ω[1])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (-J[:, :, 1, :] .* ξ3x1[:, :, 1, :])[:] .* Mf,
            (-J[:, :, 1, :] .* ξ3x2[:, :, 1, :])[:] .* Mf,
            (-J[:, :, 1, :] .* ξ3x3[:, :, 1, :])[:] .* Mf,
        ]
        d, f = 3, 6
        Mf = kron(1, ω[2], ω[1])
        @test [
            (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
            (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
        ] ≈ [
            (J[:, :, Nq[d], :] .* ξ3x1[:, :, Nq[d], :])[:] .* Mf,
            (J[:, :, Nq[d], :] .* ξ3x2[:, :, Nq[d], :])[:] .* Mf,
            (J[:, :, Nq[d], :] .* ξ3x3[:, :, Nq[d], :])[:] .* Mf,
        ]
    end
    #}}}

    # Constant preserving test
    #{{{
    for FT in (Float32, Float64), N in ((5, 5, 5), (3, 4, 5), (4, 4, 5))
        Nq = N .+ 1
        Np = prod(Nq)
        Nfp = div.(Np, Nq)

        dim = length(N)
        nface = 2dim

        # Create element operators for each polynomial order
        ξω = ntuple(j -> Elements.lglpoints(FT, N[j]), dim)
        ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
        D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

        f(ξ1, ξ2, ξ3) = @.( (
            ξ2 + ξ1 * ξ3 - (ξ1^2 * ξ2^2 * ξ3^2) / 4,
            ξ3 - ((ξ1 * ξ2 * ξ3 + 1) / 2)^3 + 1,
            ξ1 + ((ξ1 + 1) / 2)^6 * ((ξ2 + 1) / 2)^6 * ((ξ3 + 1) / 2)^6,
        ))

        e2c = Array{FT, 3}(undef, dim, 8, 1)
        e2c[:, :, 1] = [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            -1 -1 -1 -1 1 1 1 1
        ]

        nelem = size(e2c, 3)

        (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, f)
        vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)

        (Cx1, Cx2, Cx3) = (zeros(FT, Nq...), zeros(FT, Nq...), zeros(FT, Nq...))

        J =
            (@view vgeo[:, :, :, _M, :]) ./
            reshape(kron(reverse(ω)...), Nq..., 1)
        ξ1x1 = @view vgeo[:, :, :, _ξ1x1, :]
        ξ1x2 = @view vgeo[:, :, :, _ξ1x2, :]
        ξ1x3 = @view vgeo[:, :, :, _ξ1x3, :]
        ξ2x1 = @view vgeo[:, :, :, _ξ2x1, :]
        ξ2x2 = @view vgeo[:, :, :, _ξ2x2, :]
        ξ2x3 = @view vgeo[:, :, :, _ξ2x3, :]
        ξ3x1 = @view vgeo[:, :, :, _ξ3x1, :]
        ξ3x2 = @view vgeo[:, :, :, _ξ3x2, :]
        ξ3x3 = @view vgeo[:, :, :, _ξ3x3, :]

        e = 1
        for k in 1:Nq[3]
            for j in 1:Nq[2]
                Cx1[:, j, k] += D[1] * (J[:, j, k, e] .* ξ1x1[:, j, k, e])
                Cx2[:, j, k] += D[1] * (J[:, j, k, e] .* ξ1x2[:, j, k, e])
                Cx3[:, j, k] += D[1] * (J[:, j, k, e] .* ξ1x3[:, j, k, e])
            end
        end

        for k in 1:Nq[3]
            for i in 1:Nq[1]
                Cx1[i, :, k] += D[2] * (J[i, :, k, e] .* ξ2x1[i, :, k, e])
                Cx2[i, :, k] += D[2] * (J[i, :, k, e] .* ξ2x2[i, :, k, e])
                Cx3[i, :, k] += D[2] * (J[i, :, k, e] .* ξ2x3[i, :, k, e])
            end
        end

        for j in 1:Nq[2]
            for i in 1:Nq[1]
                Cx1[i, j, :] += D[3] * (J[i, j, :, e] .* ξ3x1[i, j, :, e])
                Cx2[i, j, :] += D[3] * (J[i, j, :, e] .* ξ3x2[i, j, :, e])
                Cx3[i, j, :] += D[3] * (J[i, j, :, e] .* ξ3x3[i, j, :, e])
            end
        end
        @test maximum(abs.(Cx1)) ≤ 300 * eps(FT)
        @test maximum(abs.(Cx2)) ≤ 300 * eps(FT)
        @test maximum(abs.(Cx3)) ≤ 300 * eps(FT)
    end
    #}}}

    #N = 0 test
    #{{{
    let
        for FT in (Float32, Float64)
            N = (4, 4, 0)
            Nq = N .+ 1
            Np = prod(Nq)
            Nfp = div.(Np, Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(
                j ->
                    Nq[j] == 1 ? Elements.glpoints(FT, N[j]) :
                    Elements.lglpoints(FT, N[j]),
                dim,
            )
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            fx1(ξ1, ξ2, ξ3) = ξ1 + (1 + ξ1)^2 * (1 + ξ2)^2 + ξ3 / 10
            fx1ξ1(ξ1, ξ2, ξ3) = 1 + 2 * (1 + ξ1) * (1 + ξ2)^2
            fx1ξ2(ξ1, ξ2, ξ3) = (1 + ξ1)^2 * 2 * (1 + ξ2)
            fx1ξ3(ξ1, ξ2, ξ3) = 1 / 10

            fx2(ξ1, ξ2, ξ3) = ξ2 - (1 + ξ1)^2 + (2 + ξ3) / 2
            fx2ξ1(ξ1, ξ2, ξ3) = -2 * (1 + ξ1)
            fx2ξ2(ξ1, ξ2, ξ3) = 1
            fx2ξ3(ξ1, ξ2, ξ3) = 1 / 2

            fx3(ξ1, ξ2, ξ3) = ξ3 + (1 + ξ1)^2 * (1 + ξ2)^2 / 10
            fx3ξ1(ξ1, ξ2, ξ3) = 2 * (1 + ξ1) * (1 + ξ2)^2 / 10
            fx3ξ2(ξ1, ξ2, ξ3) = (1 + ξ1)^2 * 2 * (1 + ξ2) / 10
            fx3ξ3(ξ1, ξ2, ξ3) = 1

            e2c = Array{FT, 3}(undef, 3, 8, 1)
            e2c[:, :, 1] = [
                -1 +1 -1 +1 -1 +1 -1 +1
                -1 -1 +1 +1 -1 -1 +1 +1
                -1 -1 -1 -1 +1 +1 +1 +1
            ]
            nelem = size(e2c, 3)

            # Create the metrics
            (x1, x2, x3, x1ξ1, x1ξ2, x1ξ3, x2ξ1, x2ξ2, x2ξ3, x3ξ1, x3ξ2, x3ξ3) =
                let
                    ξ1 = zeros(FT, Nq..., nelem)
                    ξ2 = zeros(FT, Nq..., nelem)
                    ξ3 = zeros(FT, Nq..., nelem)
                    Metrics.creategrid!(ξ1, ξ2, ξ3, e2c, ξ...)
                    (
                        fx1.(ξ1, ξ2, ξ3),
                        fx2.(ξ1, ξ2, ξ3),
                        fx3.(ξ1, ξ2, ξ3),
                        fx1ξ1.(ξ1, ξ2, ξ3),
                        fx1ξ2.(ξ1, ξ2, ξ3),
                        fx1ξ3.(ξ1, ξ2, ξ3),
                        fx2ξ1.(ξ1, ξ2, ξ3),
                        fx2ξ2.(ξ1, ξ2, ξ3),
                        fx2ξ3.(ξ1, ξ2, ξ3),
                        fx3ξ1.(ξ1, ξ2, ξ3),
                        fx3ξ2.(ξ1, ξ2, ξ3),
                        fx3ξ3.(ξ1, ξ2, ξ3),
                    )
                end
            J = @.(
                x1ξ1 * (x2ξ2 * x3ξ3 - x3ξ2 * x2ξ3) +
                x2ξ1 * (x3ξ2 * x1ξ3 - x1ξ2 * x3ξ3) +
                x3ξ1 * (x1ξ2 * x2ξ3 - x2ξ2 * x1ξ3)
            )
            ξ1x1 = (x2ξ2 .* x3ξ3 - x2ξ3 .* x3ξ2) ./ J
            ξ1x2 = (x3ξ2 .* x1ξ3 - x3ξ3 .* x1ξ2) ./ J
            ξ1x3 = (x1ξ2 .* x2ξ3 - x1ξ3 .* x2ξ2) ./ J
            ξ2x1 = (x2ξ3 .* x3ξ1 - x2ξ1 .* x3ξ3) ./ J
            ξ2x2 = (x3ξ3 .* x1ξ1 - x3ξ1 .* x1ξ3) ./ J
            ξ2x3 = (x1ξ3 .* x2ξ1 - x1ξ1 .* x2ξ3) ./ J
            ξ3x1 = (x2ξ1 .* x3ξ2 - x2ξ2 .* x3ξ1) ./ J
            ξ3x2 = (x3ξ1 .* x1ξ2 - x3ξ2 .* x1ξ1) ./ J
            ξ3x3 = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1) ./ J

            M = J .* reshape(kron(reverse(ω)...), Nq..., 1)

            meshwarp(ξ1, ξ2, ξ3) =
                (fx1(ξ1, ξ2, ξ3), fx2(ξ1, ξ2, ξ3), fx3(ξ1, ξ2, ξ3))
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)
            vgeo = reshape(vgeo, Nq..., _nvgeo, nelem)

            @test x1 ≈ vgeo[:, :, :, _x1, :]
            @test x2 ≈ vgeo[:, :, :, _x2, :]
            @test x3 ≈ vgeo[:, :, :, _x3, :]

            @test M ≈ vgeo[:, :, :, _M, :]

            @test (@view vgeo[:, :, :, _ξ1x1, :]) ≈ ξ1x1
            @test (@view vgeo[:, :, :, _ξ1x2, :]) ≈ ξ1x2
            @test (@view vgeo[:, :, :, _ξ1x3, :]) ≈ ξ1x3

            @test (@view vgeo[:, :, :, _ξ2x1, :]) ≈ ξ2x1
            @test (@view vgeo[:, :, :, _ξ2x2, :]) ≈ ξ2x2
            @test (@view vgeo[:, :, :, _ξ2x3, :]) ≈ ξ2x3

            @test (@view vgeo[:, :, :, _ξ3x1, :]) ≈ ξ3x1
            @test (@view vgeo[:, :, :, _ξ3x2, :]) ≈ ξ3x2
            @test (@view vgeo[:, :, :, _ξ3x3, :]) ≈ ξ3x3

            # check the normals?
            sM = @view sgeo[_sM, :, :, :]
            n1 = @view sgeo[_n1, :, :, :]
            n2 = @view sgeo[_n2, :, :, :]
            n3 = @view sgeo[_n3, :, :, :]
            @test all(
                hypot.(
                    n1[1:Nfp[1], 1:2, :],
                    n2[1:Nfp[1], 1:2, :],
                    n3[1:Nfp[1], 1:2, :],
                ) .≈ 1,
            )
            @test all(
                hypot.(
                    n1[1:Nfp[2], 3:4, :],
                    n2[1:Nfp[2], 3:4, :],
                    n3[1:Nfp[2], 3:4, :],
                ) .≈ 1,
            )
            @test all(
                hypot.(
                    n1[1:Nfp[3], 5:6, :],
                    n2[1:Nfp[3], 5:6, :],
                    n3[1:Nfp[3], 5:6, :],
                ) .≈ 1,
            )

            d, f = 1, 1
            Mf = kron(1, ω[3], ω[2])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (-J[1, :, :, :] .* ξ1x1[1, :, :, :])[:] .* Mf,
                (-J[1, :, :, :] .* ξ1x2[1, :, :, :])[:] .* Mf,
                (-J[1, :, :, :] .* ξ1x3[1, :, :, :])[:] .* Mf,
            ]
            d, f = 1, 2
            Mf = kron(1, ω[3], ω[2])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (J[Nq[d], :, :, :] .* ξ1x1[Nq[d], :, :, :])[:] .* Mf,
                (J[Nq[d], :, :, :] .* ξ1x2[Nq[d], :, :, :])[:] .* Mf,
                (J[Nq[d], :, :, :] .* ξ1x3[Nq[d], :, :, :])[:] .* Mf,
            ]
            d, f = 2, 3
            Mf = kron(1, ω[3], ω[1])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (-J[:, 1, :, :] .* ξ2x1[:, 1, :, :])[:] .* Mf,
                (-J[:, 1, :, :] .* ξ2x2[:, 1, :, :])[:] .* Mf,
                (-J[:, 1, :, :] .* ξ2x3[:, 1, :, :])[:] .* Mf,
            ]
            d, f = 2, 4
            Mf = kron(1, ω[3], ω[1])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (J[:, Nq[d], :, :] .* ξ2x1[:, Nq[d], :, :])[:] .* Mf,
                (J[:, Nq[d], :, :] .* ξ2x2[:, Nq[d], :, :])[:] .* Mf,
                (J[:, Nq[d], :, :] .* ξ2x3[:, Nq[d], :, :])[:] .* Mf,
            ]

            # for these faces we need the N = 1 metrics
            (x1ξ1, x1ξ2, x1ξ3, x2ξ1, x2ξ2, x2ξ3, x3ξ1, x3ξ2, x3ξ3) = let
                @assert Nq[1] != 1 && Nq[2] != 1 && Nq[3] == 1
                Nq_N1 = max.(2, Nq)
                ξ1 = zeros(FT, Nq_N1..., nelem)
                ξ2 = zeros(FT, Nq_N1..., nelem)
                ξ3 = zeros(FT, Nq_N1..., nelem)
                Metrics.creategrid!(
                    ξ1,
                    ξ2,
                    ξ3,
                    e2c,
                    ξ[1],
                    ξ[2],
                    Elements.lglpoints(FT, 1)[1],
                )
                (
                    fx1ξ1.(ξ1, ξ2, ξ3),
                    fx1ξ2.(ξ1, ξ2, ξ3),
                    fx1ξ3.(ξ1, ξ2, ξ3),
                    fx2ξ1.(ξ1, ξ2, ξ3),
                    fx2ξ2.(ξ1, ξ2, ξ3),
                    fx2ξ3.(ξ1, ξ2, ξ3),
                    fx3ξ1.(ξ1, ξ2, ξ3),
                    fx3ξ2.(ξ1, ξ2, ξ3),
                    fx3ξ3.(ξ1, ξ2, ξ3),
                )
            end
            J = @.(
                x1ξ1 * (x2ξ2 * x3ξ3 - x3ξ2 * x2ξ3) +
                x2ξ1 * (x3ξ2 * x1ξ3 - x1ξ2 * x3ξ3) +
                x3ξ1 * (x1ξ2 * x2ξ3 - x2ξ2 * x1ξ3)
            )
            ξ3x1 = (x2ξ1 .* x3ξ2 - x2ξ2 .* x3ξ1) ./ J
            ξ3x2 = (x3ξ1 .* x1ξ2 - x3ξ2 .* x1ξ1) ./ J
            ξ3x3 = (x1ξ1 .* x2ξ2 - x1ξ2 .* x2ξ1) ./ J

            d, f = 3, 5
            Mf = kron(1, ω[2], ω[1])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (-J[:, :, 1, :] .* ξ3x1[:, :, 1, :])[:] .* Mf,
                (-J[:, :, 1, :] .* ξ3x2[:, :, 1, :])[:] .* Mf,
                (-J[:, :, 1, :] .* ξ3x3[:, :, 1, :])[:] .* Mf,
            ]

            d, f = 3, 6
            Mf = kron(1, ω[2], ω[1])
            @test [
                (sM[1:Nfp[d], f, :] .* n1[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n2[1:Nfp[d], f, :])[:],
                (sM[1:Nfp[d], f, :] .* n3[1:Nfp[d], f, :])[:],
            ] ≈ [
                (J[:, :, 2, :] .* ξ3x1[:, :, 2, :])[:] .* Mf,
                (J[:, :, 2, :] .* ξ3x2[:, :, 2, :])[:] .* Mf,
                (J[:, :, 2, :] .* ξ3x3[:, :, 2, :])[:] .* Mf,
            ]
        end
    end
    #}}}

    # Constant preserving test for N = 0
    #{{{
    let
        for FT in (Float64, Float32),
            N in ((4, 4, 0), (0, 0, 2), (0, 3, 4), (2, 0, 3))

            Nq = N .+ 1

            Np = prod(Nq)
            Nfp = div.(Np, Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(
                j ->
                    Nq[j] == 1 ? Elements.glpoints(FT, N[j]) :
                    Elements.lglpoints(FT, N[j]),
                dim,
            )
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            rng = MersenneTwister(777)
            fx1(ξ1, ξ2, ξ3) = ξ1 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10
            fx2(ξ1, ξ2, ξ3) = ξ2 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10
            fx3(ξ1, ξ2, ξ3) = ξ3 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10

            e2c = Array{FT, 3}(undef, 3, 8, 1)
            e2c[:, :, 1] = [
                -1 +1 -1 +1 -1 +1 -1 +1
                -1 -1 +1 +1 -1 -1 +1 +1
                -1 -1 -1 -1 +1 +1 +1 +1
            ]
            nelem = size(e2c, 3)

            meshwarp(ξ1, ξ2, ξ3) =
                (fx1(ξ1, ξ2, ξ3), fx2(ξ1, ξ2, ξ3), fx2(ξ1, ξ2, ξ3))
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)

            M = vgeo[:, _M, :]
            ξ1x1 = vgeo[:, _ξ1x1, :]
            ξ2x1 = vgeo[:, _ξ2x1, :]
            ξ3x1 = vgeo[:, _ξ3x1, :]
            ξ1x2 = vgeo[:, _ξ1x2, :]
            ξ2x2 = vgeo[:, _ξ2x2, :]
            ξ3x2 = vgeo[:, _ξ3x2, :]
            ξ1x3 = vgeo[:, _ξ1x3, :]
            ξ2x3 = vgeo[:, _ξ2x3, :]
            ξ3x3 = vgeo[:, _ξ3x3, :]

            M = vgeo[:, _M, :]
            ξ1x1 = vgeo[:, _ξ1x1, :]
            ξ2x1 = vgeo[:, _ξ2x1, :]
            ξ1x2 = vgeo[:, _ξ1x2, :]
            ξ2x2 = vgeo[:, _ξ2x2, :]

            I1 = Matrix(I, Nq[1], Nq[1])
            I2 = Matrix(I, Nq[2], Nq[2])
            I3 = Matrix(I, Nq[3], Nq[3])
            D1 = kron(I3, I2, D[1])
            D2 = kron(I3, D[2], I1)
            D3 = kron(D[3], I2, I1)

            # Face interpolation operators
            L = (
                kron(I3, I2, I1[1, :]'),
                kron(I3, I2, I1[Nq[1], :]'),
                kron(I3, I2[1, :]', I1),
                kron(I3, I2[Nq[2], :]', I1),
                kron(I3[1, :]', I2, I1),
                kron(I3[Nq[3], :]', I2, I1),
            )
            sM = ntuple(f -> sgeo[_sM, 1:Nfp[cld(f, 2)], f, :], nface)
            n1 = ntuple(f -> sgeo[_n1, 1:Nfp[cld(f, 2)], f, :], nface)
            n2 = ntuple(f -> sgeo[_n2, 1:Nfp[cld(f, 2)], f, :], nface)
            n3 = ntuple(f -> sgeo[_n3, 1:Nfp[cld(f, 2)], f, :], nface)

            # If constant preserving then:
            #   \sum_{j} = D' * M * ξjxk = \sum_{f} L_f' * sM_f * n1_f
            @test D1' * (M .* ξ1x1) + D2' * (M .* ξ2x1) + D3' * (M .* ξ3x1) ≈
                  mapreduce((L, sM, n1) -> L' * (sM .* n1), +, L, sM, n1)
            @test D1' * (M .* ξ1x2) + D2' * (M .* ξ2x2) + D3' * (M .* ξ3x2) ≈
                  mapreduce((L, sM, n2) -> L' * (sM .* n2), +, L, sM, n2)
            @test D1' * (M .* ξ1x3) + D2' * (M .* ξ2x3) + D3' * (M .* ξ3x3) ≈
                  mapreduce((L, sM, n3) -> L' * (sM .* n3), +, L, sM, n3)
        end
    end
    #}}}

    # Constant preserving test with all N = 0
    #{{{
    let
        for FT in (Float64, Float32)
            N = (0, 0, 0)
            Nq = N .+ 1

            Np = prod(Nq)
            Nfp = div.(Np, Nq)

            dim = length(N)
            nface = 2dim

            # Create element operators for each polynomial order
            ξω = ntuple(
                j ->
                    Nq[j] == 1 ? Elements.glpoints(FT, N[j]) :
                    Elements.lglpoints(FT, N[j]),
                dim,
            )
            ξ, ω = ntuple(j -> map(x -> x[j], ξω), 2)
            D = ntuple(j -> Elements.spectralderivative(ξ[j]), dim)

            rng = MersenneTwister(777)
            fx1(ξ1, ξ2, ξ3) = ξ1 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10
            fx2(ξ1, ξ2, ξ3) = ξ2 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10
            fx3(ξ1, ξ2, ξ3) = ξ3 + (ξ1 * ξ2 * ξ3 * rand(rng) + rand(rng)) / 10

            e2c = Array{FT, 3}(undef, 3, 8, 1)
            e2c[:, :, 1] = [
                -1 +1 -1 +1 -1 +1 -1 +1
                -1 -1 +1 +1 -1 -1 +1 +1
                -1 -1 -1 -1 +1 +1 +1 +1
            ]
            nelem = size(e2c, 3)

            meshwarp(ξ1, ξ2, ξ3) =
                (fx1(ξ1, ξ2, ξ3), fx2(ξ1, ξ2, ξ3), fx2(ξ1, ξ2, ξ3))
            (vgeo, sgeo, _) = Grids.computegeometry(e2c, D, ξ, ω, meshwarp)

            M = vgeo[:, _M, :]
            ξ1x1 = vgeo[:, _ξ1x1, :]
            ξ2x1 = vgeo[:, _ξ2x1, :]
            ξ3x1 = vgeo[:, _ξ3x1, :]
            ξ1x2 = vgeo[:, _ξ1x2, :]
            ξ2x2 = vgeo[:, _ξ2x2, :]
            ξ3x2 = vgeo[:, _ξ3x2, :]
            ξ1x3 = vgeo[:, _ξ1x3, :]
            ξ2x3 = vgeo[:, _ξ2x3, :]
            ξ3x3 = vgeo[:, _ξ3x3, :]

            M = vgeo[:, _M, :]
            ξ1x1 = vgeo[:, _ξ1x1, :]
            ξ2x1 = vgeo[:, _ξ2x1, :]
            ξ1x2 = vgeo[:, _ξ1x2, :]
            ξ2x2 = vgeo[:, _ξ2x2, :]

            I1 = Matrix(I, Nq[1], Nq[1])
            I2 = Matrix(I, Nq[2], Nq[2])
            I3 = Matrix(I, Nq[3], Nq[3])
            D1 = kron(I3, I2, D[1])
            D2 = kron(I3, D[2], I1)
            D3 = kron(D[3], I2, I1)

            # Face interpolation operators
            L = (
                kron(I3, I2, I1[1, :]'),
                kron(I3, I2, I1[Nq[1], :]'),
                kron(I3, I2[1, :]', I1),
                kron(I3, I2[Nq[2], :]', I1),
                kron(I3[1, :]', I2, I1),
                kron(I3[Nq[3], :]', I2, I1),
            )
            sM = ntuple(f -> sgeo[_sM, 1:Nfp[cld(f, 2)], f, :], nface)
            n1 = ntuple(f -> sgeo[_n1, 1:Nfp[cld(f, 2)], f, :], nface)
            n2 = ntuple(f -> sgeo[_n2, 1:Nfp[cld(f, 2)], f, :], nface)
            n3 = ntuple(f -> sgeo[_n3, 1:Nfp[cld(f, 2)], f, :], nface)

            # If constant preserving then \sum_{f} L_f' * sM_f * n1_f ≈ 0
            @test abs(mapreduce(
                (L, sM, n1) -> L' * (sM .* n1),
                +,
                L,
                sM,
                n1,
            )[1]) < 10 * eps(FT)
            @test abs(mapreduce(
                (L, sM, n2) -> L' * (sM .* n2),
                +,
                L,
                sM,
                n2,
            )[1]) < 10 * eps(FT)
            @test abs(mapreduce(
                (L, sM, n3) -> L' * (sM .* n3),
                +,
                L,
                sM,
                n3,
            )[1]) < 10 * eps(FT)
        end
    end
    #}}}
end

ext: jl

max_forks_repo_name: jm-c/CLIMA


text:
"""
    ProductNode{T, C} <: AbstractProductNode

Data node representing a Cartesian product of several spaces each represented by subtree stored in iterable of type `T`. May store metadata of type `C`.

See also: [`AbstractProductNode`](@ref), [`AbstractMillNode`](@ref), [`ProductModel`](@ref).
"""
struct ProductNode{T, C} <: AbstractProductNode
    data::T
    metadata::C

    function ProductNode{T, C}(data::T, metadata::C) where {T, C}
        @assert(length(data) >= 1 && all(x -> nobs(x) == nobs(data[1]), data),
                "All subtrees must have an equal amount of instances!")
        new{T, C}(data, metadata)
    end
end

"""
    ProductNode(ds, m=nothing)

Construct a new [`ProductNode`](@ref) with data `ds`, and metadata `m`. `ds` should be an iterable (preferably `Tuple` or `NamedTuple`) and all its elements must contain the same number of observations.

# Examples
```jldoctest
julia> ProductNode((ArrayNode(zeros(2, 2)), ArrayNode(Flux.onehotbatch([1, 2], 1:2))))
ProductNode 	# 2 obs, 16 bytes
  ├── ArrayNode(2×2 Array with Float64 elements) 	# 2 obs, 80 bytes
  └── ArrayNode(2×2 OneHotArray with Bool elements) 	# 2 obs, 64 bytes

julia> ProductNode((x1 = ArrayNode(NGramMatrix(["Hello", "world"])),
                    x2 = BagNode(ArrayNode([1 2; 3 4]), [1:3, 4:4])))
ProductNode 	# 2 obs, 48 bytes
  ├── x1: ArrayNode(2053×2 NGramMatrix with Int64 elements) 	# 2 obs, 146 bytes
  └── x2: BagNode 	# 2 obs, 96 bytes
            └── ArrayNode(2×2 Array with Int64 elements) 	# 2 obs, 80 bytes

julia> ProductNode((ArrayNode([1 2; 3 4]), ArrayNode([1; 3])))
ERROR: AssertionError: All subtrees must have an equal amount of instances!
[...]
```

See also: [`AbstractProductNode`](@ref), [`AbstractMillNode`](@ref), [`ProductModel`](@ref).
"""
ProductNode(ds::T) where {T} = ProductNode{T, Nothing}(ds, nothing)
ProductNode(ds::T, m::C) where {T, C} = ProductNode{T, C}(ds, m)

Flux.@functor ProductNode
mapdata(f, x::ProductNode) = ProductNode(map(i -> mapdata(f, i), x.data), x.metadata)

dropmeta(x::ProductNode) = ProductNode(map(dropmeta, x.data))

Base.getindex(x::ProductNode, i::Symbol) = x.data[i]
Base.keys(x::ProductNode) = keys(x.data)

_check_idxs(as::Vector{<:Union{Vector, Tuple}}) = all(isequal(length(as[1])), (length(a) for a in as))
_check_idxs(as::Vector{<:NamedTuple{K}}) where K = true
_check_idxs(as) = false

_cattrees(as::Vector{<:Vector}) = [reduce(catobs, [a[i] for a in as]) for i in eachindex(as[1])]
_cattrees(as::Vector{T}) where T <: Tuple = T(reduce(catobs, [a[i] for a in as]) for i in eachindex(as[1]))
_cattrees(as::Vector{T}) where T <: NamedTuple = T(reduce(catobs, [a[i] for a in as]) for i in keys(as[1]))
@generated function _cattrees(xs::Vector{NamedTuple{K, T}}) where {K, T}
    es = map(K, T.parameters) do k, t
        quote $k = reduce(catobs, $t[x.$k for x in xs]) end
    end
    quote
        $(es...)
        return NamedTuple{$K}(tuple($(K...)))
    end
end

function reduce(::typeof(catobs), as::Vector{<:ProductNode})
    ds = data.(as)
    if !_check_idxs(ds)
        ArgumentError("Trying to `catobs` `ProductNode`s with different subtrees") |> throw
    end
    ProductNode(_cattrees(ds), reduce(catobs, metadata.(as)))
end

Base.getindex(x::ProductNode, i::VecOrRange{<:Int}) = ProductNode(subset(data(x), i), subset(metadata(x), i))
Base.getindex(x::ProductNode{<:Vector}, i::VecOrRange{<:Int}) = ProductNode(map(x -> getindex(x, i), data(x)), subset(metadata(x), i))

Base.hash(e::ProductNode, h::UInt) = hash((e.data, e.metadata), h)
(e1::ProductNode == e2::ProductNode) = e1.data == e2.data && e1.metadata == e2.metadata
Base.isequal(e1::ProductNode, e2::ProductNode) = isequal(e1.data, e2.data) && isequal(e1.metadata, e2.metadata)


text:
abstract type AbstractProgressiveHedgingExecution end
# Execution API #
# ------------------------------------------------------------
initialize_subproblems!(ph::AbstractProgressiveHedging, scenarioproblems::AbstractScenarioProblems, penaltyterm::AbstractPenaltyterm) = initialize_subproblems!(ph, ph.execution, scenarioproblems, penaltyterm)
finish_initilization!(ph::AbstractProgressiveHedging, penalty::AbstractFloat) = finish_initilization!(ph.execution, penalty)
restore_subproblems!(ph::AbstractProgressiveHedging) = restore_subproblems!(ph, ph.execution)
iterate!(ph::AbstractProgressiveHedging) = iterate!(ph, ph.execution)
start_workers!(ph::AbstractProgressiveHedging) = start_workers!(ph, ph.execution)
close_workers!(ph::AbstractProgressiveHedging) = close_workers!(ph, ph.execution)
resolve_subproblems!(ph::AbstractProgressiveHedging) = resolve_subproblems!(ph, ph.execution)
update_iterate!(ph::AbstractProgressiveHedging) = update_iterate!(ph, ph.execution)
update_subproblems!(ph::AbstractProgressiveHedging) = update_subproblems!(ph, ph.execution)
update_dual_gap!(ph::AbstractProgressiveHedging) = update_dual_gap!(ph, ph.execution)
calculate_objective_value(ph::AbstractProgressiveHedging) = calculate_objective_value(ph, ph.execution)
# ------------------------------------------------------------
include("common.jl")
include("serial.jl")
include("channels.jl")
include("distributed.jl")
include("synchronous.jl")
include("asynchronous.jl")


text:
# Triangle neighbors
function neighbors(grid::Grid{Tria,Periodic}, cell::Cell{Tria})
    return [north(grid, cell), south(grid, cell), east(grid, cell), west(grid, cell)]
end

function neighbors(grid::Grid{Tria,NonPeriodic}, cell::Cell{Tria})
    nrows, ncols = grid.nrows, grid.ncols
    row, col = cell.row, cell.col
    upright = iseven(row + col)
    nbs = []
    !(row == 1) && !upright && push!(nbs, north(grid, cell))
    !(row == nrows) && upright && push!(nbs, south(grid, cell))
    !(col == 1) && push!(nbs, west(grid, cell))
    !(col == ncols) && push!(nbs, east(grid, cell))
    return nbs
end

# Rectangle neighbors
function neighbors(grid::Grid{Rect,Periodic}, cell::Cell{Rect})
    return [north(grid, cell), south(grid, cell), east(grid, cell), west(grid, cell)]
end

function neighbors(grid::Grid{Rect,NonPeriodic}, cell::Cell{Rect})
    nrows, ncols = grid.nrows, grid.ncols
    row, col = cell.row, cell.col
    nbs = []
    !(row == 1) && push!(nbs, north(grid, cell))
    !(row == nrows) && push!(nbs, south(grid, cell))
    !(col == 1) && push!(nbs, west(grid, cell))
    !(col == ncols) && push!(nbs, east(grid, cell))
    return nbs
end

# Hexagon neighbors
function neighbors(grid::Grid{Hexa,Periodic}, cell::Cell{Hexa})
    return [
        north(grid, cell),
        south(grid, cell),
        northwest(grid, cell),
        northeast(grid, cell),
        southwest(grid, cell),
        southeast(grid, cell),
    ]
end

function neighbors(grid::Grid{Hexa,NonPeriodic}, cell::Cell{Hexa})
    nrows, ncols = grid.nrows, grid.ncols
    row, col = cell.row, cell.col
    nbs = []
    !(row == 1) && push!(nbs, north(grid, cell))
    !(row == nrows) && push!(nbs, south(grid, cell))
    !(col == ncols) && !(row == nrows) && push!(nbs, southeast(grid, cell))
    !(col == 1) && !(row == nrows) && push!(nbs, southwest(grid, cell))
    !(col == ncols) && !(row == 1) && push!(nbs, northeast(grid, cell))
    !(col == 1) && !(row == 1) && push!(nbs, northwest(grid, cell))
    return nbs
end


text:
adapters = sort(parse.(Int, readlines("a10.txt")))

diffs = adapters[2:end] .- adapters[1:end-1]

diffcounts = reduce(function (a,v) a[v] += 1; a end, diffs; init = zeros(Int, 3))
diffcounts[adapters[1]] += 1
diffcounts[3] += 1

part1 = diffcounts[1] * diffcounts[3]
println("$part1")

configurations = zeros(Int, adapters[end] + 3)
configurations[end] = 1
for adapter in reverse(adapters)
  configurations[adapter] = sum(configurations[adapter+1:adapter+3])
end

part2 = sum(configurations[1:3])
println("$part2")


text:
using ImageViewer
using ImageViewer.Colors
using Test

@testset "ImageViewer.jl" begin
    x1 = [0.1 0.2]
    x2 = [0.1 0.2; 0.3 0.4]
    x3 = [0.1 0.2 0.3; 0.4 0.5 0.6]
    x4 = [0.1 0.2; 0.3 0.4; 0.5 0.6]
    x5 = [0.1, 0.2]

    @testset "size(x) = $(size(x))" for x in [x1, x2, x3, x4]
        img = Gray.(x);
        img_flipped = Gray.(x');
        @test image(x) == img_flipped
        @test image(x; flip = false) == img
        @test image(x; flip = true) == img_flipped
    end

    @test_throws MethodError image(x5)
end


text:
using Pumas, DataFrames, LinearAlgebra, Plots

model = @model begin
  @param begin
    θ ∈ VectorDomain(4)
    Ω ∈ PSDDomain(3)
    σ_prop ∈ RealDomain(init=0.1)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @covariates isPM Wt

  @pre begin
    TVCL = isPM == 1 ? θ[1] : θ[4]
    CL = θ[1]*(Wt/70)^0.75*exp(η[1])
    V = θ[2]*(Wt/70)^0.75*exp(η[2])
    Ka = θ[3]*exp(η[3])
  end

  @dynamics begin
    Depot'   = -Ka*Depot
    Central' =  Ka*Depot - Central*CL/V
  end

  @vars begin
    conc = Central/V
  end

  @derived begin
    dv ~ @.Normal(conc,sqrt(conc^2*σ_prop+ eps()))
  end

end

fixeffs = (
  θ = [0.4,20,1.1,2],
  Ω = diagm(0 => [0.04,0.04,0.04]),
  σ_prop = 0.04
  )

ev = DosageRegimen(100, time=0)
first(ev.data)

s1 = Subject(id=1,evs=ev,cvs=(isPM=0, Wt=70))
for fn in fieldnames(Subject)
           x = getproperty(s1, fn)
           if !isa(x, Nothing)
               println(fn)
               println(x)
           end
end

s1.events

s1.covariates

obs = simobs(model,s1,fixeffs,obstimes=0:0.1:120)
plot(obs)

s2 = Subject(id=2,evs=ev,cvs=(isPM=1,Wt=70))

twosubjs =  Population([s1,s2])

twosubjs[1]

twosubjs[2]

obs = simobs(model,twosubjs,fixeffs,obstimes=0:0.1:120)

plot(obs)

choose_covariates() = (isPM = rand([1, 0]),
                    Wt = rand(55:80))

cvs = [ choose_covariates() for i in 1:10 ]
DataFrame(cvs)

pop_with_covariates = Population(map(i -> Subject(id=i,evs=ev,cvs=choose_covariates()),1:10))

obs = simobs(model,pop_with_covariates,fixeffs,obstimes=0:0.1:120);

plot(obs)

md =  DosageRegimen(100,ii=24,addl=6)

s3 = Subject(id=3,evs=md, cvs=(isPM=0,Wt=70))

obs = simobs(model, s3, fixeffs,obstimes=0:0.1:240)
plot(obs)

ldmd = DosageRegimen([500,100],cmt=1, time=[0,24], addl=[0,6],ii=[0,24])

s4 = Subject(id=4, evs=ldmd, cvs=(isPM=0,Wt=70))
obs = simobs(model, s4, fixeffs,obstimes=0:0.1:120)
plot(obs, ylims=(0,50))

e1 = DosageRegimen(500,cmt=1, time=0, addl=0,ii=0)
e2 = DosageRegimen(100,cmt=1, time=24, addl=6,ii=24)
evs = DosageRegimen(e1,e2)
obs = simobs(model, s4, fixeffs,obstimes=0:0.1:120)
plot(obs, ylims=(0,50))

e1 = DosageRegimen(100, ii=24, addl=6)
e2 = DosageRegimen(50,  ii=12, addl=13)
e3 = DosageRegimen(200, ii=24, addl=2)

pop1 = Population(map(i -> Subject(id=i,evs=e1,cvs=choose_covariates()),1:5))
pop2 = Population(map(i -> Subject(id=i,evs=e2,cvs=choose_covariates()),6:8))
pop3 = Population(map(i -> Subject(id=i,evs=e3,cvs=choose_covariates()),9:10))
pop = Population(vcat(pop1,pop2,pop3))

obs = simobs(model,pop,fixeffs,obstimes=0:0.1:120)
plot(obs)

inf = DosageRegimen(100, rate=3, cmt=2)

s5 = Subject(id=5, evs=inf, cvs=(isPM=0,Wt=70))
obs = simobs(model, s5, fixeffs, obstimes=0:0.1:120)
plot(obs)

using PumasTutorials
PumasTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module ProfileView

using Compose, Lazy

include("javascript.jl")
include("data.jl")

githuburl(file, line) = "https://github.com/JuliaLang/julia/tree/$(Base.GIT_VERSION_INFO.commit)/base/$file#L$line"

function basepath(file)
  path = joinpath(JULIA_HOME,"..","share","julia","base",file) |> normpath
  isfile(path) || (path = nothing)
  return path
end

maprange(x1, x2, y1, y2, p) = (p-x1)/(x2-x1)*(y2-y1)+y1

fixedscale(node::ProfileTree) = ones(length(node.children))
# widthscale(node::ProfileTree) = childwidths(node)
widthscale(node::ProfileTree) = map(w -> maprange(0, 1, 1/5, 1, w), childwidths(node))

function fileattribute(li)
  if isabspath(li.file)
    svgattribute("data-file", "$(li.file):$(li.line)")
  else
    path = basepath(li.file)
    svgattribute("data-file", "$path:$(li.line)")
  end
end

function render_(tree::ProfileTree; childscale = fixedscale, count = 0)
  widths = childwidths(tree)
  offsets = cumsum([0, widths[1:end-1]...])
  scale = childscale(tree)
  li = tree.data.line
  compose(context(),
          (context(), rectangle(),
           JS.framedata(li, tree.data.count/count),
           JS.frametooltip,
           svgclass("file-link"),
           fileattribute(li)),
          [compose(context(offsets[i], 1, widths[i], scale[i]),
                   render_(tree.children[i], childscale=childscale, count=count))
           for i = 1:length(tree.children)]...)
end

maxheight(node::ProfileTree; childscale = fixedscale) =
  isleaf(node) ? 1 :
    1 + maximum(childscale(node) .*
                  map(node->maxheight(node, childscale=childscale),
                      node.children))

render(tree::ProfileTree; childscale = widthscale) =
  compose(context(),
          (context(0,0,1,1/maxheight(tree, childscale = childscale)),
           render_(tree,childscale = childscale, count = tree.data.count),
           svgclass("tree")),
          (context(), rectangle(), svgclass("background")),
          JS.mapzoom, JS.mapdrag, JS.nonscalingstroke, JS.tooltip, JS.settooltip,
          svgclass("root"))

function Base.writemime(io::IO, ::MIME"text/html", tree::ProfileTree)
  write(io, """

  """)
  draw(SVGJS(io, 5inch, 3inch, false), render(tree))
  write(io, """

  """)
end

end