text:
module GeometricalOptics

export zplane, spherical, aspherical, zernike, bundle, bundle_as_array, opticalstack,
rms_spot, splot, rac, trace, bigtrace, bigtrace_to_bundle, bundle_as_array_big,
trace_extract_ray, trace_extract_terminus

using Plots
using PlotlyBase
Plots.plotly()
Plots.PlotlyBackend()
using ForwardDiff
using LinearAlgebra
using ZernikePolynomials

"""
    zplane(x, y, coeffs)

Composite type representing a plane z(x, y) = zpos. Argument coeffs of type `Vector{T} where T<:Real` has
entries `coeffs = [zpos]`.

# Examples
```julia-repl
julia> zplane(0., 2., [3.])
3.0
```
"""
function zplane(x::T where T<:Real,y::T where T<:Real,coeffs::Vector{T} where T<:Real)
    return coeffs[1]
end

"""
    zernike(x, y, coeffs)

Composite type representing OSA/ANSI type Zernike polynomial surface as a function of lateral coordinate `x, y`. Argument `coeffs` of type `Vector{T} where T<:Real` has
entries `coeffs = [rad, C0, C1, C2, ...]`. The normalization radius is `rad` and the corresponding Zernike surface is given
by the contributions of the Zⱼ with arguments over the cartesian coordinates normalized to the unit disk, i.e.:

``z(x, y) = ∑ⱼCⱼZⱼ(r, θ)`` where ``r ≡ √( x² + y² ) / rad`` and ``θ ≡ atan(y / x)``

__NOTE this function evaluates to `NaN` outside of the normalization radius.__

# Examples
```julia-repl
julia> zernike(0., 2., [2., 0., 1.])
2.0
```
"""
function zernike(x::T where T<:Real,y::T where T<:Real,coeffs::Vector{T} where T<:Real)
    normRad=coeffs[1]
    if (x/normRad)^2 + (y/normRad)^2 >1.0
        return NaN
    else
        z=0.0
        for i=1:(length(coeffs)-1)
            Z=Zernike(i-1,index=:OSA,coord=:cartesian)
            z=z+coeffs[i+1]*Z(x/normRad,y/normRad)
        end
        return z
    end
end

"""
    spherical(x, y, coeffs)

Composite type representing spherical surface with apex at defined point. Argument `coeffs` of type `Vector{T} where T<:Real` has
entries `coeffs = [zpos, signed_radius]`. Here `zpos` is the apex and `signed_radius` has magnitude equal to the radius of the sphere
and sign which sets the direction of curvature of the sphere.

# Examples
```julia-repl
julia> spherical(0., 2., [1., -10.])
0.7979589711327115
```
"""
function spherical(x::T where T<:Real,y::T where T<:Real,coeffs::Vector{T} where T<:Real)
    z=coeffs[1]+coeffs[2]-sign(coeffs[2])*(sqrt(coeffs[2]^2-x^2-y^2))
    return z
end

"""
    aspherical(x, y, coeffs)

Composite type representing aspherical surface as a function of lateral coordinate `x,y`. Argument `coeffs` of type `Vector{T} where T<:Real` has
entries `coeffs = [zpos, R, ``κ``, A4, A6, ...]` where the surface is

``z(x,y) = zpos + r² / R / (1 + √(1 - (1 + κ)r² / R²)) + A₄r⁴ + A₆r⁶ + ...`` where ``r² ≡ x² + y²``

# Examples
```julia-repl
julia> aspherical(0., 2., [1., -10., -4., 0.001])
1.0616649185805458
```
"""
function aspherical(x::T where T<:Real,y::T where T<:Real,coeffs::Vector{T} where T<:Real)
    rsq=x^2+y^2
    Rsq=coeffs[2]^2
    zb=coeffs[1]+rsq/coeffs[2]/(1+sqrt(1-(1+coeffs[3])*rsq/coeffs[2]^2))
    for j=1:length(coeffs)-3
        zb=zb+coeffs[j+3]*rsq^(2*(1+j))
    end
    return zb
end

struct impz{coeffs,surfDef}
    coeffs::coeffs
    surfDef::surfDef
end
function (iz::impz)(r)
    zret=iz.surfDef(r[1],r[2],iz.coeffs)
    return (r[3]-zret)^2
end

struct normalVec{impSurfIn}
    impSurfIn::impSurfIn
end
function (nv::normalVec)(R)
    gs=R -> ForwardDiff.gradient(nv.impSurfIn,R)
    gs(R) ./ sqrt(dot(gs(R),gs(R)))
end

struct refrSf{impSurf}
    impSurf::impSurf
end
function (rf::refrSf)(R0,D0,n0,n1;isReorient=true)
    nvf=normalVec(rf.impSurf)
    N=nvf(R0)
    if n1!=-n0
        if isReorient
            sg=sign(dot(N,D0))
            N=sg*N
        end
        if (n1^2-n0^2+n0^2*dot(N,D0)^2)<0; @error "TIR violates requested refraction for current ray"; end;
        k = sqrt(n1^2-n0^2+n0^2*dot(N,D0)^2)-n0*dot(N,D0)
    else
        k=sqrt(n1^2-n0^2+n0^2*dot(N,D0)^2)-n0*dot(N,D0)
    end
    (k*N+n0*D0)/n1
end

struct getIntersect_preInit{impSurfIn,sInit}
    impSurfIn::impSurfIn
    sInit::sInit
end
function (gi::getIntersect_preInit)(R,D)
    # for start position R and direction cosine vector D, implicit surface impSurfIn, and coefficient coeff,
    f(s)=gi.impSurfIn(R.+D.*s)
    gf=s->ForwardDiff.derivative(f,s)
    sg=gi.sInit
    tol=0.000001
    adj=tol+1.0
    while abs(adj)>tol
        adj=f(sg)/gf(sg)
        sg=sg-adj
    end
    sg
end

"""
    bundle(x, y, angx, angy, zpos)

Composite type representing the origination of a bundle of rays. Input zpos is type `T<:Real`
and is z-plane start position for the rays.
Other four arguments are typed in one of the following two patterns:
1. x and y are Type `AbstractArray{T,1} where T<:Real`, angx and angy are Type `T where T<:Real`
2. x and y are Type `T where T<:Real`, angx and angy are Type `AbstractArray{T,1} where T<:Real`}

When a variable is declared this type with specific parameters and called
as a function with no arguments, it returns a tuple (position, direction) each
entry of which is an output array of vectors corresponding to the input information.

# Examples
```julia-repl
julia> eval_bundle = bundle([0., 1.], [0.,1.], 0., 0., 0.)
bundle([0.0, 1.0], [0.0, 1.0], 0.0, 0.0, 0.0)

julia> eval_bundle()
([[0.0, 0.0, 0.0] [0.0, 1.0, 0.0]; [1.0, 0.0, 0.0] [1.0, 1.0, 0.0]], [[0.0, 0.0, 1.0] [0.0, 0.0, 1.0]; [0.0, 0.0, 1.0] [0.0, 0.0, 1.0]])
```
"""
struct bundle
    x::Union{AbstractArray{T,1} where T<:Real,Real}
    y::Union{AbstractArray{T,1} where T<:Real,Real}
    angx::Union{AbstractArray{T,1} where T<:Real,Real}# these are direction tangents.
    angy::Union{AbstractArray{T,1} where T<:Real,Real}
    zpos:: T where T<: Real
end
function (cb::bundle)()
    if typeof(cb.x)<:AbstractArray && typeof(cb.y)<:AbstractArray && typeof(cb.angx)<:Real && typeof(cb.angy)<:Real
        R0=[[cb.x[i].+eps(),cb.y[j].+eps(),cb.zpos.+eps()]  for i in 1:length(cb.x),j in 1:length(cb.y)]
        D0=[[sin(cb.angx),sin(cb.angy),sqrt(1-sin(cb.angx)^2-sin(cb.angy)^2)] for i in 1:length(cb.x),j in 1:length(cb.y)]
    elseif typeof(cb.angx)<:AbstractArray && typeof(cb.angy)<:AbstractArray && typeof(cb.x)<:Real && typeof(cb.y)<:Real
        R0=[[cb.x.+eps(),cb.y.+eps(),cb.zpos.+eps()]  for i in 1:length(cb.angx),j in 1:length(cb.angy)]
        D0=[[sin(cb.angx[i]),sin(cb.angy[j]),sqrt(1-sin(cb.angx[i])^2-sin(cb.angy[j])^2)] for i in 1:length(cb.angx),j in 1:length(cb.angy)]
    else
        error("First four attributes only allowed in one of two configurations for bundle: \n
            1. x and y are Type{AbstractArray{T,1} where T<:Real}, angx and angy are Type{T} where T<:Real \n
            2. x and y are Type{T} where T<:Real, angx and angy are Type{AbstractArray{T,1} where T<:Real} ")
    end
    return R0,D0
end

"""
    bundle_as_array(x, y, angx, angy, zpos)

Composite type representing the origination of a set of rays. First four arguments are
Type `AbstractMatrix{T} where T<:Real`; last argument is a scalar for the z position.

When a variable is declared this type with specific parameters and called
as a function with no arguments, it returns a tuple (position, direction) each
entry of which is an output array of vectors corresponding to the input information.

# Examples
```julia-repl
julia> eval_bundle = bundle_as_array([0.5 1.], [0.5 1.], [0. 0.], [0. 0.], 0.)
bundle_as_array{Float64}([0.5 1.0], [0.5 1.0], [0.0 0.0], [0.0 0.0], 0.0)

julia> eval_bundle()
([[0.5000000000000002, 0.5000000000000002, 2.220446049250313e-16] [1.0000000000000002, 1.0000000000000002, 2.220446049250313e-16]], [[0.0, 0.0, 1.0] [0.0, 0.0, 1.0]])
```
"""
struct bundle_as_array{T<: Real}
    x::AbstractArray{T,2}
    y::AbstractArray{T,2}
    angx::AbstractArray{T,2}# these are direction tangents.
    angy::AbstractArray{T,2}
    zpos:: Union{AbstractArray{T,2},Real}
end
function (cb::bundle_as_array)()

    if length(size(cb.zpos))==2
        R0=[[cb.x[i,j].+eps(),cb.y[i,j].+eps(),cb.zpos[i,j].+eps()] for i in 1:size(cb.x)[1],j in 1:size(cb.y)[2]]
    else
        R0=[[cb.x[i,j].+eps(),cb.y[i,j].+eps(),cb.zpos.+eps()] for i in 1:size(cb.x)[1],j in 1:size(cb.y)[2]]
    end
    D0=[[sin(cb.angx[i,j]),sin(cb.angy[i,j]),sqrt(1-sin(cb.angx[i,j])^2-sin(cb.angy[i,j])^2)] for i in 1:size(cb.x)[1],j in 1:size(cb.y)[2]]
    return R0,D0
end

"""
    bundle_as_array_big(x, y, z, Dx, Dy, Dz)

Composite type representing the origination of a set of rays. All arguments are Type `AbstractArray{T,2} where T<:Real`.

"""
struct bundle_as_array_big{T<: Real}
    x::AbstractArray{T,2}
    y::AbstractArray{T,2}
    z::AbstractArray{T,2}
    Dx::AbstractArray{T,2}# these are direction tangents.
    Dy::AbstractArray{T,2}
    Dz::AbstractArray{T,2}
end
function (cb::bundle_as_array_big)()
    R0=[[cb.x[i,j].+eps(),cb.y[i,j].+eps(),cb.z[i,j].+eps()] for i in 1:size(cb.x)[1],j in 1:size(cb.y)[2]]
    D0=[[cb.Dx[i,j],cb.Dy[i,j],cb.Dz[i,j]] for i in 1:size(cb.x)[1],j in 1:size(cb.y)[2]]
    return R0,D0
end

"""
    trace(xr, yr, zr, nsr)

Composite type holding position and accumulated optical path for parsing. Each is Type `Matrix{Vector{T}} where T<:Real`,
vector holding information along propagation and the matrix mapping to the chosen input bundle subscripting.
"""
struct trace{T<:Real}
    xr::Matrix{Vector{T}}
    yr::Matrix{Vector{T}}
    zr::Matrix{Vector{T}}
    nsr::Matrix{Vector{T}}
end

"""
    bigtrace(xr, yr, zr, Dxr, Dyr, Dzr, nsr)

Composite type holding position, direction, and accumulated optical path for parsing. Each is Type `Matrix{Vector{T}} where T<:Real`,
vector holding information along propagation and the matrix mapping to the chosen input bundle subscripting.

"""
struct bigtrace{T<:Real}
    xr::Matrix{Vector{T}}
    yr::Matrix{Vector{T}}
    zr::Matrix{Vector{T}}
    Dxr::Matrix{Vector{T}}
    Dyr::Matrix{Vector{T}}
    Dzr::Matrix{Vector{T}}
    nsr::Matrix{Vector{T}}
end

"""
    bigtrace_to_bundle(tracething::bigtrace, pos_prop::Int)

Takes a `bigtrace` object and integer position `pos_prop` and returns a bundle object.

"""
function bigtrace_to_bundle(bT::bigtrace, pos_prop::Int)
    x=[(bT.xr[i,j][pos_prop] +eps()) for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    y=[(bT.yr[i,j][pos_prop] +eps()) for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    z=[(bT.zr[i,j][pos_prop] +eps()) for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dx=[bT.Dxr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dy=[bT.Dyr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dz=[bT.Dzr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    return bundle_as_array_big(x,y,z,Dx,Dy,Dz)
end

"""
    trace_extract_terminus(tracething::Union{trace, bigtrace}, pos_prop::Int; coord="all")

Takes a `trace` or `bigtrace` object and integer position pos_prop and returns matrices of coordinates according to the bundle.

Coordinates are `x,y,z`, direction cosines are `Dx,Dy,Dz` and accumulated optical path is `ns`.

Optional argument coord can be set "all" to return all information as Type `Matrix`.  This will be passed out as
- x, y, z, ns for `trace` object
- x, y, z, Dx, Dy, Dz, ns for `bigtrace` object

Or optional argument can be set to a different string to only pass out one `Matrix` of information.
- "x", "y", "z", "ns" for `trace` object
- "x", "y", "z", "Dx", "Dy", "Dz", "ns" for `bigtrace` object

"""
function trace_extract_terminus(bT::bigtrace, pos_prop;coord="all")
    x=[bT.xr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    y=[bT.yr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    z=[bT.zr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dx=[bT.Dxr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dy=[bT.Dyr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    Dz=[bT.Dzr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    ns=[bT.nsr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    if coord=="all"
        return x,y,z,Dx,Dy,Dz
    elseif coord=="x"
        return x
    elseif coord=="y"
        return y
    elseif coord=="z"
        return z
    elseif coord=="Dx"
        return Dx
    elseif coord=="Dy"
        return Dy
    elseif coord=="Dz"
        return Dz
    elseif coord=="ns"
        return ns
    end
end
function trace_extract_terminus(bT::trace, pos_prop;coord="all")
    x=[bT.xr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    y=[bT.yr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    z=[bT.zr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    ns=[bT.nsr[i,j][pos_prop] for i in 1:size(bT.xr)[1], j in 1:size(bT.xr)[2]]
    if coord=="all"
        return x,y,z
    elseif coord=="x"
        return x
    elseif coord=="y"
        return y
    elseif coord=="z"
        return z
    elseif coord=="ns"
        return ns
    end
end

"""
    trace_extract_ray(tracething::Union{trace, bigtrace}, i::Int, j::Int; coord = "all")

Takes a `trace` or `bigtrace` object and integers i,j and returns the ray at position in the traced bundle array i,j.

Coordinates are `x, y, z`, direction cosines are `Dx, Dy, Dz` and accumulated optical path is `ns`.

Optional argument coord can be set "all" to return all information as Matrices.  This will be passed out as
- x, y, z, ns for `trace` object
- x, y, z, Dx, Dy, Dz, ns for `bigtrace` object

Or optional argument can be set to a different string to only pass out one matrix of information.
- "x", "y", "z", "ns" for `trace` object
- "x", "y", "z", "Dx", "Dy", "Dz", "ns" for `bigtrace` object

"""
function trace_extract_ray(bT::bigtrace,i::Int,j::Int;coord="all")
    x=bT.xr[i,j]
    y=bT.yr[i,j]
    z=bT.zr[i,j]
    Dx=bT.Dxr[i,j]
    Dy=bT.Dyr[i,j]
    Dz=bT.Dzr[i,j]
    ns=bT.nsr[i,j]
    if coord=="all"
        return x,y,z,Dx,Dy,Dz
    elseif coord=="x"
        return x
    elseif coord=="y"
        return y
    elseif coord=="z"
        return z
    elseif coord=="Dx"
        return Dx
    elseif coord=="Dy"
        return Dy
    elseif coord=="Dz"
        return Dz
    elseif coord=="ns"
        return ns
    end
end
function trace_extract_ray(bT::trace,i,j;coord="all")
    x=bT.xr[i,j]
    y=bT.yr[i,j]
    z=bT.zr[i,j]
    ns=bT.nsr[i,j]
    if coord=="all"
        return x,y,z
    elseif coord=="x"
        return x
    elseif coord=="y"
        return y
    elseif coord=="z"
        return z
    elseif coord=="ns"
        return ns
    end
end

struct diffSf{impSurf,volPhase}
    impSurf::impSurf
    volPhase::volPhase # volPhase(x,y,z) is the definition here. grating period is defined as the deriv phase fn divided by 2pi
end
function (rf::diffSf)(R0,D0,m,lambda;isReorient=true)
    nvf=normalVec(rf.impSurf)
    N=nvf(R0)
    vpc(R)=rf.volPhase(R[1],R[2],R[3])
    gradPhase(R)=ForwardDiff.gradient(vpc,R)
    kvec=cross(cross(N,gradPhase(R0)),N)
    D1=D0+kvec*lambda/2/pi*m
    D1./sqrt(dot(D1,D1))
end

function propSingleGRIN2(R,D,indexRef,indexRefGrad;dt=0.1)
    # Appl.Opt. Sharma et al 1982
    n = indexRef(R)
    ∇n = indexRefGrad(R)
    TT = n.*D
    A = dt.*∇n .*n
    RB = R .+ dt./2. .*TT.+dt./8 .*A
    nB = indexRef(RB)
    ∇nB = indexRefGrad(RB)
    B = dt.*∇nB .*nB
    RC = R .+ dt .* TT .+ dt./ 2. .* B
    nC = indexRef(RC)
    ∇nC = indexRefGrad(RC)
    C = dt.*∇nC .*nC
    Tnew = TT .+ 1. ./ 6. .* (A .+ 4. .* B .+ C)
    Ro=R .+ dt*(TT .+ 1. ./ 6. .* (A.+ 2. .* B))
    Dro=Tnew./n
    return Ro,Dro
end

"""
    opticalstack_instance = opticalstack( ... )

    ... = opticalstack_instance( ... )

For surface reflection and refraction use this method to define an instance of the `opticalstack` type

    opticalstack(coeffslist, surfslist, nlist)

For gradient index (and surface reflection and refraction) use this method to define an instance of the `opticalstack` type

    opticalstack(coeffslist, surfslist, nlist, ncoeffslist, grindt)

For grating/hologram/metasurface (and surface reflection and refraction) use this method to define an instance of the `opticalstack` type

    opticalstack(coeffslist, surfslist, nlist, diffractiveslist, diffcoeffslist, diffmlist, wavelength)

To use all features together use this method to define an instance of the `opticalstack` type

    opticalstack(coeffslist, surfslist, nlist, diffractiveslist, diffcoeffslist, diffmlist, wavelength, ncoeffslist, grindt)

Arguments to create an instance of `opticalstack` are

- `coeffslist` is a `Vector` set of coefficients for a corresponding surface; each set of coefficients is itself also a `Vector`.
- `surfslist` is a `Vector` of functions callable on `(x, y, coeffs)` to return a z-zurface position.
- `nlist` is a `Vector` of index values or index functions corresponding to the medium before the corresponding surface.
    - Vector values T<:Real work universally
    - Vector input of mixed type, `Real` and functions `n` where `n(x, y, z, coeffs) = ...`, only works for declarations of `opticalstack` which define `ncoeffslist` and `grindt` as well.

- `ncoeffslist::Vector{Vector{T}} where T<:Real` is a `Vector` set of coefficients for a corresponding gradient index; each set of coefficients is itself also a `Vector`.
- `grindt::Real` is the Sharma ray increment in the gradient index medium.
    - This step parameter scales the physical quantity ``OP/n²`` (optical path length and index of refraction) at every step, and raytracing is by the algorithm of Sharma et al. 1982.

-  `diffractiveslist` is a `Vector` of zeros and functions. Where a vector entry is not zero, the entry must be a valid phase function, e.g. ``ϕ`` where ``ϕ(x, y, z, coeffs) = ...`` sampled by the surface at a corresponding position to make the diffractive.
- `diffcoeffslist::Vector{Vector{T}} where T<:Real` is a `Vector` set of coefficients for a corresponding diffractive; each set of coefficients is itself also a `Vector`.
- `diffmlist::Vector{Int}` must have the integer corresponding to the diffraction order in the vector position corresponding to the position of the diffractive surface.
- `wavelength::Real` is the wavelength in system units simulated for the definition of the diffractive.

When an instance is declared, `opticalstack_instance = opticalstack( ... )`, the declared instance is called on a bundle instance `bundlething` and optional arguments as

    opticalstack_instance(bundlething; rend = nothing, color = :red, issurfplot = true, numsurfpoints = 20, halfdomain = nothing, xdom = nothing, ydom = nothing, plobj = nothing, isbigtrace = false, plotplane_at = 0., init_prop_guesses=nothing)

- When `rend` is `nothing` and `isbigtrace = false` the output is a `trace` type.
- When `rend` is `nothing` and `isbigtrace = true` the output is a `bigtrace` type.
- When `rend` is not `nothing`, `rend` must be `"YZ"` or `"XZ"` or `"3Dcirc"` or `"3Dsq"`. The output will be a plot object compatible with Plots.jl
- When `rend` is not `nothing` and `issurfplot = false` the surfaces being modeled will not be plotted.

- The domain plot can be scaled globally by entering a vector of desired plot extent into `xdom` or `ydom`. This is preferred over using `numsurfpoints` and `halfdomain.`
- `plotplane_at` is the value at which a 2D plot is sliced in the non-plotted dimension.
- `init_prop_guesses::Vector{T} where T<:Real` is an option to declare a first guess in numerical propagation preceding surface in the corresponding vector entry.
    - Defaults to zero.
    - Not recommended to use a nonzero value when traversing a gradient index function.

__Note__ that the data scaling of the output of gradient index analysis is not only with the number of rays but also with the resolution of propagation.  Assuming
that propagation occurs over some position difference ``Δz``, this data scaling is APPROXIMATELY ``Δz``/`dt`. As expected, `dt` trades precision for computation time.

__Note__ for gradient index surfaces the output ray coordinates are stacked according to propagation needed for the "longest" gradient index ray. Therefore a returned `trace` or `bigtrace` will not be easily indexed according to its surface, but all propagation positions will properly align when entering and exiting the gradient index.
Shorter rays will replicate their terminated values until the longest ray terminates. Then all rays will exit the gradient index together at the appropriate surface.

__Note__ for diffractive surfaces, the optical path contribution of the corresponding surface is that due to its preceding propagation PLUS that due to the diffraction phase function.

__Note__ that the index of refraction must be equal on both sides of the diffractive.  To work with cases where the diffractive sits on a substrate of different index of refraction, two surfaces should be collocated in sequence, one for reflection/refraction and one for diffraction.

# Examples
```julia-repl
julia> opticalstack_test = opticalstack([[1., 6.], [2., -6.], [4., 8.]], [spherical, spherical, zernike], [1., 1.5, 1.])
opticalstack{Float64}([[1.0, 6.0], [2.0, -6.0], [4.0, 8.0]], UnionAll[spherical, spherical, zernike], [1.0, 1.5, 1.0])

julia> bundle_test = bundle([0.5, 1.], [0.5, 1.], 0., 0., 0.)
bundle([0.5, 1.0], [0.5, 1.0], 0.0, 0.0, 0.0)

julia> tracething = opticalstack_test(bundle_test)
trace{Float64}([[0.5, 0.5, 0.4742890889441142, -0.024645519233983726] ... ; ...])
```

"""
struct opticalstack{T<:Real}
    coeffslist::Vector{Vector{T}}
    surfslist::Array{S,1} where S<:Any #surfslist::Array{UnionAll,1}
    nlist::Vector{S} where S<:Any
    diffractiveslist::Union{Vector{S} where S<:Any,Nothing}
    diffcoeffslist::Union{Vector{Vector{T}},Nothing}
    diffmlist::Union{Vector{Int},Nothing}
    wavelength::Union{T,Nothing}
    ncoeffslist::Union{Vector{Vector{T}},Nothing}
    grindt::Union{T,Nothing}
end
opticalstack(
    coeffslist::Vector{Vector{T}} where T<:Real,
    surfslist::Array{S,1} where S<:Any,
    nlist::Vector{S} where S<:Any)=opticalstack(
        coeffslist,
        surfslist,
        nlist,
        nothing,nothing,nothing,nothing,nothing,nothing)
opticalstack(
    coeffslist::Vector{Vector{T}} where T<:Real,
    surfslist::Array{S,1} where S<:Any,
    nlist::Vector{S} where S<:Any,
    diffractiveslist::Union{Vector{S} where S<:Any,Nothing},
    diffcoeffslist::Union{Vector{Vector{T}},Nothing} where T<:Real,
    diffmlist::Union{Vector{Int},Nothing},
    wavelength::Union{T,Nothing} where T<:Real) = opticalstack(
        coeffslist,
        surfslist,
        nlist,
        diffractiveslist,
        diffcoeffslist,
        diffmlist,
        wavelength,
        nothing,nothing)
opticalstack(
    coeffslist::Vector{Vector{T}} where T<:Real,
    surfslist::Array{S,1} where S<:Any,
    nlist::Vector{S} where S<:Any,
    ncoeffslist::Union{Vector{Vector{T}},Nothing} where T<:Real,
    grindt::Union{T,Nothing} where T<:Real) = opticalstack(
        coeffslist,
        surfslist,
        nlist,
        nothing,nothing,nothing,nothing,
        ncoeffslist,
        grindt)
function (optSta::opticalstack)(bundleEv;rend=nothing,color=:red,issurfplot=true,numsurfpoints=20,halfdomain=nothing,xdom=nothing,ydom=nothing,plobj=nothing,isbigtrace=false,plotplane_at=0.,init_prop_guesses=nothing)
    R0list,D0list=bundleEv()
    Rmaster=[R0list]
    Dmaster=[D0list]
    nsmaster=[zeros(size(D0list))]
    nlistApp=reduce(vcat,[optSta.nlist,optSta.nlist[end]]) # add final index to list (equivalent to making sure image surface is immersed.)

    if optSta.diffractiveslist===nothing
        diffList=zeros(length(optSta.coeffslist))
    else
        diffList=optSta.diffractiveslist
    end
    for i in 1:length(optSta.coeffslist)
        impz1=impz(optSta.coeffslist[i],optSta.surfslist[i])
        if init_prop_guesses===nothing
            gi1=getIntersect_preInit(impz1,0.)
        else
            gi1=getIntersect_preInit(impz1,init_prop_guesses[i])
        end

        if !(typeof(nlistApp[i])<:Real)
            nu(R)=nlistApp[i](R[1],R[2],R[3],optSta.ncoeffslist[i])
            ∇nu(R)=ForwardDiff.gradient(nu,R)
            fntraR(R,D)=propSingleGRIN2(R,D,nu,∇nu;dt=optSta.grindt)[1]
            fntraD(R,D)=propSingleGRIN2(R,D,nu,∇nu;dt=optSta.grindt)[2]
            sc=gi1.(Rmaster[end],Dmaster[end])
            while any(sc.>=0)
                Rn=fntraR.(Rmaster[end],Dmaster[end])
                Dn=fntraD.(Rmaster[end],Dmaster[end])
                sc=gi1.(Rn,Dn)
                Rn=Rn.*0.5.*(sign.(sc).+1.) .+ Rmaster[end].*0.5.*(sign.(-sc).+1.)
                Dn=Dn.*0.5.*(sign.(sc).+1.) .+ Dmaster[end].*0.5.*(sign.(-sc).+1.)
                nsmaster=reduce(vcat,[nsmaster,[nsmaster[i].+optSta.grindt.*nu.(Rmaster[end]).^2]])
                Rmaster=reduce(vcat,[Rmaster,[Rn]])
                Dmaster=reduce(vcat,[Dmaster,[Dn]])
            end
        end
        s=gi1.(Rmaster[end],Dmaster[end])
        R1=Rmaster[end].+s.*Dmaster[end]

        if diffList[i]!=0.
            if nlistApp[i]!=nlistApp[i+1]; @error "refractive indices must be equal on either side of a diffractive surface function";end
            if nlistApp[i]!=1.0; @info "diffractive is not immersed in n=1, be careful that phase function definition compensates for this"; end
            diffUse(x,y,z)=diffList[i](x,y,z,optSta.diffcoeffslist[i])
            ds1=diffSf(impz1,diffUse)
            D1=ds1.(R1,Dmaster[end],optSta.diffmlist[i],optSta.wavelength)
        else
            rs1=refrSf(impz1)
            if typeof(nlistApp[i+1])<:Real
                n2=nlistApp[i+1]
            else
                if length(optSta.ncoeffslist)==i; grC=optSta.ncoeffslist[i]; else; grC=optSta.ncoeffslist[i+1]; end
                nR2(R)=nlistApp[i+1](R[1],R[2],R[3],grC)
                n2=nR2.(R1)
            end
            if typeof(nlistApp[i])<:Real
                n1=nlistApp[i]
            else
                nR1(R)=nlistApp[i](R[1],R[2],R[3],optSta.ncoeffslist[i])
                n1=nR1.(R1)
            end
            D1=rs1.(R1,Dmaster[end],n1,n2)
        end
        if diffList[i]!=0.
            diffUseR(R)=diffUse(R[1],R[2],R[3])
            nsmaster=reduce(vcat,[nsmaster,[nsmaster[end].+nlistApp[i].*(s).+ 1. ./ 2 ./ pi .* diffUseR.(Rmaster[end])]])
        else
            nsmaster=reduce(vcat,[nsmaster,[nsmaster[end].+n1.*(s)]]) # this syntax is more diff-able
        end
        Rmaster=reduce(vcat,[Rmaster,[R1]])
        Dmaster=reduce(vcat,[Dmaster,[D1]])
    end
    Rm=Rmaster
    nsm=nsmaster
    xr=[[Rm[i][j,k][1] for i in 1:length(Rm)] for j in 1:size(Rm[1])[1],k in 1:size(Rm[1])[2]]
    yr=[[Rm[i][j,k][2] for i in 1:length(Rm)] for j in 1:size(Rm[1])[1],k in 1:size(Rm[1])[2]]
    zr=[[Rm[i][j,k][3] for i in 1:length(Rm)] for j in 1:size(Rm[1])[1],k in 1:size(Rm[1])[2]]
    if isbigtrace
        Dm=Dmaster
        Dxr=[[Dm[i][j,k][1] for i in 1:length(Rm)] for j in 1:size(Dm[1])[1],k in 1:size(Dm[1])[2]]
        Dyr=[[Dm[i][j,k][2] for i in 1:length(Rm)] for j in 1:size(Dm[1])[1],k in 1:size(Dm[1])[2]]
        Dzr=[[Dm[i][j,k][3] for i in 1:length(Rm)] for j in 1:size(Dm[1])[1],k in 1:size(Dm[1])[2]]
    end
    nsr=[[nsm[i][j,k] for i in 1:length(nsm)] for j in 1:size(nsm[1])[1],k in 1:size(nsm[1])[2]]

    halfdomainl=0
    if halfdomain===nothing
        if rend=="3Dcirc"
            halfdomainv_mx=maximum.(Rm).*1.41
            halfdomainv_mn=minimum.(Rm).*1.41
        else
            halfdomainv_mx=maximum.(Rm)
            halfdomainv_mn=minimum.(Rm)
        end
        halfdomainl=[max(halfdomainv_mx[j][1],halfdomainv_mx[j][2],halfdomainv_mn[j][1],halfdomainv_mn[j][2]) for j in 1:length(halfdomainv_mx)]
    end

    if rend!==nothing
        if plobj===nothing
            plobj=Plots.plot()
        end
        if rend=="YZ"
            if issurfplot
                for ij in 1:length(optSta.surfslist)
                    zmprof(y)=optSta.surfslist[ij](plotplane_at,y,optSta.coeffslist[ij])
                    if halfdomainl!=0; halfdomain=halfdomainl[ij+1]; end
                    res=2*halfdomain/(numsurfpoints-1)
                    if ydom===nothing; ydom=-halfdomain:res:halfdomain; end
                    Plots.plot!(plobj,zmprof.(ydom),ydom,color=:black,aspect_ratio=:equal)
                end
            end
            Plots.plot!(plobj,zr[:],yr[:],color=color,legend=:none,xlabel="z",ylabel="y")
        elseif rend=="XZ"
            if issurfplot
                for ij in 1:length(optSta.surfslist)
                    zmprof(x)=optSta.surfslist[ij](x,plotplane_at,optSta.coeffslist[ij])
                    if halfdomainl!=0; halfdomain=halfdomainl[ij+1]; end
                    res=2*halfdomain/(numsurfpoints-1)
                    if xdom===nothing; xdom=-halfdomain:res:halfdomain; end
                    Plots.plot!(plobj,zmprof.(xdom),xdom,color=:black,aspect_ratio=:equal)
                end
            end
            Plots.plot!(plobj,zr[:],xr[:],color=color,legend=:none,xlabel="z",ylabel="x")
        elseif rend=="3Dsq"
            if issurfplot
                for ij in 1:length(optSta.surfslist)
                    if halfdomainl!=0; halfdomain=halfdomainl[ij+1]; end
                    res=2*halfdomain/(numsurfpoints-1)
                    if xdom===nothing; xdom=-halfdomain:res:halfdomain; end
                    if ydom===nothing; ydom=-halfdomain:res:halfdomain; end
                    x=xdom'.*ones(length(ydom))
                    y=ones(length(xdom))'.*ydom
                    zmprof(x,y)=optSta.surfslist[ij](x,y,optSta.coeffslist[ij])
                    Plots.surface!(plobj,xdom,ydom,zmprof.(x,y),color=:blue,alpha=0.7,legend=:none,colorbar=nothing)
                end
            end
            Plots.plot!(plobj,xr[:],yr[:],zr[:],color=color,legend=:none,xlabel="x",ylabel="y",zlabel="z")
        elseif rend=="3Dcirc"

            if issurfplot
                for ij in 1:length(optSta.surfslist)
                    if halfdomainl!=0; halfdomain=halfdomainl[ij+1]; end
                    res=2*halfdomain/(numsurfpoints-1)
                    if xdom===nothing; xdom=-halfdomain:res:halfdomain; end
                    if ydom===nothing; ydom=-halfdomain:res:halfdomain; end
                    x=xdom'.*ones(length(ydom))
                    y=ones(length(xdom))'.*ydom
                    zmprof(x,y)=optSta.surfslist[ij](x,y,optSta.coeffslist[ij])
                    zp=zmprof.(x,y)
                    rad=halfdomain
                    zp[x.^2 .+ y.^2 .>rad.^2].=NaN
                    Plots.surface!(plobj,xdom,ydom,zp,color=:blue,alpha=0.7,legend=:none,colorbar=nothing)
                end
            end
            Plots.plot!(plobj,xr[:],yr[:],zr[:],color=color,legend=:none,xlabel="x",ylabel="y",zlabel="z")
        end
        return plobj
    else
        if !isbigtrace
            return trace(xr,yr,zr,nsr)
        else
            return bigtrace(xr,yr,zr,Dxr,Dyr,Dzr,nsr)
        end
    end
end

"""
    rms_spot(tracething::Union{trace, bigtrace}; isunsqrt = false, pos = nothing, iszincluded = true)

Evaluates the rms spot size for a trace intersecting a surface.

Pass in the result of a trace, and calcuate the rms spot based on the ray values for the subscript `pos`. If
pos===nothing then the final values of the trace are used (i.e. the final surface is taken to be the evaluation surface).
Option `isunsqrt` removes the square root from the calculation (returns the sum of squared error instead). The z coordinate variation is included in the calculation when `iszincluded == true`.

# Examples
```julia-repl
julia> opticalstack_test = opticalstack([[1., 6.], [2., -6.], [4., 8.]], [spherical, spherical, zernike], [1., 1.5, 1.])
opticalstack{Float64}([[1.0, 6.0], [2.0, -6.0], [4.0, 8.0]], UnionAll[spherical, spherical, zernike], [1.0, 1.5, 1.0])

julia> bundle_test = bundle([0.5, 1.], [0.5, 1.], 0., 0., 0.)
bundle([0.5, 1.0], [0.5, 1.0], 0.0, 0.0, 0.0)

julia> tracething = opticalstack_test(bundle_test)
trace{Float64}([[0.5, 0.5, 0.4742890889441142, -0.024645519233983726] [0.5, 0.5, 0.47749522247588994, -0.04245690108918665]; [1.0, 1.0, 0.9549904449517799, -0.0849138021783733] [1.0, 1.0, 0.9616736859140215, -0.1253539680837904]], [[0.5, 0.5, 0.4742890889441142, -0.024645519233983726] [1.0, 1.0, 0.9549904449517799, -0.0849138021783733]; [0.5, 0.5, 0.47749522247588994, -0.04245690108918665] [1.0, 1.0, 0.9616736859140215, -0.1253539680837904]], [[0.0, 1.0418113625438212, 1.9623895629300963, 7.999999280260749] [0.0, 1.105086411777687, 1.904234300100513, 7.999999273328101]; [0.0, 1.105086411777687, 1.904234300100513, 7.999999273328101] [0.0, 1.169047547709123, 1.8438309411610754, 7.999999266127459]], [[0.0, 1.0418113625438212, 2.9446610481344417, 9.023361779810676] [0.0, 1.105086411777687, 2.8587256850369682, 9.064376469832027]; [0.0, 1.105086411777687, 2.8587256850369682, 9.064376469832027] [0.0, 1.169047547709123, 2.7690064459109145, 9.114214613215701]])

julia> rms_spot(tracething)
0.1125476814135891
```
"""
function rms_spot(tracething;isunsqrt=false,pos=nothing,iszincluded=true)
    # extract the last position of ray information x,y, then compute sse for that data
    xru=tracething.xr[:]
    yru=tracething.yr[:]
    zru=tracething.zr[:]
    if pos===nothing
        xe=[xru[j][end] for j =1:length(xru)]
        ye=[yru[j][end] for j =1:length(yru)]
        ze=[zru[j][end] for j =1:length(zru)]
    else
        xe=[xru[j][pos] for j =1:length(xru)]
        ye=[yru[j][pos] for j =1:length(yru)]
        ze=[zru[j][pos] for j =1:length(zru)]
    end
    xe=xe.-sum(xe)/length(xe)
    ye=ye.-sum(ye)/length(ye)
    ze=ze.-sum(ze)/length(ze)
    if iszincluded
        mse=sum(xe.^2 .+ ye.^2 .+ ze.^2)/length(xe)
    else
        mse=sum(xe.^2 .+ ye.^2)/length(xe)
    end
    if isunsqrt
        return mse
    else
        return sqrt(mse)
    end
end

"""
    splot(tracething::Union{trace, bigtrace}; pos=nothing, pl_splot = nothing, color = :red, markersize = 2)

Plots the spot of a particular trace or bigtrace object at position `pos`. If `pos === nothing`, the last position in the trace is
used. Can optionally pass in a handle to existing plot pl_splot, and change the marker size from default `markersize = 2` and the plot color from default `color = :red` to any
of the conventional options supported by Plots.jl.

"""
function splot(tracething;pos=nothing,pl_splot=nothing,color=:red,markersize=2)

    if pl_splot===nothing
        pl_splot=Plots.plot()
    end
    xru=tracething.xr[:]
    yru=tracething.yr[:]
    if pos===nothing
        xe=[xru[j][end] for j =1:length(xru)]
        ye=[yru[j][end] for j =1:length(yru)]
    else
        xe=[xru[j][pos] for j =1:length(xru)]
        ye=[yru[j][pos] for j =1:length(yru)]
    end
    if size(tracething.xr)[1]==1
        Plots.scatter!(pl_splot,xe,ye,legend=:none,xlims=[-0.1,0.1],xlabel="x",ylabel="y",color=color,markersize=markersize)
    elseif size(tracething.xr)[2]==1
        Plots.scatter!(pl_splot,xe,ye,legend=:none,ylims=[-0.1,0.1],xlabel="x",ylabel="y",color=color,markersize=markersize)
    else
        Plots.scatter!(pl_splot,xe,ye,legend=:none,aspect_ratio=:equal,xlabel="x",ylabel="y",color=color,markersize=markersize)
    end
    return pl_splot
end

"""
    rac(tracething::Union{trace, bigtrace}, pos_in, pos_out; pl_rac = nothing, color = :red, axistype = "y", pmajoraxis = 2)

Plots the chosen physical coordinate `axistype=y` or `x` at an end position `pos_out` of a particular `trace` or `bigtrace` object as a function of the same coordinate at input position `pos_in`. If `pos===nothing`, the last position in the trace is
used. The coordinate plots are stitched together in lines based on the chosen matrix subscript in the `trace` or `bigtrace`.  This matrix subscript is set by the optional argument `pmajoraxis = 2` or `1`. Can optionally pass in a handle to existing plot pl_rac, and change the plot color from default `color=:red` to any
of the conventional options supported by Plots.jl.

NOTE that when tracing a ray fan (a bundle with one dimension of the array having length 1), the `pmajoraxis` must be set to the correct number to span the nonunitary length of the array. Else the plot may not appear to show anything.

"""
function rac(tracething,pos_in,pos_out;pl_rac=nothing,color=:red,axistype="y",pmajoraxis=2)

    if pl_rac===nothing
        pl_rac=Plots.plot()
    end
    if axistype=="x"
        ru=tracething.xr
    elseif axistype=="y"
        ru=tracething.yr
    end

    re_i=[ru[i,j][pos_in] for i =1:size(ru)[1],j =1:size(ru)[2]]
    re_o=[ru[i,j][pos_out] for i =1:size(ru)[1],j =1:size(ru)[2]]

    if pmajoraxis==1
        Plots.plot!(pl_rac,re_i,re_o,legend=:none,xlabel=axistype*" pos. "*string(pos_in),ylabel=axistype*" pos. "*string(pos_out),color=color)
    elseif pmajoraxis==2
        Plots.plot!(pl_rac,re_i',re_o',legend=:none,xlabel=axistype*" pos. "*string(pos_in),ylabel=axistype*" pos. "*string(pos_out),color=color)
    end
    return pl_rac
end

end # module


text:
"""
	getcurrentdate()

Returns the current date as a string conforming to ISO8601 _basic_ format.

This is used to generate filenames in a cross-platform compatible way.
"""
function getcurrentdate()
	return Dates.format(Dates.now(), "yyyymmddHHMMSS")
end

ext: jl


text:
using Random, Tar, Printf, BenchmarkTools
import TranscodingStreams: TranscodingStream
import CodecZlib: GzipCompressor

function fill_tree!(root::String, depth::Int, max_dirs::Int, num_files::Int, file_size::Int)
    if depth <= 0
        return
    end

    # First, create files
    for file_idx in 1:num_files
        fname = randstring(rand(5:20))
        open(joinpath(root, fname), "w") do io
            write(io, rand(UInt8, file_size))
        end
    end

    # Next, create directories and recurse
    for dir_idx in 1:max_dirs
        dname = joinpath(root, randstring(rand(5:20)))
        mkdir(dname)
        fill_tree!(dname, depth-1, max_dirs, num_files, file_size)
    end
end

function many_small_files(dir, depth=3, max_dirs=3, num_files=1000, file_size=1000)
    fill_tree!(dir, depth, max_dirs, num_files, file_size)
end
function some_medium_files(dir, depth=2, max_dirs=3, num_files=20, file_size=5*1000*1000)
    fill_tree!(dir, depth, max_dirs, num_files, file_size)
end
function few_large_files(dir, depth=1, max_dirs=2, num_files=3, file_size=100*1000*1000)
    fill_tree!(dir, depth, max_dirs, num_files, file_size)
end

function time_tar(genfunc::Function; kwargs...)
    mktempdir() do dir
        t_gen = @elapsed genfunc(dir; kwargs...)
        mktempdir() do outdir
            # First, test performance of bundling up just tar
            tarjl_out = joinpath(outdir, "out_tarjl.tar")
            tar_out = joinpath(outdir, "out_tar.tar")
            t_create_tarjl = @benchmark Tar.create($(dir), $(tarjl_out))
            t_create_tar = cd(dir) do
                @benchmark run($(`tar -cf $(tar_out) .`))
            end

            s_tarjl = filesize(tarjl_out)
            s_tar = filesize(tar_out)

            t_list_tarjl = @benchmark Tar.list($(tarjl_out))
            t_list_tar = @benchmark Tar.list($(tar_out))

            # Next, test performance of tar + gzip in all three configurations:
            # Tar.jl piped into gzip
            # Tar.jl chained into GzipCompressor
            # tar with -z
            tarjl_codec_out = joinpath(outdir, "out_tarjl_codec.tar.gz")
            tarjl_gzip_out = joinpath(outdir, "out_tarjl_gzip.tar.gz")
            tar_gzip_out = joinpath(outdir, "out_tar_gzip.tar.gz")
            t_create_tarjl_codec = @benchmark open($(tarjl_codec_out), "w") do io
                Tar.create($(dir), TranscodingStream(GzipCompressor(), io))
            end

            t_create_tarjl_gzip = @benchmark begin
                gzip_proc = open(pipeline(`gzip - `, stdout=$(tarjl_gzip_out)); write=true)
                Tar.create($(dir), gzip_proc)
                close(gzip_proc.in)
            end
            t_create_tar_gzip = cd(dir) do
                @benchmark run($(`tar -czf $(tar_gzip_out) .`))
            end

            s_tarjl_codec = filesize(tarjl_codec_out)
            s_tarjl_gzip = filesize(tarjl_gzip_out)
            s_tar_gzip = filesize(tar_gzip_out)

            # Convert to time in ms
            ms(t) = minimum(t).time/1e6
            @info(@sprintf("%s: %.1fms to generate", string(genfunc), t_gen*1e3),
                ms(t_create_tarjl),
                ms(t_create_tar),
                s_tarjl,
                s_tar,
                ms(t_list_tarjl),
                ms(t_list_tar),
                "GZIP_COMPRESSION",
                ms(t_create_tarjl_codec),
                ms(t_create_tarjl_gzip),
                ms(t_create_tar_gzip),
                s_tarjl_codec,
                s_tarjl_gzip,
                s_tar_gzip,
            )
        end
    end
end

Random.seed!(0)
time_tar(many_small_files)
time_tar(some_medium_files)
time_tar(few_large_files)


text:
using Test: @testset, @test, @test_throws, @test_broken
using CBindingGen

@testset "CBindingGen" begin
	@testset "convert headers" begin
		cvts = convert_header(_ -> false, "test.h", args = ["-I."])
		@test isempty(cvts)

		cvts = convert_header(_ -> true, "test.h", args = ["-I."])
		@test !isempty(cvts)
	end

	cvts = convert_header(_ -> true, "test.h", args = ["-I."])

	i = 1
	check(expr::String, exports::String...) = check(false, expr, exports...)
	function check(broken::Bool, expr::String, exports::String...)
		cvt = cvts[i]
		i += 1

		broken ? (@test_broken cvt.expr == expr) : (@test cvt.expr == expr)
		@test length(cvt.comments) == length(exports)
		for e in exports
			@test e in keys(cvt.comments)
		end
	end

	@testset "converted" begin
		@testset "simple typedefs" begin
			check("𝐣𝐥.@ctypedef T1 𝐣𝐥.Cvoid", "T1")
			check("𝐣𝐥.@ctypedef T2 T1", "T2")

			check("𝐣𝐥.@ctypedef T3 𝐣𝐥.Cchar", "T3")
			check("𝐣𝐥.@ctypedef T4 𝐣𝐥.Cshort", "T4")
			check("𝐣𝐥.@ctypedef T5 𝐣𝐥.Cint", "T5")
			check("𝐣𝐥.@ctypedef T6 𝐣𝐥.Clong", "T6")
			check("𝐣𝐥.@ctypedef T7 𝐣𝐥.Clonglong", "T7")

			check("𝐣𝐥.@ctypedef T8 𝐣𝐥.Cchar", "T8")
			check("𝐣𝐥.@ctypedef T9 𝐣𝐥.Cshort", "T9")
			check("𝐣𝐥.@ctypedef T10 𝐣𝐥.Cint", "T10")
			check("𝐣𝐥.@ctypedef T11 𝐣𝐥.Clong", "T11")
			check("𝐣𝐥.@ctypedef T12 𝐣𝐥.Clonglong", "T12")

			check("𝐣𝐥.@ctypedef T13 𝐣𝐥.Cuchar", "T13")
			check("𝐣𝐥.@ctypedef T14 𝐣𝐥.Cushort", "T14")
			check("𝐣𝐥.@ctypedef T15 𝐣𝐥.Cuint", "T15")
			check("𝐣𝐥.@ctypedef T16 𝐣𝐥.Culong", "T16")
			check("𝐣𝐥.@ctypedef T17 𝐣𝐥.Culonglong", "T17")

			check("𝐣𝐥.@ctypedef T18 𝐣𝐥.Cfloat", "T18")
			check("𝐣𝐥.@ctypedef T19 𝐣𝐥.Cdouble", "T19")
			check("𝐣𝐥.@ctypedef T20 𝐣𝐥.Clongdouble", "T20")

			check("𝐣𝐥.@ctypedef T21 𝐣𝐥.Complex{𝐣𝐥.Cfloat}", "T21")
			check("𝐣𝐥.@ctypedef T22 𝐣𝐥.Complex{𝐣𝐥.Cdouble}", "T22")
			check("𝐣𝐥.@ctypedef T23 𝐣𝐥.Complex{𝐣𝐥.Clongdouble}", "T23")
		end

		@testset "enums" begin
			check("""
				𝐣𝐥.@cenum E1 {
					E1_1 = 𝐣𝐥.Cuint(0)
					E1_2 = 𝐣𝐥.Cuint(1)
					E1_3 = 𝐣𝐥.Cuint(2)
				}""", "E1", "E1_1", "E1_2", "E1_3")

			check("""
				𝐣𝐥.@cenum E2 {
					E2_1 = 𝐣𝐥.Culonglong(1)
					E2_2 = 𝐣𝐥.Culonglong(1)
					E2_3 = 𝐣𝐥.Culonglong(11)
					E2_4 = 𝐣𝐥.Culonglong(0)
					E2_5 = 𝐣𝐥.Culonglong(4294967295)
					E2_6 = 𝐣𝐥.Culonglong(281474976710655)
				}""", "E2", "E2_1", "E2_2", "E2_3", "E2_4", "E2_5", "E2_6")
		end

		@testset "structs" begin
			check("""
				𝐣𝐥.@cstruct S1 { }""", "S1")

			check("""
				𝐣𝐥.@cstruct S2 {
					i::𝐣𝐥.Cint
				}""", "S2")

			check("""
				𝐣𝐥.@cstruct S3 {
					s::𝐣𝐥.@cstruct S1
				}""", "S3")

			check("""
				𝐣𝐥.@cstruct S4 {
					𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
						s::𝐣𝐥.@cstruct S3
					}
				}""", "S4")

			check("""
				𝐣𝐥.@cstruct S5 {
					s::𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
						s::𝐣𝐥.@cstruct S3
					}
				}""", "S5")

			check("""
				𝐣𝐥.@cstruct S6 {
					(s, s1)::𝐣𝐥.@cstruct S7 {
						i::𝐣𝐥.Cint
						s::𝐣𝐥.@cstruct S3
					}
					i::𝐣𝐥.Cint
					s2::𝐣𝐥.@cstruct S7
				}""", "S6", "S7")

			check("""
				𝐣𝐥.@cstruct S8 {
					s::𝐣𝐥.@cstruct S7
				}""", "S8")

			check("""
				𝐣𝐥.@cstruct S9 {
					(s1, s2)::𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					}
				}""", "S9")

			check("""
				𝐣𝐥.@cstruct S10 {
					𝐣𝐥.@cstruct S11 {
						i::𝐣𝐥.Cint
					}
				}""", "S10", "S11")

			check("""
				𝐣𝐥.@cstruct S12 {
					(s1, s2::𝐣𝐥.Ptr{_}, s3::𝐣𝐥.Ptr{𝐣𝐥.Ptr{_}}, s4::𝐣𝐥.Ptr{𝐣𝐥.Ptr{𝐣𝐥.Ptr{_}}})::𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					}
				}""", "S12")

			check("""
				𝐣𝐥.@cstruct S13 {
					s::𝐣𝐥.Ptr{𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					}}
				}""", "S13")

			check("""
			𝐣𝐥.@cstruct S14 {
				(s1, s2::𝐣𝐥.Ptr{_}, s3::_[2], s4::(_[2])[5], s5::(𝐣𝐥.Ptr{_}[2])[5], s6::(𝐣𝐥.Ptr{_}[2])[])::𝐣𝐥.@cstruct {
					i::𝐣𝐥.Cint
				}
			}""", "S14")

			check("""
				𝐣𝐥.@cstruct S15 {
					s::(𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					})[2]
				}""", "S15")

			check("""
				𝐣𝐥.@cstruct S16 {
					i::𝐣𝐥.Cint
					u::𝐣𝐥.Cuint
				}""", "S16")

			check("""
				𝐣𝐥.@cstruct S17 {
					c::𝐣𝐥.Cchar
					i::𝐣𝐥.Cint
				}""", "S17")

			check("""
				𝐣𝐥.@cstruct S18 {
					c::𝐣𝐥.Cchar
					i::𝐣𝐥.Cint
				} __packed__""", "S18")

			check("""
				𝐣𝐥.@cstruct S19 {
					𝐣𝐥.@cstruct {
						c::𝐣𝐥.Cchar
					} __packed__
					𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					} __packed__
				}""", "S19")

			check("""
				𝐣𝐥.@cstruct S20 {
					𝐣𝐥.@cstruct {
						c::𝐣𝐥.Cchar
					}
					𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					}
				} __packed__""", "S20")

			check("""
				𝐣𝐥.@cstruct S21 {
					f::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Cvoid, 𝐣𝐥.Tuple{}, 𝐣𝐥.CDECL}}
				}""", "S21")

			check("""
				𝐣𝐥.@cstruct S22 {
					f::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Ptr{𝐣𝐥.Cfunction{(𝐣𝐥.@cunion {
						i::𝐣𝐥.Cint
						f::𝐣𝐥.Cfloat
						𝐣𝐥.@cstruct {
							j::𝐣𝐥.Cint
						}
					}), 𝐣𝐥.Tuple{𝐣𝐥.Cint, 𝐣𝐥.Cfloat}, 𝐣𝐥.CDECL}}, 𝐣𝐥.Tuple{(𝐣𝐥.@cstruct S22), 𝐣𝐥.Cint}, 𝐣𝐥.CDECL}}
				}""", "S22")

			check("""
				𝐣𝐥.@cstruct S23 {
					f::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Cint, 𝐣𝐥.Tuple{𝐣𝐥.Ptr{𝐣𝐥.Cconst(𝐣𝐥.Cchar)}, 𝐣𝐥.Vararg}, 𝐣𝐥.CDECL}}
				}""", "S23")

			check("""
				𝐣𝐥.@cstruct S24""", "S24")

			check("""
				𝐣𝐥.@cstruct _S24""", "_S24")

			check("""
				𝐣𝐥.@cstruct S25 {
					s1::𝐣𝐥.Ptr{𝐣𝐥.@cstruct S24}
					s2::𝐣𝐥.Ptr{𝐣𝐥.@cstruct _S24}
					s3::𝐣𝐥.Ptr{𝐣𝐥.@cstruct S25}
				}""", "S25")

			check("""
				𝐣𝐥.@cstruct S24 {
					s1::𝐣𝐥.Ptr{𝐣𝐥.@cstruct S24}
					s2::𝐣𝐥.Ptr{𝐣𝐥.@cstruct _S24}
					s3::𝐣𝐥.@cstruct S25
				}""", "S24")

			check("""
				𝐣𝐥.@cstruct S26 {
					(e1, e2::𝐣𝐥.Ptr{_})::𝐣𝐥.@cenum E3 {
						E3_1 = 𝐣𝐥.Cuint(0)
					}
					e3::𝐣𝐥.@cenum E3
					(e4, e5::𝐣𝐥.Ptr{_})::𝐣𝐥.@cenum E1
				}""", "S26", "E3", "E3_1")
		end

		@testset "typedefs" begin
			check("""
				𝐣𝐥.@ctypedef S27 𝐣𝐥.@cstruct S27_ {
					i::𝐣𝐥.Cint
				}""", "S27", "S27_")

			check("""
				𝐣𝐥.@ctypedef S28 𝐣𝐥.@cstruct {
					i::𝐣𝐥.Cint
				}""", "S28")

			check("""
				𝐣𝐥.@cstruct S29""", "S29")

			check("""
				𝐣𝐥.@ctypedef S29_ 𝐣𝐥.@cstruct S29""", "S29_")

			check("""
				𝐣𝐥.@ctypedef S29_ptr 𝐣𝐥.Ptr{𝐣𝐥.@cstruct S29}""", "S29_ptr")

			check("""
				𝐣𝐥.@cstruct S29 {
					s1::𝐣𝐥.Ptr{S29_}
					s2::S29_ptr
				}""", "S29")

			check("""
				𝐣𝐥.@cstruct S30""", "S30")

			check("""
				𝐣𝐥.@ctypedef F1 𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Cvoid, 𝐣𝐥.Tuple{𝐣𝐥.@cstruct S30}, 𝐣𝐥.CDECL}}""", "F1")

			check("""
				𝐣𝐥.@cstruct S30 {
					(i, j, k)::𝐣𝐥.Cint
					f::F1
				}""", "S30")

			check("""
				𝐣𝐥.@ctypedef F2 𝐣𝐥.Ptr{𝐣𝐥.Cfunction{(𝐣𝐥.@cstruct F2ret {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}), 𝐣𝐥.Tuple{𝐣𝐥.@cunion {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}}, 𝐣𝐥.CDECL}}""", "F2", "F2ret")

			check(true, """
				𝐣𝐥.@ctypedef F3 𝐣𝐥.Ptr{𝐣𝐥.Cfunction{(𝐣𝐥.@cstruct F3ret {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}), 𝐣𝐥.Tuple{(𝐣𝐥.@cunion {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}[1]), 𝐣𝐥.Ptr{F3ret}}, 𝐣𝐥.CDECL}}[3]""", "F3", "F3ret")

			check("""
				𝐣𝐥.@cstruct _S31""", "_S31")

			check("""
				𝐣𝐥.@ctypedef S31 𝐣𝐥.@cstruct _S31""", "S31")

			check("""
				𝐣𝐥.@cstruct S32 {
					s::𝐣𝐥.Ptr{S31}
				}""", "S32")

			check("""
				𝐣𝐥.@ctypedef S33 𝐣𝐥.@cstruct S33 {
					i::𝐣𝐥.Cint
				}""", "S33")

			check("""
				𝐣𝐥.@ctypedef S34 𝐣𝐥.@cstruct {
					i::𝐣𝐥.Cconst(𝐣𝐥.Cint)
					p::𝐣𝐥.Ptr{𝐣𝐥.Cconst(𝐣𝐥.Cint)}
					cp::𝐣𝐥.Cconst(𝐣𝐥.Ptr{𝐣𝐥.Cconst(𝐣𝐥.Cint)})
					s::𝐣𝐥.Ptr{𝐣𝐥.Cconst(S33)}
					cs::𝐣𝐥.Cconst(𝐣𝐥.Ptr{𝐣𝐥.Cconst(S33)})
				}""", "S34")

			check("""
				𝐣𝐥.@ctypedef S35 𝐣𝐥.@cstruct {
					(cp::𝐣𝐥.Ptr{𝐣𝐥.Cconst(_)}, p::𝐣𝐥.Ptr{𝐣𝐥.Cconst(_)}, s::𝐣𝐥.Cconst(_))::𝐣𝐥.@cstruct {
						i::𝐣𝐥.Cint
					}
				}""", "S35")

			check("""
				𝐣𝐥.@ctypedef S36 𝐣𝐥.@cstruct {
					length::𝐣𝐥.Cint
					data::𝐣𝐥.Cchar[]
				}""", "S36")

			check("""
				𝐣𝐥.@cstruct S37 {
					(x, y::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{_, 𝐣𝐥.Tuple{}, 𝐣𝐥.CDECL}}, z::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{_, 𝐣𝐥.Tuple{𝐣𝐥.Cint}, 𝐣𝐥.CDECL}})::𝐣𝐥.Cint
				}""", "S37")

			check("""
				𝐣𝐥.@ctypedef E4 𝐣𝐥.@cenum E4_ {
					E4_1 = 𝐣𝐥.Cuint(0)
				}""", "E4", "E4_", "E4_1")

			check("""
				𝐣𝐥.@ctypedef E5 𝐣𝐥.@cenum {
					E5_1 = 𝐣𝐥.Cuint(0)
				}""", "E5", "E5_1")

			check("""
				𝐣𝐥.@ctypedef E6 𝐣𝐥.@cenum {
					E6_0 = 𝐣𝐥.Cuchar(0)
					E6_1 = 𝐣𝐥.Cuchar(1)
				} __packed__""", "E6", "E6_0", "E6_1")

			check("""
				𝐣𝐥.@ctypedef T24 𝐣𝐥.Cint[2]""", "T24")

			check("""
				𝐣𝐥.@ctypedef T25 (𝐣𝐥.Cint[2])[5]""", "T25")

			check("""
				𝐣𝐥.@ctypedef T26 𝐣𝐥.Ptr{𝐣𝐥.Cconst(𝐣𝐥.Cint)}[2]""", "T26")

			check("""
				𝐣𝐥.@ctypedef T27 (𝐣𝐥.@cstruct {
					i::𝐣𝐥.Cint
				})[3]""", "T27")
		end

		@testset "variables" begin
			check("""
				𝐣𝐥.@cextern g1::𝐣𝐥.Cconst(𝐣𝐥.Cint)""", "g1")

			check("""
				𝐣𝐥.@cextern g2::𝐣𝐥.Cconst(𝐣𝐥.@cstruct S24)""", "g2")

			check("""
				𝐣𝐥.@cextern g3::𝐣𝐥.Cconst(𝐣𝐥.Ptr{𝐣𝐥.@cstruct S24})""", "g3")

			check("""
				𝐣𝐥.@cextern g4::𝐣𝐥.Cint""", "g4")

			check("""
				𝐣𝐥.@cextern g5::𝐣𝐥.@cstruct S24""", "g5")

			check("""
				𝐣𝐥.@cextern g6::𝐣𝐥.Ptr{𝐣𝐥.@cstruct S24}""", "g6")
		end

		@testset "functions" begin
			check("""
				𝐣𝐥.@cextern f1()::𝐣𝐥.Cvoid""", "f1")

			check("""
				𝐣𝐥.@cextern f2()::𝐣𝐥.Cvoid""", "f2")

			check("""
				𝐣𝐥.@cextern f3()::𝐣𝐥.Cint""", "f3")

			check("""
				𝐣𝐥.@cextern f4(i::𝐣𝐥.Cint)::𝐣𝐥.Cvoid""", "f4")

			check("""
				𝐣𝐥.@cextern f5(i::𝐣𝐥.Cint, var"?vararg?"...)::𝐣𝐥.Cint""", "f5")

			check("""
				𝐣𝐥.@cextern f6(p::𝐣𝐥.Ptr{𝐣𝐥.Cconst(𝐣𝐥.Cchar)})::𝐣𝐥.Cvoid""", "f6")

			check("""
				𝐣𝐥.@cextern f7(i1::𝐣𝐥.Cint, i2::𝐣𝐥.Cint, var"?3?"::𝐣𝐥.Cint, var"?4?"::𝐣𝐥.Cint, i5::𝐣𝐥.Cint)::𝐣𝐥.Cvoid""", "f7")

			check("""
				𝐣𝐥.@cextern f8(i1::𝐣𝐥.Cint, var"?2?"::𝐣𝐥.Cint, _2::𝐣𝐥.Cint)::𝐣𝐥.Cvoid""", "f8")

			check("""
				𝐣𝐥.@cextern f9(arg::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Cint, 𝐣𝐥.Tuple{𝐣𝐥.Cint, 𝐣𝐥.Cint}, 𝐣𝐥.CDECL}})::𝐣𝐥.Cvoid""", "f9")

			check("""
				𝐣𝐥.@cextern f10()::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{𝐣𝐥.Cint, 𝐣𝐥.Tuple{𝐣𝐥.Cint, 𝐣𝐥.Cint}, 𝐣𝐥.CDECL}}""", "f10")

			check("""
				𝐣𝐥.@cextern f11::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{(𝐣𝐥.@cstruct {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}), 𝐣𝐥.Tuple{𝐣𝐥.@cunion {
					i::𝐣𝐥.Cint
					f::𝐣𝐥.Cfloat
				}}, 𝐣𝐥.CDECL}}""", "f11")

			check(true, """
				𝐣𝐥.@cextern (g8, f12::𝐣𝐥.Ptr{𝐣𝐥.Cfunction{_, 𝐣𝐥.Tuple{(𝐣𝐥.@cunion {
					i::𝐣𝐥.Cint
					j::𝐣𝐥.Cfloat
				}[2]), (𝐣𝐥.@cunion N1 {
					i::𝐣𝐥.Cint
					j::𝐣𝐥.Cfloat
				}), 𝐣𝐥.Ptr{SG1}}, 𝐣𝐥.CDECL}}, f13(u::(𝐣𝐥.@cunion {
					i::𝐣𝐥.Cint
					j::𝐣𝐥.Cfloat
				}), n::(𝐣𝐥.@cunion N2 {
					i::𝐣𝐥.Cint
					j::𝐣𝐥.Cfloat
				}), s::𝐣𝐥.Ptr{𝐣𝐥.@cstruct SG1})::_)::𝐣𝐥.@cstruct SG1 {
					i::𝐣𝐥.Cint
				}""", "SG1", "N1", "N2", "f12", "f13", "g8")
		end

		@testset "pre-processor directives" begin
			check("""
				macro DEFINE_CHAR() return quote ('\\'') end end
				const DEFINE_CHAR = (@DEFINE_CHAR)""", "@DEFINE_CHAR", "DEFINE_CHAR")

			check("""
				macro DEFINE_STRING() return quote ("\\"") end end
				const DEFINE_STRING = (@DEFINE_STRING)""", "@DEFINE_STRING", "DEFINE_STRING")

			check("""
				macro DEFINE_INT() return quote (1234) end end
				const DEFINE_INT = (@DEFINE_INT)""", "@DEFINE_INT", "DEFINE_INT")

			check("""
				macro DEFINE_OxINT() return quote (0x0000abcd) end end
				const DEFINE_OxINT = (@DEFINE_OxINT)""", "@DEFINE_OxINT", "DEFINE_OxINT")

			check("""
				macro DEFINE_OINT() return quote (0x000001ff) end end
				const DEFINE_OINT = (@DEFINE_OINT)""", "@DEFINE_OINT", "DEFINE_OINT")

			check("""
				macro DEFINE_INTu() return quote (0x000010e1) end end
				const DEFINE_INTu = (@DEFINE_INTu)""", "@DEFINE_INTu", "DEFINE_INTu")

			check("""
				macro DEFINE_FLOAT() return quote (0.1234f0) end end
				const DEFINE_FLOAT = (@DEFINE_FLOAT)""", "@DEFINE_FLOAT", "DEFINE_FLOAT")

			check("""
				macro DEFINE_DOUBLE() return quote (0.1234) end end
				const DEFINE_DOUBLE = (@DEFINE_DOUBLE)""", "@DEFINE_DOUBLE", "DEFINE_DOUBLE")

			check("""
				macro DEFINE_DOT_FLOAT() return quote (0.1234f0) end end
				const DEFINE_DOT_FLOAT = (@DEFINE_DOT_FLOAT)""", "@DEFINE_DOT_FLOAT", "DEFINE_DOT_FLOAT")

			check("""
				macro DEFINE_DOT_DOUBLE() return quote (0.1234) end end
				const DEFINE_DOT_DOUBLE = (@DEFINE_DOT_DOUBLE)""", "@DEFINE_DOT_DOUBLE", "DEFINE_DOT_DOUBLE")

			check("""
				macro DEFINE_FLOAT_DOT() return quote (1234.0f0) end end
				const DEFINE_FLOAT_DOT = (@DEFINE_FLOAT_DOT)""", "@DEFINE_FLOAT_DOT", "DEFINE_FLOAT_DOT")

			check("""
				macro DEFINE_DOUBLE_DOT() return quote (1234.0) end end
				const DEFINE_DOUBLE_DOT = (@DEFINE_DOUBLE_DOT)""", "@DEFINE_DOUBLE_DOT", "DEFINE_DOUBLE_DOT")

			check("""
				macro DEFINE_EXP_FLOAT() return quote (-2140.0f0) end end
				const DEFINE_EXP_FLOAT = (@DEFINE_EXP_FLOAT)""", "@DEFINE_EXP_FLOAT", "DEFINE_EXP_FLOAT")

			check("""
				macro DEFINE_EXP_DOUBLE() return quote (-2140.0) end end
				const DEFINE_EXP_DOUBLE = (@DEFINE_EXP_DOUBLE)""", "@DEFINE_EXP_DOUBLE", "DEFINE_EXP_DOUBLE")

			check("""
				macro DEFINE_INT_TYPE() return quote (\$(𝐣𝐥.Cint)) end end
				const DEFINE_INT_TYPE = (@DEFINE_INT_TYPE)""", "@DEFINE_INT_TYPE", "DEFINE_INT_TYPE")

			check("""
				macro DEFINE_UINT_TYPE() return quote (\$(𝐣𝐥.Cuint)) end end
				const DEFINE_UINT_TYPE = (@DEFINE_UINT_TYPE)""", "@DEFINE_UINT_TYPE", "DEFINE_UINT_TYPE")

			check("""
				macro DEFINE_DEFINED() return quote (@DEFINE_INT_TYPE) end end""", "@DEFINE_DEFINED")

			check("""
				macro E6_0() return quote 0 end end""", "@E6_0")
		end
	end

	@testset "generated" begin
		io = IOBuffer()
		@test (generate(io, nothing => cvts) ; true)
	end
end


text:
### A Pluto.jl notebook ###
# v0.16.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 123c3f17-4a4a-478d-ae0e-5ec7f9d4ad44
begin
    let
        using Pkg
		Pkg.activate(mktempdir())
        Pkg.Registry.update()
        Pkg.add("PlutoUI")
        Pkg.add("Tar")
        Pkg.add("MLDataPattern")
        Pkg.add("Glob")
        Pkg.add("NIfTI")
        Pkg.add("DataAugmentation")
        Pkg.add("CairoMakie")
        Pkg.add("ImageCore")
        Pkg.add("DataLoaders")
        # Pkg.add("CUDA")
        Pkg.add("FastAI")
        Pkg.add("StaticArrays")
    end

    using PlutoUI
    using Tar
    using MLDataPattern
    using Glob
    using NIfTI
    using DataAugmentation
    using DataAugmentation: OneHot, Image
    using CairoMakie
    using ImageCore
    using DataLoaders
    # using CUDA
    using FastAI
    using StaticArrays
end

# ╔═╡ d781f284-16c1-4505-9f5a-dbd90c0334d2
TableOfContents()

# ╔═╡ 6e2a1536-e659-4366-93de-ac7d2e7a27a4
md"""
## Load data
Part of the [Medical Decathlon Dataset](http://medicaldecathlon.com/)
"""

# ╔═╡ 855407b1-b8ce-465f-8531-1af5d014abe9
data_dir = raw"/Users/daleblack/Google Drive/Datasets/Task02_Heart"

# ╔═╡ 4cefe3ca-d7f1-4b2a-b63c-e39d734d2514
function loadfn_label(p)
    a = NIfTI.niread(string(p)).raw
    convert_a = convert(Array{UInt8}, a)
    convert_a = convert_a .+ 1
    return convert_a
end

# ╔═╡ 838a17df-7479-42c2-a2c5-ce626a4016fe
function loadfn_image(p)
    a = NIfTI.niread(string(p)).raw
    convert_a = convert(Array{Float32}, a)
    convert_a = convert_a / max(convert_a...)
    return convert_a
end

# ╔═╡ 11d5a9b4-fdc3-48e4-a6ad-f7dec5fabb8b
begin
    niftidata_image(dir) = mapobs(loadfn_image, Glob.glob("*.nii*", dir))
    niftidata_label(dir) =  mapobs(loadfn_label, Glob.glob("*.nii*", dir))
    data = (
        niftidata_image(joinpath(data_dir, "imagesTr")),
        niftidata_label(joinpath(data_dir, "labelsTr")),
    )
end

# ╔═╡ 778cc0f9-127c-4d4b-a7de-906dfcc29cae
train_files, val_files = MLDataPattern.splitobs(data, 0.8)

# ╔═╡ 349d843a-4a5f-44d7-9371-38c140b9972d
md"""
## Create learning method
"""

# ╔═╡ 019e666e-e2e4-4e0b-a225-b346c7c70939
struct ImageSegmentationSimple <: DLPipelines.LearningMethod
    imagesize
end

# ╔═╡ f7274fa9-8231-44fd-8d00-1c7ab7fc855c
image_size = (112, 112, 96)

# ╔═╡ 9ac18928-9fe4-46ed-ab9c-916791739157
method = ImageSegmentationSimple(image_size)

# ╔═╡ edf2b37a-2775-44c0-8d2e-5e2350b454c4
md"""
### Set up `encode` pipelines
"""

# ╔═╡ 7cf0cc6b-8ff9-4198-9646-9d0787e1013d
begin
  function DLPipelines.encode(
          method::ImageSegmentationSimple,
          context::Training,
          (image, target)::Union{Tuple, NamedTuple}
          )

      tfm_proj = RandomResizeCrop(method.imagesize)
      tfm_im = DataAugmentation.compose(
			ImageToTensor(),
			NormalizeIntensity()
          )
      tfm_mask = OneHot()

      items = Image(Gray.(image)), MaskMulti(target)
      item_im, item_mask = apply(tfm_proj, (items))

      return itemdata(apply(tfm_im, item_im)), itemdata(apply(tfm_mask, item_mask))
  end

  function DLPipelines.encode(
          method::ImageSegmentationSimple,
          context::Validation,
          (image, target)::Union{Tuple, NamedTuple}
          )

      tfm_proj = CenterResizeCrop(method.imagesize)
      tfm_im = DataAugmentation.compose(
          ImageToTensor(),
          NormalizeIntensity()
          )
      tfm_mask = OneHot()

      items = Image(Gray.(image)), MaskMulti(target)
      item_im, item_mask = apply(tfm_proj, (items))

      return itemdata(apply(tfm_im, item_im)), itemdata(apply(tfm_mask, item_mask))
  end
end

# ╔═╡ f2a92ff7-0f94-44d9-aba7-4b7db9fa4a56
begin
	methoddata_train = DLPipelines.MethodDataset(train_files, method, Training())
	methoddata_valid = DLPipelines.MethodDataset(val_files, method, Validation())
end

# ╔═╡ b8b18728-cb5a-445a-809c-986e3965aad3
let
    x, y = MLDataPattern.getobs(methoddata_valid, 1)
    @assert size(x) == (image_size..., 1)
    @assert size(y) == (image_size..., 2)
end

# ╔═╡ bfe051e5-66a3-4b4b-b446-48195d8f3868
md"""
## Visualize
"""

# ╔═╡ 4b7bff28-ec60-4a9d-8b38-648ab871ed16
begin
    x, y = MLDataPattern.getobs(methoddata_valid, 3)
end;

# ╔═╡ 6d1e9ca8-da31-4143-b46f-44b459a2cfc3
size(x)

# ╔═╡ 792a709c-6740-4fc0-b2a3-43ad2ae8035a
size(y)

# ╔═╡ 08be9516-a44c-4fe5-ac46-f586600d586a
@bind b PlutoUI.Slider(1:size(x)[3], default=50, show_value=true)

# ╔═╡ e2670a05-2c09-4bd2-b40c-0cc46fef2344
heatmap(x[:, :, b, 1], colormap=:grays)

# ╔═╡ 427ff0c3-d773-4099-a483-bb8f7d4b8ba1
heatmap(y[:, :, b, 2], colormap=:grays)

# ╔═╡ 616aa780-cbcf-4bf2-b6c5-a984f2530482
md"""
## Dataloader
"""

# ╔═╡ bda7309e-ae97-4ac0-96b6-36a776e9215e
begin
    train_loader = DataLoaders.DataLoader(methoddata_train, 2)
    val_loader = DataLoaders.DataLoader(methoddata_valid, 4)
end

# ╔═╡ 3da0e60e-0196-4538-ab08-49f21b46679a
train_loader

# ╔═╡ e70816b8-4597-4a54-b4f7-735880df6132
val_loader

# ╔═╡ ff962aea-2501-42f5-90bc-72f29deb42af
md"""
## Model
"""

# ╔═╡ 5987b24a-a1da-4d13-89d1-c854de2c3ba0
begin
    # 3D layer utilities
    conv = (stride, in, out) -> Flux.Conv((3, 3, 3), in=>out, stride=stride, pad=Flux.SamePad())
    tran = (stride, in, out) -> Flux.ConvTranspose((3, 3, 3), in=>out, stride=stride, pad=Flux.SamePad())

    conv1 = (in, out) -> Flux.Chain(conv(1, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
    conv2 = (in, out) -> Flux.Chain(conv(2, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
    tran2 = (in, out) -> Flux.Chain(tran(2, in, out), Flux.BatchNorm(out), x -> leakyrelu.(x))
end

# ╔═╡ f451e763-d7d0-4de7-a1f3-76776a194022
begin
    function unet3D(in_chs, lbl_chs)
        # Contracting layers
        l1 = Flux.Chain(conv1(in_chs, 4))
        l2 = Flux.Chain(l1, conv1(4, 4), conv2(4, 16))
        l3 = Flux.Chain(l2, conv1(16, 16), conv2(16, 32))
        l4 = Flux.Chain(l3, conv1(32, 32), conv2(32, 64))
        l5 = Flux.Chain(l4, conv1(64, 64), conv2(64, 128))

        # Expanding layers
        l6 = Flux.Chain(l5, tran2(128, 64), conv1(64, 64))
        l7 = Flux.Chain(Flux.Parallel(+, l6, l4), tran2(64, 32), conv1(32, 32))
        l8 = Flux.Chain(Flux.Parallel(+, l7, l3), tran2(32, 16), conv1(16, 16))
        l9 = Flux.Chain(Flux.Parallel(+, l8, l2), tran2(16, 4), conv1(4, 4))
        l10 = Flux.Chain(l9, conv1(4, lbl_chs))
    end
end

# ╔═╡ ad342428-920a-4414-a4cc-cab281a681dc
model = unet3D(1, 2);

# ╔═╡ d073cbac-18c8-4682-a508-307a619e84fc
md"""
## Helper functions
"""

# ╔═╡ 5726b8bd-e1e7-44eb-8872-a4a8d26be0f9
function dice_metric(ŷ, y)
    dice = 2 * sum(ŷ .& y) / (sum(ŷ) + sum(y))
    return dice
end

# ╔═╡ 7b974fc5-7fc3-45e7-bf6c-48ea4f8eff16
function as_discrete(array, logit_threshold)
    array = array .>= logit_threshold
    return array
end

# ╔═╡ b87d764e-dd31-48b7-bcf1-0538ef5c5e6a
md"""
## Loss functions
"""

# ╔═╡ 6dee15db-5cf9-4bfb-ab42-48f1c2777c04
function dice_loss(ŷ, y)
    ϵ = 1e-5
    return loss = 1 - ((2 * sum(ŷ .* y) + ϵ) / (sum(ŷ .* ŷ) + sum(y .* y) + ϵ))
end

# ╔═╡ bdc1ca32-4847-440f-8a44-98bf7e822803
function hd_loss(ŷ, y, ŷ_dtm, y_dtm)
    M = (ŷ .- y) .^ 2 .* (ŷ_dtm .^ 2 .+ y_dtm .^ 2)
    return loss = mean(M)
end

# ╔═╡ aafb8a3f-6e6c-4475-92db-550eb6999741
md"""
## Training
"""

# ╔═╡ 3d0f1399-54c3-45b6-840a-8e287513bbe7
ps = Flux.params(model);

# ╔═╡ 06061dc5-eae9-47cb-99aa-d521c5cd37dd
loss_function = dice_loss

# ╔═╡ 3a13a9a9-380a-4613-b137-99fdef8cb92f
optimizer = Flux.ADAM(0.01)

# ╔═╡ 36501d5d-d29a-4281-9e75-cebc97bb68cd
begin
  max_epochs = 2
  val_interval = 1
  epoch_loss_values = []
  val_epoch_loss_values = []
  dice_metric_values = []
end

# ╔═╡ ad1bf410-f187-4811-b752-09252e1b3a32
# begin
#   for epoch in 1:max_epochs
#       step = 0
#       @show epoch

#       # Loop through training data
#       for (xs, ys) in train_loader
#           step += 1
#           @show step

#           gs = Flux.gradient(ps) do
#               ŷs = model(xs)
#               loss = loss_function(ŷs[:, :, :, 2, :], ys[:, :, :, 2, :])
#               return loss
#           end
#           Flux.update!(optimizer, ps, gs)
#       end

#       # Loop through validation data
#       if (epoch + 1) % val_interval == 0
#           val_step = 0
#           for (val_xs, val_ys) in val_loader
#               val_step += 1
#               @show val_step

#               local val_ŷs = model(val_xs)
#               local val_loss = loss_function(val_ŷs[:, :, :, 2, :], val_ys[:, :, :, 2, :])
#               val_ŷs, val_ys = as_discrete(val_ŷs, 0.5), as_discrete(val_ys, 0.5)
#           end
#       end
#   end
# end

# ╔═╡ Cell order:
# ╠═123c3f17-4a4a-478d-ae0e-5ec7f9d4ad44
# ╠═d781f284-16c1-4505-9f5a-dbd90c0334d2
# ╟─6e2a1536-e659-4366-93de-ac7d2e7a27a4
# ╠═855407b1-b8ce-465f-8531-1af5d014abe9
# ╠═4cefe3ca-d7f1-4b2a-b63c-e39d734d2514
# ╠═838a17df-7479-42c2-a2c5-ce626a4016fe
# ╠═11d5a9b4-fdc3-48e4-a6ad-f7dec5fabb8b
# ╠═778cc0f9-127c-4d4b-a7de-906dfcc29cae
# ╟─349d843a-4a5f-44d7-9371-38c140b9972d
# ╠═019e666e-e2e4-4e0b-a225-b346c7c70939
# ╠═f7274fa9-8231-44fd-8d00-1c7ab7fc855c
# ╠═9ac18928-9fe4-46ed-ab9c-916791739157
# ╟─edf2b37a-2775-44c0-8d2e-5e2350b454c4
# ╠═7cf0cc6b-8ff9-4198-9646-9d0787e1013d
# ╠═f2a92ff7-0f94-44d9-aba7-4b7db9fa4a56
# ╠═b8b18728-cb5a-445a-809c-986e3965aad3
# ╟─bfe051e5-66a3-4b4b-b446-48195d8f3868
# ╠═4b7bff28-ec60-4a9d-8b38-648ab871ed16
# ╠═6d1e9ca8-da31-4143-b46f-44b459a2cfc3
# ╠═792a709c-6740-4fc0-b2a3-43ad2ae8035a
# ╟─08be9516-a44c-4fe5-ac46-f586600d586a
# ╠═e2670a05-2c09-4bd2-b40c-0cc46fef2344
# ╠═427ff0c3-d773-4099-a483-bb8f7d4b8ba1
# ╟─616aa780-cbcf-4bf2-b6c5-a984f2530482
# ╠═bda7309e-ae97-4ac0-96b6-36a776e9215e
# ╠═3da0e60e-0196-4538-ab08-49f21b46679a
# ╠═e70816b8-4597-4a54-b4f7-735880df6132
# ╟─ff962aea-2501-42f5-90bc-72f29deb42af
# ╠═5987b24a-a1da-4d13-89d1-c854de2c3ba0
# ╠═f451e763-d7d0-4de7-a1f3-76776a194022
# ╠═ad342428-920a-4414-a4cc-cab281a681dc
# ╟─d073cbac-18c8-4682-a508-307a619e84fc
# ╠═5726b8bd-e1e7-44eb-8872-a4a8d26be0f9
# ╠═7b974fc5-7fc3-45e7-bf6c-48ea4f8eff16
# ╟─b87d764e-dd31-48b7-bcf1-0538ef5c5e6a
# ╠═6dee15db-5cf9-4bfb-ab42-48f1c2777c04
# ╠═bdc1ca32-4847-440f-8a44-98bf7e822803
# ╟─aafb8a3f-6e6c-4475-92db-550eb6999741
# ╠═3d0f1399-54c3-45b6-840a-8e287513bbe7
# ╠═06061dc5-eae9-47cb-99aa-d521c5cd37dd
# ╠═3a13a9a9-380a-4613-b137-99fdef8cb92f
# ╠═36501d5d-d29a-4281-9e75-cebc97bb68cd
# ╠═ad1bf410-f187-4811-b752-09252e1b3a32


text:
export set_dynamics!

function set_dynamics!( model::Model, sysf::Function; mm::MassMatrix=I, dx::CallbackArg=:autodiff, du::CallbackArg=:autodiff, dp::CallbackArg=:autodiff )
    # Create/store the mass matrix in the model
    if isa( mm, UniformScaling )
        model.massMatrix = convert( Matrix{Float64}, mm(model.Nx) )
    else
        model.massMatrix = mm
    end

    model.callbacks.ffct_cb = sysf

    if isa( dx, Function )
        model.callbacks.dfdx_vec_cb = dx
    elseif dx == :autodiff

    else
        throw( ArgumentError( "Unknown option $(String(dx)) for dx" ) )
    end

    if isa( du, Function )
        model.callbacks.dfdu_vec_cb = du
    elseif du == :autodiff

    else
        throw( ArgumentError( "Unknown option $(String(du)) for du" ) )
    end

    if isa( dp, Function )
        model.callbacks.dfdp_vec_cb = dp
    elseif dp == :autodiff

    else
        throw( ArgumentError( "Unknown option $(String(dp)) for dp" ) )
    end
end


text:
# mask operator

export joMask
"""
    julia> op = joMask(n,idx;[DDT=...,][RDT=...,][name=...])

Mask operator with index array

# Signature

    function joMask(n::Integer,idx::Vector{VDT};
        DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joMask") where {VDT<:Integer}

# Arguments

- `n`: size of square operator
- `idx`: vector of true indecies
- keywords
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

    A=joMask(3,[1,3])

examples with DDT/RDT

    A=joMask(3,[1,3]; DDT=Float32)
    A=joMask(3,[1,3]; DDT=Float32,RDT=Float64)

"""
function joMask(n::Integer,idx::Vector{VDT};
    DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joMask") where {VDT<:Integer}

    l=length(idx)
    n>=l || throw(joLinearFunctionException("joMask length(idx) must be <= n"))
    n>=max(idx...) || throw(joLinearFunctionException("joMask max(idx) must be <= n"))
    mask=falses(n)
    mask[idx].=true
    return joLinearFunctionFwd(n,n,
        v1->jo_convert(RDT,mask.*v1,false),
        v2->jo_convert(DDT,mask.*v2,false),
        v3->jo_convert(DDT,mask.*v3,false),
        v4->jo_convert(RDT,mask.*v4,false),
        DDT,RDT;
        name=name,
        fMVok=true
        )
end
"""
    julia> op = joMask(mask;[makecopy=...,][DDT=...,][RDT=...,][name=...])

Mask operator with BitArray mask

# Signature

    function joMask(mask::BitArray{1};
        makecopy::Bool=true,DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joMask")

# Arguments
- `mask`: BitArray mask of true indecies
- keywords
    - `makecopy`: copy mask array
    - `DDT`: domain data type
    - `RDT`: range data type
    - `name`: custom name

# Examples

define BitArray mask and operator

    mask=falses(3)
    mask[[1,3]]=true

    A=joMask(mask)

examples with RDT/DDT

    A=joMask(mask; DDT=Float32)
    A=joMask(mask; DDT=Float32,RDT=Float64)

"""
function joMask(mask::BitArray{1};
    makecopy::Bool=true,DDT::DataType=joFloat,RDT::DataType=DDT,name::String="joMask")

    n=length(mask)
    mymask= makecopy ? Base.deepcopy(mask) : mask
    return joLinearFunctionFwd(n,n,
        v1->jo_convert(RDT,mymask.*v1,false),
        v2->jo_convert(DDT,mymask.*v2,false),
        v3->jo_convert(DDT,mymask.*v3,false),
        v4->jo_convert(RDT,mymask.*v4,false),
        DDT,RDT;
        name=name,
        fMVok=true
        )
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, Baseball

makedocs(;
    modules=[Baseball],
    format=Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/benelsen/Baseball.jl/blob/{commit}{path}#L{line}",
    sitename="Baseball.jl",
    authors="Ben Elsen",
    # assets=[],
)

deploydocs(;
    repo="github.com/benelsen/Baseball.jl",
)


text:
# ---
# title: 1643. Kth Smallest Instructions
# id: problem1643
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/kth-smallest-instructions/description/>
# hidden: true
# ---
#
# Bob is standing at cell `(0, 0)`, and he wants to reach `destination`: `(row,
# column)`. He can only travel **right** and **down**. You are going to help Bob
# by providing **instructions** for him to reach `destination`.
#
# The **instructions** are represented as a string, where each character is
# either:
#
#   * `'H'`, meaning move horizontally (go **right** ), or
#   * `'V'`, meaning move vertically (go **down** ).
#
# Multiple **instructions** will lead Bob to `destination`. For example, if
# `destination` is `(2, 3)`, both `"HHHVV"` and `"HVHVH"` are valid
# **instructions**.
#
# However, Bob is very picky. Bob has a lucky number `k`, and he wants the `kth`
# **lexicographically smallest instructions** that will lead him to
# `destination`. `k` is **1-indexed**.
#
# Given an integer array `destination` and an integer `k`, return _the_`kth`
# _**lexicographically smallest instructions**  that will take Bob to
# _`destination`.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/12/ex1.png)
#
#
#
#     Input: destination = [2,3], k = 1
#     Output: "HHHVV"
#     Explanation: All the instructions that reach (2, 3) in lexicographic order are as follows:
#     ["HHHVV", "HHVHV", "HHVVH", "HVHHV", "HVHVH", "HVVHH", "VHHHV", "VHHVH", "VHVHH", "VVHHH"].
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/10/12/ex2.png)**
#
#
#
#     Input: destination = [2,3], k = 2
#     Output: "HHVHV"
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/10/12/ex3.png)**
#
#
#
#     Input: destination = [2,3], k = 3
#     Output: "HHVVH"
#
#
#
#
# **Constraints:**
#
#   * `destination.length == 2`
#   * `1 <= row, column <= 15`
#   * `1 <= k <= nCr(row + column, row)`, where `nCr(a, b)` denotes `a` choose `b`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test, ConvolutionalOperatorLearning
using Random

# Load reference implementation
include(joinpath(@__DIR__,"Reference.jl"))

# Random gaussian images (2D)
@testset "Random gaussian images (2D): vectorized x" begin
	rng = MersenneTwister(1)
	x = [randn(rng,16,16) for _ in 1:4] # vector form

	# 3 x 3 filters
	R = (3,3)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:8
		@testset "3 x 3 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)

			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 4 x 4 filters in list form
	R = (4,4)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:5:15
		@testset "4 x 4 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			h, htrace, objtrace, Hdifftrace =
				CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)

			@test ConvolutionalOperatorLearning._filtermatrix(h)[1]          				== refH
			@test [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in htrace]   == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 2 x 3 filters: Reference test does not support rectangular filters, so
	# test that cost function is decending
	R = (2,3)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 100, 0.0
	for p in 0:2:4
		@testset "2 x 3 filters (p = $p)" begin
			H, Htrace, objtrace, Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test all(diff(objtrace) .<= 100*eps(objtrace[1]))

			H2 = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == H2
		end
	end

	# 4 x 3 filters in vectorized form
	R = (4,3)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 100, 0.0
	for p in 0:2:4
		@testset "4 x 3 filters in vectorized form (p = $p)" begin
			h, htrace, objtrace, Hdifftrace =
				CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)
			@test all(diff(objtrace) .<= 100*eps(objtrace[1]))

			H2 = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test ConvolutionalOperatorLearning._filtermatrix(h)[1] == H2
		end
	end

	# 4 x 1 filters - verify error is thrown for 1d filters with 2d signals
	R = (4,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "1d filters with 2d sigmals (p = $p)" begin
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test_throws AssertionError CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol)
		end
	end
end

# Random overlaid squares (2D); all X's are now in matrix form
@testset "Random overlaid squares (2D): matrix X" begin
	rng = MersenneTwister(1)
	square((x,y),w,n) = [y-w<i<y+w && x-w<j<x+w ? 1.0 : 0.0 for i in 1:n, j in 1:n]
	randsquare(rng,wrange,n,T) = square(
		rand(rng,1+first(wrange):n-first(wrange),2),
		rand(rng,wrange), n)
	x = [sum(randsquare(rng,2:3,16,1) for _ in 1:5) for _ in 1:4]
	X = cat(x..., dims=ndims(x[1])+1) # in matrix form

	# 3 x 3 filters
	R = (3,3)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-2, 4, 0.0
	for p in 0:2:8
		@testset "3 x 3 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 4 x 4 filters in list form
	R = (4,4)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-2, 4, 0.0
	for p in 0:4:15
		@testset "4 x 4 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			h, htrace, objtrace, Hdifftrace =
				CAOL(X,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)

			@test ConvolutionalOperatorLearning._filtermatrix(h)[1]          				== refH
			@test [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in htrace]   == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 4 x 2 filters
	R = (4,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-2, 100, 0.0
	for p in 0:2:6
		@testset "4 x 2 filters (p = $p)" begin
			H, Htrace, objtrace, Hdifftrace =
				CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test all(diff(objtrace) .<= 100*eps(objtrace[1]))

			H2 = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == H2
		end
	end

	# 3 x 4 filters in vectorized form
	R = (3,4)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 100, 0.0
	for p in 0:2:4
		@testset "3 x 4 filters in vectorized form (p = $p)" begin
			h, htrace, objtrace, Hdifftrace =
				CAOL(X,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)
			@test all(diff(objtrace) .<= 100*eps(objtrace[1]))

			H2 = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test ConvolutionalOperatorLearning._filtermatrix(h)[1] == H2
		end
	end

	# 4 x 1 filters - verify error is thrown for 1d filters with 2d signals
	R = (4,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "1d filters with 2d signals (p = $p)" begin
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test_throws AssertionError CAOL(X,λ,h0,maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(X,λ,h0,maxiters=iters,tol=tol)
		end
	end

	# 2 x 2 x 2 filters - verify error is thrown for 3d filters with 2d signals
	R = (2,2,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "2 x 3 filters (p = $p)" begin
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R))
			@test_throws AssertionError CAOL(X,λ,h0,trace=true)
			@test_throws AssertionError CAOL(X,λ,h0,maxiters=iters,tol=tol)
		end
	end
end

# Random gaussian signals (1D)
@testset "Random gaussian signals (1D)" begin
	rng = MersenneTwister(1)
	x = [randn(rng,16) for _ in 1:4]
	X = cat(x..., dims=ndims(x[1])+1) # in matrix form

	# length 3 filters
	R = (3,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:1:2
		@testset "3 x 1 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)

			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# length 3 filters with data in matrix form
	R = (3,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:1:2
		@testset "3 x 1 filters, data in matrix form (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)

			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 6 x 1 filters
	R = (6,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "6 x 1 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 6 x 1 filters in vectorized form
	R = (6,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "6 x 1 vectorized filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			h, htrace, objtrace, Hdifftrace =
				CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)

			@test ConvolutionalOperatorLearning._filtermatrix(h)[1]          				== refH
			@test [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in htrace]   == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 6 x 1 filters in vectorized form and data in matrix form
	R = (6,)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "6 x 1 vectorized filters with data in matrix form (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			h, htrace, objtrace, Hdifftrace =
				CAOL(X,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)

			@test ConvolutionalOperatorLearning._filtermatrix(h)[1]          				== refH
			@test [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in htrace]   == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 2 x 2 filters - verify error is thrown for 2d filters with 1d signals
	R = (2,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "2d filters with 1d signal (p = $p)" begin
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)

			@test_throws AssertionError CAOL(x,λ,h0,maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,h0,maxiters=iters,tol=tol)
			@test_throws AssertionError CAOL(X,λ,h0,maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(X,λ,h0,maxiters=iters,tol=tol)
		end
	end
end

# Random gaussian signals (3D)
@testset "Random gaussian signals (3D)" begin
	rng = MersenneTwister(1)
	x = [randn(rng,16,16,16) for _ in 1:4]

	# 2x3x4 filters in matrix form
	R = (2,3,4)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:5:20
		@testset "2x3x4 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			H, Htrace, objtrace, Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test H          == refH
			@test Htrace     == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 2 x 2 x 2 filters in list form
	R = (2,2,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "4 x 4 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)
			h, htrace, objtrace, Hdifftrace =
				CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)

			@test ConvolutionalOperatorLearning._filtermatrix(h)[1]          				== refH
			@test [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in htrace]   == refHtrace
			@test objtrace   == refobjtrace
			@test Hdifftrace == refHdifftrace

			H = CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end

	# 2 x 2 filters - verify error is thrown for 2d filters with 3d signals
	R = (2,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "2 x 3 filters (p = $p)" begin
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
		end
	end

	# 2x2x2x2 filters - verify error is thrown for 4d filters with 3d signals
	R = (2,2,2,2)
	H0 = generatefilters(:DCT,R,form=:matrix)
	λ, iters, tol = 1e-4, 4, 0.0
	for p in 0:2:4
		@testset "2 x 3 filters (p = $p)" begin
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)
			@test_throws AssertionError CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
		end
	end
end

# need to test:
# + termination condition (make sure all agree)
@testset "Termination condition test" begin
	rng = MersenneTwister(1)
	square((x,y),w,n) = [y-w<i<y+w && x-w<j<x+w ? 1.0 : 0.0 for i in 1:n, j in 1:n]
	randsquare(rng,wrange,n,T) = square(
		rand(rng,1+first(wrange):n-first(wrange),2),
		rand(rng,wrange), n)
	x = [sum(randsquare(rng,2:3,16,1) for _ in 1:5) for _ in 1:4]
	X = cat(x..., dims=ndims(x[1])+1) # in matrix form

	# 3 x 3 filters
	R = (3,3)
	H0 = generatefilters(:DCT,R,form=:matrix)
	h0 = generatefilters(:DCT,R,form=:list)
	λ, iters, tol = 1e-2, 10000, 1e-8
	for p in 0:2:8
		@testset "3 x 3 filters (p = $p)" begin
			refH, refHtrace, refobjtrace, refHdifftrace =
				Reference.CAOL(x,λ,H0[:,p+1:size(H0,2)],R,iters,tol)

			Hx, Hx_Htrace, Hx_objtrace, Hx_Hdifftrace =
				CAOL(x,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)

			HX, HX_Htrace, HX_objtrace, HX_Hdifftrace =
				CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol,trace=true)

			hx, hx_Htrace, hx_objtrace, hx_Hdifftrace =
				CAOL(x,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)
			hx = ConvolutionalOperatorLearning._filtermatrix(hx)[1]
			hx_Htrace = [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in hx_Htrace]

			hX, hX_Htrace, hX_objtrace, hX_Hdifftrace =
				CAOL(X,λ,h0[p+1:end],maxiters=iters,tol=tol,trace=true)
			hX = ConvolutionalOperatorLearning._filtermatrix(hX)[1]
			hX_Htrace = [ConvolutionalOperatorLearning._filtermatrix(ht)[1] for ht in hX_Htrace]

			@test length(hX_Hdifftrace) < iters # otherwise not testing convergence!

			@test refH 		    == Hx 			 == HX            == hX 		   == hx
			@test refHtrace     == Hx_Htrace     == HX_Htrace     == hX_Htrace     == hx_Htrace
			@test refobjtrace   == Hx_objtrace   == HX_objtrace   == hX_objtrace   == hx_objtrace
			@test refHdifftrace == Hx_Hdifftrace == HX_Hdifftrace == hX_Hdifftrace == hx_Hdifftrace

			H = CAOL(X,λ,(H0[:,p+1:size(H0,2)],R),maxiters=iters,tol=tol)
			@test H == refH
		end
	end
end


text:
using DataDeps

register(DataDep(
  "statlog-project-landsat-satellite",
  "http://archive.ics.uci.edu/ml/datasets/Statlog+Project",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/statlog/satimage/sat.trn",
  "e896dc88a960fa2404160fc4c3cb3dc53fcf4afd80ba920bf2d261bd42d12613",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=37,
        feature_indices=1:36,
        delim=' ',
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using HDF5
using Test
using DSGE

path = dirname(@__FILE__)

file = "$path/../reference/solve.h5"
TTT_expected = h5read(file, "TTT")
CCC_expected = h5read(file, "CCC")
RRR_expected = h5read(file, "RRR")

m = AnSchorfheide()
TTT, RRR, CCC = solve(m)

@testset "Check state-space system matches reference" begin
    @test @test_matrix_approx_eq TTT_expected TTT
    @test @test_matrix_approx_eq RRR_expected RRR
    @test @test_matrix_approx_eq CCC_expected CCC
end

m = Model1002()
m <= Setting(:alternative_policy, taylor93())
TTT_alt93, RRR_alt93, CCC_alt93 = solve(m, apply_altpolicy = true)
file = "$path/../reference/solve_alt.h5"
TTT_alt_expected93 = h5read(file, "TTT_alt93")
CCC_alt_expected93 = h5read(file, "CCC_alt93")
RRR_alt_expected93 = h5read(file, "RRR_alt93")

@testset "Check state-space system with altpolicy matches reference for Taylor93" begin
    @test @test_matrix_approx_eq TTT_alt_expected93 TTT_alt93
    @test @test_matrix_approx_eq RRR_alt_expected93 RRR_alt93
    @test @test_matrix_approx_eq CCC_alt_expected93 CCC_alt93
end

m = Model1002()
m <= Setting(:alternative_policy, taylor99())
TTT_alt99, RRR_alt99, CCC_alt99 = solve(m, apply_altpolicy = true)
file = "$path/../reference/solve_alt.h5"
TTT_alt_expected99 = h5read(file, "TTT_alt99")
CCC_alt_expected99 = h5read(file, "CCC_alt99")
RRR_alt_expected99 = h5read(file, "RRR_alt99")

@testset "Check state-space system with altpolicy matches reference for Taylor99" begin
    @test @test_matrix_approx_eq TTT_alt_expected99 TTT_alt99
    @test @test_matrix_approx_eq RRR_alt_expected99 RRR_alt99
    @test @test_matrix_approx_eq CCC_alt_expected99 CCC_alt99
end

nothing


text:
module DiffEqsWrappersTests

using Test

@testset "DiffEqWrappers" begin include("DiffEqsTests.jl") end

end # module


text:
using Test
include("../src/info.jl")

@testset "Mutual Information" begin
    let as = [1,1,1,1,2,2,2,2]
        bs = [2,2,2,2,1,1,1,1]
        @test abs(mutualinfo(as, bs) - 1.0) < 1e-6
        @test abs(mutualinfo(bs, as) - 1.0) < 1e-6
    end

    let as = [1,1,2,2,2,2,1,1,1]
        bs = [2,2,1,1,1,1,2,2,2]
        @test abs(mutualinfo(as, bs) - 0.991076) < 1e-6
        @test abs(mutualinfo(bs, as) - 0.991076) < 1e-6
    end

    let as = [2,2,1,2,1,2,2,2,1]
        bs = [2,2,1,1,1,2,1,2,2]
        @test abs(mutualinfo(as, bs) - 0.072780) < 1e-6
        @test abs(mutualinfo(bs, as) - 0.072780) < 1e-6
    end

    let as = [1,1,1,1,1,1,1,1,1]
        bs = [2,2,2,1,1,1,2,2,2]
        @test abs(mutualinfo(as, bs) - 0.0) < 1e-6
        @test abs(mutualinfo(bs, as) - 0.0) < 1e-6
    end

    let as = [2,2,2,2,1,1,1,1,2]
        bs = [2,2,2,1,1,1,2,2,2]
        @test abs(mutualinfo(as, bs) - 0.072780) < 1e-6
        @test abs(mutualinfo(bs, as) - 0.072780) < 1e-6
    end

    let as = [2,2,1,1,2,2,1,1,2]
        bs = [2,2,2,1,1,1,2,2,2]
        @test abs(mutualinfo(as, bs) - 0.018311) < 1e-6
        @test abs(mutualinfo(bs, as) - 0.018311) < 1e-6
    end

    let xs = rand(1:2, 60), ys = rand(1:2, 60)
        @time for _ in 1:810000
            mutualinfo(xs, ys)
        end
    end
end

@testset "Transfer Entropy" begin
    let as = [2,2,2,1,1], bs = [2,2,1,1,2]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 2/3) < 1e-6
    end

    let as = [1,1,2,2,2,1,1,1,1,2]
        bs = [2,2,1,1,1,1,1,1,2,2]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 0.106844) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 1/2) < 1e-6
    end

    let as = [1,2,1,2,1,1,2,2,1,1]
        bs = [1,1,2,1,2,2,2,1,2,2]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 1/4) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 0.344361) < 1e-6
    end

    let as = [1,1,1,2,2], bs = [1,1,2,2,1]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 2/3) < 1e-6
    end

    let as = [2,2,1,1,1,2,2,2,2,1]
        bs = [1,1,2,2,2,2,2,2,1,1]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 0.106844) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 1/2) < 1e-6
    end

    let as = [2,1,2,1,2,2,1,1,2,2]
        bs = [2,2,1,2,1,1,1,2,1,1]
        @test abs(transferentropy(as, as, 2) - 0.0) < 1e-6
        @test abs(transferentropy(bs, bs, 2) - 0.0) < 1e-6
        @test abs(transferentropy(as, bs, 2) - 1/4) < 1e-6
        @test abs(transferentropy(bs, as, 2) - 0.344361) < 1e-6
    end

    let xs = rand(1:2, 60), ys = rand(1:2, 60)
        @time for _ in 1:810000
            transferentropy(xs, ys, 2)
        end
        @time for _ in 1:810000
            transferentropy(xs, ys, 3)
        end
    end
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
using SupremumInfimum
using Test

@testset "supremum-infimum   " begin include("supinf.jl") end
@testset "check ordering type" begin include("check.jl") end
@testset "integer limits"      begin include("limits.jl") end