text:
## Hermite
@register0 Hermite AbstractCCOP0
export Hermite

"""
    Hermite

The [Hermite](https://en.wikipedia.org/wiki/Hermite_polynomials) polynomials have two versions the physicists (`Hermite`  or  `H`) and the probablalists (`ChebyshevHermite` or  `Hₑ`). They are  related through  `Hᵢ(x) =  2^(i/2) Hₑᵢ(√2 x)`.

The Hermite   polynomials have weight  function `w(x)=exp(-x^2/2)` and domain the real line.

## Examples
```jldoctest
julia> using Polynomials,  SpecialPolynomials

julia> x = variable(Polynomial{Rational{Int}})
Polynomial(x)

julia> [basis(Hermite, i)(x) for i in 0:5]
6-element Array{Polynomial{Float64},1}:
 Polynomial(1.0)
 Polynomial(2.0*x)
 Polynomial(-2.0 + 4.0*x^2)
 Polynomial(-12.0*x + 8.0*x^3)
 Polynomial(12.0 - 48.0*x^2 + 16.0*x^4)
 Polynomial(120.0*x - 160.0*x^3 + 32.0*x^5)

julia> [basis(ChebyshevHermite, i)(x) for i in 0:5]
6-element Array{Polynomial{Float64},1}:
 Polynomial(1.0)
 Polynomial(1.0*x)
 Polynomial(-1.0 + 1.0*x^2)
 Polynomial(-3.0*x + 1.0*x^3)
 Polynomial(3.0 - 6.0*x^2 + 1.0*x^4)
 Polynomial(15.0*x - 10.0*x^3 + 1.0*x^5)
```

!!! note
    The Hermite family needs help, as the computed values for `Bn`,and,`Cn` are  both 0.
"""
Hermite

basis_symbol(::Type{<:Hermite}) = "H"
Polynomials.domain(::Type{<:Hermite}) = Polynomials.Interval(-Inf, Inf)

abcde(::Type{<:Hermite})  = NamedTuple{(:a,:b,:c,:d,:e)}((1,0,0,-2,0))

k0(P::Type{<:Hermite}) = one(eltype(P))
function k1k0(P::Type{<:Hermite}, k::Int)
    val = 2*one(eltype(P))
    val
end

norm2(::Type{<:Hermite}, n) = sqrt(pi) * 2^n * gamma(n+1)
weight_function(::Type{<: Hermite})  = x -> exp(-x^2)
generating_function(::Type{<:Hermite}) = (t, x)  -> exp(2*t*x - t^2)
function classical_hypergeometric(::Type{<:Hermite}, n, x)
    as = iseven(n) ? (-n ÷ 2, -(n-1)/2) : (-n/2, -(n-1)÷2)
    bs = ()
    (2x)^n * pFq(as, bs, -inv(x)^2)
end

gauss_nodes_weights(p::Type{P}, n) where {P <: Hermite} =
    FastGaussQuadrature.gausshermite(n)

## Overrides
# Use override here, as we get  An,Bn,Cn =  1, 0,0, otherwise
An(P::Type{<:Hermite}, n::Int) = 2*one(eltype(P))
Bn(P::Type{<:Hermite}, n::Int) = zero(eltype(P))
Cn(P::Type{<:Hermite}, n::Int) = 2*n*one(eltype(P))

β̃n(P::Type{<:Hermite},  n::Int) = zero(eltype(P))
γ̃n(P::Type{<:Hermite},  n::Int) = zero(eltype(P))
b̂̃n(P::Type{<:Hermite}, ::Val{N}) where {N} = zero(eltype(P)) #where {M} = error("Don't call me")#zero(S)
ĉ̃n(P::Type{<:Hermite}, ::Val{N}) where {N}  = zero(eltype(P)) #where {M} = error("Don't call me")#zero(S)

## https://arxiv.org/pdf/1901.01648.pdf. Connection formula (14)
##  x^n  = n! sum(H_{n-2j}/ (2^j(n-2j)!j!) j = 0:floor(n/2))

Base.convert(P::Type{<:Hermite}, q::Polynomial) = connection(P,q)
function Base.iterate(o::Connection{P, Q}, state=nothing) where
    {P<:Hermite,
     Q<:Polynomials.StandardBasisPolynomial}

    k, n = o.k, o.n

    if state == nothing
        i = k
        j = 0
        i > n && return nothing
        val = __hermite_lambda(i,k)
    elseif state[1] + 2 > n # terminate
        return nothing
    else
        j,val1 = state
        #val1 *= (i+1)*(i+2)/4/(j+1/4)
        j += 1
        i = k + 2j
        val = __hermite_lambda(i,k)
    end

    return(i, val), (j, val)
end

function __hermite_lambda(n,k)
    val = gamma(1+n)/2^n
    val /= gamma(1+k)
    val /= gamma(1 + (n-k)/2)
    val
end

# 2x more performant
function Polynomials.integrate(p::P, C::Number=0) where {T,P<:Hermite{T}}
    # int H_n = 1/(n+1) H_{n+1}
    R = eltype(one(T)/1)
    d = degree(p)
    qs = zeros(R, d+2)

    for i in 0:d
        qs[i+2] = p[i]/(2(i+1))
    end

    q = ⟒(P)(qs, p.var)
    q = q - q(0) + R(C)

    return q
end

##
## --------------------------------------------------
##

@register0 MonicHermite AbstractCCOP0
export MonicHermite
ϟ(::Type{<:MonicHermite}) = Hermite
ϟ(::Type{<:MonicHermite{T}}) where {T} = Hermite{T}
@register_monic(MonicHermite)

pqr_start(::Type{<:MonicHermite}) = 0
pqr_symmetry(::Type{<:MonicHermite}) = true
pqr_weight(P::Type{<:MonicHermite}, n, x, dπx) = (one(P)*2)*exp(-x^2)/(dπx*dπx)

##
## --------------------------------------------------
##
@register0 ChebyshevHermite AbstractCCOP0
export ChebyshevHermite
"""
    ChebyshevHermite

Type for the Probabalist's  [Hermite](https://en.wikipedia.org/wiki/Hermite_polynomials) polynomials.

"""
ChebyshevHermite

basis_symbol(::Type{<:ChebyshevHermite}) = "Hₑ"
Polynomials.domain(::Type{<:ChebyshevHermite}) = Polynomials.Interval(-Inf, Inf)

# https://arxiv.org/pdf/1901.01648.pdf eqn 17
abcde(::Type{<:ChebyshevHermite})  = NamedTuple{(:a,:b,:c,:d,:e)}((0,0,1,-1,0))

kn(P::Type{<:ChebyshevHermite}, n::Int) =  one(eltype(P))
k1k0(P::Type{<:ChebyshevHermite}, k::Int)  = one(eltype(P))
k1k_1(P::Type{<:ChebyshevHermite}, k::Int) =  one(eltype(P))
ismonic(::Type{<:ChebyshevHermite}) = true

weight_function(P::Type{ChebyshevHermite})  = x -> exp(-(one(eltype(P))*x^2)/2)
generating_function(::Type{<:ChebyshevHermite}) = (t, x)  -> exp(t*x - t^2/2)
function classical_hypergeometric(::Type{<:ChebyshevHermite}, n, x)
    2^(-n/2)*classical_hypergeometric(Hermite, n,  x/sqrt(2))
end

function  gauss_nodes_weights(P::Type{<:ChebyshevHermite}, n)  where {α,β}
    xs, ws = glaser_liu_rokhlin_gauss_nodes(basis(ChebyshevHermite, n))
    xs, ws
end

pqr_start(::Type{<:ChebyshevHermite}) = 0
pqr_symmetry(::Type{<:ChebyshevHermite}) = true
pqr_weight(P::Type{<:ChebyshevHermite}, n, x, dπx) =  sqrt(2)*gamma(1+n)*sqrt(pi)/(dπx)^2

## Overrides
#An(P::Type{<:ChebyshevHermite}, n::Int) = one(eltype(P))
#Bn(P::Type{<:ChebyshevHermite}, n::Int) = zero(eltype(P))
#Cn(P::Type{<:ChebyshevHermite}, n::Int) = -one(eltype(P))*n

##
## --------------------------------------------------
##

##  Multiplication
## TODO: work  on  types
AbstractHermite{T} = Union{Hermite{T}, ChebyshevHermite{T}}

# Default is slow. Instead
# directly use a linearization formula from
# https://arxiv.org/pdf/1901.01648.pdf
# and  for Hermite, Hn(x) =2^(n/2)Hₑ_n(sqrt(2) x)
# so Hm⋅Hn = ∑ C_{m,n,j} 2^j H_{m+n-2j}
function ⊗(p::P, q::Q) where {T,S, P<:AbstractHermite{T}, Q<:AbstractHermite{S}}
    R = eltype(one(promote_type(T,S))/1)
    N,M = degree(p), degree(q)
    N == -1 && return zero(⟒(P){R}, q.var)
    M == -1 && return zero(⟒(P){R}, p.var)
    N == 0  && return q * p[0]
    M == 0  && return p * q[0]
    q.var != p.var && throw(ArgumentError("Variable names must   match"))

    as =  zeros(R, 1+N+M)
    for i  in eachindex(p)
        pᵢ = p[i]
        iszero(pᵢ) && continue
        for  j in eachindex(q)
            qⱼ = q[j]
            pᵢqⱼ = pᵢ*qⱼ
            iszero(qⱼ) &&  continue
            λᵤ = one(R)
            for  k in 0:min(i,j)
                as[1 + i  +  j - 2k] += pᵢqⱼ * λᵤ
                λᵤ *= (i-k) * (j-k) / (k+1) * hermite_α(P)
            end
        end
    end
    ⟒(P){R}(as, p.var)
end

hermite_α(P::Type{<:ChebyshevHermite}) = 1
hermite_α(P::Type{<:Hermite}) = 2

#⊗(p::Hermite, q::Hermite) = linearization_product(p, q)
#⊗(p::ChebyshevHermite, q::ChebyshevHermite) = linearization_product(p, q)

# function Base.iterate(o::Linearization{P,R}, state =  nothing) where {P <: AbstractHermite,R}

#     l, m, n = o.l, o.m, o.n
#     l  > m + n && return nothing

#     if state == nothing

#         #k =  0
#         j = 0
#         p = min(l, n)
#         q = l - p
#         val = linearization_α(P, R, j, l, p, q)

#     else
#         # we work with l + 2j as there is a parity consideration
#         j, p, q, val  = state
#         s = l + j # l +2j + p + q = l +2j + l = 2l + 2j, so s=l+j
#         if p == 0 ||  q == m || q  > s
#             j += 1
#             l + 2j > n+m && return nothing #  l + 2j  is too  big now
#             p = min(l + j, n)  # p <= s
#             q = l + 2j - p
#             q > s + 1 && return nothing
#             val = linearization_α(P, R, j, l, p, q)
#         else
#             p -= 1
#             q += 1

#             λ = q/(p+1)*(s-(q-1))/(s-p)
#             val *= λ
#         end

#     end

#     return (p,q,val), (j, p, q, val)
# end

# #  https://arxiv.org/pdf/1901.01648.pdf equation (74)
# function linearization_α(P::Type{<:AbstractHermite}, R, j, l, n, m)
#     s = l + j # pass in j to avoid s = divrem(l + m + n,  2) call
#     val = one(R)
#     val *= gamma(1+m)*gamma(1+n)/gamma(1+s-m)/gamma(1+s-n)/gamma(1+s-l)
#     val *= linearization_λ(P, l, m,n)
# end
# linearization_λ(::Type{<:Hermite}, l, m, n) = 2.0^((m+n-l)/2)
# linearization_λ(::Type{<:ChebyshevHermite}, l, m, n) = 1


text:
using FactCheck

include("returns.jl")
include("methods.jl")
include("const.jl")
include("helpers.jl")

exitstatus()

ext: jl


text:
using Base.Test
using ODLCommonTools
import ODLCommonTools.sview
using SummationByParts  # SBP operators

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))

using PDESolver
#using Base.Test
using PdePumiInterface  # common mesh interface - pumi
using EllipticEquationMod
using LinearSolvers
using NonlinearSolvers   # non-linear solvers
using OptimizationInterface
using ArrayViews
import ArrayViews.view
using EllipticEquationMod
using Utils
using MPI
using Input
using PETSc2
using Jacobian

#------------------------------------------------------------------------------
# define tests and tags

using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const EllipticTests = TestList()

include("test_conv_rate.jl")

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Elliptic tests -----" begin
  runTestSystem(EllipticTests, solvePDE, ARGS)
end

#------------------------------------------------------------------------------
# cleanup


text:
"""
$(TYPEDEF)
An abstract type that contains delay stochastic simulation algorithms:
- DelayDirect
-
"""
abstract type AbstractDelayAggregatorAlgorithm  end
"""
$(TYPEDEF)

"""
struct DelayRejection <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)

"""
struct DelayDirect <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)
A modifed Composition-Rejection Direct Method with delays (DelayDirectCR), implementation combining features from the original article and from the code in  `DiffEqJump` package: DirectCR :
*A constant-time kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks*,
by A. Slepoy, A.P. Thompson and S.J. Plimpton, J. Chem. Phys, 128, 205101 (2008).
and
*Efficient Formulations for Exact Stochastic Simulation of Chemical Systems*,
by S. Mauch and M. Stalzer, ACM Trans. Comp. Biol. and Bioinf., 8, No. 1, 27-35 (2010).
"""
struct DelayDirectCR <: AbstractDelayAggregatorAlgorithm end
"""
$(TYPEDEF)

"""
struct DelayMNRM <: AbstractDelayAggregatorAlgorithm end


text:
using Base: @deprecate

# BEGIN TimeZones 1.0 deprecations
# END TimeZones 1.0 deprecations


text:
@inline function _to_cartesian(a, i::CanonicalInt)
    @inbounds(CartesianIndices(ntuple(dim -> indices(a, dim), Val(ndims(a))))[i])
end
@inline function _to_linear(a, i::Tuple{CanonicalInt,Vararg{CanonicalInt}})
    _strides2int(offsets(a), size_to_strides(size(a), static(1)), i) + static(1)
end

"""
    is_splat_index(::Type{T}) -> StaticBool

Returns `static(true)` if `T` is a type that splats across multiple dimensions.
"""
is_splat_index(@nospecialize(x)) = is_splat_index(typeof(x))
is_splat_index(::Type{T}) where {T} = static(false)
_is_splat(::Type{I}, i::StaticInt) where {I} = is_splat_index(_get_tuple(I, i))

"""
    ndims_index(::Type{I}) -> StaticInt

Returns the number of dimension that an instance of `I` maps to when indexing. For example,
`CartesianIndex{3}` maps to 3 dimensions. If this method is not explicitly defined, then `1`
is returned.

"""
ndims_index(@nospecialize(i)) = ndims_index(typeof(i))
ndims_index(::Type{I}) where {I} = static(1)
ndims_index(::Type{<:AbstractCartesianIndex{N}}) where {N} = static(N)
ndims_index(::Type{<:AbstractArray{T}}) where {T} = ndims_index(T)
ndims_index(::Type{<:AbstractArray{Bool,N}}) where {N} = static(N)
ndims_index(::Type{<:LogicalIndex{<:Any,<:AbstractArray{Bool,N}}}) where {N} = static(N)
_ndims_index(::Type{I}, i::StaticInt) where {I} = ndims_index(_get_tuple(I, i))

"""
    to_indices(A, I::Tuple) -> Tuple

Converts the tuple of indexing arguments, `I`, into an appropriate form for indexing into `A`.
Typically, each index should be an `Int`, `StaticInt`, a collection with values of `Int`, or a collection with values of `CartesianIndex`
This is accomplished in three steps after the initial call to `to_indices`:

# Extended help

This implementation differs from that of `Base.to_indices` in the following ways:

*  `to_indices(A, I)` never results in recursive processing of `I` through
  `to_indices(A, axes(A), I)`. This is avoided through the use of an internal `@generated`
  method that aligns calls of `to_indices` and `to_index` based on the return values of
  `ndims_index`. This is beneficial because the compiler currently does not optimize away
  the increased time spent recursing through
    each additional argument that needs converting. For example:
    ```julia
    julia> x = rand(4,4,4,4,4,4,4,4,4,4);

    julia> inds1 = (1, 2, 1, 2, 1, 2, 1, 2, 1, 2);

    julia> inds2 = (1, CartesianIndex(1, 2), 1, CartesianIndex(1, 2), 1, CartesianIndex(1, 2), 1);

    julia> inds3 = (fill(true, 4, 4), 2, fill(true, 4, 4), 2, 1, fill(true, 4, 4), 1);

    julia> @btime Base.to_indices(\$x, \$inds2)
    1.105 μs (12 allocations: 672 bytes)
    (1, 1, 2, 1, 1, 2, 1, 1, 2, 1)

    julia> @btime ArrayInterface.to_indices(\$x, \$inds2)
    0.041 ns (0 allocations: 0 bytes)
    (1, 1, 2, 1, 1, 2, 1, 1, 2, 1)

    julia> @btime Base.to_indices(\$x, \$inds3);
    340.629 ns (14 allocations: 768 bytes)

    julia> @btime ArrayInterface.to_indices(\$x, \$inds3);
    11.614 ns (0 allocations: 0 bytes)

    ```
* Recursing through `to_indices(A, axes, I::Tuple{I1,Vararg{Any}})` is intended to provide
  context for processing `I1`. However, this doesn't tell use how many dimensions are
  consumed by what is in `Vararg{Any}`. Using `ndims_index` to directly align the axes of
  `A` with each value in `I` ensures that a `CartesiaIndex{3}` at the tail of `I` isn't
  incorrectly assumed to only consume one dimension.
* `Base.to_indices` may fail to infer the returned type. This is the case for `inds2` and
  `inds3` in the first bullet on Julia 1.6.4.
* Specializing by dispatch through method definitions like this:
  `to_indices(::ArrayType, ::Tuple{AxisType,Vararg{Any}}, ::Tuple{::IndexType,Vararg{Any}})`
  require an excessive number of hand written methods to avoid ambiguities. Furthermore, if
  `AxisType` is wrapping another axis that should have unique behavior, then unique parametric
  types need to also be explicitly defined.
* `to_index(axes(A, dim), index)` is called, as opposed to `Base.to_index(A, index)`. The
  `IndexStyle` of the resulting axis is used to allow indirect dispatch on nested axis types
  within `to_index`.
"""
to_indices(A, ::Tuple{}) = ()
@inline function to_indices(a::A, inds::I) where {A,I}
    _to_indices(
        a,
        inds,
        IndexStyle(A),
        static(ndims(A)),
        eachop(_ndims_index, nstatic(Val(known_length(I))), I),
        eachop(_is_splat, nstatic(Val(known_length(I))), I)
    )
end
@generated function _to_indices(A, inds::I, ::S, ::StaticInt{N}, ::NDI, ::IS) where {I,S,N,NDI,IS}
    cnt = zeros(Int, known_length(NDI))
    splat_position = 0
    remaining = N
    for i in 1:known_length(NDI)
        ndi = known(NDI.parameters[i])
        splat = known(IS.parameters[i])
        if splat && splat_position === 0
            splat_position = i
        else
            remaining -= ndi
            cnt[i] = ndi
        end
    end
    if splat_position !== 0
        cnt[splat_position] = max(0, remaining)
    else
        # if there are additional trailing dimensions not consumed by the index then we have
        # to assume it's linear indexing or that these are trailing dimensions.
        cnt[end] += max(0, remaining)
    end

    t = Expr(:tuple)
    dim = 0
    for i in 1:known_length(NDI)
        if i === known_length(NDI) && S <: IndexLinear
            ICall = :LinearIndices
        else
            ICall = :CartesianIndices
        end
        c = cnt[i]
        iexpr = :(@inbounds(getfield(inds, $i))::$(I.parameters[i]))
        if dim === N
            push!(t.args, :(to_index($(ICall)(()), $iexpr)))
        elseif c === 1
            dim += 1
            push!(t.args, :(to_index(@inbounds(getfield(axs, $dim)), $iexpr)))
        else
            subaxs = Expr(:tuple)
            for _ in 1:c
                if dim < N
                    dim += 1
                    push!(subaxs.args, :(@inbounds(getfield(axs, $dim))))
                end
            end
            push!(t.args, :(to_index($(ICall)($subaxs), $iexpr)))
        end
    end
    Expr(:block,
        Expr(:meta, :inline),
        Expr(:(=), :axs, :(lazy_axes(A))),
        :(_flatten_tuples($t))
    )
end
@generated function _flatten_tuples(inds::I) where {I}
    t = Expr(:tuple)
    for i in 1:known_length(I)
        p = I.parameters[i]
        if p <: Tuple
            for j in 1:known_length(p)
                push!(t.args, :(@inbounds(getfield(getfield(inds, $i), $j))))
            end
        else
            push!(t.args, :(@inbounds(getfield(inds, $i))))
        end
    end
    t
end

"""
    to_index([::IndexStyle, ]axis, arg) -> index

Convert the argument `arg` that was originally passed to `getindex` for the dimension
corresponding to `axis` into a form for native indexing (`Int`, Vector{Int}, etc.). New
axis types with unique behavior should use an `IndexStyle` trait:

```julia
to_index(axis::MyAxisType, arg) = to_index(IndexStyle(axis), axis, arg)
to_index(::MyIndexStyle, axis, arg) = ...
```
"""
to_index(x, i::Slice) = i
to_index(x, ::Colon) = indices(x)
# logical indexing
to_index(x, i::AbstractArray{Bool}) = LogicalIndex(i)
to_index(x::LinearIndices, i::AbstractArray{Bool}) = LogicalIndex{Int}(i)
# cartesian indexing
@inline to_index(x, i::CartesianIndices{0}) = i
@inline to_index(x, i::CartesianIndices) = axes(i)
@inline to_index(x, i::CartesianIndex) = Tuple(i)
@inline to_index(x, i::NDIndex) = Tuple(i)
@inline to_index(x, i::AbstractArray{<:AbstractCartesianIndex}) = i
# integer indexing
to_index(x, i::AbstractArray{<:Integer}) = i
to_index(x, @nospecialize(i::StaticInt)) = i
to_index(x, i::Integer) = Int(i)
@inline to_index(x, i) = to_index(IndexStyle(x), x, i)
function to_index(S::IndexStyle, x, i)
    throw(ArgumentError(
        "invalid index: $S does not support indices of type $(typeof(i)) for instances of type $(typeof(x))."
    ))
end

"""
    unsafe_reconstruct(A, data; kwargs...)

Reconstruct `A` given the values in `data`. New methods using `unsafe_reconstruct`
should only dispatch on `A`.
"""
function unsafe_reconstruct(axis::OneTo, data; kwargs...)
    if axis === data
        return axis
    else
        return OneTo(data)
    end
end
function unsafe_reconstruct(axis::UnitRange, data; kwargs...)
    if axis === data
        return axis
    else
        return UnitRange(first(data), last(data))
    end
end
function unsafe_reconstruct(axis::OptionallyStaticUnitRange, data; kwargs...)
    if axis === data
        return axis
    else
        return OptionallyStaticUnitRange(static_first(data), static_last(data))
    end
end
function unsafe_reconstruct(A::AbstractUnitRange, data; kwargs...)
    return static_first(data):static_last(data)
end

"""
    to_axes(A, inds) -> Tuple

Construct new axes given the corresponding `inds` constructed after
`to_indices(A, args) -> inds`. This method iterates through each pair of axes and
indices calling [`to_axis`](@ref).
"""
@inline function to_axes(A, inds::Tuple)
    if ndims(A) === 1
        return (to_axis(axes(A, 1), first(inds)),)
    elseif length(inds) === 1
        return (to_axis(eachindex(IndexLinear(), A), first(inds)),)
    else
        return to_axes(A, axes(A), inds)
    end
end
# drop this dimension
to_axes(A, a::Tuple, i::Tuple{<:Integer,Vararg{Any}}) = to_axes(A, tail(a), tail(i))
to_axes(A, a::Tuple, i::Tuple{I,Vararg{Any}}) where {I} = _to_axes(ndims_index(I), A, a, i)
function _to_axes(::StaticInt{1}, A, axs::Tuple, inds::Tuple)
    return (to_axis(first(axs), first(inds)), to_axes(A, tail(axs), tail(inds))...)
end
@propagate_inbounds function _to_axes(::StaticInt{N}, A, axs::Tuple, inds::Tuple) where {N}
    axes_front, axes_tail = Base.IteratorsMD.split(axs, Val(N))
    if IndexStyle(A) === IndexLinear()
        axis = to_axis(LinearIndices(axes_front), getfield(inds, 1))
    else
        axis = to_axis(CartesianIndices(axes_front), getfield(inds, 1))
    end
    return (axis, to_axes(A, axes_tail, tail(inds))...)
end
to_axes(A, ::Tuple{Ax,Vararg{Any}}, ::Tuple{}) where {Ax} = ()
to_axes(A, ::Tuple{}, ::Tuple{}) = ()

"""
    to_axis(old_axis, index) -> new_axis

Construct an `new_axis` for a newly constructed array that corresponds to the
previously executed `to_index(old_axis, arg) -> index`. `to_axis` assumes that
`index` has already been confirmed to be in bounds. The underlying indices of
`new_axis` begins at one and extends the length of `index` (i.e., one-based indexing).
"""
@inline function to_axis(axis, inds)
    if !can_change_size(axis) &&
       (known_length(inds) !== nothing && known_length(axis) === known_length(inds))
        return axis
    else
        return to_axis(IndexStyle(axis), axis, inds)
    end
end

# don't need to check size b/c slice means it's the entire axis
@inline function to_axis(axis, inds::Slice)
    if can_change_size(axis)
        return copy(axis)
    else
        return axis
    end
end
to_axis(S::IndexLinear, axis, inds) = StaticInt(1):static_length(inds)

"""
    ArrayInterface.getindex(A, args...)

Retrieve the value(s) stored at the given key or index within a collection. Creating
another instance of `ArrayInterface.getindex` should only be done by overloading `A`.
Changing indexing based on a given argument from `args` should be done through,
[`to_index`](@ref), or [`to_axis`](@ref).
"""
function getindex(A, args...)
    inds = to_indices(A, args)
    @boundscheck checkbounds(A, inds...)
    unsafe_getindex(A, inds...)
end
function getindex(A; kwargs...)
    inds = to_indices(A, order_named_inds(dimnames(A), values(kwargs)))
    @boundscheck checkbounds(A, inds...)
    unsafe_getindex(A, inds...)
end
@propagate_inbounds getindex(x::Tuple, i::Int) = getfield(x, i)
@propagate_inbounds getindex(x::Tuple, ::StaticInt{i}) where {i} = getfield(x, i)

## unsafe_getindex ##
function unsafe_getindex(a::A) where {A}
    parent_type(A) <: A && throw(MethodError(unsafe_getindex, (A,)))
    unsafe_getindex(parent(a))
end

# TODO Need to manage index transformations between nested layers of arrays
function unsafe_getindex(a::A, i::CanonicalInt) where {A}
    if IndexStyle(A) === IndexLinear()
        parent_type(A) <: A && throw(MethodError(unsafe_getindex, (A, i)))
        return unsafe_getindex(parent(a), i)
    else
        return unsafe_getindex(a, _to_cartesian(a, i)...)
    end
end
function unsafe_getindex(a::A, i::CanonicalInt, ii::Vararg{CanonicalInt}) where {A}
    if IndexStyle(A) === IndexLinear()
        return unsafe_getindex(a, _to_linear(a, (i, ii...)))
    else
        parent_type(A) <: A && throw(MethodError(unsafe_getindex, (A, i)))
        return unsafe_getindex(parent(a), i, ii...)
    end
end

unsafe_getindex(a, i::Vararg{Any}) = unsafe_get_collection(a, i)

unsafe_getindex(A::Array) = Base.arrayref(false, A, 1)
unsafe_getindex(A::Array, i::CanonicalInt) = Base.arrayref(false, A, Int(i))
@inline function unsafe_getindex(A::Array, i::CanonicalInt, ii::Vararg{CanonicalInt})
    unsafe_getindex(A, _to_linear(A, (i, ii...)))
end

unsafe_getindex(A::LinearIndices, i::CanonicalInt) = Int(i)
unsafe_getindex(A::CartesianIndices, i::CanonicalInt, ii::Vararg{CanonicalInt}) = CartesianIndex(i, ii...)
unsafe_getindex(A::CartesianIndices, i::CanonicalInt) = @inbounds(A[i])

unsafe_getindex(A::ReshapedArray, i::CanonicalInt) = @inbounds(parent(A)[i])
function unsafe_getindex(A::ReshapedArray, i::CanonicalInt, ii::Vararg{CanonicalInt})
    @inbounds(parent(A)[_to_linear(A, (i, ii...))])
end

unsafe_getindex(A::SubArray, i::CanonicalInt) = @inbounds(A[i])
unsafe_getindex(A::SubArray, i::CanonicalInt, ii::Vararg{CanonicalInt}) = @inbounds(A[i, ii...])

# This is based on Base._unsafe_getindex from https://github.com/JuliaLang/julia/blob/c5ede45829bf8eb09f2145bfd6f089459d77b2b1/base/multidimensional.jl#L755.
#=
    unsafe_get_collection(A, inds)

Returns a collection of `A` given `inds`. `inds` is assumed to have been bounds-checked.
=#
function unsafe_get_collection(A, inds)
    axs = to_axes(A, inds)
    dest = similar(A, axs)
    if map(Base.unsafe_length, axes(dest)) == map(Base.unsafe_length, axs)
        Base._unsafe_getindex!(dest, A, inds...)
    else
        Base.throw_checksize_error(dest, axs)
    end
    return dest
end
_ints2range(x::Integer) = x:x
_ints2range(x::AbstractRange) = x
@inline function unsafe_get_collection(A::CartesianIndices{N}, inds) where {N}
    if (length(inds) === 1 && N > 1) || stride_preserving_index(typeof(inds)) === False()
        return Base._getindex(IndexStyle(A), A, inds...)
    else
        return CartesianIndices(to_axes(A, _ints2range.(inds)))
    end
end
@inline function unsafe_get_collection(A::LinearIndices{N}, inds) where {N}
    if length(inds) === 1 && isone(_ndims_index(typeof(inds), static(1)))
        return @inbounds(eachindex(A)[first(inds)])
    elseif stride_preserving_index(typeof(inds)) === True()
        return LinearIndices(to_axes(A, _ints2range.(inds)))
    else
        return Base._getindex(IndexStyle(A), A, inds...)
    end
end

"""
    ArrayInterface.setindex!(A, args...)

Store the given values at the given key or index within a collection.
"""
@propagate_inbounds function setindex!(A, val, args...)
    can_setindex(A) || error("Instance of type $(typeof(A)) are not mutable and cannot change elements after construction.")
    inds = to_indices(A, args)
    @boundscheck checkbounds(A, inds...)
    unsafe_setindex!(A, val, inds...)
end
@propagate_inbounds function setindex!(A, val; kwargs...)
    can_setindex(A) || error("Instance of type $(typeof(A)) are not mutable and cannot change elements after construction.")
    inds = to_indices(A, order_named_inds(dimnames(A), values(kwargs)))
    @boundscheck checkbounds(A, inds...)
    unsafe_setindex!(A, val, inds...)
end

## unsafe_setindex! ##
function unsafe_setindex!(a::A, v) where {A}
    parent_type(A) <: A && throw(MethodError(unsafe_setindex!, (A, v)))
    return unsafe_setindex!(parent(a), v)
end
# TODO Need to manage index transformations between nested layers of arrays
function unsafe_setindex!(a::A, v, i::CanonicalInt) where {A}
    if IndexStyle(A) === IndexLinear()
        parent_type(A) <: A && throw(MethodError(unsafe_setindex!, (A, v, i)))
        return unsafe_setindex!(parent(a), v, i)
    else
        return unsafe_setindex!(a, v, _to_cartesian(a, i)...)
    end
end
function unsafe_setindex!(a::A, v, i::CanonicalInt, ii::Vararg{CanonicalInt}) where {A}
    if IndexStyle(A) === IndexLinear()
        return unsafe_setindex!(a, v, _to_linear(a, (i, ii...)))
    else
        parent_type(A) <: A && throw(MethodError(unsafe_setindex!, (A, v, i, ii...)))
        return unsafe_setindex!(parent(a), v, i, ii...)
    end
end

function unsafe_setindex!(A::Array{T}, v) where {T}
    Base.arrayset(false, A, convert(T, v)::T, 1)
end
function unsafe_setindex!(A::Array{T}, v, i::CanonicalInt) where {T}
    return Base.arrayset(false, A, convert(T, v)::T, Int(i))
end

unsafe_setindex!(a, v, i::Vararg{Any}) = unsafe_set_collection!(a, v, i)

# This is based on Base._unsafe_setindex!.
#=
    unsafe_set_collection!(A, val, inds)

Sets `inds` of `A` to `val`. `inds` is assumed to have been bounds-checked.
=#
unsafe_set_collection!(A, v, i) = Base._unsafe_setindex!(IndexStyle(A), A, v, i...)


text:
struct Plsr
    T::Matrix{Float64}
    P::Matrix{Float64}
    R::Matrix{Float64}
    W::Matrix{Float64}
    C::Matrix{Float64}
    TT::Vector{Float64}
    xmeans::Vector{Float64}
    ymeans::Vector{Float64}
    weights::Vector{Float64}
    ## For consistency with plsrannar
    U::Union{Array{Float64}, Nothing}
end

struct Pcr
    fm_pca
    T::Matrix{Float64}
    R::Matrix{Float64}
    C::Matrix{Float64}
    xmeans::Vector{Float64}
    ymeans::Vector{Float64}
    weights::Vector{Float64}
end

"""
    plskern(X, Y, weights = ones(size(X, 1)); nlv)
Partial Least Squares Regression (PLSR) with the
"Improved kernel algorithm #1" (Dayal & McGegor, 1997).
* `X` : X-data.
* `Y` : Y-data.
* `weights` : Weights of the observations.
* `nlv` : Nb. latent variables (LVs) to compute.

For the row-weighting in PLS algorithms (`weights`), see in particular Schaal et al. 2002,
Siccard & Sabatier 2006, Kim et al. 2011 and Lesnoff et al. 2020.

Vector `weights` is internally normalized to sum to 1.

`X` and `Y` are internally centered. The model is computed with an intercept.

## References

Dayal, B.S., MacGregor, J.F., 1997. Improved PLS algorithms.
Journal of Chemometrics 11, 73-85.

Kim, S., Kano, M., Nakagawa, H., Hasebe, S., 2011. Estimation of active
pharmaceutical ingredients content using locally weighted partial
least squares and statistical wavelength selection. Int. J. Pharm., 421, 269-274.

Lesnoff, M., Metz, M., Roger, J.M., 2020. Comparison of locally weighted
PLS strategies for regression and discrimination on agronomic NIR Data.
Journal of Chemometrics. e3209.
https://onlinelibrary.wiley.com/doi/abs/10.1002/cem.3209

Schaal, S., Atkeson, C., Vijayamakumar, S. 2002. Scalable techniques
from nonparametric statistics for the real time robot learning.
Applied Intell., 17, 49-60.

Sicard, E. Sabatier, R., 2006. Theoretical framework for local PLS1 regression
and application to a rainfall data set. Comput. Stat. Data Anal., 51, 1393-1410.
"""
function plskern(X, Y, weights = ones(size(X, 1)); nlv)
    plskern!(copy(X), copy(Y), weights; nlv = nlv)
end

function plskern!(X, Y, weights = ones(size(X, 1)); nlv)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    n = size(X, 1)
    p = size(X, 2)
    q = size(Y, 2)
    nlv = min(nlv, n, p)
    weights = mweights(weights)
    xmeans = colmeans(X, weights)
    ymeans = colmeans(Y, weights)
    center!(X, xmeans)
    center!(Y, ymeans)
    D = Diagonal(weights)
    XtY = X' * (D * Y)                   # = Xd' * Y = X' * D * Y  (Xd = D * X   Very costly!!)
    # Pre-allocation
    T = similar(X, n, nlv)
    P = similar(X, p, nlv)
    W = copy(P)
    R = copy(P)
    C = similar(X, q, nlv)
    TT = similar(X, nlv)
    t   = similar(X, n)
    dt  = similar(X, n)
    zp  = similar(X, p)
    w   = similar(X, p)
    r   = similar(X, p)
    c   = similar(X, q)
    tmp = similar(XtY)
    # End
    @inbounds for a = 1:nlv
        if q == 1
            w .= vcol(XtY, 1)
        else
            tmp .= XtY
            u = svd!(tmp).V           # = svd(XtY').U = eigen(XtY' * XtY).vectors[with ordering]
            mul!(w, XtY, vcol(u, 1))
        end
        w ./= sqrt(dot(w, w))         # w .= w ./ sqrt(dot(w, w))
        r .= w
        if a > 1
            @inbounds for j = 1:(a - 1)
                r .-= dot(w, vcol(P, j)) .* vcol(R, j)    # r .= r .- ...
            end
        end
        mul!(t, X, r)                 # t = X * r
        dt .= weights .* t            # dt = D * t
        tt = dot(t, dt)               # tt = t' * dt = t' * D * t
        mul!(c, XtY', r)
        c ./= tt                      # c = XtY' * r / tt
        mul!(zp, X', dt)              # zp = (D * X)' * t = X' * (D * t)
        XtY .-= mul!(tmp, zp, c')     # XtY = XtY - zp * c' ; deflation of the kernel matrix
        P[:, a] .= zp ./ tt           # (cont') ==> the metric applied to covariance is applied outside the loop,
        T[:, a] .= t                  # (cont') conversely to other algorithms such as nipals
        W[:, a] .= w
        R[:, a] .= r
        C[:, a] .= c
        TT[a] = tt
     end
     Plsr(T, P, R, W, C, TT, xmeans, ymeans, weights, nothing)
end

"""
    summary(object::Plsr, X)
Summarize the maximal (i.e. with maximal nb. LVs) fitted model.
* `object` : The fitted model.
* `X` : The X-data that was used to fit the model.
"""
function Base.summary(object::Plsr, X)
    n, nlv = size(object.T)
    X = center(X, object.xmeans)
    # Could be center! but changes x
    # If too heavy ==> Makes summary!
    sstot = sum(object.weights' * (X.^2))
    tt = object.TT
    tt_adj = vec(sum(object.P.^2, dims = 1)) .* tt
    pvar = tt_adj / sstot
    cumpvar = cumsum(pvar)
    xvar = tt_adj / n
    explvar = DataFrame(nlv = 1:nlv, var = xvar, pvar = pvar, cumpvar = cumpvar)
    (explvar = explvar,)
end

"""
    transform(object::Plsr, X; nlv = nothing)
Compute LVs ("scores" T) from a fitted model and a matrix X.
* `object` : The maximal fitted model.
* `X` : Matrix (m, p) for which LVs are computed.
* `nlv` : Nb. LVs to consider. If nothing, it is the maximum nb. LVs.
"""
function transform(object::Plsr, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    T = center(X, object.xmeans) * vcol(object.R, 1:nlv)
    # Could be center! but changes x
    # If too heavy ==> Makes summary!
    T
end

"""
    coef(object::Union{Plsr, Pcr}; nlv = nothing)
Compute the X b-coefficients of a model fitted with `nlv` LVs.
* `object` : The maximal fitted model.
* `nlv` : Nb. LVs to consider. If nothing, it is the maximum nb. LVs.

If X is (n, p) and Y is (n, q), the returned object `B` is a matrix (p, q).
If `nlv` = 0, `B` is a matrix of zeros.
The returned object `int` is the intercept.
"""
function coef(object::Union{Plsr, Pcr}; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    beta = object.C[:, 1:nlv]'
    B = vcol(object.R, 1:nlv) * beta
    int = object.ymeans' .- object.xmeans' * B
    (B = B, int = int)
end

"""
    predict(object::Union{Plsr, Pcr}, X; nlv = nothing)
Compute Y-predictions from a fitted model.
* `object` : The fitted model.
* `X` : X-data for which predictions are computed.
* `nlv` : Nb. LVs, or collection of nb. LVs, to consider.
    If nothing, it is the maximum nb. LVs.
"""
function predict(object::Union{Plsr, Pcr}, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = (max(minimum(nlv), 0):min(maximum(nlv), a))
    le_nlv = length(nlv)
    pred = list(le_nlv, Matrix{Float64})
    @inbounds for i = 1:le_nlv
        z = coef(object; nlv = nlv[i])
        pred[i] = z.int .+ X * z.B
    end
    le_nlv == 1 ? pred = pred[1] : nothing
    (pred = pred,)
end


text:
using EconPDEs

struct WangWangYangModel
    μ::Float64
    σ::Float64
    r::Float64
    ρ::Float64
    γ::Float64
    ψ::Float64
    wmax::Float64
end

function WangWangYangModel(;μ = 0.015, σ = 0.1, r = 0.035, ρ = 0.04, γ = 3, ψ = 1.1, wmax = 5000.0)
    WangWangYangModel(μ, σ, r, ρ, γ, ψ, wmax)
end

function initialize_stategrid(m::WangWangYangModel; n = 500)
    OrderedDict(:w => range(0.0, stop = sqrt(m.wmax), length = n).^2)
end

function initialize_y(m::WangWangYangModel, stategrid)
    OrderedDict(:p => 1 .+ stategrid[:w])
end

function (m::WangWangYangModel)(state::NamedTuple, y::NamedTuple)
    μ = m.μ ;  σ = m.σ ;  r = m.r ;  ρ = m.ρ ;  γ = m.γ ;  ψ = m.ψ ; wmax = m.wmax
    w = state.w
    p, pw, pww = y.p, y.pw, y.pww
    c = (r + ψ * (ρ - r)) * p * pw^(-ψ)
    μw = (r - μ + σ^2) * w + 1 - c
    # One only needs a ghost node if μw <= 0 (since w^2p_ww = 0). In this case, we obtain a formula for pw so that c <= 1
    if w ≈ 0.0 && μw <= 0.0
       pw = ((r + ψ * (ρ - r)) * p)^(1 / ψ)
       c = (r + ψ * (ρ - r)) * p * pw^(-ψ)
       μw = (r - μ + σ^2) * w + 1 - c
    end
    # At the top, I use the solution of the unconstrainted, i.e. pw = 1 (I could also do reflecting boundary but less elegant)
    pt = (((r + ψ * (ρ - r)) * pw^(1 - ψ) - ψ * ρ) / (ψ - 1) + μ - γ * σ^2 / 2) * p + ((r - μ + γ * σ^2) * w + 1) * pw + σ^2 * w^2 / 2  * (pww - γ * pw^2 / p)
    μw = 1 + (r - μ + σ^2) * w - c
    return (pt,), (μw,), (w = w, p = p, pw = pw, pww = pww, μw = μw, c = c)
end

m = WangWangYangModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
y, result, distance = pdesolve(m, stategrid, y0; bc = OrderedDict(:pw => (1.0, 1.0)))


text:
@testset "Shape" begin
    import NaiveGAflux: ΔShape, ShapeAdd, ShapeMul, ShapeDiv, ShapeFlatten, fshape, revert, combine, combine, filter_noops, ShapeTrace, shapetrace, Δshapes, squashshapes, orderΔshapes, ndimsout, ndimsin

    @testset "ΔShapes" begin

        @test fshape(tuple(), (1,2,3,4)) == (1,2,3,4)

        @testset "ShapeAdd" begin
            @test fshape(ShapeAdd(1,2), (3,4)) == (4,6)
            @test fshape(combine(ShapeAdd(1,2), ShapeAdd(-2,-1)), (3,4)) == (2,5)
            @test combine(ShapeAdd(1,2), ShapeAdd(3,4)) == tuple(ShapeAdd(4,6))
            @test revert(ShapeAdd(1,2)) == ShapeAdd(-1,-2)
            @test filter_noops(ShapeAdd(0,1,2,3)) == tuple(ShapeAdd(0,1,2,3))
            @test filter_noops(ShapeAdd(0,0)) == tuple()
        end

        @testset "ShapeMul" begin
            @test fshape(ShapeMul(1,2), (3,4)) == (3,8)
            @test fshape(combine(ShapeMul(1,2), ShapeMul(2,3)), (3,4)) == (6,24)
            @test combine(ShapeMul(1,2), ShapeMul(3,4)) == tuple(ShapeMul(3,8))
            @test revert(ShapeMul(1,2)) == ShapeDiv(1,2)
            @test filter_noops(ShapeMul(0,1,2,3)) == tuple(ShapeMul(0,1,2,3))
            @test filter_noops(ShapeMul(1,1)) == tuple()
        end

        @testset "ShapeDiv" begin
            @test fshape(ShapeDiv(1,2), (3,4)) == (3,2)
            @test fshape(combine(ShapeDiv(1,2), ShapeDiv(2,3)), (3,4)) == (2,1)
            @test combine(ShapeDiv(1,2), ShapeDiv(3,4)) == tuple(ShapeDiv(3,8))
            @test revert(ShapeDiv(1,2)) == ShapeMul(1,2)
            @test filter_noops(ShapeDiv(0,1,2,3)) == tuple(ShapeDiv(0,1,2,3))
            @test filter_noops(ShapeDiv(1,1)) == tuple()
        end

        @testset "ShapeFlatten" begin
            @test fshape(ShapeFlatten(), (1,2)) == fshape(ShapeFlatten(), (1,2,3)) == tuple()
            @test ShapeFlatten() |> revert |> revert == ShapeFlatten()
            @test ShapeFlatten{2}() |> revert |> revert == ShapeFlatten{2}()
            @test combine(ShapeFlatten{2}(), ShapeFlatten{2}()) == tuple(ShapeFlatten{2}(), ShapeFlatten{2}())
            @test filter_noops(ShapeFlatten()) == tuple(ShapeFlatten())
        end

        @testset "ndims" begin
            @test ndimsin(ShapeDiv(1,2,3)) == ndimsout(ShapeDiv(1,2,3)) == 3
            @test ndimsin((ShapeAdd(2,2), ShapeDiv(1,2,3))) == 2
            @test ndimsout((ShapeAdd(2,2), ShapeDiv(1,2,3))) == 3
        end

        @testset "Combine $s1 and $s2" for (s1,s2,exp) in (
            (ShapeDiv(2,2), ShapeMul(2,2), (ShapeDiv(2,2), ShapeMul(2,2))),
            (ShapeMul(2,2), ShapeDiv(2,2), tuple(ShapeMul(1,1))),
            (ShapeDiv(3,3), ShapeMul(2,2), (ShapeDiv(3,3), ShapeMul(2,2))),
            (ShapeMul(3,3), ShapeDiv(2,2), (ShapeMul(3,3), ShapeDiv(2,2))),
            (ShapeDiv(2,2), ShapeMul(3,3), (ShapeDiv(2,2), ShapeMul(3,3))),
            (ShapeMul(2,2), ShapeDiv(3,3), (ShapeMul(2,2), ShapeDiv(3,3))),
            (ShapeDiv(4,4), ShapeMul(2,2), (ShapeDiv(4,4), ShapeMul(2,2))),
            (ShapeMul(4,4), ShapeDiv(2,2), tuple(ShapeMul(2,2))),
            (ShapeDiv(2,2), ShapeMul(4,4), tuple(ShapeDiv(2,2), ShapeMul(4,4))),
            (ShapeMul(2,2), ShapeDiv(4,4), tuple(ShapeMul(2,2), ShapeDiv(4,4)))
            )
            act = combine(s1,s2)
            @test act == exp
            @testset "Same fshape with insize $insize" for insize in 1:5:100
                is = (insize,insize+1)
                @test fshape((s1,s2), is) == fshape(act, is) == fshape(exp, is)
            end
        end

        @testset "Swap shapes" begin
            import NaiveGAflux:  swapΔshape
            @test swapΔshape(ShapeAdd(1,2), ShapeAdd(2,3)) == (ShapeAdd(2,3), ShapeAdd(1,2))

            sasm = ShapeAdd(1,2), ShapeMul(2,3)
            @test swapΔshape(sasm...) == (ShapeMul(2,3), ShapeAdd(2,6))
            @test fshape(swapΔshape(sasm...), (11,12)) == fshape(sasm, (11,12)) == (24,42)
            @testset "swap $(typeof.(sasm)) size $s" for s in 1:10
                @test fshape(swapΔshape(sasm...), (s,s)) == fshape(sasm, (s,s))
            end

            @test swapΔshape(ShapeMul(2,3), ShapeAdd(1,2)) == (ShapeMul(2,3), ShapeAdd(1,2))
            @test swapΔshape(ShapeMul(2,3), ShapeAdd(4,12)) == (ShapeAdd(2, 4), ShapeMul(2,3))

            sdsa = ShapeDiv(2,3), ShapeAdd(1,2)
            @test swapΔshape(sdsa...) == (ShapeAdd(2,6), ShapeDiv(2,3))
            @test fshape(swapΔshape(sdsa...), (13,14)) == fshape(sdsa, (13,14)) == (8,7)
            @testset "swap $(typeof.(sdsa)) size $s" for s in 1:10
                @test fshape(swapΔshape(sdsa...), (s,s)) == fshape(sdsa, (s,s))
            end

            @test swapΔshape(ShapeAdd(1,2), ShapeDiv(2,3)) == (ShapeAdd(1,2), ShapeDiv(2,3))
            @test swapΔshape(ShapeAdd(4,12), ShapeDiv(2,3)) == (ShapeDiv(2,3), ShapeAdd(2, 4))

            @test swapΔshape(ShapeMul(2,4), ShapeDiv(4,8)) == (ShapeMul(2,4), ShapeDiv(4,8))

            @test swapΔshape(ShapeAdd(2,4), ShapeFlatten()) == (ShapeAdd(2,4), ShapeFlatten())
        end

        @testset "Order shapes" begin
            @test orderΔshapes(tuple(ShapeAdd(1,2))) == tuple(ShapeAdd(1,2))
            @test orderΔshapes((ShapeAdd(1,2), ShapeMul(2,3), ShapeAdd(4,6))) == (ShapeAdd(2,2), ShapeAdd(1,2), ShapeMul(2,3))
            @test orderΔshapes((ShapeAdd(1,2), ShapeMul(2,3), ShapeAdd(4,6), ShapeFlatten())) == (ShapeAdd(1,2), ShapeAdd(2,2), ShapeMul(2,3), ShapeFlatten())

            @testset "Mixed ShapeAdds" begin
                s = (ShapeAdd(1,2), ShapeDiv(2,3), ShapeMul(2,2), ShapeAdd(4,4), ShapeAdd(3,5), ShapeAdd(6,8), ShapeMul(2,3), ShapeMul(1,1), ShapeAdd(3,7))
                os =  orderΔshapes(s)
                @test os == (ShapeAdd{2}((6, 12)), ShapeAdd{2}((1, 2)), ShapeAdd{2}((4, 6)), ShapeDiv{2}((2, 3)), ShapeMul{2}((2, 2)), ShapeAdd{2}((3, 5)), ShapeMul{2}((2, 3)), ShapeMul{2}((1, 1)), ShapeAdd{2}((3, 7)))

                @testset "same shape with input size $insize" for insize in 1:10
                    @test fshape(s, (insize,insize)) == fshape(os, (insize,insize))
                end
            end

            @testset "Single ShapeAdd with ShapeMul and ShapeDiv" begin
                s = (ShapeMul(2,2), ShapeDiv(12,12), ShapeAdd(2,3), ShapeMul(2,2), ShapeMul(3,3))
                os = orderΔshapes(s)
                @test os == (ShapeMul(2,2), ShapeDiv(12,12), ShapeMul(3,3), ShapeMul(2,2), ShapeAdd(12,18))
                @testset "same shape with input size $insize" for insize in 1:5:100
                    @test fshape(s, (insize,insize)) == fshape(os, (insize,insize))
                end
            end
        end

        @testset "squash shapes" begin
            @test squashshapes(ShapeAdd(1,2)) == tuple(ShapeAdd(1,2))
            @test squashshapes(ShapeAdd(1,2), ShapeAdd(3,4)) == tuple(ShapeAdd(4,6))
            @test squashshapes(ShapeAdd(1,2), ShapeAdd(3,4), ShapeAdd(5,6)) == tuple(ShapeAdd(9,12))
            @test squashshapes(ShapeAdd(1,2), ShapeAdd(3,4), ShapeAdd(-4,-6)) == tuple()
            @test squashshapes(ShapeAdd(1,2), ShapeMul(3,4)) == (ShapeAdd(1,2), ShapeMul(3,4))
            @test squashshapes(ShapeAdd(1,2), ShapeAdd(-1,1), ShapeFlatten()) == (ShapeAdd(0, 3), ShapeFlatten())

            as = (ShapeAdd(1,2), ShapeMul(3,4))
            sa = revert(as)
            @test squashshapes((as..., sa...)) == tuple()
            @test squashshapes(ShapeAdd(4,2), ShapeDiv(2,2), ShapeMul(2,2), ShapeAdd(-4,-2)) == (ShapeDiv(2,2), ShapeMul(2,2))

            @testset "Squash mix add first" begin
                s = (ShapeAdd(2,3), ShapeMul(2,2), ShapeMul(2,2), ShapeMul(3,3), ShapeDiv(12,12))
                sq = squashshapes(s)
                @test sq == tuple(ShapeAdd(2,3))
                @testset "squashed shape with input size $insize" for insize in 1:5:100
                    @test fshape(s, (insize, insize)) == fshape(sq, (insize,insize))
                end
            end

            @testset "Squash mix add last" begin
                s = (ShapeMul(2,2), ShapeMul(2,2), ShapeMul(3,3), ShapeDiv(12,12), ShapeAdd(2,3))
                sq = squashshapes(s)
                @test sq == tuple(ShapeAdd(2,3))
                @testset "squashed shape with input size $insize" for insize in 1:5:100
                    @test fshape(s, (insize, insize)) == fshape(sq, (insize,insize))
                end
            end

            @testset "Squash mix add mid" begin
                s = (ShapeMul(2,2), ShapeMul(2,2), ShapeAdd(2,3), ShapeMul(3,3), ShapeDiv(12,12))
                sq = squashshapes(s)
                @test sq == (ShapeMul(12, 12), ShapeAdd(6, 9), ShapeDiv(12, 12))
                @testset "squashed shape with input size $insize" for insize in 1:5:100
                    @test fshape(s, (insize, insize)) == fshape(sq, (insize,insize))
                end
            end

            @testset "Squash mix add and div mid" begin
                s = (ShapeMul(2,2), ShapeMul(2,2), ShapeDiv(12,12), ShapeAdd(2,3), ShapeMul(3,3))
                sq = squashshapes(s)
                @test sq == (ShapeMul(4, 4), ShapeDiv(12,12), ShapeMul(3,3), ShapeAdd(6, 9))
                @testset "squashed shape with input size $insize" for insize in 1:5:100
                    @test fshape(s, (insize, insize)) == fshape(sq, (insize,insize))
                end
            end
        end

        @testset "Diff shapes" begin
            import NaiveGAflux: Δshapediff

            @testset "Primitives" begin
                @test Δshapediff(ShapeAdd(2,2), ShapeAdd(2,2)) == tuple()
                @test Δshapediff(ShapeAdd(2,3), ShapeAdd(3,4)) == tuple(ShapeAdd(-1,-1))

                @test Δshapediff(ShapeMul(2,2), ShapeMul(2,2)) == tuple()
                @test Δshapediff(ShapeMul(2,3), ShapeMul(3,4)) == (ShapeDiv(3,4), ShapeMul(2,3))

                @test Δshapediff(ShapeDiv(2,2), ShapeDiv(2,2)) == tuple()
                @test Δshapediff(ShapeDiv(2,3), ShapeDiv(3,4)) == (ShapeMul(3,4), ShapeDiv(2,3))
                @test Δshapediff(ShapeDiv(2,3), ShapeDiv(4,6)) == tuple(ShapeMul(2,2))

                @test Δshapediff(ShapeFlatten(), ShapeFlatten()) == tuple()
                @test Δshapediff(ShapeFlatten{2}(), ShapeFlatten{3}()) == tuple(ShapeFlatten{2}(), ShapeFlatten{3}())
            end

            @testset "Tuples" begin
                @testset "Basic" begin
                    s1 = (ShapeAdd(2,2), ShapeDiv(2,2), ShapeMul(4,2))
                    @test Δshapediff(s1, s1) == tuple()

                    @test Δshapediff(s1, s1[1:2]) == tuple(ShapeMul(4,2))
                    @test Δshapediff(s1, s1[2:3]) == tuple(ShapeAdd(4,2))
                    @test Δshapediff(s1, s1[[1,3]]) == (ShapeDiv(8, 4), ShapeMul(4, 2))

                    @test Δshapediff(s1[1:2], s1) == tuple(ShapeDiv(4,2))
                    @test Δshapediff(s1[2:3], s1) == tuple(ShapeAdd(-4,-2))
                    @test Δshapediff(s1[[1,3]], s1) == (ShapeDiv(4, 2), ShapeMul(8, 4))

                    s2 = (ShapeDiv(4,4), ShapeAdd(3,5), ShapeMul(8,4), ShapeDiv(2,2))
                    @test Δshapediff(s1,s2) == (ShapeDiv(4, 2), ShapeMul(8, 4), ShapeAdd(-20, -18))
                end
            end
        end
    end

    @testset "ShapeTrace" begin
        import NaiveGAflux: GlobalPool
        import NaiveNASflux: named
        iv(N=2) = convinputvertex("in", 1, N)
        bv(in, name) = fluxvertex(name, BatchNorm(nout(in)), in)
        pv(in, name; ks=(3,3), stride=ntuple(i->1, length(ks)), kwargs...) = fluxvertex(name, MaxPool(ks; stride=stride, kwargs...), in)
        cv(in, name; ks=(3,3), kwargs...) = fluxvertex(name, Conv(ks, nout(in) => nout(in); kwargs...), in)
        gp(in, name) = invariantvertex(GlobalPool{MaxPool}(), in, traitdecoration=named(name))

        @testset "Trace $vf" for vf in (cv,pv)
            vi = iv()
            v1 = vf(vi, "v1"; ks=(1,1))
            v2 = vf(v1, "v2"; ks=(3,3))
            v3 = vf(v2, "v3"; ks=(2,4), pad=(5,7))
            v4 = vf(v3, "v4"; ks=(3,3), stride=(1, 2))
            v5 = vf(v4, "v5"; ks=(4,5), pad=(2,3,4,5), stride=(3,4))

            @test filter_noops(Δshapes(v1)) == tuple()
            @testset "shape for $(name(vn))" for vn in (v2,v3,v4,v5)
                s = ShapeTrace(vn).trace
                @test fshape(s, (10,9)) == size(vn(ones(Float32, 10,9, nout(vi), 1)))[1:2]
            end

            @testset "shapetrace graph" begin
                g = CompGraph(vi, v5)
                sg = g(ShapeTrace(vi)).trace
                sv = shapetrace(v5; trfun = v -> ShapeTrace(v)).trace
                @test fshape(sg, (30,31)) == fshape(sv, (30,31))== size(g(ones(Float32, 30,31, nout(vi), 1)))[1:2]
            end
        end

        @testset "Trace $vf merge" for vf in (cv, pv)
            vi = iv()
            v1 = vf(vi, "v1"; ks=(2,5))
            va1 = vf(v1, "va1"; ks=(1,1), stride=(2,2))
            va2 = vf(va1, "va2"; ks=(2,2))
            vb1 = vf(v1, "vb1"; ks=(3,3))
            vb2 = vf(vb1, "vb2";ks=(1,1), stride=(2,2))
            mv = concat("concat", va2,vb2)
            v2 = vf(mv, "v2"; ks=(3,3))

            tr = shapetrace(v2)

            @test fshape(squashshapes(tr), (13,17)) == size(CompGraph(vi,v2)(ones(Float32, 13,17,1,1)))[1:2]
        end

        @testset "Conv dilation" begin
            vi = iv()
            v1 = cv(vi, "v1"; ks=(1,1), dilation=2)
            v2 = cv(v1, "v2"; ks=(2,2), dilation=2)
            v3 = cv(v2, "v3"; ks=(3,4), dilation=(4,5), pad=(1,2,3,4))
            v4 = cv(v3, "v4"; ks=(1,2), dilation=(2,3), pad=(4,5), stride=(6,7))

            @test filter_noops(Δshapes(v1)) == tuple()
            @testset "shape for $(name(vn))" for vn in (v2,v3,v4)
                s = ShapeTrace(vn).trace
                @test fshape(s, (20,19)) == size(vn(ones(Float32, 20,19, nout(vi), 1)))[1:2]
            end

            @testset "shapetrace graph" begin
                g = CompGraph(vi, v4)
                sg = g(ShapeTrace(vi)).trace
                sv = shapetrace(v4; trfun = v -> ShapeTrace(v)).trace
                @test fshape(sg, (30,31)) == fshape(sv, (30,31))== size(g(ones(Float32, 30,31, nout(vi), 1)))[1:2]
            end
        end

        @testset "1D  $vf" for vf in (cv,pv)
            vi = iv(1)
            v1 = vf(vi, "v1"; ks=(1,), pad=1)
            v2 = vf(v1, "v2", ks=(3,), stride=(2,), pad=(1,1))

            @testset "shape for $(name(vn))" for vn in (v1, v2)
                s = ShapeTrace(vn).trace
                @test fshape(s, (7,)) == size(vn(ones(Float32, 7, nout(vi), 1)))[1:1]
            end
        end

        @testset "3D  $vf" for vf in (cv,pv)
            vi = iv(3)
            v1 = vf(vi, "v1"; ks=(1,2,3), pad=1)
            v2 = vf(v1, "v2", ks=(3,3,3), stride=(2,1,3), pad=(1,2,3,4,5,6))

            @testset "shape for $(name(vn))" for vn in (v1, v2)
                s = ShapeTrace(vn).trace
                @test fshape(s, (17,11,20)) == size(vn(ones(Float32, 17, 11, 20, nout(vi), 1)))[1:3]
            end
        end

        @testset "GlobalPool" begin
            v1 = gp(iv(3), "v1")
            @test Δshapes(v1) == tuple(ShapeFlatten{2}())

            s = ShapeTrace(v1).trace
            @test fshape(s, (11,12)) == tuple()
        end

        @testset "Flux.flatten" begin
            v1 = invariantvertex(Flux.flatten, iv(3); traitdecoration = named("v1"))
            @test Δshapes(v1) == tuple(ShapeFlatten{2}())

            s = ShapeTrace(v1).trace
            @test fshape(s, (11,12)) == tuple()
        end

        @testset "Squash misaligned branches" begin
            vi = iv()
            v1 = pv(vi, "v1"; ks=(2,2))
            va1 = pv(v1, "va1"; ks=(1,1), stride=(2,2))
            va2 = pv(va1, "va2"; ks=(3,3))
            vb1 = pv(v1, "vb1"; ks=(3,3), pad=(2,2))
            vb2 = pv(vb1, "vb2"; ks=(9,9))
            vb3 = pv(vb2, "vb3"; ks=(3,3))
            av = "add" >> va2 + vb3
            v2 = pv(av, "v2"; ks=(2,3))

            g = CompGraph(vi ,v2)
            @test size(g(ones(Float32, 13,14, nout(vi), 1)))[1:2] == (3,3)

            trv2 = shapetrace(v2)
            sv2 = squashshapes(trv2)

            trva2 = shapetrace(va2)
            sva2 = squashshapes(trva2)

            trvb3 = shapetrace(vb3)
            svb3 = squashshapes(trvb3)

            @test sv2[1][1] == sva2
            @test sv2[1][2] == svb3
            @test sv2[2] == filter_noops(Δshapes(v2))

            @test fshape((sva2...,sv2[2]...), (13,14)) == (3,3)
            @test fshape((svb3...,sv2[2]...), (13,14)) == (3,3)

            @test fshape((sva2...,sv2[2]...), (16,17)) == (5,4)
            @test fshape((svb3...,sv2[2]...), (16,17)) == (6,6)
        end

        @testset "Squash nested branches" begin

            # The fact that I could use the shapetrace to align this monstrosity is perhaps a sign that it is not completely useless...
            vi = iv()
            v1 = pv(vi, "v1"; ks=(3,3))

            va1 = pv(v1, "va1"; ks=(2,2), stride=(2,2))
            va2 = pv(va1, "va2"; ks=(5,5))

            vaa1 = pv(va2, "vaa1"; ks=(1,1), stride=(3,3))
            vaa2 = pv(vaa1, "vaa2"; ks=(3,3))
            vab1 = pv(va2, "vab1"; ks=(7,7), stride=(3,3))

            mva1 = "mva1" >> vaa2 + vab1
            va3 = pv(mva1, "va3"; ks=(5,5))

            vb1 = pv(v1, "vb1"; ks=(4,4), stride=(3,3))
            vb2 = pv(vb1, "vb2"; ks=(3,3), stride=(2,2))
            vb3 = pv(vb2, "vb3"; ks=(2,2))
            vb4 = pv(vb3, "vb4"; ks=(6,6))

            mv1 = "mv1" >> va3 + vb4
            v2 = pv(mv1, "v2"; ks=(3,3))

            @test squashshapes(shapetrace(v2)) == squashshapes(shapetrace(v2, v1)) == (ShapeAdd(-59, -59), ShapeDiv(6, 6))
        end

        @testset "Squash with different start types" begin
            import NaiveGAflux: allΔshapetypes
            vi = iv()

            va1 = pv(vi, "va1"; ks=(3,3), stride=(2,2), pad=(1,1))
            va2 = pv(va1, "va2"; ks = (3,3))

            vb1 = pv(vi, "vb1"; ks=(5,5))
            vb2 = pv(vb1,"vb1"; ks=(3,3), stride=(2,2), pad=(1,1))
            mv = "add" >> va2 + vb2

            st = shapetrace(mv);
            # All vertices insert ΔShapes in the same order, so we need to do this to make them have different order
            st = NaiveGAflux.@set st.trace[1][1].trace = filter_noops(st.trace[1][1].trace)

            @test allΔshapetypes(st) == [ShapeDiv{2}, ShapeAdd{2}]

            @test squashshapes(st) == (ShapeDiv(2, 2), ShapeAdd(-2,-2))
        end
    end
end


text:
module CHIP8

using OffsetArrays

include("sprites.jl")           # Defines sprites for all hex digits
include("utils.jl")
include("keyboard.jl")
include("tests.jl")

export Chip8, parseopcode

const memsize = 4096            # 4KiB of memory
const program_start = 0x200 + 1 # Starting address of programs in memory (+1 because of Julia 1 indexing ...)
const disp_columns = 64         # Number of horizontal pixels
const disp_rows = 32            # Number of vertical pixels

"Holds the state of a CHIP8 system"
mutable struct Chip8
  # Registers
  V::AbstractVector{UInt8}      # General purpose registers
  I::UInt16                     # Address register
  DT::UInt8                     # Timer register
  ST::UInt8                     # Sound register
  PC::UInt16                    # Program counter
  # Stack pointer is only needed if the stack is managed manually
  # PC::UInt16                  # Stack pointer

  # Memory
  mem::AbstractVector{UInt8}    # Main memory
  disp::AbstractMatrix{Bool}    # Display buffer

  # Stack
  stack::AbstractVector{UInt16}

  # Keyboard
  keys::AbstractVector{Bool}    # Models which of the 16 keys are pressed at the moment
end

function Chip8()
  ch8 = Chip8(OffsetVector(zeros(UInt8, 16), 0:15),                # General purpose registers
              0, 0, 0,                                             # Address, timer, sound register
              program_start,                                       # Program counter
              OffsetVector(zeros(UInt8, memsize), 0:memsize-1),    # Main memory
              OffsetMatrix(zeros(Bool, 64, 32), 0:63, 0:32),       # Display buffer
              Vector{UInt16}(),                                    # Stack
              OffsetVector(zeros(Bool, 16), 0:15),                 # Keys
              )
  # Write predefined digit sprites to memory (starting at 0x100)
  ch8.mem[0x101:0x101+length(sprites)-1] = sprites

  return ch8
end

include("graphics.jl")

"""
`decode(b::Integer) = reverse(digits(b, base = 2, pad = sizeof(b)))`

Returns the digits of b as an array, where the most significant bit has the lowest index.
Mainly used for decoding lines of sprites saved in memory as bytes.
"""
decode(b::Integer) = reverse(digits(b, base = 2, pad = sizeof(b)))

# Define our instructions

struct Instruction
  asm::String
  opcode::String
  description::String
  f::Function
end

#TODO: Find a good solution for the 1-index problem in Julia ...

instructions = [
  Instruction("CLS", "00E0", "Clear Screen",
              function (c)
              c.disp .= false
              end)
  Instruction("RET", "00EE", "Return from subroutine",
              function (c)
              c.PC = pop(c.stack)
              end)
  Instruction("SYS addr", "0nnn", "Call system subroutine",
              function (c, nnn)
              error("Tried to call system subroutine on address $nnn. Not possible in emulation!")
              end)
  Instruction("JP addr", "1nnn", "Jump to address",
              function (c, nnn)
              c.PC = nnn - 2    # -2 because the address is going to be advanced by 2 bytes later
              end)
  Instruction("CALL addr", "2nnn", "Call subroutine at address",
              function (c, nnn)
              push!(c.stack, c.PC)
              c.PC = nnn - 2
              end)
  Instruction("SE Vx, val", "3xkk", "Jump over next instruction if (Vx == val,)",
              function (c, x, kk)
              (c.V[x] == kk) && (c.PC += 2)
              end)
  Instruction("SNE Vx, val", "4xkk", "Jump over next instruction if (Vx != val)",
              function (c, x, kk)
              (c.V[x] != kk) && (c.PC += 2)
              end)
  Instruction("SE Vx, Vy", "5xy0", "Jump over next instruction if (Vx == Vy)",
              function (c, x, y)
              (c.V[x] == c.V[y]) && (c.PC += 2)
              end)
  Instruction("LD Vx, val", "6xkk", "Vx = val",
              function (c, x, kk)
              c.V[x] == kk
              end)
  Instruction("ADD Vx, val", "7xkk", "Vx = Vx + val",
              function (c, x, kk)
              c.V[x] += kk
              end)
  Instruction("LD Vx, Vy", "8xy0", "Vx = Vy",
              function (c, x, y)
              c.V[x] = c.V[y]
              end)
  Instruction("OR Vx, Vy", "8xy1", "Vx = Vx or Vy",
              function (c, x, y)
              c.V[x] = c.V[x] | c.V[y]
              end)
  Instruction("AND Vx, Vy", "8xy2", "Vx = Vx and Vy",
              function (c, x, y)
              c.V[x] = c.V[x] & c.V[y]
              end)
  Instruction("XOR Vx, Vy", "8xy3", "Vx = Vx xor Vy",
              function (c, x, y)
              c.V[x] = c.V[x] ⊻ c.V[y]
              end)
  Instruction("ADD Vx, Vy", "8xy4", "Vx = Vx + Vy",
              function (c, x, y)
              c.V[16] = (UInt(c.V[x]) + UInt(x.V[y])) > 255 # Set carry flag
              c.V[x] += c.V[y]
              end)
  Instruction("SUB Vx, Vy", "8xy5", "Vx = Vx - Vy",
              function (c, x, y)
              c.V[16] = c.V[x] > x.V[y] # Set carry flag
              c.V[x] -= c.V[y]
              end)
  Instruction("SHR Vx {, Vy}", "8xy6", "Vx = Vx >> 1",
              function (c, x, y)
              c.V[16] = c.V[x] & 0x01
              c.V[x] >>=  1
              end)
  Instruction("SUBN Vx, Vy", "8xy7", "Vx = Vy - Vx",
              function (c, x, y)
              c.V[16] = c.V[y] > c.V[x]
              c.V[x] = c.V[y] - c.V[x]
              end)
  Instruction("SHL Vx {, Vy}", "8xyE", "Vx = Vx << 1",
              function (c, x, y)
              c.V[16] = c.V[x] & 0x01 != 0x01
              c.V[x] <<= 1
              end)
  Instruction("SNE Vx, Vy", "9xy0", "Jump over next instruction if (Vx != Vy)",
              function (c, x, y)
              if c.V[x] != c.V[y]
              end
              end)
  Instruction("LD I, addr", "Annn", "I = addr",
              function (c, nnn)
              c.I = nnn
              end)
  Instruction("JP V0, addr", "Bnnn", "Jump to address: V0 + addr",
              function (c, nnn)
              c.I = c.V[1] + nnn
              end)
  Instruction("RND Vx, kk", "Cxkk", "Vx = rand() and kk",
              function (c, x, kk)
              x.V[x] = rand(UInt8) & kk
              # Do we have to set V[16] here?
              end)
  Instruction("DRW Vx, Vy, n", "Dxyn", "Display n rows of sprites beginning at I to the position Vx, Vy",
              function (c, x, y, n)
              c.V[16] = 0       # Will contain 1 if there was a collision
              for line in 1:n
                ypos = (c.V[y] + line) % disp_rows
                spriteline = decode(c.mem[c.I + line - 1]) # Decode one line (one byte) of sprite
                for pixel in eachindex(spriteline)
                  xpos = (c.V[x] + pixel) % disp_columns
                  # Check if there has been a collision while drawing the pixel
                  (spriteline[pixel] == 1) && (c.disp[xpos, ypos] == 1) && (c.V[16] = 1)
                  # Draw pixel of sprite. This is done using xor.
                  # i.e. drawing a white pixel on a white pixel will make it black again
                  c.disp[xpos, ypos] ⊻= spriteline[pixel]
                  end
                end
              end)
  Instruction("SKP Vx", "Ex9E", "Skip next instruction if key Vx is pressed",
              function (c, x)

              end)
  Instruction("SKNP Vx", "ExA1", "Skip next instruction if key Vx is not pressed",
              function (c, x)

              end)
  Instruction("LD Vx, DT", "Fx07", "Vx = DT",
              function (c, x)

              end)
  Instruction("LD Vx, K", "Fx0A", "Wait for key press and put pressed key into Vx",
              function (c, x)

              end)
  Instruction("LD DT, Vx", "Fx15", "DT = Vx",
              function (c, x)

              end)
  Instruction("LD ST, Vx", "Fx18", "ST = Vx",
              function (c, x)

              end)
  Instruction("ADD I, Vx", "Fx1E", "I = I + Vx",
              function (c, x)

              end)
  Instruction("LD F, Vx", "Fx29", "I = address to sprite representing number in Vx",
              function (c, x)

              end)
  Instruction("LD B, Vx", "Fx33", "BCD of Vx is saved to I, I+1, and I+2",
              function (c, x)

              end)
  Instruction("LD [I], Vx", "Fx55", "V0 to Vx are safed to memory beginning at I",
              function (c, x)

              end)
  Instruction("LD Vx, [I]", "Fx65", "V0 to Vx are read from memory beginning at I",
              function (c, x)

              end)
]

"Parses the opcode description and returns mask, pattern to match and mask, type and name of parameters"
function parseopcode(op::String)
  mask = 0x0000
  match = 0x0000
  parameters = Dict{Symbol, Int}()

  for c in op
    mask <<= 4
    match <<= 4
    if isxdigit(c)  # We have a valid hex digit
      mask |= 0xF
      match |= parse(UInt8, c, base = 16)
    else            # We have something else (i.e. a parameter)
      parameters[Symbol(c)] = get(parameters, Symbol(c), 0) + 1
      # Shift mask and match one nibble to the left to set next hex digit
    end
  end

  # Convert number of occurences of parameter into correct type
  parameters = Dict(s => c < 3 ? UInt8 : UInt16 for (s,c) in parameters)

  return mask, match, parameters
end

end # module


text:
import Pandas

@testset "summarystats" begin
    rng = MersenneTwister(42)
    nchains, ndraws = 4, 10
    idata = convert_to_inference_data(Dict(
        "a" => randn(rng, nchains, ndraws),
        "b" => randn(rng, nchains, ndraws, 3, 4),
    ))

    @test summarystats(idata) isa Pandas.DataFrame
    @test summarystats(idata; fmt = "wide") isa Pandas.DataFrame
    @test summarystats(idata; fmt = "long") isa Pandas.DataFrame
    s = summarystats(idata)
    @test "a" ∈ Pandas.index(s)
    @test "b" ∉ Pandas.index(s)
    @test "b[0,0]" ∉ Pandas.index(s)
    @test "b[1,1]" ∈ Pandas.index(s)
    @test "b[0,0]" ∈ Pandas.index(summarystats(idata; index_origin = 0))

    s2 = summarystats(idata; fmt = "xarray")
    @test s2 isa ArviZ.Dataset
end

@testset "ArviZ.summary" begin
    rng = MersenneTwister(42)
    nchains, ndraws = 4, 10
    data = Dict(
        "a" => randn(rng, nchains, ndraws),
        "b" => randn(rng, nchains, ndraws, 3, 4),
    )

    @test ArviZ.summary(data) isa Pandas.DataFrame
end


text:
#--------------------------------------------------------------------------
# t-test
#--------------------------------------------------------------------------
struct TTReturn
    title::String
    colnms::Vector{String}
    array::Array
    levels::Vector
    μ0::Real
    t::Float64
    dof::Float64 # for UnequalVarianceTTest
    p_left::Float64
    p_both::Float64
    p_right::Float64
    paired::Bool
    welch::Bool
end
#pvalue(x::TTReturn; tail=:both) = pvalue(TDist(x.dof), x.t; tail=tail)

"""
    ttest(df::DataFrame,varname::Symbol;byvar::Union{Nothing,Symbol}=nothing,sig=95,welch=false)
    ttest(df::DataFrame,varname1::Symbol,varname2::Symbol;paired::Bool=false,welch=false,sig=95)
    ttest(df::DataFrame,varname::Symbol,μ0::Real;sig=95)
    ttest(var1::AbstractVector,var2::AbstractVector;paired::Bood=false,welch::Bool=false,sig=95)

Performs one-way and two-way t tests and produces a summary table with t statistic,
degree of freedom, and associated p-values. For paired t tests, specify `paired = true`.
For two samples with unequal variances, use `welch = true` option.

### Example 1: One sample t test comparing a variable to a value (μ under H₀)

```
julia> t = ttest(auto, :price, 25000)
One-sample t test
| Variable |      N |      Mean |        SD |       SE |    95% LB |    95% UB |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
|    price |     74 | 6165.2568 | 2949.4959 | 342.8719 | 5481.9140 | 6848.5995 |
diff = mean(price)
H₀: diff = 25000
t = -54.9323 (df = 73)

Hₐ: diff < 25000                Hₐ: diff != 25000               Hₐ: diff > 25000
Pr(T < t) = 0.0000           Pr(|T| > |t|) = 0.0000           Pr(T > t) = 1.0000
```

### Example 2: Two sample t test by groups (equal variances)

```
julia> ttest(auto, :price, byvar = :foreign)
Two-sample t test with equal variances
|  foreign |      N |      Mean |        SD |       SE |    95% LB |    95% UB |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
|        0 |     52 | 6072.4231 | 3097.1043 | 429.4911 | 5210.1837 | 6934.6624 |
|        1 |     22 | 6384.6818 | 2621.9151 | 558.9942 | 5222.1898 | 7547.1738 |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
| combined |     74 | 6165.2568 | 2949.4959 | 342.8719 | 5481.9140 | 6848.5995 |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
|     diff |        | -312.2587 |           | 754.4488 | -1816.225 | 1191.7075 |

diff = mean(0) - mean(1)
H₀: diff = 0
t = -0.4139 (df = 72)

Hₐ: diff < 0                      Hₐ: diff != 0                     Hₐ: diff > 0
Pr(T < t) = 0.3401           Pr(|T| > |t|) = 0.6802           Pr(T > t) = 0.6599
```

### Example 3: Two sample t test with unequal variances

```
julia> ttest(auto,:price,by=:foreign, welch=true)
Two-sample t test with unequal variances
|  foreign |      N |      Mean |        SD |       SE |    95% LB |    95% UB |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
|        0 |     52 | 6072.4231 | 3097.1043 | 429.4911 | 5210.1837 | 6934.6624 |
|        1 |     22 | 6384.6818 | 2621.9151 | 558.9942 | 5222.1898 | 7547.1738 |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
| combined |     74 | 6165.2568 | 2949.4959 | 342.8719 | 5481.9140 | 6848.5995 |
├──────────┼────────┼───────────┼───────────┼──────────┼───────────┼───────────┤
|     diff |        | -312.2587 |           | 704.9376 | -1730.856 | 1106.3386 |

diff = mean(0) - mean(1)
H₀: diff = 0
t = -0.4430 (df = 46.4471)

Hₐ: diff < 0                      Hₐ: diff != 0                     Hₐ: diff > 0
Pr(T < t) = 0.3299           Pr(|T| > |t|) = 0.6599           Pr(T > t) = 0.6701
```

### Example 4: Two sample paired t test comparing two variables

```
julia> ttest(fuel, :mpg1, :mpg2, paired = true)
Paired t test
| Variable |      N |     Mean |       SD |       SE |   95% LB |   95% UB |
├──────────┼────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
|     mpg1 |     12 | 21.00000 |  2.73030 |  0.78817 | 19.26525 | 22.73475 |
|     mpg2 |     12 | 22.75000 |  3.25087 |  0.93845 | 20.68449 | 24.81551 |
├──────────┼────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
|     diff |     12 | -1.75000 |  2.70101 |  0.77971 | -3.46614 | -0.03386 |

diff = mean(mpg1) - mean(mpg2)
H₀: diff = 0
t = -2.2444 (df = 11)

Hₐ: diff < 0                   Hₐ: diff != 0                    Hₐ: diff > 0
Pr(T < t) = 0.0232         Pr(|T| > |t|) = 0.0463         Pr(T > t) = 0.9768
```

### Example 5: Two sample unpaired t test comparing two variables

```
julia> ttest(fuel, :mpg1, :mpg2)
Two-sample t test with equal variances
| Variable |      N |     Mean |       SD |       SE |   95% LB |   95% UB |
├──────────┼────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
|     mpg1 |     12 | 21.00000 |  2.73030 |  0.78817 | 19.26525 | 22.73475 |
|     mpg2 |     12 | 22.75000 |  3.25087 |  0.93845 | 20.68449 | 24.81551 |
├──────────┼────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
| combined |     24 | 21.87500 |  3.06895 |  0.62645 | 20.57909 | 23.17091 |
├──────────┼────────┼──────────┼──────────┼──────────┼──────────┼──────────┤
|     diff |        | -1.75000 |          |  1.22552 | -4.29157 |  0.79157 |

diff = mean(mpg1) - mean(mpg2)
H₀: diff = 0
t = -1.4280 (df = 22)

Hₐ: diff < 0                   Hₐ: diff != 0                    Hₐ: diff > 0
Pr(T < t) = 0.0837         Pr(|T| > |t|) = 0.1673         Pr(T > t) = 0.9163
```

"""
function ttest(df::DataFrame,varname::Symbol; by::Union{Nothing,Symbol} = nothing, sig = 95, welch = false)
    if by == nothing
        error("`by` is required.")
    end

    dft = df[completecases(df[[varname,by]]),[varname,by]]

    lev = sort(collect(Set(dft[by])))
    if length(lev) != 2
        error(by," must have two levels; it has ",length(lev), " levels.")
    end

    val = Vector(2)
    val[1] = Vector(dft[dft[by] .== lev[1],varname])
    val[2] = Vector(dft[dft[by] .== lev[2],varname])
    if length(val[1]) == 0 || length(val[2]) == 0
        error("On or both groups have zero observations")
    end
    tt = ttest(val[1],val[2],paired=false,welch=welch,sig=sig,levels=lev)
    tt.colnms[1] = string(by)
    return tt
end
function ttest(df::DataFrame, var1::Symbol, var2::Symbol; sig = 95, paired = false, welch = false)

    if paired == true
        ba = completecases(df[[var1,var2]])
        x = Vector(df[ba,var1])
        y = Vector(df[ba,var2])
    else
        x = Vector(dropna(df[var1]))
        y = Vector(dropna(df[var2]))
    end

    length(x) > 0 && length(y) > 0 || error("One or both variables are empty")

    return ttest(x,y,paired=paired,welch=welch,levels=[var1,var2])
end
function ttest(x::AbstractVector,y::AbstractVector; paired::Bool=false,welch::Bool=false,sig=95,levels=Any[:x,:y])

    paired && length(x) != length(y) && error("Paired t test requires equal lengths in input vectors")

    if paired == true
        title = "Paired t test"
        tt = OneSampleTTest(x, y, 0.0)
    elseif welch == false
        title = "Two-sample t test with equal variances"
        tt = EqualVarianceTTest(x, y)
    else
        title = "Two-sample t test with unequal variances"
        tt = UnequalVarianceTTest(x, y)
    end

    # compute standard errors and confidence intervals
    N = Vector{Int64}(4)
    MEAN = Vector{Float64}(4)
    SD = Vector{Float64}(4)
    SE = Vector{Float64}(4)
    LB = Vector{Float64}(4)
    UB = Vector{Float64}(4)

    val = [x,y]
    for i = 1:2
        N[i] = length(val[i])
        MEAN[i] = mean(val[i])
        SD[i] = std(val[i])
        SE[i] = SD[i] / sqrt(N[i])
        LB[i],UB[i] = StatsBase.confint(OneSampleTTest(val[i],0),1-sig/100)
    end

    # combined
    N[3] = N[1] + N[2]
    MEAN[3] = mean(vcat(val[1],val[2]))
    SD[3] = std(vcat(val[1],val[2]))
    SE[3] = SD[3] / sqrt(N[3])
    LB[3],UB[3] = StatsBase.confint(OneSampleTTest(vcat(val[1],val[2]),0),1 .- sig/100)

    # diff
    N[4] = paired ? N[1] : 0 # not to be used
    MEAN[4] = tt.xbar
    SD[4] = paired ? tt.stderr*sqrt(N[4]) : 0.0 # not to be used
    SE[4] = tt.stderr
    LB[4],UB[4] = StatsBase.confint(tt,1 .- sig/100)

    return TTReturn(title,
        ["Variable", "N", "Mean", "SD", "SE", string(sig,"% LB"), string(sig,"% UB")],
        [vcat(levels,"combined","diff"),N,MEAN,SD,SE,LB,UB],
        levels,
        tt.μ0,
        tt.t,
        tt.df,
        pvalue(tt, tail = :left),
        pvalue(tt),
        pvalue(tt, tail = :right),
        paired,
        welch)
end
function ttest(df::DataFrame, varname::Symbol, μ0::Real; sig = 95)
    v = Vector(df[completecases(df[[varname]]),varname])

    tt = ttest(v,μ0,sig=sig)
    tt.array[1][1] = string(varname)
    return tt
end
function ttest(var::AbstractVector,μ0::Real = 0; sig = 95)

    N = length(var)
    if N == 0
        error("Vector is empty")
    end
    MEAN = mean(var)
    SD = std(var)

    tt = OneSampleTTest(var, μ0)
    title = "One-sample t test"

    # compute standard errors and confidence intervals
    SE = SD / sqrt(N)
    (LB, UB) = StatsBase.confint(OneSampleTTest(var,0),1-sig/100)

    return TTReturn(title,
        ["Variable", "N", "Mean", "SD", "SE", string(sig,"% LB"), string(sig,"% UB")],
        [["x"],[N],[MEAN],[SD],[SE],[LB],[UB]],
        [],
        tt.μ0,
        tt.t,
        tt.df,
        pvalue(tt, tail = :left),
        pvalue(tt),
        pvalue(tt, tail = :right),
        false,
        false)
end


text:
const default_initializer = ArrayInitializer(t -> t<:Union{Tuple, NamedTuple, Pair}, (T, sz) -> similar(arrayof(T), sz))

"""
    collect_columns(itr)

Collect an iterable as a `Columns` object if it iterates `Tuples` or `NamedTuples`, as a normal
`Array` otherwise.

# Examples

    s = [(1,2), (3,4)]
    collect_columns(s)

    s2 = Iterators.filter(isodd, 1:8)
    collect_columns(s2)
"""
collect_columns(args...) = vec(collect_structarray(args..., initializer = default_initializer))
collect_columns(s::StructVector) = s
collect_empty_columns(itr) = collect_empty_structarray(itr, initializer = default_initializer)

grow_to_columns!(args...) = grow_to_structarray!(args...)
collect_to_columns!(args...) = collect_to_structarray!(args...)

_first(x) = first(x)
_first(p::Pair) = p.first => first(p.second)

function collect_columns_flattened(itr::T) where {T}
    elem = iterate(itr)
    if (elem === nothing)
        S = Core.Compiler.return_type(_first∘first, Tuple{T})
        return default_initializer(S, (0,))
    end
    el, st = elem
    collect_columns_flattened(itr, el, st)
end

function collect_columns_flattened(itr, el, st)
    fr = iterate(el)
    while fr === nothing
        elem = iterate(itr, st)
        elem === nothing && return collect_empty_columns(el)
        el, st = elem
        fr = iterate(el)
    end
    dest = collect_columns(el, fr)
    collect_columns_flattened!(dest, itr, el, st)
end

function collect_columns_flattened!(dest, itr, el, st)
    while true
        elem = iterate(itr, st)
        elem === nothing && break
        el, st = elem
        dest = grow_to_columns!(dest, el)
    end
    return dest
end

function collect_columns_flattened(itr, el::Pair, st)
    fr = iterate(el.second)
    while fr === nothing
        elem = iterate(itr, st)
        elem === nothing && return collect_empty_columns(el.first => i for i in el.second)
        el, st = elem
        fr = iterate(el.second)
    end
    dest_data = collect_columns(el.second, fr)
    dest_key = collect_columns(el.first for i in dest_data)
    init = Columns{Pair{eltype(dest_key), eltype(dest_data)}}((dest_key, dest_data))
    collect_columns_flattened!(init, itr, el, st)
end

function collect_columns_flattened!(dest::Columns{<:Pair}, itr, el::Pair, st)
    dest_key, dest_data = columns(dest)
    while true
        elem = iterate(itr, st)
        elem === nothing && break
        el, st = elem
        n = length(dest_data)
        dest_data = grow_to_columns!(dest_data, el.second)
        dest_key = grow_to_columns!(dest_key, el.first for i in (n+1):length(dest_data))
    end
    return Columns{Pair{eltype(dest_key), eltype(dest_data)}}((dest_key, dest_data))
end


text:
using TimeModels, DSP

"""
    kalman_filter(s::SSR; kwargs...)

Use a kalman filter to estimate signal amplitude

#### Arguments

* 'freq_of_interest': Frequencies to analyse (modulation rate)
* 'model_type': Type of model to use (single channel ASSR)
* 'ID': Value to store as ID (empty "")
* 'results_key': Where to store results in processing dictionary ("statistics")
* 'proc_noise_cov': Process noise covariance (1e-10)
* 'obs_noise_cov': Observation noise covariance (covariance of signal)
* 'error_cov': Error covariance, measure of accuracy of state estimate (1.0)
* 'x0': Initial state estimates

"""
function TimeModels.kalman_filter(s::SSR;
            freq_of_interest::Union{Real, AbstractArray} = modulationrate(s),
            model_type::Function = acoustic_model, reduction_method::Function = mean,
            ID::AbstractString = "", results_key::AbstractString = "statistics",
            proc_noise_cov = 1e-10, obs_noise_cov = cov(s.data[find(~isnan(s.data[:, 1])), :]),
            error_cov = 1.0, x0=[0.0, 0.0], kwargs...)

    Logging.info("Running single channel Kalman filter on SSR data for $(size)")

    # Convert data to 64 bit type
    s.data = convert(Array{Float64, 2}, s.data)
    obs_noise_cov = convert(Array{typeof(proc_noise_cov)}, obs_noise_cov)

    for freq in freq_of_interest
        for elec in channelnames(s)

            elec_idx = findfirst(channelnames(s), elec)

            debug("Processing channel $elec in index $elec_idx")

            model = model_type(freq, samplingrate(s), proc_noise_cov, obs_noise_cov[elec_idx, elec_idx], error_cov, x0)

            f = kalman_smooth(s.data[:, elec_idx], model)

            amp = reduction_method(model_amplitude(f, s))
            noi = std(model_amplitude_time(f, s))
            pha = reduction_method(DSP.unwrap(model_phase(f), range = pi))

            result = DataFrame( ID                  = vec([ID]),
                                Channel             = vec([elec]),
                                ModulationRate      = modulationrate(s),
                                AnalysisType        = vec(["Kalman"]),
                                AnalysisFrequency   = freq,
                                SignalAmplitude     = amp,
                                SignalPhase         = pha,
                                NoiseAmplitude      = noi,
                                SNRdB               = 10 * log10((amp^2) / (noi^2)),
                                StateVariable1      = f.smoothed[end, 1],
                                StateVariable2      = f.smoothed[end, 2],
                                ProcessCov          = proc_noise_cov,
                                ObservationCov      = obs_noise_cov[elec_idx, elec_idx],
                                ErrorCov            = error_cov)

            result = add_dataframe_static_rows(result, kwargs)

            if haskey(s.processing, results_key)
                s.processing[results_key] = vcat(s.processing[results_key], result)
            else
                s.processing[results_key] = result
            end
        end
    end

    return s
end

"""
Return the amplitude for the SSR state space model
"""
function model_amplitude{T <: Number}(filte::Union{KalmanSmoothed{T}, TimeModels.KalmanSmoothedMinimal{T}}, s::SSR)

    amp = filte.smoothed[find(~isnan(s.data[:, 1])), :]
    amp = sqrt(amp[:, 1].^2 .+ amp[:, 2].^2)
end

function model_amplitude{T <: Number}(filte::Union{KalmanSmoothed{T}, TimeModels.KalmanSmoothedMinimal{T}})

    amp = filte.smoothed
    amp = sqrt(amp[:, 1].^2 .+ amp[:, 2].^2)
end

function model_amplitude_time{T <: Number}(filte::Union{KalmanSmoothed{T}, TimeModels.KalmanSmoothedMinimal{T}}, s::SSR)

    amp = filte.smoothed[find(~isnan(s.data[:, 1])), :]
    amp = sqrt(amp[:, 1].^2 .+ amp[:, 2].^2)
end

function model_amplitude_time{T <: Number}(filte::Union{KalmanSmoothed{T}, TimeModels.KalmanSmoothedMinimal{T}})

    amp = filte.smoothed
    amp = sqrt(amp[:, 1].^2 .+ amp[:, 2].^2)
end

function model_phase{T <: Number}(filte::Union{KalmanSmoothed{T}, TimeModels.KalmanSmoothedMinimal{T}})

    pha = filte.smoothed
    pha = atan(pha[:, 2] ./ pha[:, 1])
end

function model_phase{T <: Number}(filte::KalmanFiltered{T})

    pha = filte.filtered
    pha = atan(pha[:, 2] ./ pha[:, 1])
end

"""
Build a basic SSR model with single sinusoid at the modulation rate.

Returns a state space model.

#### Parameters

* 'modulation_rate': Modulation rate of sinusoid to match
* 'sample_rate': Sample rate of signal
* 'v': Process noise covariance (1e-10)
* 'w': Observation noise covariance (covariance of signal)
* 'p': Error covariance, measure of accuracy of state estimate (1.0)
* 'x0': Initial state estimates
* 'num_sensors': Number of sensors

"""
function acoustic_model{T}(modulation_rate::T, sample_rate::T, v::T, w::T, p::T, x0::Vector{T};
            num_sensors::Int = 1, kwargs...)

    num_states = length(x0)                 # Number of state variables
    V = v * eye(num_states)                 # Process noise covariance
    W = w * eye(num_sensors)                # Observation noise covariance
    P0 = p * eye(num_states)                # Initial noice covariance

    acoustic_model(num_sensors, modulation_rate, sample_rate, V, W, P0, x0; kwargs...)
end

function acoustic_model{T}(num_sensors::Int, modulation_rate::T, sample_rate::T,
            V::Array{T, 2}, W::Array{T, 2}, P0::Array{T, 2}, x0::Array{T, 1}; kwargs...)

    Δt = 1 / sample_rate
    ω = 2 * pi * modulation_rate
    num_states = length(x0)

    ## Process transition
    F = eye(num_states)

    # G = Array(Function, num_sensors, num_states)
    #
    # ## Observation
    # for sensor in 1:num_sensors
    #     G[sensor, 1] = k ->  cos(ω * Δt * k)
    #     G[sensor, 2] = k -> -sin(ω * Δt * k)
    # end

    function G(k::Int)
        [cos(ω * Δt * k); -sin(ω * Δt * k)]'
    end

    StateSpaceModel(_->F, _->V, G, _->W, x0, P0)
end

"""
Build a SSR model with single sinusoid at the modulation rate and artifacts at the carrier rate.

Returns a state space model.

"""
function artifact_model{T}(modulation_rate::T, sample_rate::T, v::Vector{T}, w::T, p::Vector{T}, x0::Vector{T}, fa::T; kwargs...)

    # Put in check that fa was passed in and v2

    Δt = 1 / sample_rate
    ω = 2 * pi * modulation_rate
    n = length(x0)

    ## Process transition and noise covariance

    F = eye(n)
    V = diagm(v)  # We expect little variation in the state matrix

    ## Observation and noise covariance

    function G1(k);  cos(ω * Δt * k); end
    function G2(k); -sin(ω * Δt * k); end
    function G3(k);  triangle_artifact(Δt * k, fa; kwargs...); end
    G = reshape([G1, G2, G3], 1, n)
    W = diagm([w])   # Noise Covariance

    ## Inital guesses at state and error covariance
    # x0 set in function
    P0 = diagm(p)

    StateSpaceModel(F, V, G, W, x0, P0)
end

# Internal function to create an artifact model.
function triangle_artifact{T <: AbstractFloat}(t::AbstractArray{T, 1}, fa::Number; kwargs...)

    info("Creating a trianglar artifact model")

    a = Array(T, size(t))

    for i in 1:length(t)
        a[i] = triangle_artifact(t[i], fa; kwargs...)
    end
    return a
end

# Internal function to create an artifact model.
function triangle_artifact{T <: AbstractFloat}(t::T, fa::Number; adj = 0.0, offset = 0.0, tol = 0.0001)

    dis1 = mod(-(t + offset), 1 / (fa+adj))
    dis2 = mod( (t + offset), 1 / (fa+adj))
    dis = min(dis1, dis2)

    if dis < tol
        x = 1.0 - abs(dis / tol)
    else
        x = 0.0
    end

    return x
end

"""
    set_nans(a)

Set samples to nan that are likely to be artifacts.

Set values to be excluded from analysis to NaN.
Usually in ASSR analysis you would remove the epochs, but for Kalman filter the timing would then be lost.
Instead we set the value to NaN and the filter does not updates its state estimates over these values, but does propogate the internal system.
The artifact ± `boundary` are set to NaN. The first and last 5 samples are not altered.

Currently only works on a single channel.
"""
function set_nans(a::SSR; boundary::Int = round(Int, 8192/2), perc=0.1)

    @assert size(a.data, 2) == 1 "Set NaNs only works for single channel SSR"

    info("Removing artifacts (NaN method)")
    debug("Input data consists of $(round(sum(isnan(a.data)) / length(a.data), 3))% NaNs")

    # We want to maintain the first and last 5 samples for Kalman filtering
    # Save them here and replace them at the end in case they were removed
    start_data = a.data[1:5, :]
    end_data = a.data[end-5:end, :]

    while (sum(isnan(a.data)) / length(a.data)) <= perc

        idx = findfirst(a.data, maximum(a.data))

        debug("Removing a segment of length $(2 * boundary) samples centered at sample $idx")

        new_min = maximum([1, idx - boundary])
        new_max = minimum([length(a.data), idx + boundary])

        a.data[new_min:new_max] = NaN
    end

    a.data[1:5, :] = start_data
    a.data[end-5:end, :] = end_data

    debug("Output data consists of $(round(sum(isnan(a.data)) / length(a.data), 3))% NaNs")

    return a
end


text:
# Autogenerated wrapper script for Antic_jll for aarch64-linux-musl
export libantic

using FLINT_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Antic")
JLLWrappers.@declare_library_product(libantic, "libantic.so.0")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libantic,
        "lib/libantic.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module PartitionedArrays

using SparseArrays
using SparseMatricesCSR
using LinearAlgebra
using Printf
import MPI
import IterativeSolvers
import Distances

export AbstractBackend
export prun, prun_debug
export AbstractPData
export SequentialData
export MPIData
export num_parts
export i_am_main
export map_parts
export get_part_ids
export get_backend
export MAIN
export map_main
export get_part
export get_main_part
export gather!
export gather_all!
export gather
export gather_all
export scatter
export scatter!
export emit
export reduce_main
export reduce_all
export iscan
export iscan_main
export iscan_all
export xscan
export xscan_main
export xscan_all
export async_exchange!
export exchange!
export exchange
export async_exchange
export Table
export get_data
export get_ptrs
export length_to_ptrs!
export rewind_ptrs!
export discover_parts_snd
export AbstractIndexSet
export IndexSet
export ExtendedIndexRange
export IndexRange
export num_gids
export num_lids
export num_oids
export num_hids
export oids_are_equal
export hids_are_equal
export lids_are_equal
export touched_hids
export get_lid_to_gid
export get_lid_to_part
export get_oid_to_lid
export get_hid_to_lid
export get_lid_to_ohid
export get_gid_to_lid
export Exchanger
export empty_exchanger
export allocate_rcv_buffer
export allocate_snd_buffer
export PRange
export WithGhost
export with_ghost
export NoGhost
export no_ghost
export add_gid!
export add_gids
export add_gids!
export to_lids!
export to_gids!
export PCartesianIndices
export PVector
export local_view
export global_view
export async_assemble!
export assemble!
export PSparseMatrix
export matrix_exchanger
export nzindex
export nziterator
export compresscoo
export indextype
export PTimer
export tic!
export toc!
export print_timer

export SequentialBackend
export sequential

export MPIBackend
export mpi

include("Helpers.jl")

include("Interfaces.jl")

include("SparseUtils.jl")

include("IndexSets.jl")

include("PTimers.jl")

include("SequentialBackend.jl")

include("MPIBackend.jl")

end


text:
function rodrigues!(Rmat::Array{Float64,2}, rvec::Vector{Float64})
  th = norm(rvec)
  rotv = rvec./th
  cth = cos(th)
  rotv1c = rotv.*(1-cth)

  Rmat[1,1] = 0.0
  Rmat[1,2] = -rotv[3]
  Rmat[1,3] = rotv[2]
  Rmat[2,1] = rotv[3]
  Rmat[2,2] = 0
  Rmat[2,3] = -rotv[1]
  Rmat[3,1] = -rotv[2]
  Rmat[3,2] = rotv[1]
  Rmat[3,3] = 0
  Rmat .*= sin(th)

  Rmat[1,1] += rotv1c[1]*rotv[1]
  Rmat[1,2] += rotv1c[1]*rotv[2]
  Rmat[1,3] += rotv1c[1]*rotv[3]
  Rmat[2,1] += rotv1c[2]*rotv[1]
  Rmat[2,2] += rotv1c[2]*rotv[2]
  Rmat[2,3] += rotv1c[2]*rotv[3]
  Rmat[3,1] += rotv1c[3]*rotv[1]
  Rmat[3,2] += rotv1c[3]*rotv[2]
  Rmat[3,3] += rotv1c[3]*rotv[3]

  Rmat[1,1] += cth
  Rmat[2,2] += cth
  Rmat[3,3] += cth
  nothing
end

function buildtagdict(cTt,
                      Ql::Rotations.Quat,
                      tvec, # ::Union{Vector{Float64},StaticArrays.SArray{Tuple{3},<:Real,1,3}}
                      tagsize::Float64,
                      bTcl )
  #
  # @show tvec, q
  cTt = tvec ∘ LinearMap(Ql)
  # @show bTcl, cTt
  bTt = bTcl ∘ cTt
  bP2t = getTagPP2(bTt)

  onetag = Dict{Symbol, Any}()
  onetag[:pos] = tvec.translation[:]
  onetag[:quat] = Float64[Ql.w; Ql.x; Ql.y; Ql.z]
  onetag[:tagxy] = Float64[tagsize; tagsize]
  onetag[:bearing] = [atan(-tvec.translation[1], tvec.translation[3]);]
  onetag[:range] = [norm(tvec.translation[[1;3]]);]
  # onetag[:tRYc] = convert(RotXYZ, q).theta2
  onetag[:bP2t] = bP2t
  onetag
end
# function buildtagdict(cTt,
#                       q::Rotations.Quat,
#                       tvec::CoordinateTransformations.Translation{T},
#                       tagsize::Float64,
#                       bTcl ) where T
#   buildtagdict(cTt, q, tvec, tagsize, bTcl )
# end

function getAprilTagTransform(tag::AprilTag,
                              camK::Array{Float64,2},
                              k1::Float64,
                              k2::Float64,
                              tagsize::Float64 )
  #
  imgPts = zeros(4,2,1)
  for i in 1:4, j in 1:2
    imgPts[i,j,1] = tag.p[i][j]
  end

  tsz = tagsize/2.0
  objPts = zeros(4,3,1)
  objPts[1,1,1] = -tsz; objPts[1,2,1] = -tsz;
  objPts[2,1,1] =  tsz; objPts[2,2,1] = -tsz;
  objPts[3,1,1] =  tsz; objPts[3,2,1] =  tsz;
  objPts[4,1,1] = -tsz; objPts[4,2,1] =  tsz;

  distCoeffs = zeros(5)
  distCoeffs[1] = k1
  distCoeffs[2] = k2

  # Python OpenCV
  # ret, rvec, tvec = cv2.solvePnP(objPts, imgPts, camK, distCoeffs)
  # Rmat = zeros(3,3)
  # rodrigues!(Rmat,rvec[:])
  # q = convert(Quat, RotMatrix{3}(Rmat))

  cTt = AprilTags.homographytopose(tag.H, camK[1,1], camK[2,2], camK[1,3], camK[2,3], taglength=tagsize)

  q = Quat(cTt[1:3,1:3])
  return q, Translation(SVector(cTt[1:3,4]...)), camK
end
# objPts.push_back(cv::Point3f(-s,-s, 0));
# objPts.push_back(cv::Point3f( s,-s, 0));
# objPts.push_back(cv::Point3f( s, s, 0));
# objPts.push_back(cv::Point3f(-s, s, 0));

# objectPoints = np.random[:random]((4,3,1))
# imagePoints = np.random[:random]((4,2,1))
# cameraMatrix = np.eye(3)

function getAprilTagTransform(tag::AprilTag;
                              fx::Float64=100.0,
                              fy::Float64=100.0,
                              cx::Float64=320.0,
                              cy::Float64=240.0,
                              k1::Float64=0.0,
                              k2::Float64=0.0,
                              tagsize::Float64=0.172 )
  camK = [fx 0 cx; 0 fy cy; 0 0 1]
  getAprilTagTransform(tag,
                       camK,
                       k1,
                       k2,
                       tagsize )
end

function prepTagBag(TAGS; iterposes::Int=99999999)
  tvec = Translation(0.0,0,0)
  q = Quat(1.0,0,0,0)
  # package tag detections for all keyframes in a tag_bag
  tag_bag = Dict{Int, Any}()
  for psid in 0:minimum([(length(TAGS)-1), iterposes])
    tag_det = Dict{Int, Any}()
    for tag in TAGS[psid+1]
      q, tvec = getAprilTagTransform(tag, camK, k1, k2, tagsize)
      cTt = tvec ∘ CoordinateTransformations.LinearMap(q)
      tag_det[tag.id] = buildtagdict(cTt, q, tvec, tagsize, bTc)
    end
    tag_bag[psid] = tag_det
  end
  return tag_bag
end

# get Pose2Pose2 tag orientation transform
function getTagPP2(bTt)
  # @show bTt
  cVz = LinearMap(Quat(bTt.linear))([1.0;0;0])
  wYt = atan(cVz[2],cVz[1])
  Translation(bTt.translation[1],bTt.translation[2],0) ∘ LinearMap(RotZ(wYt))
end

## DETECT APRILTAGS FROM IMAGE DATA
function detectTagsInImgs(datafolder, imgfolder, resultsdir, camidxs; iterposes::Int=999999999)

    # prep keyframe image data
    camlookup = prepCamLookup(camidxs)

    # detect tags and extract pose transform
    IMGS, TAGS = detectTagsViaCamLookup(camlookup, joinpath(datafolder,imgfolder), resultsdir)

    # prep dictionary with all tag detections and poses
    tag_bag = prepTagBag(TAGS, iterposes=iterposes)

    # save the tag detections for later comparison
    fid=open(resultsdir*"/tags/pose_tags.csv","w")
    for pose in sort(collect(keys(tag_bag)))
      println(fid, "$pose, $(collect(keys(tag_bag[pose])))")
    end
    close(fid)
    return tag_bag
end

#


text:
abstract type AbstractDocument end

"""
    Document(id::AbstractString; score=nothing, payload=nothing, kwargs...) -> Document

A RediSearch Document that is returned from a search query.

# Attributes
- `id::String`: Redis Key id that the document is stored in.
- `score::Union{Float64,Nothing}`: Match Score for the document based on the search
criteria.
- `payload::Union{Nothing,AbstractDict}`: Payload returned buy the query.
- `kwargs`: Additional redis fields.
"""
mutable struct Document <: AbstractDocument
    id::String
    score::Union{Float64,Nothing}
    payload::Union{Nothing,Dict}
end

function Document(
    id::AbstractString;
    score::Union{Float64,Nothing},
    payload::Union{Nothing,AbstractDict}=Dict(),
    kwargs...
)
    for (k, v) in kwargs
        payload[k] = v
    end

    return Document(id, score, payload)
end