text:
##################
# AbstractConfig #
##################

abstract type AbstractConfig end

Base.show(io::IO, cfg::AbstractConfig) = print(io, typeof(cfg).name)

##################
# GradientConfig #
##################

struct GradientConfig{I} <: AbstractConfig
    input::I
    tape::InstructionTape
    # disable default outer constructor
    GradientConfig{I}(input, tape) where {I} = new{I}(input, tape)
end

# "private" convienence constructor
_GradientConfig(input::I, tape::InstructionTape) where {I} = GradientConfig{I}(input, tape)

"""
    ReverseDiff.GradientConfig(input, tp::InstructionTape = InstructionTape())

Return a `GradientConfig` instance containing the preallocated tape and work buffers used
by the `ReverseDiff.gradient`/`ReverseDiff.gradient!` methods.

Note that `input` is only used for type and shape information; it is not stored or modified
in any way. It is assumed that the element type of `input` is same as the element type of
the target function's output.

See `ReverseDiff.gradient` for a description of acceptable types for `input`.
"""
GradientConfig(input::AbstractArray{T}, tp::InstructionTape = InstructionTape()) where {T} = GradientConfig(input, T, tp)

GradientConfig(input::Tuple, tp::InstructionTape = InstructionTape()) = GradientConfig(input, eltype(first(input)), tp)

"""
    ReverseDiff.GradientConfig(input, ::Type{D}, tp::InstructionTape = InstructionTape())

Like `GradientConfig(input, tp)`, except the provided type `D` is assumed to be the element
type of the target function's output.
"""
function GradientConfig(input::Tuple, ::Type{D}, tp::InstructionTape = InstructionTape()) where D
    return _GradientConfig(map(x -> track(similar(x), D, tp), input), tp)
end

function GradientConfig(input::AbstractArray, ::Type{D}, tp::InstructionTape = InstructionTape()) where D
    return _GradientConfig(track(similar(input), D, tp), tp)
end

##################
# JacobianConfig #
##################

struct JacobianConfig{I,O} <: AbstractConfig
    input::I
    output::O
    tape::InstructionTape
    # disable default outer constructor
    JacobianConfig{I,O}(input, output, tape) where {I,O} = new{I,O}(input, output, tape)
end

# "private" convienence constructor
_JacobianConfig(input::I, output::O, tape::InstructionTape) where {I,O} = JacobianConfig{I,O}(input, output, tape)

"""
    ReverseDiff.JacobianConfig(input, tp::InstructionTape = InstructionTape())

Return a `JacobianConfig` instance containing the preallocated tape and work buffers used
by the `ReverseDiff.jacobian`/`ReverseDiff.jacobian!` methods.

Note that `input` is only used for type and shape information; it is not stored or modified
in any way. It is assumed that the element type of `input` is same as the element type of
the target function's output.

See `ReverseDiff.jacobian` for a description of acceptable types for `input`.

    ReverseDiff.JacobianConfig(input, ::Type{D}, tp::InstructionTape = InstructionTape())

Like `JacobianConfig(input, tp)`, except the provided type `D` is assumed to be the element
type of the target function's output.
"""
function JacobianConfig(args...)
    gcfg = GradientConfig(args...)
    return _JacobianConfig(gcfg.input, nothing, gcfg.tape)
end

"""
    ReverseDiff.JacobianConfig(output::AbstractArray, input, tp::InstructionTape = InstructionTape())

Return a `JacobianConfig` instance containing the preallocated tape and work buffers used
by the `ReverseDiff.jacobian`/`ReverseDiff.jacobian!` methods. This method assumes the
target function has the form `f!(output, input)`

Note that `input` and `output` are only used for type and shape information; they are not
stored or modified in any way.

See `ReverseDiff.jacobian` for a description of acceptable types for `input`.
"""
function JacobianConfig(output::AbstractArray{D}, input::Tuple, tp::InstructionTape = InstructionTape()) where D
    cfg_input = map(x -> track(similar(x), D, tp), input)
    cfg_output = track!(similar(output, TrackedReal{D,D,Nothing}), output, tp)
    return _JacobianConfig(cfg_input, cfg_output, tp)
end

# we dispatch on V<:Real here because InstructionTape is actually also an AbstractArray
function JacobianConfig(output::AbstractArray{D}, input::AbstractArray{V}, tp::InstructionTape = InstructionTape()) where {D,V<:Real}
    cfg_input = track(similar(input), D, tp)
    cfg_output = track!(similar(output, TrackedReal{D,D,Nothing}), output, tp)
    return _JacobianConfig(cfg_input, cfg_output, tp)
end

"""
    ReverseDiff.JacobianConfig(result::DiffResults.DiffResult, input, tp::InstructionTape = InstructionTape())

A convenience method for `JacobianConfig(DiffResults.value(result), input, tp)`.
"""
JacobianConfig(result::DiffResult, input, tp::InstructionTape) = JacobianConfig(DiffResults.value(result), input, tp)

#################
# HessianConfig #
#################

struct HessianConfig{G<:GradientConfig,J<:JacobianConfig} <: AbstractConfig
    gradient_config::G
    jacobian_config::J
end

"""
    ReverseDiff.HessianConfig(input::AbstractArray, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape())

Return a `HessianConfig` instance containing the preallocated tape and work buffers used
by the `ReverseDiff.hessian`/`ReverseDiff.hessian!` methods. `gtp` is the tape used for
the inner gradient calculation, while `jtp` is used for outer Jacobian calculation.

Note that `input` is only used for type and shape information; it is not stored or modified
in any way. It is assumed that the element type of `input` is same as the element type of
the target function's output.
"""
function HessianConfig(input::AbstractArray, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape())
    return HessianConfig(input, eltype(input), gtp, jtp)
end

"""
    ReverseDiff.HessianConfig(input::AbstractArray, ::Type{D}, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape())

Like `HessianConfig(input, tp)`, except the provided type `D` is assumed to be the element
type of the target function's output.
"""
function HessianConfig(input::AbstractArray, ::Type{D}, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape()) where D
    jcfg = JacobianConfig(input, D, jtp)
    gcfg = GradientConfig(jcfg.input, gtp)
    return HessianConfig(gcfg, jcfg)
end

"""
    ReverseDiff.HessianConfig(result::DiffResults.DiffResult, input::AbstractArray, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape())

Like `HessianConfig(input, tp)`, but utilize `result` along with `input` to construct work
buffers.

Note that `result` and `input` are only used for type and shape information; they are not
stored or modified in any way.
"""
function HessianConfig(result::DiffResult, input::AbstractArray, gtp::InstructionTape = InstructionTape(), jtp::InstructionTape = InstructionTape())
    jcfg = JacobianConfig(DiffResults.gradient(result), input, jtp)
    gcfg = GradientConfig(jcfg.input, gtp)
    return HessianConfig(gcfg, jcfg)
end


text:
using Roots
using PrettyTables

fock_matrix_element(f, S, u, v, tmp) = dot(u, S, mul!(tmp, f, v))

krylov_wrapper(kw::KrylovWrapper) = kw
krylov_wrapper(okw::OrthogonalKrylovWrapper) = okw.A

function analyze_symmetry_orbitals(fock, P::AbstractVecOrMat{T},
                                   kws, ϵ=1e-3; verbosity=0) where T
    S = fock.S

    m = length(fock.equations)
    data = zeros(T, sum(length(s)*(length(s)+1)÷2 for s in fock.symmetries), 3)
    labels = Vector{String}()
    tmp = zeros(T, size(P,1))
    ii = 1
    should_rotate = Vector{Tuple{Int,Int}}()
    for sym ∈ fock.symmetries
        for i ∈ sym
            vPi = view(P, :, i)
            A = krylov_wrapper(kws[i])
            for j ∈ sym
                j < i && continue
                vPj = view(P, :, j)
                data[ii,1] = dot(vPi, S, vPj)
                afb = data[ii,2] = fock_matrix_element(A, S, vPi, vPj, tmp)
                data[ii,3] = fock_matrix_element(A, S, vPj, vPi, tmp)
                push!(labels, "$i – $j")
                ii += 1
                i ≠ j && abs(afb) > ϵ && push!(should_rotate, (i,j))
            end
        end
    end
    if verbosity > 0
        pretty_table(hcat(labels, data), header=["i – j", "⟨i|j⟩", "⟨i|𝔣|j⟩", "⟨j|𝔣|i⟩"],
                     hlines=[1], vlines=[1])
    end

    should_rotate
end

function rotate!(u::AbstractVector, v::AbstractVector, η, tmp::AbstractVector)
    a = inv(√(1+η^2))
    b = η*a
    copyto!(tmp, u)
    BLAS.axpy!(b, v, tmp)
    BLAS.axpy!(-b, u, v)
    copyto!(u, tmp)
    u,v
end

function rotate!(P::AbstractVecOrMat, fock::Fock, A,
                 i::Integer, j::Integer, ϵ; verbosity=0)
    S = fock.S

    a = view(P,:,i)
    b = view(P,:,j)
    tmp = similar(a)

    fba = fock_matrix_element(A, S, b, a, tmp)
    faa = fock_matrix_element(A, S, a, a, tmp)
    fbb = fock_matrix_element(A, S, b, b, tmp)

    g = (faa-fbb)
    g′ = fba
    verbosity > 0 && @show i,j,g,g′,ϵ
    abs(g′) < ϵ && return
    η̃ = g/2g′
    # TODO: Derive this properly
    η̂ = -(η̃ + √(η̃^2 + 1))
    abs(η̂) > 2 && return
    verbosity > 0 && @show i,j,η̂

    rotate!(a, b, η̂, tmp)
end

rotate!(P::AbstractVecOrMat, fock::Fock, okws::Vector,
        i::Integer, j::Integer, ϵ;
        kwargs...) =
            rotate!(P, fock, krylov_wrapper(okws[i]),
                    i, j, ϵ; kwargs...)

function rayleigh_ritz!(P, fock, kws; verbosity=0)
    verbosity > 0 && @info "Rayleigh–Ritz rotation"
    T = eltype(P)
    S = fock.S
    tmp = zeros(T, size(P,1))
    for sym ∈ fock.symmetries
        ns = length(sym)
        if verbosity > 2
            println(repeat("-", 100))
            println("Symmetry: ", sym)
        end
        f = zeros(T, ns, ns)
        Ps = P[:,sym]
        for (ii,i) ∈ enumerate(sym)
            vPi = view(Ps, :, ii)
            A = krylov_wrapper(kws[i])
            for (jj,j) ∈ enumerate(sym)
                j < i && continue
                vPj = view(Ps, :, jj)
                f[ii,jj] = fock_matrix_element(A, S, vPi, vPj, tmp)
                f[jj,ii] = fock_matrix_element(A, S, vPj, vPi, tmp)
            end
        end
        ee = eigen(f)
        mul!(view(P,:,sym)', ee.vectors', Ps')
        if verbosity > 2
            println("Fock submatrix:")
            display(f)
            println("Energies:")
            display(ee.values')
            display(ee.vectors)
            println()
        end
    end
    update!(fock)
    true
end

function rotate_orbitals!(P, fock, kws;
                          verbosity=0,
                          rotate_orbitals=true, rotϵ=1e-3,
                          rotation_method=:pairwise,
                          kwargs...)
    did_rotate = if rotate_orbitals
        if rotation_method == :pairwise
            should_rotate=analyze_symmetry_orbitals(fock, P, kws, rotϵ,
                                                    verbosity=verbosity)
            for (i,j) in should_rotate
                rotate!(P, fock, kws, i, j, rotϵ, verbosity=verbosity)
            end
            !isempty(should_rotate)
        elseif rotation_method == :eigen
            rayleigh_ritz!(P, fock, kws; verbosity=verbosity)
        else
            throw(ArgumentError("Unknown orbital rotation method $(rotation_method)"))
        end
    else
        false
    end
    did_rotate && verbosity > 2 &&
        analyze_symmetry_orbitals(fock, P, kws, verbosity=verbosity)
    did_rotate
end

ext: jl


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      PositionState, Product1pBasis
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf, test_fio,
                       println_slim

##
@info("Testing B1pMultiplier")

module M
    import ACE
    struct TestMult{TF} <: ACE.B1pMultiplier{Float64}
        f::TF
    end

    ACE._inner_evaluate(mult::TestMult, X) = mult.f(X)
end

Bsel = SimpleSparseBasis(3, 5)

RnYlm = ACE.Utils.RnYlm_1pbasis()
ACE.init1pspec!(RnYlm, Bsel)

@info("some basic tests")

const mult_val = 1.234
mult1 = M.TestMult(X -> mult_val)

B1p = mult1 * RnYlm
ACE.init1pspec!(B1p, Bsel)

println_slim(@test length(B1p) == length(RnYlm))

nX = 10
Xs = rand(PositionState{Float64}, RnYlm.bases[1], nX)
cfg = ACEConfig(Xs)

A1 = evaluate(RnYlm, cfg)
A2 = evaluate(B1p, cfg)
println_slim(@test( A1 * mult_val ≈ A2 ))

##

@info("test against manual summation")

_f = X -> exp(- ACE.sumsq(X.rr .- mult_val))
mult2 = M.TestMult(_f)
B1p2 = mult2 * RnYlm
ACE.init1pspec!(B1p2, Bsel)

A1 = sum( evaluate(RnYlm, X) * _f(X) for X in Xs )
A2 = evaluate(B1p2, cfg)
println_slim(@test A1 ≈ A2)


text:
include("../AF.jl")

module AFTest

using AF
using Base.Test

function testOnAllBackends(f, title)
    for b in getSupportedBackends()
        println("DEGIN: $title on $b")
        af = ArrayFire{b}()
        scope!(af) do this
            f(af)
        end
        println("PASSED")
    end
end

include("device.jl")
include("array.jl")
include("binary.jl")
include("unary.jl")
include("create.jl")
include("vectorAlgos.jl")
include("modify.jl")
include("reduction.jl")

end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
import Base: ==, ⋅, ×, +, -, *, /, norm, normalize, zero, ∩
import Images: colorview, RGBA
import ImageView: imshow

struct Vector3
    x::Float64
    y::Float64
    z::Float64
end
Vector3() = Vector3(0, 0, 0)

==(𝐚::Vector3, 𝐛::Vector3) = 𝐚.x == 𝐛.x && 𝐚.y == 𝐛.y && 𝐚.z == 𝐛.z

⋅(𝐚::Vector3, 𝐛::Vector3) = 𝐚.x * 𝐛.x + 𝐚.y * 𝐛.y + 𝐚.z * 𝐛.z
×(𝐚::Vector3, 𝐛::Vector3) = Vector3(
    𝐚.y * 𝐛.z - 𝐚.z * 𝐛.y,
    𝐚.z * 𝐛.x - 𝐚.x * 𝐛.z,
    𝐚.x * 𝐛.y - 𝐚.y * 𝐛.x
)
+(𝐚::Vector3, 𝐛::Vector3) = Vector3(𝐚.x + 𝐛.x, 𝐚.y + 𝐛.y, 𝐚.z + 𝐛.z)
-(𝐚::Vector3, 𝐛::Vector3) = Vector3(𝐚.x - 𝐛.x, 𝐚.y - 𝐛.y, 𝐚.z - 𝐛.z)
*(𝐚::Vector3, μ::Real) = Vector3(𝐚.x * μ, 𝐚.y * μ, 𝐚.z * μ)
/(𝐚::Vector3, μ::Real) = Vector3(𝐚.x / μ, 𝐚.y / μ, 𝐚.z / μ)
*(μ::Real, 𝐚::Vector3) = 𝐚 * μ
-(𝐯::Vector3) = Vector3(-𝐯.x, -𝐯.y, -𝐯.z)

norm²(𝐯::Vector3) = 𝐯 ⋅ 𝐯
norm(𝐯::Vector3) = √norm²(𝐯)
normalize(𝐯::Vector3) = 𝐯 / norm(𝐯)
zero(::Type{Vector3}) = Vector3(0, 0, 0)

∥(𝐚::Vector3, 𝐛::Vector3) = 𝐚 × 𝐛 == zero(Vector3)
⟂(𝐚::Vector3, 𝐛::Vector3) = 𝐚 ⋅ 𝐛 == 0

struct Ray3
    origin::Vector3
    direction::Vector3
end
get_point(ray::Ray3, t::Real) = ray.origin + ray.direction * t

const Canvas = Array{Float64, 3}
Canvas(width, height) = zeros(4, width, height)
width(c::Canvas) = size(c)[2]
height(c::Canvas) = size(c)[3]

imshow(c::Canvas, args...) = imshow(colorview(RGBA, c), args...)

abstract type Color end

struct ColorRGB <: Color
    r::Float64
    g::Float64
    b::Float64
end
const Black = ColorRGB(0, 0, 0)
const White = ColorRGB(1, 1, 1)
const Red = ColorRGB(1, 0, 0)
const Green = ColorRGB(0, 1, 0)
const Blue = ColorRGB(0, 0, 1)

+(a::ColorRGB, b::ColorRGB) = ColorRGB(a.r + b.r, a.g + b.g, a.b + b.b)
*(c::ColorRGB, μ::Real) = ColorRGB(c.r * μ, c.g * μ, c.b * μ)
*(μ::Real, c::ColorRGB) = c * μ
/(c::ColorRGB, μ::Real) = c * (1 / μ)

modulate(a::ColorRGB, b::ColorRGB) = ColorRGB(a.r * b.r, a.g * b.g, a.b * b.b)


text:
export NFFTOp
import Base.adjoint
import LinearOperators.FuncOrNothing

mutable struct NFFTOp{T,F1<:FuncOrNothing,F2<:FuncOrNothing} <:
                      AbstractLinearOperator{T,Function,F1,F2}
  nrow :: Int
  ncol :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod :: Function
  tprod :: F1
  ctprod :: F2
end

"""
    NFFTOp(shape::Tuple, tr::Trajectory; nodes=nothing, kargs...)

generates a `NFFTOp` which evaluates the MRI Fourier signal encoding operator using the NFFT.

# Arguments:
* `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
* `tr::Trajectory`          - Trajectory with the kspace nodes to sample
* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
* (`kargs`)                   - additional keyword arguments
"""
function NFFTOp(shape::Tuple, tr::Trajectory; nodes=nothing, kargs...)
  nodes==nothing ? nodes=kspaceNodes(tr) : nothing
  plan = NFFTPlan(nodes, shape, 3, 1.25)

  function produ(x::Vector{T}) where T<:Union{Real,Complex}
    y = nfft(plan,reshape(x[:],shape))
    return vec(y)
  end

  function ctprodu(y::Vector{T}) where T<:Union{Real,Complex}
    x = nfft_adjoint(plan, y[:])
    return vec(x)
  end

  return NFFTOp{ComplexF64,Nothing,Function}(size(nodes,2), prod(shape), false, false
            , produ
            , nothing
            , ctprodu)
end

function adjoint(op::NFFTOp{T}) where T
  return LinearOperator{T,Function,Nothing,Function}(op.ncol, op.nrow, op.symmetric, op.hermitian,
                        op.ctprod, nothing, op.prod)
end


text:
###################################################################################################
# COMMON
###################################################################################################
"""
    `get_all_∇logZs(n,θs)`

    Function that returns gradient vectors for all parameter vectors θ passed.

    # Examples
    ```julia-repl
    julia> n=[5]; θs=[[0.0,0.0],[0.0,0.0]];
    julia> CpelTdm.get_all_∇logZs(n,θs)
    2-element Array{Array{Float64,1},1}:
     [0.0, 0.0]
     [0.0, 0.0]
    ```
"""
function get_all_∇logZs(n::Vector{Int64},θs::Vector{Vector{Float64}})::Vector{Vector{Float64}}

    # Generate all gradients
    ∇logZs = []
    @inbounds for s=1:length(θs)
        push!(∇logZs,get_∇logZ(n,θs[s]))
    end

    # Return vector of gradients
    return ∇logZs

end
###################################################################################################
# UNMATCHED TESTS
###################################################################################################
"""
    `comp_unmat_stat_mml(n,∇logZ1s,∇logZ2s)`

    Function that computes MML difference between groups (unmatched case).

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[0.0,0.0],[0.0,0.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_stat_mml(n,∇logZ1s,∇logZ2s)
    -0.48867352
    ```
"""
function comp_unmat_stat_mml(n::Vector{Int64},∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Float64

    # Compute mean for group 1
    mml1 = 0.0
    @inbounds for ∇logZ in ∇logZ1s
        mml1 += comp_mml(n,∇logZ)
    end
    mml1 /= length(∇logZ1s)

    # Compute mean for group 2
    mml2 = 0.0
    @inbounds for ∇logZ in ∇logZ2s
        mml2 += comp_mml(n,∇logZ)
    end
    mml2 /= length(∇logZ2s)

    # Return Tmml
    return mml1-mml2

end
"""
    `comp_unmat_perm_stat_mml(n,∇logZ1s,∇logZ2s,perm_ids)`

    Function that produces a permutation statistic for Tmml in unmatched case.

    # Examples
    ```julia-repl
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5); perm_ids=[1,3,5,7,9];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_perm_stat_mml(n,∇logZ1s,∇logZ2s,perm_ids)
    0.15564650799999996
    ```
"""
function comp_unmat_perm_stat_mml(n::Vector{Int64},∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}},perm_ids::Vector{Int64})::Float64

    # Get vectors for each group
    ∇logZ2sp = vcat(∇logZ1s,∇logZ2s)
    ∇logZ1sp = ∇logZ2sp[perm_ids]
    deleteat!(∇logZ2sp,perm_ids)

    # Return
    return comp_unmat_stat_mml(n,∇logZ1sp,∇logZ2sp)

end
"""
    `comp_unmat_stat_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)`

    Function that computes NME difference between groups (unmatched case).

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[0.0,0.0],[0.0,0.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_stat_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    0.91986634
    ```
"""
function comp_unmat_stat_nme(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Float64

    # Compute mean for group 1
    nme1 = 0.0
    @inbounds for s=1:length(θ1s)
        nme1 += comp_nme(n,θ1s[s],∇logZ1s[s])
    end
    nme1 /= length(θ1s)

    # Compute mean for group 2
    nme2 = 0.0
    @inbounds for s=1:length(θ2s)
        nme2 += comp_nme(n,θ2s[s],∇logZ2s[s])
    end
    nme2 /= length(θ2s)

    # Return difference in nme
    return nme1-nme2

end
"""
    `comp_unmat_perm_stat_nme(n,θ1s,θ2s,perm_ids)`

    Function that produces a permutation statistic for Tnme in unmatched case.

    # Examples
    ```julia-repl
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5); perm_ids=[1,3,5,7,9];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_perm_stat_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s,perm_ids)
    0.0
    ```
"""
function comp_unmat_perm_stat_nme(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}},perm_ids::Vector{Int64})::Float64

    # Get vectors for each group
    θ2sp = vcat(θ1s,θ2s)
    θ1sp = θ2sp[perm_ids]
    deleteat!(θ2sp,perm_ids)

    # Get vectors for each group
    ∇logZ2sp = vcat(∇logZ1s,∇logZ2s)
    ∇logZ1sp = ∇logZ2sp[perm_ids]
    deleteat!(∇logZ2sp,perm_ids)

    # Return
    return comp_unmat_stat_nme(n,θ1sp,θ2sp,∇logZ1sp,∇logZ2sp)

end
"""
    `comp_unmat_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s)`

    Function that computes test statistic Tgsd in unmatched case.

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[-1.0,1.0],[-1.0,1.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    0.7888652058295635
    ```
"""
function comp_unmat_stat_pdm(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Float64

    ## Compute stat
    cmds = []
    @inbounds for s1=1:length(θ1s)
        @inbounds for s2=1:length(θ2s)
            push!(cmds,comp_cmd(n,θ1s[s1],θ2s[s2],∇logZ1s[s1],∇logZ2s[s2]))
        end
    end

    # Return sum of GJSD
    return sum(cmds)/(length(θ1s)*length(θ2s))

end
"""
    `comp_unmat_perm_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s,perm_ids)`

    Function that produces a permutation statistic for tpdm in unmatched case.

    # Examples
    ```julia-repl
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5); perm_ids=[1,3,5,7,9];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_unmat_perm_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s,perm_ids)
    0.1691594717040136
    ```
"""
function comp_unmat_perm_stat_pdm(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}},perm_ids::Vector{Int64})::Float64

    # Get vectors for each group
    θ2sp = vcat(θ1s,θ2s)
    θ1sp = θ2sp[perm_ids]
    deleteat!(θ2sp,perm_ids)

    # Get vectors for each group
    ∇logZ2sp = vcat(∇logZ1s,∇logZ2s)
    ∇logZ1sp = ∇logZ2sp[perm_ids]
    deleteat!(∇logZ2sp,perm_ids)

    # Return
    return comp_unmat_stat_pdm(n,θ1sp,θ2sp,∇logZ1sp,∇logZ2sp)

end
"""
    `unmat_tests(n,θ1s,θ2s)`

    Function that performs hypothesis testing in unmatched samples group comparison.

    # Examples
    ```julia-repl
    julia> using Random; Random.seed!(1234);
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.unmat_tests(n,θ1s,θ2s)
    ((0.7782325400000001, 0.007936507936507936), (0.0, 1.0), (0.3253066763538721, 0.007936507936507936))
    julia> n=[10]; θ1s=fill([0.5,0.5],8); θ2s=fill([-0.5,0.5],8);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.unmat_tests(n,θ1s,θ2s)
    ((0.77823254, 0.0020833333333333333), (0.0, 1.0), (0.325306676353872, 0.0020833333333333333))
    julia> n=[10]; θ1s=fill([0.0,0.0],8); θ2s=fill([0.0,2.5],8);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.unmat_tests(n,θ1s,θ2s)
    ((0.0, 1.0), (0.84782978, 0.0010309278350515464), (0.43786215630140907, 0.0010309278350515464))
    ```
"""
function unmat_tests(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}};Lmax::Int64=1000)::NTuple{3,NTuple{2,Float64}}

    # Compute gradients
    ∇logZ1s = get_all_∇logZs(n,θ1s)
    ∇logZ2s = get_all_∇logZs(n,θ2s)

    # Compute observed stats
    tmml_obs = comp_unmat_stat_mml(n,∇logZ1s,∇logZ2s)
    tnme_obs = comp_unmat_stat_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    tpdm_obs = comp_unmat_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s)

    # Compute number of possible randomizations
    L = binomial(length(θ1s)+length(θ2s),length(θ1s))
    exact = L<Lmax

    # Create iteratable object with all combinations
    comb_iter = combinations(1:(length(θ1s)+length(θ2s)),length(θ1s))

    # Get group label combinations to use
    comb_iter_used = []
    if exact
        # Use all group assignments
        comb_iter_used = comb_iter
    else
        # Use Lmax group assignments
        ind_subset = rand(1:L,Lmax)
        @inbounds for (ind,comb) in enumerate(comb_iter)
            (ind in ind_subset) && push!(comb_iter_used,comb)
        end
    end

    # Use method for random permutation
    tmml_perms = map(x->comp_unmat_perm_stat_mml(n,∇logZ1s,∇logZ2s,x),comb_iter_used)
    tnme_perms = map(x->comp_unmat_perm_stat_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s,x),comb_iter_used)
    tpdm_perms = map(x->comp_unmat_perm_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s,x),comb_iter_used)

    # Compute p-values two-sided test
    if exact
        # Compute exact p-value
        tmml_pval = sum(abs.(tmml_perms).>=abs(tmml_obs))/length(tmml_perms)
        tnme_pval = sum(abs.(tnme_perms).>=abs(tnme_obs))/length(tnme_perms)
        tpdm_pval = sum(tpdm_perms.>=tpdm_obs)/length(tpdm_perms)
    else
        # Compute p-value in MC setting
        tmml_pval = (1.0+sum(abs.(tmml_perms).>=abs(tmml_obs)))/(1.0+length(tmml_perms))
        tnme_pval = (1.0+sum(abs.(tnme_perms).>=abs(tnme_obs)))/(1.0+length(tnme_perms))
        tpdm_pval = (1.0+sum(tpdm_perms.>=tpdm_obs))/(1.0+length(tpdm_perms))
    end

    # Return stat-pval pairs
    return (tmml_obs,tmml_pval),(tnme_obs,tnme_pval),(tpdm_obs,tpdm_pval)

end
###################################################################################################
# MATCHED TESTS
###################################################################################################
"""
    `comp_mat_diff_mml(n,∇logZ1s,∇logZ2s)`

    Function that computes MML differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[0.0,0.0],[0.0,0.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_mat_diff_mml(n,∇logZ1s,∇logZ2s)
    2-element Array{Float64,1}:
     -0.48867352
     -0.48867352
    ```
"""
function comp_mat_diff_mml(n::Vector{Int64},∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Vector{Float64}

    # Compute mean differences
    diffs = []
    @inbounds for s=1:length(∇logZ1s)
        push!(diffs,comp_mml(n,∇logZ1s[s])-comp_mml(n,∇logZ2s[s]))
    end

    # Return vector of differences
    return diffs

end
"""
    `comp_mat_diff_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)`

    Function that computes NME differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[0.0,0.0],[0.0,0.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_mat_diff_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    2-element Array{Float64,1}:
     0.91986634
     0.91986634
    ```
"""
function comp_mat_diff_nme(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Vector{Float64}

    # Compute entropy differences
    diffs = []
    @inbounds for s=1:length(θ1s)
        push!(diffs,comp_nme(n,θ1s[s],∇logZ1s[s])-comp_nme(n,θ2s[s],∇logZ2s[s]))
    end

    # Return vector of differences
    return diffs

end
"""
    `comp_mat_j_stat(diffs,j)`

    Function that computes permutation statistic for j-th sign assigment given vector of differences
    between matched pairs (matched case).

    # Examples
    ```julia-repl
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5);
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> mml_diffs = CpelTdm.comp_mat_diff_mml(n,∇logZ1s,∇logZ2s);
    julia> CpelTdm.comp_mat_j_stat(mml_diffs,0)
    -0.7782325400000001
    ```
"""
function comp_mat_j_stat(diffs::Vector{Float64},j::Int64)::Float64

    # Init diff in permutation
    diffs_perm = 0.0

    # Change sign if pertinent
    @inbounds for i=1:length(diffs)
        diffs_perm += Bool(j & 1) ? diffs[i] : -diffs[i]
        j >>= 1
    end

    # Return permutation statistic
    return diffs_perm/length(diffs)

end
"""
    `comp_mat_perm_stats(n,θ1s,θ2s)`

    Function that computes permutation statistics given a vector of differences between pairs (matched case).

    # Examples
    ```julia-repl
    julia> n=[10]; θ1s=fill([0.5,0.5],2); θ2s=fill([-0.5,0.5],2);
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> mml_diffs = CpelTdm.comp_mat_diff_mml(n,∇logZ1s,∇logZ2s);
    julia> js = collect(0:2^length(θ1s)-1);
    julia> CpelTdm.comp_mat_perm_stats(mml_diffs,js)
    4-element Array{Float64,1}:
     -0.7782325400000001
     0.0
     0.0
     0.7782325400000001
    ```
"""
function comp_mat_perm_stats(diffs::Vector{Float64},js::Vector{Int64})::Vector{Float64}

    # Return all possible signed sums
    return [comp_mat_j_stat(diffs,j) for j in js]

end
"""
    `comp_mat_stat_pdm(n,θ1s,θ2s)`

    Function that computes GJSD between pairs (matched case).

    # Examples
    ```julia-repl
    julia> n=[5]; θ1s=[[-1.0,1.0],[-1.0,1.0]]; θ2s=[[1.0,1.0],[1.0,1.0]];
    julia> ∇logZ1s = [CpelTdm.get_∇logZ(n,θ1) for θ1 in θ1s];
    julia> ∇logZ2s = [CpelTdm.get_∇logZ(n,θ2) for θ2 in θ2s];
    julia> CpelTdm.comp_mat_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    0.7888652058295635
    ```
"""
function comp_mat_stat_pdm(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}},
    ∇logZ1s::Vector{Vector{Float64}},∇logZ2s::Vector{Vector{Float64}})::Float64

    # Compute stat
    cmds = 0.0
    @inbounds for s=1:length(θ1s)
        cmds += comp_cmd(n,θ1s[s],θ2s[s],∇logZ1s[s],∇logZ2s[s])
    end

    # Return mean GJSD
    return cmds/length(θ1s)

end
"""
    `mat_tests(n,θ1s,θ2s)`

    Function that performs hypothesis testing in matched samples group comparison. Note that the matched test
    is only performed for Tmml and Tnme.

    # Examples
    ```julia-repl
    julia> using Random; Random.seed!(1234);
    julia> n=[10]; θ1s=fill([0.5,0.5],5); θ2s=fill([-0.5,0.5],5);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.mat_tests(n,θ1s,θ2s)
    ((0.7782325400000001, 0.0625), (0.0, 1.0), (0.3253066763538722, 0.03125))
    julia> n=[10]; θ1s=fill([0.5,0.5],8); θ2s=fill([-0.5,0.5],8);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.mat_tests(n,θ1s,θ2s)
    ((0.7782325400000002, 0.00390625), (0.0, 1.0), (0.3253066763538723, 1.0))
    julia> n=[4]; θ1s=fill([0.5,0.5],8); θ2s=fill([0.0,0.0],8);
    julia> tmml_test,tnme_test,tpdm_test = CpelTdm.mat_tests(n,θ1s,θ2s)
    ((0.7782325400000002, 0.00390625), (0.0, 1.0), (0.3253066763538723, 1.0))
    ```
"""
function mat_tests(n::Vector{Int64},θ1s::Vector{Vector{Float64}},θ2s::Vector{Vector{Float64}};Lmax::Int64=1000)::NTuple{3,NTuple{2,Float64}}

    # Compute gradients
    ∇logZ1s = get_all_∇logZs(n,θ1s)
    ∇logZ2s = get_all_∇logZs(n,θ2s)

    # Compute observed stats
    mml_diffs = comp_mat_diff_mml(n,∇logZ1s,∇logZ2s)
    nme_diffs = comp_mat_diff_nme(n,θ1s,θ2s,∇logZ1s,∇logZ2s)
    tpdm_obs = comp_mat_stat_pdm(n,θ1s,θ2s,∇logZ1s,∇logZ2s)

    # Get group label combinations to use
    exact = 2^length(θ1s)<Lmax
    js = exact ? collect(0:2^length(θ1s)-1) : rand(0:2^length(θ1s)-1,Lmax)

    # Compute permutation stats
    tmml_perms = comp_mat_perm_stats(mml_diffs,js)
    tnme_perms = comp_mat_perm_stats(nme_diffs,js)

    # Compute observed stats
    tmml_obs = sum(mml_diffs)/length(θ1s)
    tnme_obs = sum(nme_diffs)/length(θ1s)

    # Compute p-values
    if exact
        # Compute exact p-value
        tmml_pval = sum(abs.(tmml_perms).>=abs(tmml_obs))/length(tmml_perms)
        tnme_pval = sum(abs.(tnme_perms).>=abs(tnme_obs))/length(tnme_perms)
    else
        # Compute p-value in MC setting
        tmml_pval = (1.0+sum(abs.(tmml_perms).>=abs(tmml_obs)))/(1.0+length(tmml_perms))
        tnme_pval = (1.0+sum(abs.(tnme_perms).>=abs(tnme_obs)))/(1.0+length(tnme_perms))
    end

    # Return (stat,pval) pairs
    return (tmml_obs,tmml_pval),(tnme_obs,tnme_pval),(tpdm_obs,NaN)

end
###################################################################################################
# MULTIPLE HYPOTHESIS CORRECTION
###################################################################################################
"""
    `add_qvalues_to_path(PATH)`

    Function that takes in the differential output and adds BH q-value.

    # Examples
    ```julia-repl
    julia> CpelAsm.add_qvalues_to_path(path)
    ```
"""
function add_qvalues_to_path(path::String)::Nothing

    # Leave if no data
    filesize(path)>0 || return nothing

    # Get data
    all_data = readdlm(path,'\t',Any)
    qvals = fill(NaN,size(all_data)[1])

    # Multiple hypothesis testing correction
    ind = .!isnan.(all_data[:,5])
    if sum(ind)>0
        qvals[ind] = MultipleTesting.adjust(convert(Vector{Float64},all_data[ind,5]),BenjaminiHochberg())
    end

    # Append to output matrix
    all_data = hcat(all_data,qvals)

    # Write to temp output
    temp_path = path * ".tmp"
    open(temp_path,"w") do io
        writedlm(io,all_data,'\t')
    end

    # Move to original file
    mv(temp_path,path,force=true)

    # Return
    return nothing

end
"""
    `mult_hyp_corr(DIFF_PATHS)`

    Function that takes in all the differential output and adds BH q-value in each one.

    # Examples
    ```julia-repl
    julia> CpelAsm.mult_hyp_corr(diff_paths)
    ```
"""
function mult_hyp_corr(diff_paths::Vector{String})

    # Add q-values
    add_qvalues_to_path(diff_paths[1])
    add_qvalues_to_path(diff_paths[2])
    add_qvalues_to_path(diff_paths[3])

    # Return
    return nothing

end


text:
# test 2D PdePumiInterface

using PumiConfig

@testset "--- Testing PdePumiInterface --- " begin

  test_math()

  # test masked copy
  mask = [1, 3]
  a = rand(4)
  a2 = zeros(2)
  PdePumiInterface.copy_masked(a2, a, mask)
  @test ( a2 )== a[mask]

  a = rand(3, 4)
  a2 = rand(3, 2)
  PdePumiInterface.copy_masked(a2, a, mask)
  @test ( a2 )== a[:, mask]

  a = rand(3, 2, 4)
  a2 = rand(3, 2,  2)
  PdePumiInterface.copy_masked(a2, a, mask)
  @test ( a2 )== a[:, :, mask]

  a = rand(3,2, 5, 4)
  a2 = rand(3, 2, 5, 2)
  PdePumiInterface.copy_masked(a2, a, mask)
  @test ( a2 )== a[:, :, :, mask]

  a = rand(2, 3)
  a2 = rand(1, 3)
  @test_throws Exception  PdePumiInterface.copy_masked(a2, a, mask)

  test_refcounting()

  opts = Dict{Any, Any}(
    "numBC" => 1,
    "BC1" =>  [0],
    "run_type" => 4,
    "verify_coloring" => true,
    "use_edge_res" => false,
    "write_edge_vertnums" => true,
    "write_face_vertnums" => true,
    "write_boundarynums" => true,
    "write_dxidx" => true,
    "write_coords" => true,
    "write_sparsity" => true,
    "write_offsets" => true,
    "write_counts" => true,
    "write_sparsity_nodebnds" => true,
    "write_dofs" => false,
    "write_interfaces" => false,
    "write_boundaries" => false,
    "write_sharedboundaries" => false,
    "use_linear_metrics" => true,
    )

    # names of output files
    fnames = ["boundary_nums", "face_vertnums", "edge_vertnums", "sparsity_bnds", "sparsity_bnds", "sparsity_nodebnds"]

    smb_name = "tri2l.smb"
    dmg_name = ".null"
    for order = 1:4
      println("testing order ", order, " CG mesh")
      sbp = getTriSBPGamma(degree=order)
#      sbp = TriSBP{Float64}(degree=order)
      ref_verts = [-1. 1 -1; -1 -1 1]
      sbpface = TriFace{Float64}(order, sbp.cub, ref_verts.')

    mesh =  PumiMesh2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)
    @test ( mesh.numVert )== 4
    @test ( mesh.numEdge )== 5
    @test ( mesh.numFace )== mesh.numEdge
    @test ( mesh.numEl )== 2
    @test ( mesh.numEntitiesPerType )== [4, 5, 2]
    @test ( mesh.numTypePerElement )== [3, 3, 1]
    @test ( mesh.numDofPerNode )== 4
    @test ( mesh.numBoundaryFaces )== 4
    @test ( mesh.numInterfaces )== 1

    @test ( mesh.bndryfaces[1].element )== 1
    @test ( mesh.bndryfaces[1].face )== 3
    @test ( mesh.bndryfaces[2].element )== 2
    @test ( mesh.bndryfaces[2].face )== 1
    @test ( mesh.bndryfaces[3].element )== 1
    @test ( mesh.bndryfaces[3].face )== 2
    @test ( mesh.bndryfaces[4].element )== 2
    @test ( mesh.bndryfaces[4].face )== 2

  #  println("mesh.interfaces = ",  mesh.interfaces)
    @test ( mesh.interfaces[1].elementL )== 1
    @test ( mesh.interfaces[1].elementR )== 2
    @test ( mesh.interfaces[1].faceL )== 1
    @test ( mesh.interfaces[1].faceR )== 3

    @test ( mesh.order )== order
    @test ( length(mesh.bndry_funcs) )== 1
    @test ( mesh.bndry_offsets )== [1, 5]
    @test ( mesh.bndry_geo_nums[1] )== opts["BC1"]

#=
    for i=1:mesh.numBoundaryFaces
      for j=1:(sum(mesh.numNodesPerType[1:2]))
        if i == 1
          @test isapprox( mesh.bndry_normals[:, j, i], [-1.0, 1.0]) atol=1e-13
        elseif i == 2
          @test isapprox( mesh.bndry_normals[:, j, i], [1.0, -1.0]) atol=1e-13
        elseif i == 3
          @test isapprox( mesh.bndry_normals[:, j, i], [0.0, 1.0]) atol=1e-13
        elseif i == 4
          @test isapprox( mesh.bndry_normals[:, j, i], [1.0, 0.0]) atol=1e-13
        end
      end
    end
=#
#=
    for i=1:mesh.numInterfaces
      for j=1:sbp.numfacenodes
        @test isapprox( mesh.interface_normals[:, 1, j, i], [1.0, -2.0]) atol=1e-13
        @test isapprox( mesh.interface_normals[:, 2, j, i], [-2.0, 1.0]) atol=1e-13
      end
    end
=#
    # verify that dofs on a node are numbered consecutively
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        start_dof = mesh.dofs[1, j, i]
        for k=1:mesh.numDofPerNode
          @test ( mesh.dofs[k, j, i] )== start_dof + k - 1
        end
      end
    end

    @test ( mesh.color_masks[1][1] )== 1
    @test ( mesh.color_masks[1][2] )== 0
    @test ( mesh.color_masks[2][1] )== 0
    @test ( mesh.color_masks[2][2] )== 1

    @test ( mesh.neighbor_colors[1,1] )== 2
    @test ( mesh.neighbor_colors[2,1] )== 1
    @test ( mesh.neighbor_colors[1,2] )== 1
    @test ( mesh.neighbor_colors[2,2] )== 2

    @test ( mesh.neighbor_nums[1,1] )== 2
    @test ( mesh.neighbor_nums[2,1] )== 1
    @test ( mesh.neighbor_nums[2,1] )== 1
    @test ( mesh.neighbor_nums[2,2] )== 2

    @test ( mesh.pertNeighborEls[1, 1] )== 1
    @test ( mesh.pertNeighborEls[2, 1] )== 1
    @test ( mesh.pertNeighborEls[1, 2] )== 2
    @test ( mesh.pertNeighborEls[2, 2] )== 2

    @test ( mesh.color_cnt )== [1,1,0, 0]
    if order == 1
      @test ( mesh.numDof )== 16
      @test ( mesh.numNodes )== 4
      @test ( mesh.numNodesPerElement )== 3
      @test ( mesh.numNodesPerType )== [1, 0 , 0]
      @test ( mesh.typeOffsetsPerElement )== [1, 4, 4, 4]
    elseif order == 2
      @test ( mesh.numNodes )== 11
      @test ( mesh.numDof )== 44
      @test ( mesh.numNodesPerElement )== 7
      @test ( mesh.numNodesPerType )== [1, 1, 1]
      @test ( mesh.typeOffsetsPerElement )== [1, 4, 7, 8]
    elseif order == 3
      @test ( mesh.numNodes )== 20
      @test ( mesh.numDof )== 80
      @test ( mesh.numNodesPerElement )== 12
      @test ( mesh.numNodesPerType )== [1, 2, 3]
      @test ( mesh.typeOffsetsPerElement )== [1, 4, 10, 13]
      # check orientation of 3rd edge
#      println("mesh.dofs = ", mesh.dofs)
      @test ( reshape(mesh.dofs[1, 4:5, 1], 2) )== reshape(mesh.dofs[1, 8:9, 2], 2)
    elseif order == 4
      @test ( mesh.numNodes )== 31
      @test ( mesh.numDof )== 124
      @test ( mesh.numNodesPerElement )== 18
      @test ( mesh.numNodesPerType )== [1, 3, 6]
      @test ( mesh.typeOffsetsPerElement )== [1, 4, 13, 19]
#      println("mesh.dofs = ", mesh.dofs)
      # check orientation of 3rd edge
      @test ( reshape(mesh.dofs[1, 10:12, 2], 3) )==reshape(mesh.dofs[1, 4:6, 1], 3)
    end

    @test ( mesh.typeOffsetsPerElement_ )== mesh.typeOffsetsPerElement

    @test ( length(mesh.verts) )== mesh.numVert
    @test ( length(mesh.edges) )== mesh.numEdge
    @test ( length(mesh.elements) )== mesh.numEl

#    println("mesh.coords = ", mesh.coords)

    @test isapprox( mesh.jac, ones(mesh.numNodesPerElement ,2))

    fnames = ["boundary_nums", "face_vertnums", "edge_vertnums"]
    for name in fnames
      name_code = string(name, ".dat")
      name_ref = string(name, "_p", order, "true.dat")
      println("checking file ", name_code)
      println("against ", name_ref)
      data_code = readdlm(name_code)
      data_ref = readdlm(name_ref)

      for i=1:length(data_code)
        data_i = data_code[i]
        if typeof(data_i) <: Number
          @test isapprox( data_i, data_ref[i]) atol=1e-13
        else
          @test ( data_i )== data_ref[i]
        end
      end
    end

  # test vertmap
  println("testing vertmap")

  nedges_interior = 0
  for i=1:mesh.numEl
    vertnums_i = mesh.element_vertnums[:, i]
    for j=1:3
      @test  vertnums_i[j]  > 0
      @test  vertnums_i[j]  < mesh.numVert + 1
    end

    # test there are exactly 3 elements the current element shares 2 vertices
    # with
    nedges = 0
    for j=1:mesh.numEl
      vertnums_j = mesh.element_vertnums[:, j]
      if length(intersect(vertnums_i, vertnums_j)) == 2
        nedges += 1
        nedges_interior += 1
      end
    end

    @test  nedges  < 4
    @test  nedges  > 0
  end  # end loop i

  @test ( nedges_interior )== 2*(mesh.numEdge - mesh.numBoundaryFaces)

  # test getNodeXi
  if order <= 2
    fshape = apf.getFieldShape(0, order, 2)
    eshape = apf.getEntityShape(fshape, 2)  # triangle
    nodexi = PdePumiInterface.getXiCoords(order, 2)
    numnodes = size(nodexi, 2)
    for i=1:numnodes
      vals = apf.getValues(mesh.m_ptr, eshape, nodexi[:, i], numnodes)
      for j=1:numnodes
        if i == j
          @test  abs(vals[j] - 1)  < 1e-12
        else
          @test  abs(vals[j])  < 1e-12
        end  # end if else
      end   # end loop j
    end  # end loop i
  end  # end if p < 2

  # test getBoundary
  geo_bndries = opts["BC1"]
  bndries = getBoundaries(mesh, geo_bndries)
  bndries_range = mesh.bndry_offsets[1]:(mesh.bndry_offsets[2]-1)
  @test ( length(bndries) )== length(bndries_range)
  for i=1:length(bndries)
    bndry1 = bndries[i]
    bndry2 = mesh.bndryfaces[bndries_range[i]]
    @test ( bndry1.face )== bndry2.face
    @test ( bndry1.element )== bndry2.element
  end

  end  # end loop over p=1:4

  # verify adj and det are correct
  for i=1:10  # 10 random matrices
    A2 = rand(2,2)
    d = det(A2)
    @test isapprox( PdePumiInterface.det2(A2), d) atol=1e-12
    B2 = zeros(A2)
    PdePumiInterface.adjugate2(A2, B2)
    @test isapprox( B2./d, inv(A2)) atol=1e-12

    A3 = rand(3,3)
    d = det(A3)
    @test isapprox( PdePumiInterface.det3(A3), d) atol=1e-12
    B3 = zero(A3)
    PdePumiInterface.adjugate3(A3, B3)
    @test isapprox( B3/d, inv(A3)) atol=1e-12
  end

  # now test on a fully unstructured mesh
  smb_name = "vortex.smb"
  dmg_name = "vortex.dmg"

  opts = Dict{Any, Any}(
    "numBC" => 2,
    "BC1" =>  [4, 10],
    "BC2" =>  [7, 13],
    "run_type" => 4,
    "verify_coloring" => true,
    "use_edge_res" => false,
    "write_edge_vertnums" => true,
    "write_face_vertnums" => true,
    "write_boundarynums" => true,
    "write_dxidx" => true,
    "write_coords" => true,
    "write_sparsity" => true,
    "write_offsets" => true,
    "write_counts" => true,
    "write_sparsity_nodebnds" => true,
    "write_dofs" => false,
    "write_interfaces" => false,
    "write_boundaries" => false,
    "write_sharedboundaries" => false,
    "use_linear_metrics" => true,
    )

  for order = 1:4
    println("testing order ", order, " CG mesh against files")
    sbp = getTriSBPGamma(degree=order)
#    sbp = TriSBP{Float64}(degree=order)
    ref_verts = [-1. 1 -1; -1 -1 1]
    sbpface = TriFace{Float64}(order, sbp.cub, ref_verts.')

    mesh =  PumiMesh2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

    for name in fnames
      name_code = string(name, ".dat")
      name_ref = string(name, "vortex", "_p", order, "true.dat")
      println("checking file ", name_code)
      println("against reference file ", name_ref)
      data_code = readdlm(name_code)
      data_ref = readdlm(name_ref)

      for i=1:length(data_code)
        data_i = data_code[i]
        if typeof(data_i) <: Number
          @test isapprox( data_i, data_ref[i]) atol=1e-13
        else
          @test ( data_i )== data_ref[i]
        end
      end
    end

  end   # end loop order=1:4

end

@testset "----- Testing PdePumiInterfaceDG -----" begin

  order = 1
  smb_name = "tri2l.smb"
  dmg_name = ".null"

  opts = Dict{Any, Any}(
    "dmg_name" => dmg_name,
    "smb_name" => smb_name,
    "order" => order,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" =>  [0],
    "run_type" => 4,
    "verify_coloring" => true,
    "use_edge_res" => false,
    "write_edge_vertnums" => true,
    "write_face_vertnums" => true,
    "write_boundarynums" => true,
    "write_dxidx" => true,
    "write_coords" => true,
    "write_sparsity" => true,
    "write_offsets" => true,
    "write_counts" => true,
    "write_sparsity_nodebnds" => true,
    "write_dofs" => false,
    "write_interfaces" => false,
    "write_boundaries" => false,
    "write_sharedboundaries" => false,
    "exact_visualization" => true,
    "use_linear_metrics" => true,
    )
#    interp_op = [0.5 0 0; 0 0.5 0; 0 0 0.5]

    sbp = getTriSBPOmega(degree=order)
#    sbp = TriSBP{Float64}(degree=order, internal=true)
    vtx = sbp.vtx
#    interp_op = SummationByParts.buildinterpolation(sbp, vtx.')

    sbpface = TriFace{Float64}(order, sbp.cub, vtx)
    println("sbpface.numnodes = ", sbpface.numnodes)

    # make a complex mesh because it is needed later
    # create it first so the underlying Pumi mesh gets destroyed and
    # replaced by the Float64 version
    mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface, dofpernode=4)
#    mesh_c = PumiMeshDG2{Complex128, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
#    mesh =  PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

   @test ( mesh.m_ptr )!=C_NULL
   @test ( mesh.mnew_ptr )!=C_NULL
   @test ( mesh.numVert )== 4
   @test ( mesh.numEdge )== 5
   @test ( mesh.numFace )== mesh.numEdge
   @test ( mesh.numEl )== 2
   @test ( mesh.numDof )== 24
   @test ( mesh.numNodes )== 6
   @test ( mesh.numDofPerNode )== 4
   @test ( mesh.numBoundaryFaces )== 4
   @test ( mesh.numInterfaces )== 1
   @test ( mesh.numNodesPerElement )== 3
   @test ( mesh.numNodesPerType )== [0, 0, 3]
   @test ( mesh.numEntitiesPerType )== [4, 5, 2]
   @test ( mesh.numTypePerElement )== [3, 3, 1]
   @test ( mesh.typeOffsetsPerElement )== [1, 1, 1, 4]
   @test ( mesh.typeOffsetsPerElement_ )== mesh.typeOffsetsPerElement
   @test isapprox( mesh.volume, 4.0) atol=1e-12
   volume2 = PdePumiInterface.calcVolume(mesh)
   @test isapprox( volume2, mesh.volume) atol=1e-12
   @test ( mesh.dim )== 2
   @test ( mesh.isDG )== true
   @test ( mesh.coloringDistance )== 2
   @test ( mesh.numColors )== 4
   @test ( mesh.numBC )== 1
   @test ( mesh.bndry_geo_nums[1] )== opts["BC1"]
   @test ( mesh.elementNodeOffsets )== zeros(mesh.numNodesPerElement, mesh.numEl)
   @test ( mesh.typeNodeFlags[1] )== trues(3, mesh.numEl)
   tmp = trues(3, 2); tmp[1, 1] = false
   @test ( mesh.typeNodeFlags[2] )== tmp
   @test ( mesh.bndry_offsets )== [1, 5]
   iface = mesh.interfaces[1]
   @test ( iface.elementL )== 1
   @test ( iface.elementR )== 2
   @test ( iface.faceL )== 1
   @test ( iface.faceR )== 3
   @test isapprox( mesh.coords[:, :, 1], [-2/3 -2/3 1/3; 2/3 -1/3 2/3]) atol=1e-13
   @test isapprox( mesh.coords[:, :, 2], [2/3 -1/3 2/3; 1/3 -2/3 -2/3]) atol=1e-13
   @test ( sort(unique(mesh.dofs)) )== collect(1:mesh.numDof)

   @test ( mesh.color_masks[1][1] )== true
   @test ( mesh.color_masks[1][2] )== false
   @test ( mesh.color_masks[2][1] )== false
   @test ( mesh.color_masks[2][2] )== true
   @test ( mesh.neighbor_colors[1, 1] )== 2
   @test ( mesh.neighbor_colors[2, 1] )== 1
   @test ( mesh.neighbor_colors[1, 2] )== 1
   @test ( mesh.neighbor_colors[2, 2] )== 2

   @test ( mesh.neighbor_nums[1, 1] )== 2
   @test ( mesh.neighbor_nums[2, 1] )== 1
   @test ( mesh.neighbor_nums[1, 2] )== 1
   @test ( mesh.neighbor_nums[2, 2] )== 2

   @test ( mesh.pertNeighborEls[1, 1] )== 1
   @test ( mesh.pertNeighborEls[2, 1] )== 1
   @test ( mesh.pertNeighborEls[1, 2] )== 2
   @test ( mesh.pertNeighborEls[2, 2] )== 2

   @test ( mesh.color_cnt[1] )== 1
   @test ( mesh.color_cnt[2] )== 1

   @test ( mesh.coord_order )== 1
   @test ( mesh.coord_numNodesPerElement )== 3
   @test ( length(mesh.typeOffsetsPerElement) )== 4
   @test ( size(mesh.coord_xi, 1) )== 2
   @test ( size(mesh.coord_xi, 2) )== 3

   # check that adjoint variables are right size
  @test ( size(mesh.dxidx) )== size(mesh.dxidx_bar)
  @test ( size(mesh.dxidx_bndry) )== size(mesh.dxidx_bndry_bar)
  @test ( size(mesh.dxidx_face) )== size(mesh.dxidx_face_bar)
  @test ( size(mesh.dxidx_sharedface_bar) )== size(mesh.dxidx_sharedface_bar)

  @test isapprox( mesh.jac, ones(mesh.numNodesPerElement ,2))
  test_coordNumbering(mesh)

  # check if dxidx is consistent with the old way of calculating it
  dxidx2 = zeros(mesh.dxidx)
  jac2 = zeros(mesh.jac)

  SummationByParts.mappingjacobian!(sbp, mesh.coords, dxidx2, jac2)

  @test isapprox( norm(mesh.jac - jac2)/length(mesh.jac), 0.0) atol=1e-13

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      @test isapprox( norm(mesh.dxidx[:, :, j, i]  - dxidx2[:, :, j, i]), 0.0) atol=1e-13
    end
  end

  u = ones(mesh.numDofPerNode, mesh.numTypePerElement[1], mesh.numEl)
  PdePumiInterface.saveNodalSolution(mesh, u)
  writeVisFiles(mesh, "nodal_solution")

  # test accumulation at vertices
  u_volume = ones(6, mesh.numTypePerElement[1], mesh.numEl)
  u_verts = ones(6, mesh.numVert)

  PdePumiInterface.accumulateAtVerts(mesh, u_volume, u_verts)

  # check that value at each vert is the number of elements using that vert
  for i=1:mesh.numVert
    vert_i = mesh.verts[i]
    nel = apf.countAdjacent(mesh.m_ptr, vert_i, mesh.dim)
    for j=1:6
      @test ( u_verts[j, i] )== nel
    end
  end

  # test surface numbering
  testSurfaceNumbering(mesh, sbp, opts)

  # check reverse mode
  # SBP testing the correctness, these tests only verify values get to the right place

  fill!(mesh.dxidx_bar, 1.0)
  PdePumiInterface.getVertCoords_rev(mesh, sbp)

  for i=1:mesh.numEl
    @test  norm(mesh.vert_coords_bar[:, :, i])  > 0.0
  end

  # make sure the data fields are the same
  fill!(mesh_c.vert_coords, -1)
  PdePumiInterface.copy_data!(mesh_c, mesh)
  compare_meshes(mesh, mesh_c)

  # test zero_bar_arrays
  fill!(mesh.vert_coords_bar, 1.0)
  zeroBarArrays(mesh)
  @test ( vecnorm(mesh.vert_coords_bar) )== 0.0

   function test_interp(mesh::AbstractMesh{Tmsh}) where Tmsh
     sbpface = mesh.sbpface
     dxdxi_element = zeros(2, 2, mesh.numNodesPerElement, 1)
     dxdxi_face = zeros(4, sbpface.numnodes, 1)
     dxidx_face = zeros(2,2, sbpface.numnodes)
     jac_face = zeros(sbpface.numnodes)

     for i=1:mesh.numInterfaces
       el = mesh.interfaces[i].elementL
       face = mesh.interfaces[i].elementR

       for j=1:mesh.numNodesPerElement
         dxidx_hat = mesh.dxidx[:, :, j, el]
         dxidx = dxidx_hat*mesh.jac[j, el]
         dxdxi = inv(dxidx)
         dxdxi_element[:, :, j, el] = dxdxi
       end

       bndry_arr = [ Boundary(1, face)]
       dxdxi_element_rshape = reshape(dxdxi_element, 4, mesh.numNodesPerElement, 1)
       boundaryinterpolate!(sbpface, bndry_arr, dxdxi_element_rshape, dxdxi_face)
       dxdxi_face_rshape = reshape(dxdxi_face, 2, 2, sbpface.numnodes)
    for j=1:sbpface.numnodes
      dxidx = inv(dxdxi_face_rshape[:, :, j])
      jac_face[j] = det(dxidx)
      dxidx_face[:, :, j] = dxidx/det(dxidx)

      @test isapprox( jac_face[j], mesh.jac_face[j, i]) atol=1e-13
      @test isapprox( dxidx_face[:, :, j], mesh.dxidx_face[:, :, j, i]) atol=1e-13
    end

    end  # end loop over interfaces

  end  # end function

  println("testing dxidx_face")
  if size(mesh.dxidx_face, 1) > 0
    println("running the test")
    test_interp(mesh)
  end

   # check dxidx_bndry
   # for straight sided elements, dxidx is constant
   println("testing dxidx_bndry")
   if size(mesh.jac_bndry, 1) > 0
     println("running the test")
     for i=1:mesh.numBoundaryFaces
       el = mesh.bndryfaces[i].element
       dxidx_test = mesh.dxidx[:, :, 1, el]
       jac_test = mesh.jac[1, el]
       for j=1:mesh.sbpface.numnodes
         dxidx_code = mesh.dxidx_bndry[:, :, j, i]
         jac_code = mesh.jac_bndry[j, i]
         @test isapprox( dxidx_code, dxidx_test) atol=1e-13
         @test isapprox( jac_code, jac_test) atol=1e-13
       end
     end
   end

  println("testing coords bndry")

   # check coords_bndry
   vert_coords = [-1.0 -1.0; 1  1; -1 1]

   function getBndry(arr, el, face)
     for i=1:length(arr)
       bndry_i = arr[i]
       if bndry_i.element == el && bndry_i.face == face
         return bndry_i, i
       end
     end
   end

   bndry, idx = getBndry(mesh.bndryfaces, 1, 2)
   # get vertex coordinates
   v1 = vert_coords[2, :]
   v2 = vert_coords[3, :]
   diff = v2 - v1
   slope = diff[2]/diff[1]
   b = v1[2] - slope*v1[1]

   # vertify the face node coordinates are on the line
   for i=1:mesh.sbpface.numnodes
     x_i = mesh.coords_bndry[1, i, idx]
     y_i = mesh.coords_bndry[2, i, idx]
    @test isapprox( y_i, slope*x_i + b) atol=1e-13
   end

   function test_normal_orientation(mesh, ifaces::AbstractArray{I}, nrm::AbstractArray{T, 3}) where {I <: Union{Boundary, Interface}, T}
     println("-----entered test_normal_orientation-----")
     topo = mesh.topo
     numVertPerElement = mesh.numTypePerElement[1]
     numVertPerFace = numVertPerElement - 1
     println("numVertPerFace = ", numVertPerFace)
     el_verts = Array{Ptr{Void}}(numVertPerElement)
     other_vert_coords = zeros(mesh.dim)
     face_verts = Array{Ptr{Void}}(numVertPerElement - 1)
     face_vert_coords = zeros(mesh.dim, numVertPerFace)

     println("face_verts = ", topo.face_verts)
     nfaces = length(ifaces)
     for i=1:nfaces
       println("i = ", i)
       iface_i = ifaces[i]
       elnum = getElementL(iface_i)
       facenum_local = getFaceL(iface_i)

       el_i = mesh.elements[elnum]
       nverts = apf.getDownward(mesh.m_ptr, el_i, 0, el_verts)
       println("nverts = ", nverts)

       for j=1:numVertPerFace
         face_verts[j] = el_verts[topo.face_verts[j, facenum_local]]
         tmp = zeros(3)
         apf.getPoint(mesh.m_ptr, face_verts[j], 0, tmp)
         face_vert_coords[1:mesh.dim, j] = tmp[1:mesh.dim]
       end

       # get the vert not on the face
       other_vert = Ptr{Void}(0)
       for j=1:numVertPerElement
         if !(el_verts[j] in face_verts)
           other_vert = el_verts[j]
         end
       end

       tmp = zeros(3)
       apf.getPoint(mesh.m_ptr, other_vert, 0, tmp)
       other_vert_coords[1:mesh.dim] = tmp[1:mesh.dim]

       # check that the face normal is in the opposite direction as the
       # vectors from a vertex on the face to the vertex not on the face
       for j=1:mesh.numNodesPerFace
         for k=1:numVertPerFace
           r1 = other_vert_coords - face_vert_coords[:, k]

           val = dot(nrm[:, j, i], r1)
           @test  val  < 0.0
         end
       end

     end  # end loop i

     return nothing

   end  # end function test_normal_orientation

   test_normal_orientation(mesh, mesh.interfaces, mesh.nrm_face)
   test_normal_orientation(mesh, mesh.bndryfaces, mesh.nrm_bndry)

   println("testing number nodes windy")
   # check adjacency apf.reordering algorithm doesn't error out
   PdePumiInterface.numberNodesWindy(mesh, [0.0, 0.0, 0.0])

    opts["smb_name"] = "tri8l.smb"
    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
#    mesh =  PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

  # check mapping interpolation
  # should be constant within an element for straight-sided elementsa
  if size(mesh.jac_face, 1) > 0
    for i=1:mesh.numInterfaces
      iface_i = mesh.interfaces[i]
      el_i = iface_i.elementL
      dxidx_el = mesh.dxidx[:, :, 1, el_i]
      jac_el = mesh.jac[:, el_i]
      jac_face = mesh.jac_face[:, i]

      for j=1:mesh.numNodesPerFace
        dxidx_face = mesh.dxidx_face[:, :, j, i]
        for k=1:2
          for p=1:2
            @test isapprox( dxidx_face[p, k], dxidx_el[p, k]) atol=1e-13
          end
        end
        @test isapprox( jac_face[j], jac_el[j]) atol=1e-13
      end
    end  # end loop over interfaces
  end

  # test reverse mode interpolation

  test_interp_rev(mesh)
  test_coords_rev(mesh, sbp)

  # test update_coords
  println("testing update_coords")
  coords_orig = copy(mesh.vert_coords)

  # double all coordinates
  coords_i = zeros(Float64, 2, 3)
  for i=1:mesh.numEl
    for j=1:3
      coords_i[1, j] = 2*coords_orig[1, j, i]
      coords_i[2, j] = 2*coords_orig[2, j, i]
    end
    update_coords(mesh, i, coords_i)
  end

  commit_coords(mesh, sbp, opts)

#  PdePumiInterface.getCoordinatesAndMetrics(mesh, sbp)
  for i = 1:length(mesh.vert_coords)
    @test  abs(2*coords_orig[i] - mesh.vert_coords[i])  < 1e-10
  end

  # test vertmap
  println("testing vertmap")
  nedges_interior = 0
  for i=1:mesh.numEl
    vertnums_i = mesh.element_vertnums[:, i]
    for j=1:3
      @test  vertnums_i[j]  > 0
      @test  vertnums_i[j]  < mesh.numVert + 1
    end

    # test there are exactly 3 elements the current element shares 2 vertices
    # with
    nedges = 0
    for j=1:mesh.numEl
      vertnums_j = mesh.element_vertnums[:, j]
      if length(intersect(vertnums_i, vertnums_j)) == 2
        nedges += 1
        nedges_interior += 1
      end
    end

    @test  nedges  < 4
    @test  nedges  > 0
  end  # end loop i

  @test ( nedges_interior )== 2*(mesh.numEdge - mesh.numBoundaryFaces)

  # test saveSolutionToMesh interpolation
  println("testing saveSolutionToMesh interpolation")
  u_vals = zeros(mesh.numDof)
  for i=1:mesh.numDof

    # get dof, node, element
    idx = findfirst(mesh.dofs, i)
    dofidx, node, el = ind2sub(mesh.dofs, idx)

    x = mesh.coords[1, node, el]
    y = mesh.coords[2, node, el]
    order = mesh.order

    u_vals[i] = x^order + y^order + 1
  end

  saveSolutionToMesh(mesh, u_vals)
  writeVisFiles(mesh, "dg_vis_test")

  # check that the solution is interpolated exactly
  down_verts = Array{Ptr{Void}}(mesh.numTypePerElement[1])
  coords_vert = zeros(Float64, 3)
  interp_vals = zeros(mesh.numDofPerNode)
  for i=1:mesh.numEl
    order = mesh.order
    el_ptr = mesh.elements[i]
    apf.getDownward(mesh.m_ptr, el_ptr, 0, down_verts)

    for j=1:mesh.numTypePerElement[1]  # loop over vertices
      vert_j = down_verts[j]
      apf.getPoint(mesh.m_ptr, vert_j, 0, coords_vert)

      x = coords_vert[1]
      y = coords_vert[2]

      # note: this assumes mnew = m
      apf.getComponents(mesh.fnew_ptr, vert_j, 0, interp_vals)

      val_expected = x^order + y^order + 1
      for k=1:mesh.numDofPerNode
        @test  abs(interp_vals[k] - val_expected)  < 1e-12
      end
    end
  end

  # test retrieveing solution from mesh
  u_vals2 = zeros(u_vals)
  PdePumiInterface.retrieveSolutionFromMesh_interp(mesh, u_vals2)

  for i=1:length(u_vals2)
    @test abs(u_vals2[i] - u_vals[i]) < 1e-12
  end

  # just for good measure, create a new mesh
  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
#  mesh =  PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)
  mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface, dofpernode=4)
#  mesh_c =  PumiMeshDG2{Complex128, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

  compare_meshes(mesh, mesh_c)

  # test periodic
  println("testing periodic")
  @test ( mesh.numPeriodicInterfaces )== 0

  opts["smb_name"] = "tri3_px.smb"
  opts["BC1"] = [0, 2]
  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
  mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface, dofpernode=4)
#  mesh = PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

  @test ( mesh.numPeriodicInterfaces )== 3
  @test ( length(mesh.interfaces) )== 24
  @test ( mesh.numInterfaces )== 24
  @test ( mesh.numBoundaryFaces )== 6
  @test ( opts["numBC"] )== 1  # the default BC should not be created for periodic

  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    @test  iface_i.elementL  > 0
    @test  iface_i.elementL  < mesh.numEl + 1
    @test  iface_i.elementR  > 0
    @test  iface_i.elementR  < mesh.numEl + 1
    @test  iface_i.faceL  < 4
    @test  iface_i.faceR  > 0
    @test  iface_i.faceR  < 4
  end

  test_coordNumbering(mesh)

  # test curvilinear
  println("testing curvilinear")
  # check curvilinear routines reproduce linear results
  coords_face_orig = copy(mesh.coords_interface)
  nrm_face_orig = copy(mesh.nrm_face)

  PdePumiInterface.getCoordinates(mesh, sbp)
  PdePumiInterface.getFaceCoordinatesAndNormals(mesh, sbp)
  PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh, sbp)

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      @test isapprox( mesh.coords_interface[:, j, i], coords_face_orig[:, j, i]) atol=1e-12
      @test isapprox( mesh.nrm_face[:, j, i], nrm_face_orig[:, j, i]) atol=1e-12
    end
  end

  # test metrics reverse
  test_metrics_rev(mesh, mesh_c, sbp, opts)

  # a 0 - 5 square that used a sin wave to remap the nondimensionalized
  # coordinates
  opts["smb_name"] = "square_05_curve.smb"
  opts["use_linear_metrics"] = false
  opts["BC1"] = [0, 1, 2, 3]
  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
  mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface, dofpernode=4)
#  mesh =  PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

  testSurfaceNumbering(mesh, sbp, opts)
  test_coordNumbering(mesh)
  test_coord_field(mesh)
  # TODO: re-enable when Pumi's ENABLE_ZOLTAN flag is available
  #test_split(mesh)
  test_submesh_transfer(mesh)

  #TODO: check sizes of arrays

  function test_volume_curvilinear(mesh, sbp)

    volume = 0.0
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        dxidx_scaled = mesh.dxidx[:, :, j, i]
        dxidx = dxidx_scaled*mesh.jac[j, i]

        dxdxi = inv(dxidx)
        jac = det(dxdxi)
        volume += sbp.w[j]*jac
      end
    end

    println("volume = ", volume)
    @test isapprox( volume, 25.0) atol=1e-12
  end  # end function

  test_volume_curvilinear(mesh, sbp)

  # translate the mesh and verify all quantities are the same
  dxidx_orig = copy(mesh.dxidx)
  jac_orig = copy(mesh.jac)
  nrm_bndry_orig = copy(mesh.nrm_bndry)
  nrm_face_orig = copy(mesh.nrm_face)

  for i=1:mesh.numEl
    coords_i = mesh.vert_coords[:, :, i]
    coords_i += 1
    update_coords(mesh, i, coords_i)
  end

  commit_coords(mesh, sbp, opts)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.dim
        for p=1:mesh.dim
          @test isapprox( mesh.dxidx[p, k, j, i], dxidx_orig[p, k, j, i]) atol=1e-12
        end
      end

      @test isapprox( mesh.jac[j, i], jac_orig[j, i]) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_bndry[k, j, i], nrm_bndry_orig[k, j, i]) atol=1e-12
      end
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_face[k, j, i], nrm_face_orig[k, j, i]) atol=1e-12
      end
    end
  end

  # recalculate and verify all quantities are the same
  recalcCoordinatesAndMetrics(mesh, sbp, opts)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      for k=1:mesh.dim
        for p=1:mesh.dim
          @test isapprox( mesh.dxidx[p, k, j, i], dxidx_orig[p, k, j, i]) atol=1e-12
        end
      end

      @test isapprox( mesh.jac[j, i], jac_orig[j, i]) atol=1e-12
    end
  end

  for i=1:mesh.numBoundaryFaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_bndry[k, j, i], nrm_bndry_orig[k, j, i]) atol=1e-12
      end
    end
  end

  for i=1:mesh.numInterfaces
    for j=1:mesh.numNodesPerFace
      for k=1:mesh.dim
        @test isapprox( mesh.nrm_face[k, j, i], nrm_face_orig[k, j, i]) atol=1e-12
      end
    end
  end

  # test metrics reverse
  test_metrics_rev(mesh, mesh_c, sbp, opts)
  test_metrics_rev_1d(mesh_c, sbp, opts)

  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
#  mesh =  PumiMeshDG2{Float64, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)
  mesh_c = PumiMeshDG2(Complex128, sbp, opts, sbpface, dofpernode=4)
#  mesh_c =  PumiMeshDG2{Complex128, typeof(sbpface)}(dmg_name, smb_name, order, sbp, opts, sbpface, coloring_distance=2, dofpernode=4)

  compare_meshes(mesh, mesh_c)

  test_submesh()

  test_adapt_2d()

  test_ScatterData(mesh)

  # do this last since it rewrites the mesh coordinate field
  test_setPoint(mesh, opts)

  opts["smb_name"] = "meshes/tri4x4_.smb"
  opts["dmg_name"] = "meshes/tri4x4_.dmg"
  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
  test_geoNums(mesh)

  println("\nTesting coords file")
  test_coords_file(opts, sbp, sbpface)

  if HAVE_SIMMETRIX

    # load airfoil CAD mesh
    opts = Dict{Any, Any}(
    "dimensions" => 2,
    "run_type" => 5,
    "jac_type" => 2,
    "order" => 1,
    "use_DG" => true,
    "coloring_distance" => 2,
    "numBC" => 2,
    "BC1" => [40],
    "BC1_name" => "FreeStreamBC",
    "BC2" => [5],
    "BC2_name" => "noPenetrationBC",
    "smb_name" => "meshes/airfoil2_.smb",
    "dmg_name" => "meshes/airfoil2.smd",
    )

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    test_geoMapping(mesh)
    test_geoWrapping(mesh)

    # the airfoil has large curvature at the trailing edge, so its hard to find
    # a finite difference step size that works.  Use a smoother shape instead.

    # load airfoil CAD mesh
    opts = Dict{Any, Any}(
    "dimensions" => 2,
    "run_type" => 5,
    "jac_type" => 2,
    "order" => 1,
    "use_DG" => true,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" => [4, 7, 10, 13],
    "BC1_name" => "FreeStreamBC",
    "smb_name" => "meshes/UnitSquare/UnitSquare.smb",
    "dmg_name" => "meshes/UnitSquare/UnitSquare.x_t",
    )

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    @time test_geoDerivative(mesh)
    println("test_geoDerivative @time printed above")

    println("\n\nTesting geometric derivative on curved geometry")
    # load smooth CAD mesh
    opts = Dict{Any, Any}(
    "dimensions" => 2,
    "run_type" => 5,
    "jac_type" => 2,
    "order" => 1,
    "use_DG" => true,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" => [4, 7, 10, 13],
    "BC1_name" => "FreeStreamBC",
    "smb_name" => "meshes/UnitSquareCurve/UnitSquareCurve_linear.smb",
    "dmg_name" => "meshes/UnitSquareCurve/UnitSquareCurve.x_t",
    )

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    test_geoDerivative(mesh)

    test_coords_file(opts, sbp, sbpface)

    xivec = getXiCoords(mesh)
    xivec_pert = rand(length(xivec))
    for i=1:length(xivec)
      xivec[i] += 1e-3*xivec[i]
    end
    update_coordsXi(mesh, sbp, opts, xivec)
    test_geoDerivative(mesh)
    test_geoWarp(mesh, sbp, opts)
    test_update_coords(mesh, sbp, opts)

    # this following test messes up the coordinate field, so load a new mesh
    # for this test
    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    test_setPoint(mesh, opts)

    #TODO: test quadratic mesh once that works

    opts = Dict{Any, Any}(
    "dimensions" => 2,
    "run_type" => 5,
    "jac_type" => 2,
    "order" => 1,
    "use_DG" => true,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" => [4, 7, 10, 13],
    "BC1_name" => "FreeStreamBC",
    "smb_name" => "meshes/UnitSquareCurve/UnitSquareCurve_quadratic.smb",
    "dmg_name" => "meshes/UnitSquareCurve/UnitSquareCurve.x_t",
    )

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    println("testing quadratic mesh")
    test_setPoint(mesh, opts)
    test_coords_file(opts, sbp, sbpface)

    # load new mesh because test_setPoint messes up the coordinate field
    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    test_geoDerivative(mesh)

    # test mesh with geometry that has an internal edge
    println("\n\nTesting internal edge")
    opts = Dict{Any, Any}(
    "dimensions" => 2,
    "run_type" => 5,
    "jac_type" => 2,
    "order" => 1,
    "use_DG" => true,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" => [7, 56, 60, 13, 4, 10, 16],
    "BC1_name" => "FreeStreamBC",
    "smb_name" => "meshes/corner20_split/corner20_split_mesh0_.smb",
    "dmg_name" => "meshes/corner20_split/corner20_split.smd",
    )

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)
    test_interioredge(mesh)
  end

end


text:
function isequal(S::SpikeTrains, Q::SpikeTrains)
    typeof(S) == typeof(Q) || return false
    for (fp,fq) in zip(fieldnames(typeof(S)),fieldnames(typeof(Q)))
        isequal(getfield(S, fp), getfield(Q, fq)) || return false
    end
    return true
end

function hash(S::SpikeTrains, h::UInt)
    _h = h
    for fn in fieldnames(S)
        _h = hash(getfield(S, fn), _h)
    end
    return _h
end

################################################################################
#### File IO
################################################################################
"""
    savespikes(S::SpikeTrains; [filename], [dir], [old])

Save the given spike trains to disk using the `JLD` package. Returns the full path to the
saved file. Deletes file `old` if it exists in `dir`. By default `old` is `nothing` to avoid
unnecessarily deleting and rewriting files.

Default `filename` is `string(hash(S))`. Default `dir` is `Pkg.dir(Spikes)/saved`.

If `dir` doesn't exist, will use `mkpath`. If file exists, contents will be overwritten.
"""
function savespikes(S::SpikeTrains;
    filename=string(hash(S)),
    dir=joinpath(Pkg.dir("Spikes"), "saved"),
    old=nothing)

    # _dir = abspath(dir)
    if !ispath(dir)
        mkpath(dir)
    end
    if old != nothing
        _old = endswith(old, ".jld") ? old : old * ".jld"
        if isfile(joinpath(dir, _old))
            rm(joinpath(dir, _old))
        end
    end
    _fn = endswith(filename, ".jld") ? filename : filename * ".jld"
    # save(joinpath(_dir, _fn), "S", S)
    jldopen(joinpath(dir, _fn), "w") do file
        addrequire(file, Spikes)
        write(file, "S", S)
    end
    return joinpath(dir, _fn)
end

"""
    loadspikes(filename; [dir])
    loadspikes(h; [dir])

Returns the stimulus saved in `filename` (or identified by hash `h`) as saved by
`savespikes`. Default `dir` is `Pkg.dir(Spikes)/saved`

"""
function loadspikes(filename; dir=joinpath(Pkg.dir("GrayScaleStimuli"), "saved"))
    _fn = endswith(filename, ".jld") ? filename : filename * ".jld"
    _fn = basename(_fn)
    return load(joinpath(dir, _fn), "S")
end
loadspikes(h::UInt; dir=joinpath(Pkg.dir("Spikes"),"saved")) = loadspikes(string(h); dir=dir)


text:
module test_octo_pretty

using Test
using Octo.Pretty

nts = [(id=1,name="John"),(id=2,name="Tom")]
@test Pretty.table(nts) == """
|   id | name   |
| ---- | ------ |
|    1 | John   |
|    2 | Tom    |
2 rows."""

Pretty.set(colsize=30)
nts = [(id=1,case="글씨가 넘흐 길어서 짤릴 때"),
       (id=2,case="글씨에 점  하나가"),
       (id=3,case="안짤리는 넘")]
@test Pretty.table(nts) == """
|   id | case              |
| ---- | ----------------- |
|    1 | 글씨가 넘흐 길... |
|    2 | 글씨에 점  하나가 |
|    3 | 안짤리는 넘       |
3 rows."""

nts = Vector{NamedTuple{(:a,),Tuple{Int}}}()
@test Pretty.table(nts) == """
| a   |
| --- |

0 rows."""

Pretty.set(false)
buf = IOBuffer()
Base.show(buf, MIME"text/plain"(), nts)
@test String(take!(buf)) == """
NamedTuple{(:a,),Tuple{Int64}}[]"""
Pretty.set(true)

@test Pretty._regularize_text("가1", 1) == "가"

end # module test_octo_pretty


text:
# %% Ti Tests
@testset "Temporal Tests" begin

@testset "Temporal weighting" begin
    # spatiotemporal weights:
    W = zero(A)
    W[Ti(5)] .= 1.0
    res = timeagg(mean, A, W)
    @test all(res .≈ A[Ti(5)])
    @test dims(A, Lon).val == dims(res, Lon).val
    res = timeagg(std, A, W)
    @test all(x -> isapprox(x, 0; atol = 1e-8), res)

    # just time weight
    w = zeros(length(t))
    w[5] = 1.0
    res = timeagg(mean, A, w)
    @test all(res .≈ A[Ti(5)])
    @test dims(A, Lon).val == dims(res, Lon).val

    # TODO: more tests needed here, e.g. for timeagg(mean, t, a, w)
end

@testset "Advanced temporal manipulation" begin
    tyearly = Date(2000, 3, 1):Year(1):Date(2020, 3, 31)
    thourly = DateTime(2000, 3, 1):Hour(1):DateTime(2001, 4, 15)
    mdates = unique!([(year(d), month(d)) for d in tdaily])
    ydates = unique!([year(d) for d in tdaily])
    tranges = temporalrange(tdaily, Dates.month)
    yranges = temporalrange(tdaily, Dates.year)
    @testset "time sampling" begin
        @test length(tranges) == length(mdates)
        @test length(yranges) == length(ydates)

        @test temporal_sampling(t) == :monthly
        @test temporal_sampling(tdaily) == :daily
        @test temporal_sampling(tyearly) == :yearly

        # Test hourly stuff
        @test temporal_sampling(thourly) == :hourly
        @test temporal_sampling(collect(thourly)) == :hourly
        hmys = maxyearspan(thourly, :hourly)
        @test hmys < length(thourly)
        @test year(thourly[1]) == year(thourly[hmys])-1

        # Arbitrary random range
        tother = collect(thourly)
        tother[4] = DateTime(2000, 7, 2, 2, 2)
        @test temporal_sampling(tother) == :other
    end

    @testset "sametimespan" begin
        tm1 = Date(2000, 3, 15):Month(1):Date(2010, 3, 15)
        tm2 = Date(2001, 1, 1):Month(1):Date(2011, 1, 1)
        A1 = ClimArray(rand(length(tm1)), (Ti(tm1),))
        A2 = ClimArray(rand(length(tm2)), (Ti(tm2),))
        Bv = sametimespan((A1, A2))
        # dict version
        d = Dict(:A1 => A1, :A2 => A2)
        d2 = sametimespan(d)
        Bd = (d2[:A1], d2[:A2])
        for B in (Bv, Bd)
            B1, B2 = B
            @test size(B1) == size(B2)
            @test dims(B1, Ti)[1] == Date(2001, 1, 15)
            @test dims(B2, Ti)[1] == Date(2001, 1, 1)
            @test dims(B1, Ti)[end] == Date(2010, 3, 15)
            @test dims(B2, Ti)[end] == Date(2010, 3, 1)
        end
        B1, B2 = sametimespan((A1, A2); mintime=Date(2005,1,1), maxtime=Date(2008,12,30))
        @test year(gnv(dims(B1, Ti))[1]) == year(gnv(dims(B2, Ti))[1]) == 2005
        @test year(gnv(dims(B1, Ti))[end]) == year(gnv(dims(B2, Ti))[end]) == 2008
    end

    @testset "monthlyagg and co." begin
        # test yearly temporal weights (i.e., no weighting)
        X = ClimArray(rand(3,3), (Lon(1:3), Ti(tyearly[1:3])))
        W = [0, 1, 0]
        @test vec(timemean(X).data) == vec(mean(X.data; dims = 2))
        @test timemean(X, W).data == X.data[:, 2]

        C = ClimArray(zeros(length(lats), length(tdaily)), (Lat(lats), Ti(tdaily)))

        # First version: just count number of days
        for j in 1:length(tdaily)
            for i in 1:length(lats)
                C[i, j] = daysinmonth(tdaily[j])
            end
        end
        Cm = monthlyagg(C)

        @test length(unique(Array(Cm))) == 4 # there are four unique number of days
        # test that each value when rounded to an integer is an integer (for first slice only
        # all remaining slices are the same)
        for e in Cm[:, 1]
            @test round(Int, e) == e
        end
        @test step(dims(Cm, Ti).val) == Month(1)
        @test temporal_sampling(dims(Cm, Ti).val) == :monthly

        for j in 1:length(tdaily)
            for i in 1:length(lats)
                C[i, j] = daysinyear(tdaily[j])
            end
        end
        Cy = yearlyagg(C)
        @test length(unique(Array(Cy))) == 2
        for e in Cy
            @test round(Int, e) == e
        end
        @test step(dims(Cy, Ti).val) == Year(1)
        @test temporal_sampling(dims(Cy, Ti).val) == :yearly

        # Second version: test actual physics
        for j in 1:length(tdaily)
            for i in 1:length(lats)
                C[i, j] = insolation(tdaily[j], lats[i])
            end
        end
        Bz = zonalmean(B)
        Cm = monthlyagg(C)
        @test size(Cm, Ti) == size(B, Ti) + 1
        # For times besides the first and last month, they sould be the same:
        for k in 5:length(t)-5
            @test all(Cm[Ti(k)] .≈ Bz[Ti(k)])
        end

        Asea = seasonalyagg(A)
        tsea = dims(Asea, Ti).val
        @test Base.step(tsea) == Month(3)
    end

    @testset "interannual variability" begin
        C = spacemean(B)
        for y in unique(year.(t))[1:end-1]
            C[Ti(At(Date(y, 3, 15)))] = 0
        end
        dates, vals = seasonality(C)
        @test length(dates) == 12
        @test all(v -> length(v) == 20, vals)
        @test month(dates[3]) == 3
        @test all(iszero, vals[3])
        @test all(!iszero, vals[2])
    end
end

end # tempooral tests


text:
# * When we say "variable" we mean an MOI variable and
# * when we say columns, we mean a Mosek variable.
# This is because Mosek variables are numbered 1:n and are columns of the `A`
# matrix while the MOI variables numbering can have holes, i.e. when variables
# are deleted.

###############################################################################
# TASK ########################################################################
###### The `task` field should not be accessed outside this section. ##########

num_columns(task::Mosek.MSKtask) = getnumvar(task)
num_columns(m::MosekModel) = num_columns(m.task)

add_column(task::Mosek.MSKtask) = appendvars(task, 1)
add_column(m::MosekModel) = add_column(m.task)

"""
    function init_columns(task::Mosek.MSKtask, cols::ColumnIndices)

Set the bound to free, i.e. `MSK_BK_FR`, in the internal Mosek task for the
columns cols.

See [`clear_columns`](@ref) which is kind of the reverse operation.
"""
function init_columns(task::Mosek.MSKtask, cols::ColumnIndices)
    # Set each column to a free variable
    N = length(cols.values)
    bnd = zeros(Float64, N)
    putvarboundlist(task, cols.values, fill(MSK_BK_FR, N), bnd, bnd)

    if DEBUG
        for col in cols.values
            putvarname(task, col, "x$col")
        end
    end
    return
end
function init_columns(m::MosekModel, vis::Vector{MOI.VariableIndex})
    init_columns(m.task, columns(m, vis))
end

## Name #######################################################################
###############################################################################
function set_column_name(task::Mosek.MSKtask, col::ColumnIndex, name::String)
    putvarname(task, col.value, name)
end
function set_column_name(m::MosekModel, vi::MOI.VariableIndex, name::String)
    set_column_name(m.task, mosek_index(m, vi), name)
end
column_name(task::Mosek.MSKtask, col::ColumnIndex) = getvarname(task, col.value)
function column_name(m::MosekModel, vi::MOI.VariableIndex)
    column_name(m.task, mosek_index(m, vi))
end
function column_with_name(task::Mosek.MSKtask, name::String)
    asgn, col = getvarnameindex(task, name)
    if iszero(asgn)
        return nothing
    else
        return col
    end
end
column_with_name(m::MosekModel, name::String) = column_with_name(m.task, name)

"""
    function clear_columns(task::Mosek.MSKtask, cols::Vector{Int32})

Delete all coefficients for columns of indices `cols`. Note that the column
is not actually deleted since it will be reused by a new variable when one
is added.

See [`init_columns`](@ref) which is kind of the reverse operation.
"""
function clear_columns(task::Mosek.MSKtask, cols::ColumnIndices)
    N = length(cols.values)
    # Objective: Clear any non-zeros in `c` vector
    putclist(task, cols.values, zeros(Int64, N))

    # Constraints: Clear any non-zeros in columns of `A` matrix
    putacollist(task,
                cols.values,
                zeros(Int64, N),
                zeros(Int64, N),
                Int32[],
                Float64[])

    # Bounds: Fix the variable to zero to make it have very low footprint for
    #         mosek in case `MOI.optimize!` is called before a new variable is
    #         added to reuse this column.
    bnd = zeros(Float64, N)
    putvarboundlist(task, cols.values, fill(MSK_BK_FX, N), bnd, bnd)

    if DEBUG
        for col in cols.values
            # Rename deleted column to help debugging
            putvarname(task, col, "deleted$col")
        end
    end
    return
end
function clear_columns(m::MosekModel, vis::Vector{MOI.VariableIndex})
    clear_columns(m.task, columns(m, vis))
end

###############################################################################
## INDEXING ###################################################################
###############################################################################

function column(m::MosekModel, vi::MOI.VariableIndex)
    @assert variable_type(m, vi) == ScalarVariable
    return ColumnIndex(getindex(m.x_block, ref2id(vi)))
end
function columns(m::MosekModel, vis::Vector{MOI.VariableIndex})
    return ColumnIndices(Int32[
        (mosek_index(m, vi)::ColumnIndex).value for vi in vis])
end
function variable_type(m::MosekModel, vi::MOI.VariableIndex)
    if 1 ≤ vi.value ≤ length(m.x_type)
        t = m.x_type[vi.value]
        if t != Deleted
            return t
        end
    end
    throw(MOI.InvalidIndex(vi))
end
function is_scalar(m::MosekModel, vi::MOI.VariableIndex)
    return variable_type(m, vi) == ScalarVariable
end
function is_matrix(m::MosekModel, vi::MOI.VariableIndex)
    return variable_type(m, vi) == MatrixVariable
end
function decide_variable(m::MosekModel, vi::MOI.VariableIndex)
    if variable_type(m, vi) == Undecided
        allocate_variable(m, vi)
        m.x_type[vi.value] = ScalarVariable
        init_columns(m, [vi])
    end
end
function mosek_index(m::MosekModel, vi::MOI.VariableIndex)
    decide_variable(m, vi)
    if is_scalar(m, vi)
        return column(m, vi)
    else
        return m.x_sd[vi.value]
    end
end

function index_of_column(m::MosekModel, col::Int32)
    id = m.x_block.back[col]
    if iszero(id)
        return nothing
    else
        return MOI.VariableIndex(id)
    end
end

"""
    function allocate_variable(m::MosekModel)

Allocate one variable. If there is a free column (because one variable was
deleted), we use it and don't create any new column, i.e. don't call
`appendvars`. Otherwise, we create new column.

See [`clear_variable`](@ref) which is kind of the reverse operation.
"""
function allocate_variable(m::MosekModel, vi::MOI.VariableIndex)
    numvar = num_columns(m)
    alloced = ensurefree(m.x_block, 1)
    if !iszero(alloced)
        @assert isone(alloced)
        @assert length(m.x_block) == num_columns(m) + 1
        add_column(m)
        @assert length(m.x_block) == num_columns(m)
    end
    allocate_block(m.x_block, 1, vi.value)
end

## Delete #####################################################################
######### Delete variables by clearing the column. The column is reused when ##
######### new variables is added. While there exists a function in the Mosek ##
######### API to delete a column, it is costly and is best avoided.          ##

function throw_if_cannot_delete(m::MosekModel, vi::MOI.VariableIndex)
    if !allocated(m.x_block, ref2id(vi))
        @assert m.x_type[ref2id(vi)] == Deleted
        throw(MOI.InvalidIndex(vi))
    end
    if is_scalar(m, vi)
        col = column(m, vi)
        for S in [MOI.LessThan{Float64}, MOI.GreaterThan{Float64},
                  MOI.EqualTo{Float64}, MOI.Interval{Float64},
                  MOI.Integer, MOI.ZeroOne]
            if has_flag(m, vi, S)
                MOI.delete(m, MOI.ConstraintIndex{MOI.SingleVariable, S}(vi.value))
            end
        end
        # All bounds have been removed so there can only be not constraint left
        # on the variable or a `VectorOfVariable`-in-`VectorCone` constraint.
        if has_flag(m, vi, VectorCone)
            throw(MOI.DeleteNotAllowed("Cannot delete variable $vi as it is involved in a `VectorOfVariables` constraint."))
        end
    end
end

"""
    clear_variable(m::MosekModel, vi::MOI.VariableIndex)

Delete the `vi` variable (without actually doing any change to the internal
Mosek solver) to free the corresponding column for reuse by another variable.

See [`allocate_variable`](@ref) which is kind of the reverse operation.
"""
function clear_variable(m::MosekModel, vi::MOI.VariableIndex)
    m.publicnumvar -= 1
    deleteblock(m.x_block, ref2id(vi))
end

###############################################################################
# MOI #########################################################################
###############################################################################

## Add ########################################################################
###############################################################################

function MOI.add_variable(m::MosekModel)
    m.publicnumvar += 1
    id = create_block(m.x_block, 1)
    push!(m.x_type, Undecided)
    push!(m.x_constraints, 0x0)
    push!(m.x_sd, MatrixIndex(0, 0, 0))
    push!(m.variable_to_vector_constraint_id, 0)
    return MOI.VariableIndex(id)
end

function MOI.add_variables(m::MosekModel, n::Int)
    @assert n ≥ 0
    return MOI.VariableIndex[MOI.add_variable(m) for i in 1:n]
end

## Delete #####################################################################
###############################################################################
function MOI.is_valid(model::MosekModel, vi::MOI.VariableIndex)
    return allocated(model.x_block, ref2id(vi))
end
function delete_vector_of_variables_constraint(m::MosekModel, vis::Vector{MOI.VariableIndex})
    i = first(vis).value
    id = m.variable_to_vector_constraint_id[i]
    if id > 0
        S = type_cone(getconeinfo(m.task, id)[1])
        ci = MOI.ConstraintIndex{MOI.VectorOfVariables, S}(i)
        if MOI.is_valid(m, ci) && vis == MOI.get(m, MOI.ConstraintFunction(), ci).variables
            MOI.delete(m, ci)
        end
    end
end
function MOI.delete(m::MosekModel, vis::Vector{MOI.VariableIndex})
    for vi in vis
        throw_if_cannot_delete(m, vi)
    end
    delete_vector_of_variables_constraint(m, vis)
    for vi in vis
        MOI.delete(m, vi)
    end
#    clear_columns(m, vis)
#    for vi in vis
#        clear_variable(m, vi)
#    end
end
function MOI.delete(m::MosekModel, vi::MOI.VariableIndex)
    if variable_type(m, vi) == Undecided
        m.publicnumvar -= 1
    else
        throw_if_cannot_delete(m, vi)
        delete_vector_of_variables_constraint(m, [vi])
        if !iszero(m.variable_to_vector_constraint_id[vi.value])
            MOIU.throw_delete_variable_in_vov(vi)
        end
        clear_columns(m, [vi])
        clear_variable(m, vi)
    end
    m.x_type[vi.value] = Deleted
    return
end

## List #######################################################################
###############################################################################
MOI.get(m::MosekModel, ::MOI.NumberOfVariables) = m.publicnumvar
function MOI.get(m::MosekModel, ::MOI.ListOfVariableIndices)
    ids = allocatedlist(m.x_block)
    return MOI.VariableIndex[MOI.VariableIndex(vid) for vid in ids]
end
function MOI.get(m::MosekModel, ::MOI.ListOfVariableAttributesSet)
    set = MOI.AbstractVariableAttribute[]
    if m.has_variable_names
        push!(set, MOI.VariableName())
    end
    # TODO add VariablePrimalStart when get is implemented on it
    return set
end

## Name #######################################################################
###############################################################################
MOI.supports(::MosekModel, ::MOI.VariableName, ::Type{MOI.VariableIndex}) = true
function MOI.set(m::MosekModel, ::MOI.VariableName, vi::MOI.VariableIndex,
                 name::String)
    m.has_variable_names = true
    set_column_name(m, vi, name)
end
function MOI.get(m::MosekModel, ::MOI.VariableName, vi::MOI.VariableIndex)
    return column_name(m, vi)
end
function MOI.get(m::MosekModel, ::Type{MOI.VariableIndex}, name::String)
    col = column_with_name(m, name)
    if col === nothing
        return nothing
    else
        return index_of_column(m, col)
    end
end