text:
export ToeplitzOperator, HankelOperator, LaurentOperator

type ToeplitzOperator{T<:Number} <: Operator{T}
    negative::Vector{T}
    nonnegative::Vector{T}
end

ToeplitzOperator{T<:Number,Q<:Number}(V::Vector{T},W::Vector{Q}) =
    ToeplitzOperator{promote_type(T,Q)}(V,W)
ToeplitzOperator(V::AbstractVector,W::AbstractVector) =
    ToeplitzOperator(collect(V),collect(W))

function SymToeplitzOperator(V::Vector)
    W=V[2:end]
    V=copy(V)
    V[1]*=2
    ToeplitzOperator(W,V)
end

for OP in (:domainspace,:rangespace)
    @eval $OP(T::ToeplitzOperator) = ℓ⁰
end

getindex(T::ToeplitzOperator,k::Integer,j::Integer) =
    toeplitz_getindex(T.negative,T.nonnegative,k,j)

function toeplitz_getindex{T}(negative::AbstractVector{T},nonnegative::AbstractVector{T},k::Integer,j::Integer)
    if 0<k-j≤length(negative)
        negative[k-j]
    elseif 0≤j-k≤length(nonnegative)-1
        nonnegative[j-k+1]
    else
        zero(T)
    end
end

function toeplitz_getindex{T}(cfs::AbstractVector{T},k::Integer,j::Integer)
    if k==j
        2cfs[1]
    elseif 0<k-j≤length(cfs)-1
        cfs[k-j+1]
    elseif 0<j-k≤length(cfs)-1
        cfs[j-k+1]
    else
        zero(T)
    end
end

function Base.convert{T}(::Type{BandedMatrix},S::SubOperator{T,ToeplitzOperator{T},Tuple{UnitRange{Int},UnitRange{Int}}})
    ret=bzeros(S)

    kr,jr=parentindexes(S)

    neg=parent(S).negative
    pos=parent(S).nonnegative

    toeplitz_axpy!(1.0,neg,pos,kr,jr,ret)
end

bandinds(T::ToeplitzOperator)=(-length(T.negative),length(T.nonnegative)-1)

# slice of a ToeplitzOPerator is a ToeplitzOperator

function Base.getindex(T::ToeplitzOperator,kr::AbstractCount,jr::AbstractCount)
    sh=first(jr)-first(kr)
    st=step(jr)
    @assert st==step(kr)
    if sh ≥0
        ToeplitzOperator([reverse!(T.nonnegative[1:sh]);T.negative],T.nonnegative[sh+1:st:end])
    else
        ToeplitzOperator(T.negative[-sh+1:st:end],[reverse!(T.negative[1:-sh]);T.nonnegative])
    end
end

## Hankel Operator

type HankelOperator{T<:Number} <: Operator{T}
    coefficients::Vector{T}
end

for OP in (:domainspace,:rangespace)
    @eval $OP(T::HankelOperator) = ℓ⁰
end

HankelOperator(V::AbstractVector)=HankelOperator(collect(V))

HankelOperator(f::Fun)=HankelOperator(f.coefficients)

@eval Base.convert{TT}(::Type{Operator{TT}},T::HankelOperator)=HankelOperator(convert(Vector{TT},T.coefficients))

function hankel_getindex(v::AbstractVector,k::Integer,j::Integer)
   if k+j-1 ≤ length(v)
        v[k+j-1]
    else
        zero(eltype(v))
    end
end

getindex(T::HankelOperator,k::Integer,j::Integer) =
    hankel_getindex(T.coefficients,k,j)

function Base.convert{T}(::Type{BandedMatrix},S::SubOperator{T,HankelOperator{T},Tuple{UnitRange{Int},UnitRange{Int}}})
    ret=bzeros(S)

    kr,jr=parentindexes(S)
    cfs=parent(S).coefficients

    hankel_axpy!(1.0,cfs,kr,jr,ret)
end

bandinds(T::HankelOperator)=(1-length(T.coefficients),length(T.coefficients)-1)

## Laurent Operator

type LaurentOperator{T<:Number} <: Operator{T}
    negative::Vector{T}
    nonnegative::Vector{T}
end

for TYP in(:ToeplitzOperator,:LaurentOperator)
    @eval Base.convert{TT}(::Type{Operator{TT}},T::$TYP)=$TYP(convert(Vector{TT},T.negative),
                                                                            convert(Vector{TT},T.nonnegative))
end

for OP in (:domainspace,:rangespace)
    @eval $OP(T::LaurentOperator) = ℓ⁰
end

function laurent_getindex{T}(negative::AbstractVector{T},nonnegative::AbstractVector{T},k::Integer,j::Integer)
    # switch to double-infinite indices
    k=iseven(k)?-div(k,2):div(k-1,2)
    j=iseven(j)?-div(j,2):div(j-1,2)

    if 0<k-j≤length(negative)
        negative[k-j]
    elseif 0≤j-k≤length(nonnegative)-1
        nonnegative[j-k+1]
    else
        zero(T)
    end
end

getindex(T::LaurentOperator,k::Integer,j::Integer)=laurent_getindex(T.negative,T.nonnegative,k,j)

shiftbandinds(T::LaurentOperator)=-length(T.negative),length(T.nonnegative)-1
function bandinds(T::LaurentOperator)
    sbi=shiftbandinds(T)
    min(2sbi[1],-2sbi[end]),max(2sbi[end],-2sbi[1])
end

## algebra

function Base.maximum(T::ToeplitzOperator)
    if isempty(T.negative)
        maximum(T.nonnegative)
    elseif isempty(T.nonnegative)
        maximum(T.negative)
    else
        max(maximum(T.negative),maximum(T.nonnegative))
    end
end

for TYP in (:ToeplitzOperator,:LaurentOperator)
    @eval begin
        -(T::$TYP)=$TYP(-T.negative,-T.nonnegative)
        *(c::Number,T::$TYP)=$TYP(c*T.negative,c*T.nonnegative)
    end
end

-(H::HankelOperator)=HankelOperator(-H.coefficients)
*(c::Number,H::HankelOperator)=HankelOperator(c*H.coefficients)

## inv

function Base.inv(T::ToeplitzOperator)
    @assert length(T.nonnegative)==1
    ai=linsolve(T,[1.0];maxlength=100000)
    ToeplitzOperator(ai[2:end],ai[1:1])
end

function Fun(T::ToeplitzOperator)
   if length(T.nonnegative)==1
      Fun([T.nonnegative;T.negative],Taylor())
    elseif length(T.negative)==0
        Fun(T.nonnegative,Hardy{false}())
    else
        Fun(interlace(T.nonnegative,T.negative),Laurent(Circle()))
    end
end


text:
md"""
# Quick Tutorial

Check out the basic usage of [NaiveNASlib](https://github.com/DrChainsaw/NaiveNASlib.jl) for less verbose examples.

Here is a quick rundown of some common operations.

"""
@testset "Quick tutorial" begin #src
using NaiveNASflux, Flux, Test

# Create an input vertex which tells its output vertices that they can expect 2D convolutional input (i.e 4D arrays).
invertex = conv2dinputvertex("in", 3)

# Vertex type for Flux-layers is automatically inferred through [`fluxvertex`](@ref).
conv = fluxvertex(Conv((3,3), 3 => 5, pad=(1,1)), invertex)
batchnorm = fluxvertex(BatchNorm(nout(conv), relu), conv)

# Explore the graph.
@test inputs(conv) == [invertex]
@test outputs(conv) == [batchnorm]

@test nin(conv) == [3]
@test nout(conv) == 5

@test layer(conv) isa Flux.Conv
@test layer(batchnorm) isa Flux.BatchNorm

# Naming vertices is a good idea for debugging and logging purposes.
namedconv = fluxvertex("namedconv", Conv((5,5), 3=>7, pad=(2,2)), invertex)

@test name(namedconv) == "namedconv"

# Concatenate activations. Dimension is automatically inferred.
conc = concat("conc", namedconv, batchnorm)
@test nout(conc) == nout(namedconv) + nout(batchnorm)

# No problem to combine with convenience functions from NaiveNASlib.
residualconv = fluxvertex("residualconv", Conv((3,3), nout(conc) => nout(conc), pad=(1,1)), conc)
add = "add" >> conc + residualconv

@test name(add) == "add"
@test inputs(add) == [conc, residualconv]

# Computation graph for evaluation. It is basically a more general version of `Flux.Chain`.
graph = CompGraph(invertex, add)

# Access the vertices of the graph.
@test vertices(graph) == [invertex, namedconv, conv, batchnorm, conc, residualconv, add]

# `CompGraph`s can be evaluated just like any function.
x = ones(Float32, 7, 7, nout(invertex), 2)
@test size(graph(x)) == (7, 7, nout(add), 2) == (7 ,7, 12 ,2)

# Mutate number of neurons.
@test nout(add) == nout(residualconv) == nout(conv) + nout(namedconv) == 12
Δnout!(add => -3)
@test nout(add) == nout(residualconv) == nout(conv) + nout(namedconv) == 9

# Remove a layer.
@test nvertices(graph) == 7
remove!(batchnorm)
@test nvertices(graph) == 6

# Add a layer.
insert!(residualconv, v -> fluxvertex(BatchNorm(nout(v), relu), v))
@test nvertices(graph) == 7

# Change kernel size (and supply new padding).
namedconv |> KernelSizeAligned(-2,-2; pad=SamePad())

# Note: Parameters not changed yet...
@test size(NaiveNASflux.weights(layer(namedconv))) == (5, 5, 3, 7)

# ... because mutations are lazy by default so that no new parameters are created until the graph is evaluated.
@test size(graph(x)) == (7, 7, nout(add), 2) == (7, 7, 9, 2)
@test size(NaiveNASflux.weights(layer(namedconv))) == (3, 3, 3, 4)
end #src

ext: jl


text:
"""
NumDiffTools.derivative(f, x::Real, [abstol, reltol, h, t, maxiters])

Returns `df/dx` evaluted at `x`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function derivative(f, x::Real; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    fx = f(x);

    # Initialize Work Arrays
    A1 = Vector{typeof(fx)}(undef, maxiters)
    A2 = Vector{typeof(fx)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 1:maxiters

    extrapolate(h->(f(x + h) - fx)/h,             # Function to extrapolate
                A1, A2,                           # Work Arrays
                hh, k,                            # Sequence of step sizes (h)
                abstol, reltol                    # tolerances
                )
end

"""
NumDiffTools.derivative(f!, y::AbstractArray, x::Real)

Returns `df!/dx` evaluated at `x` assuming `f!` is called as `f!(y,x)`.  Result is also stored in `y`.

"""
function derivative(f!, y::AbstractArray, x::Real; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    # Initialize Work Arrays
    A1 = Vector{typeof(y)}(undef, maxiters)
    A2 = Vector{typeof(y)}(undef, maxiters)

    # Vector of step sizes
    hh = [h/t^(i - 1) for i = 1:maxiters]
    # Vector of orders
    k = 1:maxiters

    fx = similar(y)
    f!(fx,x)
    # Function to extrapolate
    function A(h, x)
        f!(y, x + h)
        return (y - fx)/h
    end
    y .= extrapolate(h->A(h,x),                        # Function to extrapolate
                     A1, A2,                           # Work Arrays
                     hh, k,                            # Sequence of step sizes (h)
                     abstol, reltol                    # tolerances
                     )
end

"""
NumDiffTools.derivative(result::AbstractArray, f, x::Real, [abstol, reltol, h, t, maxiters])

Stores `df/dx` evaluted at `x` in `result`, assuming `f` is called as `f(x)`.

Optional Arguments:\n
`abstol`: Target absolute error in derivative estimate (default `1e-12`)\n
`reltol`: Target relative error in derivative estimate (default `1e-12`)\n
`h`: Initial step size (default `0.1`)\n
`t`: Ratio of step sizes for each Richardson iteration (default `2.0001`.  Should not be an exact integer for periodic functions)\n
`maxiters`: Maximum number of Richardson iterations to perform (default `6`)\n
"""
function derivative!(result::AbstractArray, f, x::Real; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    result .= derivative(f, x; abstol = abstol, reltol = reltol, h = h, t= t, maxiters = maxiters)
    nothing
end

"""
NumDiffTools.derivative!(result::AbstractArray, f!, y::AbstractArray, x::Real)

Stores `df!/dx` evaluated at `x`  in `result` assuming `f!` is called as `f!(y,x)`.  Result is also stored in `y`.
"""
function derivative!(result::AbstractArray, f!, y::AbstractArray, x::Real; abstol::Real = 1e-12, reltol::Real = 1e-12, h::Real = 0.1, t::Real = 2.0001, maxiters::Integer = 6)
    result .= derivative(f!, y, x; abstol = abstol, reltol = reltol, h = h, t= t, maxiters = maxiters)
    nothing
end


text:
using Test
using StaticArrays, IntervalSets, LinearAlgebra

import ClimaCore: slab, Domains, Meshes, Topologies, Spaces, Fields

import ClimaCore.Geometry: Geometry
import ClimaCore.DataLayouts: IJFH

@testset "1d domain space" begin
    FT = Float64
    domain = Domains.IntervalDomain(
        Geometry.XPoint{FT}(-3) .. Geometry.XPoint{FT}(5),
        periodic = true,
    )
    mesh = Meshes.IntervalMesh(domain; nelems = 1)
    topology = Topologies.IntervalTopology(mesh)

    quad = Spaces.Quadratures.GLL{4}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace1D(topology, quad)

    coord_data = Spaces.coordinates_data(space)
    @test eltype(coord_data) == Geometry.XPoint{Float64}

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (4, 1, 1)
    coord_slab = slab(Spaces.coordinates_data(space), 1)
    @test coord_slab[1] == Geometry.XPoint{FT}(-3)
    @test coord_slab[4] == Geometry.XPoint{FT}(5)

    local_geometry_slab = slab(space.local_geometry, 1)
    dss_weights_slab = slab(space.dss_weights, 1)

    for i in 1:4
        @test Geometry.components(local_geometry_slab[i].∂x∂ξ) ≈
              @SMatrix [8 / 2]
        @test Geometry.components(local_geometry_slab[i].∂ξ∂x) ≈
              @SMatrix [2 / 8]
        @test local_geometry_slab[i].J ≈ (8 / 2)
        @test local_geometry_slab[i].WJ ≈ (8 / 2) * weights[i]
        if i in (1, 4)
            @test dss_weights_slab[i] ≈ 1 / 2
        else
            @test dss_weights_slab[i] ≈ 1
        end
    end

end

@testset "1×1 domain space" begin
    FT = Float32
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(-3) .. Geometry.XPoint{FT}(5),
        Geometry.YPoint{FT}(-2) .. Geometry.YPoint{FT}(8),
        x1periodic = true,
        x2periodic = false,
        x2boundary = (:south, :north),
    )
    mesh = Meshes.RectilinearMesh(domain, 1, 1)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{4}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (4, 4, 2, 1)
    coord_slab = slab(Spaces.coordinates_data(space), 1)
    @test coord_slab[1, 1] ≈ Geometry.XYPoint{FT}(-3.0, -2.0)
    @test coord_slab[4, 1] ≈ Geometry.XYPoint{FT}(5.0, -2.0)
    @test coord_slab[1, 4] ≈ Geometry.XYPoint{FT}(-3.0, 8.0)
    @test coord_slab[4, 4] ≈ Geometry.XYPoint{FT}(5.0, 8.0)

    local_geometry_slab = slab(space.local_geometry, 1)
    dss_weights_slab = slab(space.dss_weights, 1)

    for i in 1:4, j in 1:4
        @test Geometry.components(local_geometry_slab[i, j].∂x∂ξ) ≈
              @SMatrix [8/2 0; 0 10/2]
        @test Geometry.components(local_geometry_slab[i, j].∂ξ∂x) ≈
              @SMatrix [2/8 0; 0 2/10]
        @test local_geometry_slab[i, j].J ≈ (10 / 2) * (8 / 2)
        @test local_geometry_slab[i, j].WJ ≈
              (10 / 2) * (8 / 2) * weights[i] * weights[j]
        if i in (1, 4)
            @test dss_weights_slab[i, j] ≈ 1 / 2
        else
            @test dss_weights_slab[i, j] ≈ 1
        end
    end

    @test length(space.boundary_surface_geometries) == 2
    @test keys(space.boundary_surface_geometries) == (:south, :north)
    @test sum(parent(space.boundary_surface_geometries.north.sWJ)) ≈ 8
    @test parent(space.boundary_surface_geometries.north.normal)[1, :, 1] ≈
          [0.0, 1.0]
end

@testset "dss on 2×2 rectangular mesh (unstructured)" begin
    FT = Float64
    n1, n2 = 2, 2
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(0) .. Geometry.XPoint{FT}(4),
        Geometry.YPoint{FT}(0) .. Geometry.YPoint{FT}(4),
        x1periodic = false,
        x2periodic = false,
        x1boundary = (:west, :east),
        x2boundary = (:south, :north),
    )
    mesh = Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{4}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (4, 4, 2, 4)

    Nij = length(points)
    field = Fields.Field(IJFH{FT, Nij}(ones(Nij, Nij, 1, n1 * n2)), space)
    field_values = Fields.field_values(field)
    Spaces.horizontal_dss!(field)

    @testset "dss should not modify interior degrees of freedom of any element" begin
        result = true
        for el in 1:(n1 * n2)
            slb = slab(field_values, 1, el)
            for i in 2:(Nij - 1), j in 2:(Nij - 1)
                if slb[i, j] ≠ 1
                    result = false
                end
            end
        end
        @test result
    end
    s1 = slab(field_values, 1, 1)
    s2 = slab(field_values, 1, 2)
    s3 = slab(field_values, 1, 3)
    s4 = slab(field_values, 1, 4)

    @testset "vertex common to all (4) elements" begin
        @test (s1[Nij, Nij] == s2[1, Nij] == s3[Nij, 1] == s4[1, 1])
    end

    @testset "vertices common to (2) elements" begin
        @test s1[Nij, 1] == s2[1, 1]
        @test s1[1, Nij] == s3[1, 1]
        @test s2[Nij, Nij] == s4[Nij, 1]
        @test s3[Nij, Nij] == s4[1, Nij]
    end

    @testset "boundary faces" begin
        for fc in 2:(Nij - 1)
            @test s1[1, fc] == 1 # element 1 face 1
            @test s1[fc, 1] == 1 # element 1 face 3
            @test s2[Nij, fc] == 1 # element 2 face 2
            @test s2[fc, 1] == 1 # element 2 face 3
            @test s3[1, fc] == 1 # element 3 face 1
            @test s3[fc, Nij] == 1 # element 3 face 4
            @test s4[Nij, fc] == 1 # element 4 face 2
            @test s4[fc, Nij] == 1 # element 4 face 4
        end
    end

    @testset "interior faces" begin
        for fc in 2:(Nij - 1)
            @test (s1[Nij, fc] == s2[1, fc] == 2) # (e1, f2) == (e2, f1) == 2
            @test (s1[fc, Nij] == s3[fc, 1] == 2) # (e1, f4) == (e3, f3) == 2
            @test (s2[fc, Nij] == s4[fc, 1] == 2) # (e2, f4) == (e4, f3) == 2
            @test (s3[Nij, fc] == s4[1, fc] == 2) # (e3, f2) == (e4, f1) == 2
        end
    end
end

@testset "dss on 2×2 rectangular mesh" begin
    FT = Float64
    n1, n2 = 2, 2
    Nij = 4
    domain = Domains.RectangleDomain(
        Geometry.XPoint{FT}(0) .. Geometry.XPoint{FT}(4),
        Geometry.YPoint{FT}(0) .. Geometry.YPoint{FT}(4),
        x1periodic = false,
        x2periodic = false,
        x1boundary = (:west, :east),
        x2boundary = (:south, :north),
    )
    mesh = Meshes.RectilinearMesh(domain, n1, n2)
    grid_topology = Topologies.Topology2D(mesh)

    quad = Spaces.Quadratures.GLL{Nij}()
    points, weights = Spaces.Quadratures.quadrature_points(FT, quad)

    space = Spaces.SpectralElementSpace2D(grid_topology, quad)

    array = parent(Spaces.coordinates_data(space))
    @test size(array) == (Nij, Nij, 2, n1 * n2)

    data = zeros(Nij, Nij, 3, n1 * n2)
    data[:, :, 1, :] .= 1:Nij
    data[:, :, 2, :] .= (1:Nij)'
    data[:, :, 3, :] .= reshape(1:(n1 * n2), 1, 1, :)
    field = Fields.Field(IJFH{Tuple{FT, FT, FT}, Nij}(data), space)
    field_dss = Spaces.horizontal_dss!(copy(field))
    data_dss = parent(field_dss)

    @testset "slab 1" begin
        @test data_dss[1:(Nij - 1), 1:(Nij - 1), :, 1] ==
              data[1:(Nij - 1), 1:(Nij - 1), :, 1]
        @test data_dss[Nij, 1:(Nij - 1), :, 1] ==
              data[Nij, 1:(Nij - 1), :, 1] .+ data[1, 1:(Nij - 1), :, 2]
        @test data_dss[1:(Nij - 1), Nij, :, 1] ==
              data[1:(Nij - 1), Nij, :, 1] .+ data[1:(Nij - 1), 1, :, 3]
        @test data_dss[Nij, Nij, :, 1] ==
              data[Nij, Nij, :, 1] .+ data[1, Nij, :, 2] .+
              data[Nij, 1, :, 3] .+ data[1, 1, :, 4]
    end

    @testset "slab 2" begin
        @test data_dss[2:Nij, 1:(Nij - 1), :, 2] ==
              data[2:Nij, 1:(Nij - 1), :, 2]
        @test data_dss[1, 1:(Nij - 1), :, 2] ==
              data[Nij, 1:(Nij - 1), :, 1] .+ data[1, 1:(Nij - 1), :, 2]
        @test data_dss[2:Nij, Nij, :, 2] ==
              data[2:Nij, Nij, :, 2] .+ data[2:Nij, 1, :, 4]
        @test data_dss[1, Nij, :, 2] ==
              data[Nij, Nij, :, 1] .+ data[1, Nij, :, 2] .+
              data[Nij, 1, :, 3] .+ data[1, 1, :, 4]
    end

    @testset "slab 3" begin
        @test data_dss[1:(Nij - 1), 2:Nij, :, 3] ==
              data[1:(Nij - 1), 2:Nij, :, 3]
        @test data_dss[Nij, 2:Nij, :, 3] ==
              data[Nij, 2:Nij, :, 3] .+ data[1, 2:Nij, :, 4]
        @test data_dss[1:(Nij - 1), 1, :, 3] ==
              data[1:(Nij - 1), Nij, :, 1] .+ data[1:(Nij - 1), 1, :, 3]
        @test data_dss[Nij, 1, :, 3] ==
              data[Nij, Nij, :, 1] .+ data[1, Nij, :, 2] .+
              data[Nij, 1, :, 3] .+ data[1, 1, :, 4]
    end

    @testset "slab 3" begin
        @test data_dss[2:Nij, 2:Nij, :, 4] == data[2:Nij, 2:Nij, :, 4]
        @test data_dss[1, 2:Nij, :, 4] ==
              data[Nij, 2:Nij, :, 3] .+ data[1, 2:Nij, :, 4]
        @test data_dss[2:Nij, 1, :, 4] ==
              data[2:Nij, Nij, :, 2] .+ data[2:Nij, 1, :, 4]
        @test data_dss[1, 1, :, 4] ==
              data[Nij, Nij, :, 1] .+ data[1, Nij, :, 2] .+
              data[Nij, 1, :, 3] .+ data[1, 1, :, 4]
    end
end


text:
    """
    `RacahAlgebra.specialValue(w3j::RacahAlgebra.W3j)`
        ... attempts to find a special value for the Wigner 3j symbol w3j.
            A (istrue, rex)::Tuple{Bool, RacahExpression} is returned, and where istrue determined of whether a
            special value is returned in rex. For istrue = false, rex has no meaning.
    """
    function specialValue(w3j::RacahAlgebra.W3j)
        deltas = Kronecker[];    triangles = Triangle[];   integrals = Integral[]
        w3js = W3j[];   w6js = W6j[];   w9js = W9j[];   ylms = Ylm[];    djpqs = Djpq[]

        rexList = RacahAlgebra.symmetricForms(w3j)
        for  rex in rexList
            ww = rex.w3js[1]
            #
            #  Rule:        ( j1  j2  j3 )        J/2 [(J-2*j1)!(J-2*j2)!(J-2*j3)!] 1/2                   (J/2)!
            #  -----        (            )  = (-1)    [---------------------------]       -----------------------------------
            #               (  0   0   0 )            [          (J+1)!           ]       (J/2 - j1)! (J/2 - j2)! (J/2 - j3)!
            #
            #               with  J = j1+j2+j3  is  even; the 3j symbol is zero if J is odd.
            #
            specialW3j = W3j(ww.ja, ww.jb, ww.jc, 0, 0, 0);     J = ww.ja + ww.jb + ww.jc
            if  ww == specialW3j
                push!( deltas, Kronecker(J, Basic(:even)) )
                wa = RacahExpression( rex.summations, integrals, rex.phase + J/2,
                        rex.weight * sqrt( factorial(J - 2*ww.ja) * factorial(J - 2*ww.jb) * factorial(J - 2*ww.jb) / factorial(J+1) ) *
                                      factorial(J/2) / ( factorial(J/2-ww.ja) * factorial(J/2-ww.jb) * factorial(J/2-ww.jc) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j   j   0 )           j-m     -1/2
            #  -----        (           )   =   (-1)      [j]
            #               ( m  -m   0 )
            #
            specialW3j = W3j(ww.ja, ww.ja, 0, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.ja - ww.ma, rex.weight / sqrt(2*ww.ja+1),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j   j-1/2   1/2 )           j-m-1   (   j - m   ) 1/2
            #  -----        (                 )   =   (-1)        (-----------)
            #               ( m  -m-1/2   1/2 )                   ( 2j (2j+1) )
            #
            specialW3j = W3j(ww.ja, ww.ja-1//2, 1//2, ww.ma, -ww.ma-1//2, 1//2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.ja - ww.ma - 1,
                                      rex.weight * sqrt( (ww.ja - ww.ma)/(2*ww.ja * (2*ww.ja+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j+1    j    1 )           j-m-1  (    (j-m)(j-m+1)    ) 1/2
            #  -----        (               )   =   (-1)       (--------------------)
            #               (  m   -m-1   1 )                  ( (2j+3)(2j+2)(2j+1) )
            #
            specialW3j = W3j(ww.jb+1, ww.jb, 1, ww.ma, -ww.ma-1, 1)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma - 1,
                                      rex.weight * sqrt( (ww.jb - ww.ma)*(ww.jb - ww.ma + 1)/(2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j+1   j   1 )           j-m-1  (  2(j+m+1)(j-m+1)   ) 1/2
            #  -----        (             )   =   (-1)       (--------------------)
            #               (  m   -m   0 )                  ( (2j+3)(2j+2)(2j+1) )
            #
            specialW3j = W3j(ww.jb+1, ww.jb, 1, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma - 1,
                        rex.weight * sqrt( 2*(ww.jb + ww.ma + 1)*(ww.jb - ww.ma + 1)/(2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j    j     1 )           j-m  (  2*(j-m)(j+m+1)  ) 1/2
            #  -----        (              )   =   (-1)     (------------------)
            #               ( m   -m-1   1 )                ( (2j+2)(2j+1)(2j) )
            #
            specialW3j = W3j(ww.ja, ww.ja, 1, ww.ma, -ww.ma-1, 1)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.ja - ww.ma,
                        rex.weight * sqrt( 2*(ww.ja - ww.ma)*(ww.ja + ww.ma + 1)/(2*ww.ja+2) * (2*ww.ja+1) * (2*ww.ja) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:        ( j   j   1 )           j-m            2 m
            #  -----        (           )   =   (-1)      -----------------------
            #               ( m  -m   0 )                 ((2j+2)(2j+1)(2j))^(1/2)
            #
            specialW3j = W3j(ww.ja, ww.ja, 1, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.ja - ww.ma,
                                      rex.weight * 2*ww.ma / sqrt( (2*ww.ja+2) * (2*ww.ja+1) * (2*ww.ja) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 1      ( j+3/2    j    3/2 )           j-m+1/2   ( (j-m-1/2)(j-m+1/2)(j-m+3/2) ) 1/2
            #  -------      (                   )   =   (-1)          (-----------------------------)
            #               (   m   -m-3/2  3/2 )                     (   (2j+4)(2j+3)(2j+2)(2j+1)  )
            #
            specialW3j = W3j(ww.jb+3//2, ww.jb, 3//2, ww.ma, -ww.ma-3//2, 3//2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma + 1//2,
                                      rex.weight * sqrt( ( (ww.jb-ww.ma-1//2) * (ww.jb-ww.ma+1//2) * (ww.jb-ww.ma-3//2) )/
                                                         ( (2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 2      ( j+3/2    j    3/2 )           j-m+1/2 ( 3*(j-m+1/2)(j-m+3/2)(j+m+3/2) ) 1/2
            #  -------      (                   )   =   (-1)        (-------------------------------)
            #               (   m   -m-1/2  1/2 )                   (    (2j+4)(2j+3)(2j+2)(2j+1)   )
            #
            specialW3j = W3j(ww.jb+3//2, ww.jb, 3//2, ww.ma, -ww.ma-1//2, 1//2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma + 1//2,
                                      rex.weight * sqrt( ( 3*(ww.jb-ww.ma+1//2) * (ww.jb-ww.ma+3//2) * (ww.jb+ww.ma+3//2) )/
                                                         ( (2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 4      ( j+1/2    j    3/2 )            j-m-1/2           (        j-m+1/2         ) 1/2
            #  -------      (                   )    =   (-1)      (j+3*m+3/2) (------------------------)
            #               (   m   -m-1/2  1/2 )                              ( (2j+3)(2j+2)(2j+1)(2j) )
            #
            #
            specialW3j = W3j(ww.jb+1//2, ww.jb, 3//2, ww.ma, -ww.ma-1//2, 1//2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma - 1//2,
                                      rex.weight * (ww.jb + 3*ww.ma + 3//2) *
                                      sqrt( (ww.jb-ww.ma+1//2) / ((2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1) * 2*ww.jb) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 5      ( j+2    j   2 )           j-m  (   (j-m-1)(j-m)(j-m+1)(j-m+2)   ) 1/2
            #  -------      (              )   =   (-1)     (--------------------------------)
            #               ( m    -m-2  2 )                ( (2j+5)(2j+4)(2j+3)(2j+2)(2j+1) )
            #
            specialW3j = W3j(ww.jb+2, ww.jb, 2, ww.ma, -ww.ma-2, 2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma,
                                      rex.weight * sqrt( (ww.jb-ww.ma-1) * (ww.jb-ww.ma) * (ww.jb-ww.ma+1) *(ww.jb-ww.ma-2) /
                                                         ((2*ww.jb+5) * (2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2)  * (2*ww.jb+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 7      ( j+2   j   2 )          j-m  ( 6 (j+m+2)(j+m+1)(j-m+2)(j-m+1) ) 1/2
            #  -------      (             )  =   (-1)     (--------------------------------)
            #               ( m    -m   0 )               ( (2j+5)(2j+4)(2j+3)(2j+2)(2j+1) )
            #
            specialW3j = W3j(ww.jb+2, ww.jb, 2, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma,
                                      rex.weight * sqrt( 6* (ww.jb+ww.ma+2) * (ww.jb+ww.ma+1) * (ww.jb-ww.ma+2) * (ww.jb-ww.ma+1) /
                                                         ((2*ww.jb+5) * (2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2)  * (2*ww.jb+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 8      ( j+1    j    2 )         j-m+1    (  (j-m-1)(j-m)(j-m+1)(j+m+2)  ) 1/2
            #  -------      (               ) =   (-1)      2* (------------------------------)
            #               ( m    -m-2   2 )                  ( (2j+4)(2j+3)(2j+2)(2j+1)(2j) )
            #
            specialW3j = W3j(ww.jb+1, ww.jb, 2, ww.ma, -ww.ma-2, 2)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma + 1,
                                      rex.weight * 2* sqrt( (ww.jb-ww.ma-1) * (ww.jb-ww.ma) * (ww.jb-ww.ma+1) * (ww.jb-ww.ma+2) /
                                                            ((2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1)  * 2*ww.jb ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 9      ( j+1    j    2 )         j-m+1            (        (j-m+1)(j-m)          ) 1/2
            #  -------      (               ) =   (-1)     2*(j+2*m+2) (------------------------------)
            #               ( m    -m-1   1 )                          ( (2j+4)(2j+3)(2j+2)(2j+1)(2j) )
            #
            specialW3j = W3j(ww.jb+1, ww.jb, 2, ww.ma, -ww.ma-1, 1)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma + 1,
                                      rex.weight * 2* (ww.jb+2*ww.ma+2) * sqrt( (ww.jb-ww.ma-1) * (ww.jb-ww.ma) /
                                                           ((2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1)  * 2*ww.jb ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 10     ( j+1   j   2 )          j-m+1     (       6 (j+m+1)(j-m+1)       ) 1/2
            #  --------     (             )  =   (-1)      2*m (------------------------------)
            #               ( m    -m   0 )                    ( (2j+4)(2j+3)(2j+2)(2j+1)(2j) )
            #
            specialW3j = W3j(ww.jb+1, ww.jb, 2, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma + 1,
                                      rex.weight * 2*ww.ma * sqrt( 6* (ww.jb-ww.ma+1) * (ww.jb-ww.ma+1) /
                                                           ((2*ww.jb+4) * (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1)  * 2*ww.jb ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 13     ( j   j   2 )          j-m         2 (3*m*m - j(j+1))
            #  --------     (           )  =   (-1)     ------------------------------------
            #               ( m  -m   0 )               ((2j+3)(2j+2)(2j+1)(2j)(2j-1))^(1/2)
            #
            specialW3j = W3j(ww.jb, ww.jb, 2, ww.ma, -ww.ma, 0)
            if  ww == specialW3j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.jb - ww.ma,
                                      rex.weight * 2* (3*ww.ma*ww.ma - ww.jb*(ww.jb+1)) /
                                                   sqrt( (2*ww.jb+3) * (2*ww.jb+2) * (2*ww.jb+1)  * 2*ww.jb * (2*ww.jb-1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
        end

        return( (false, RacahExpression() ) )
    end

    """
    `RacahAlgebra.specialValue(w6j::RacahAlgebra.W6j)`
        ... attempts to find a special value for the Wigner 6-j symbol w6j.
            A (istrue, rex)::Tuple{Bool, RacahExpression} is returned, and where istrue determined of whether a
            special value is returned in rex. For istrue = false, rex has no meaning.
    """
    function specialValue(w6j::RacahAlgebra.W6j)
        deltas = Kronecker[];    triangles = Triangle[];   integrals = Integral[]
        w3js = W3j[];   w6js = W6j[];   w9js = W9j[];   ylms = Ylm[];    djpqs = Djpq[]

        rexList = RacahAlgebra.symmetricForms(w6j)
        for  rex in rexList
            ww = rex.w6js[1]
            #
            #                                            j1+j2+j3
            #  Rule:         ( j1   j2   j3 )        (-1)
            #  -----        {(              )}   =  --------------  t(j1,j2,j3)  d(j1,l2) d(j2,l1)
            #                ( l1   l2   0  )                   1/2
            #                                         [ j1, j2 ]
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, ww.d, ww.e, 0)
            if  ww == specialW6j
                push!( deltas, Kronecker(ww.a, ww.e) )
                push!( deltas, Kronecker(ww.b, ww.d) )
                push!( triangles, Triangle(ww.a, ww.b, ww.c) )
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.a + ww.b + ww.c,
                                      rex.weight / sqrt( (2*ww.a+1)*(2*ww.b+1)), deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1     j2      j3   )            j1+j2+j3    (     (j1+j3-j2)(j1+j2-j3+1)     ) 1/2
            #  -----        {(                     )}   =   (-1)            (--------------------------------)
            #                ( 1/2  j3-1/2  j2+1/2 )                        ( (2*j2+1)(2*j2+2)(2*j3)(2*j3+1) )
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1//2, ww.c -1//2,  ww.b +1//2)
            if  ww == specialW6j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.a + ww.b + ww.c,
                                      rex.weight * sqrt( (ww.a+ww.c-ww.b)*(ww.a+ww.b-ww.c+1) /
                                                         (2*ww.b+1)*(2*ww.b+2)*2*ww.c*(2*ww.c+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1     j2      j3   )            j1+j2+j3+1   (     (j1+j2+j3+2)(j2+j3-j1+1)     ) 1/2
            #  -----        {(                     )}   =   (-1)             (----------------------------------)
            #                ( 1/2  j3+1/2  j2+1/2 )                         ( (2*j2+1)(2*j2+2)(2*j3+1)(2*j3+2) )
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1//2, ww.c +1//2,  ww.b +1//2)
            if  ww == specialW6j
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.a + ww.b + ww.c + 1,
                                      rex.weight * sqrt( (ww.a+ww.b+ww.c+2)*(ww.b+ww.c-ww.a+1) /
                                                         (2*ww.b+1)*(2*ww.b+2)*(2*ww.c+1)*(2*ww.c+2) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1   j2    j3  )            S       (           S(S+1)(S-2*j1)(S-2*j1-1)           ) 1/2
            #  -----        {(                )}   =   (-1)        (----------------------------------------------)
            #                (  1  j3-1  j2-1 )                    ( (2*j2-1)(2*j2)(2*j2+1)(2*j3-1)(2*j3)(2*j3+1) )
            #
            #                with  S = j1+j2+j3
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1, ww.c -1,  ww.b -1)
            if  ww == specialW6j
                S = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + S,
                                      rex.weight * sqrt( S*(S+1)*(S-2*ww.a)*(S-2*ww.a -1) /
                                                         (2*ww.b-1)* 2*ww.b * (2*ww.b+1)* (2*ww.c-1)* 2*ww.c * (2*ww.c+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1   j2   j3  )            S   (      2(S+1)(S-2*j1)(S-2*j2)(S-2*j3+1)        ) 1/2
            #  -----        {(               )}   =   (-1)    (----------------------------------------------)
            #                (  1  j3-1  j2  )                ( (2*j2)(2*j2+1)(2*j2+2)(2*j3-1)(2*j3)(2*j3+1) )
            #
            #                with  S = j1+j2+j3
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1, ww.c -1,  ww.b)
            if  ww == specialW6j
                S = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + S,
                                      rex.weight * sqrt( 2*(S+1)*(S-2*ww.a)*(S-2*ww.b)*(S-2*ww.c +1) /
                                                         2*ww.b * (2*ww.b+1) * (2*ww.b+2) * (2*ww.c-1) * 2*ww.c * (2*ww.c+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1   j2   j3    )            S   (     (S-2*j2-1)(S-2*j2)(S-2*j3+1)(S-2*j3+2)     ) 1/2
            #  -----        {(                 )}   =   (-1)    (------------------------------------------------)
            #                (  1  j3-1  j2+1  )                ( (2*j2+1)(2*j2+2)(2*j2+3)(2*j3-1)(2*j3)(2*j3+1) )
            #
            #                with  S = j1+j2+j3
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1, ww.c -1,  ww.b +1)
            if  ww == specialW6j
                S = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + S,
                                      rex.weight * sqrt( (S-2*ww.b-1)*(S-2*ww.b)*(S-2*ww.c +1)*(S-2*ww.c +2) /
                                                         (2*ww.b+1) * (2*ww.b+2) * (2*ww.b+3) * (2*ww.c-1) * 2*ww.c * (2*ww.c+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( j1  j2  j3  )            S+1         2 [ j2(j2+1) + j3(j3+1) -  j1(j1+1) ]
            #  -----        {(             )}   =   (-1)       ----------------------------------------------  1/2
            #                (  1  j3  j2  )                   ( (2*j2)(2*j2+1)(2*j2+2)(2*j3)(2*j3+1)(2*j3+2) )
            #
            #                with  S = j1+j2+j3
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 1, ww.c,  ww.b)
            if  ww == specialW6j
                S = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + S + 1,
                                      rex.weight * 2 * ( ww.b*(ww.b+1) + ww.c*(ww.c+1) - ww.a*(ww.a+1) ) /
                                                       sqrt( 2*ww.b * (2*ww.b+1) * (2*ww.b+2) * 2*ww.c * (2*ww.c+1) * (2*ww.c+2) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 1       (  a     b      c   )            s   (         (s-1)s(s+1)(s-2*a-2)(s-2*a-1)(s-2*a)         ) 1/2
            #  -----        {(                   )}   =   (-1)    (------------------------------------------------------)
            #                ( 3/2  c-3/2  b-3/2 )                ( (2*b-2)(2*b-1)(2*b)(2*b+1)(2*c-2)(2*c-1)(2*c)(2*c+1) )
            #
            #                with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-3//2,  ww.b-3//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( (s-1)* s * (s+1) * (s-2*ww.a-2) * (s-2*ww.a-1) * (s-2*ww.a)/
                                                         (2*ww.b-2)*(2*ww.b-1)*2*ww.b*(2*ww.b+1) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 2       (  a     b      c   )            s   (       3*s(s+1)(s-2*a-1)(s-2*a)(s-2*b)(s-2*c+1)       ) 1/2
            #  -----        {(                   )}   =   (-1)    (------------------------------------------------------)
            #                ( 3/2  c-3/2  b-1/2 )                ( (2*b-1)(2*b)(2*b+1)(2*b+2)(2*c-2)(2*c-1)(2*c)(2*c+1) )
            #
            #                with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-3//2,  ww.b-1//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( 3*s * (s+1) * (s-2*ww.a-1) * (s-2*ww.a) * (s-2*ww.b) * (s-2*ww.c+1)/
                                                         ((2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 3       (  a     b      c   )            s   (   3*(s+1)(s-2*a)(s-2*b-1)(s-2*b)(s-2*c+1)(s-2*c+2)   ) 1/2
            #  -----        {(                   )}   =   (-1)    (------------------------------------------------------)
            #                ( 3/2  c-3/2  b+1/2 )                ( (2*b)(2*b+1)(2*b+2)(2*b+3)(2*c-2)(2*c-1)(2*c)(2*c+1) )
            #
            #                with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-3//2,  ww.b+1//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( 3* (s+1) * (s-2*ww.a) * (s-2*ww.b-1) * (s-2*ww.b) * (s-2*ww.c+1) * (s-2*ww.c+2)/
                                                         (2*ww.b *(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 4       (  a     b      c   )            s    (  (s-2*b-2)(s-2*b-1)(s-2*b)(s-2*c+1)(s-2*c+2)(s-2*c+3)  ) 1/2
            #  -----        {(                   )}   =   (-1)     (--------------------------------------------------------)
            #                ( 3/2  c-3/2  b+3/2 )                 ( (2*b+1)(2*b+2)(2*b+3)(2*b+4)(2*c-2)(2*c-1)(2*c)(2*c+1) )
            #
            #                with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-3//2,  ww.b+3//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( (s-2*ww.b-2) * (s-2*ww.b-1) * (s-2*ww.b) * (s-2*ww.c+1) * (s-2*ww.c+2) * (s-2*ww.c+3)/
                                                         ((2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3)*(2*ww.b+4) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule: 5       (  a     b      c   )            s                                      (                     (s+1)(s-2*a)                     ) 1/2
            #  -----        {(                   )}   =  (-1)   [2*(s-2*b)(s-2*c) - (s+2)(s-2*a-1)]  (------------------------------------------------------)
            #                ( 3/2  c-1/2  b-1/2 )                                                   ( (2*b-1)(2*b)(2*b+1)(2*b+2)(2*c-1)(2*c)(2*c+1)(2*c+2) )
            #
            #      with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-1//2,  ww.b-1//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * (2*(s-2*ww.b)*(s-2*ww.c) - (s+2)*(s-2*ww.a-1)) *   sqrt( (s+1)*(s-2*ww.a) /
                                                         ((2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2) * (2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 6       (  a     b      c   )            s                                     (                   (s-2*b)(s-2*c+1)                   ) 1/2
            #   -----        {(                   )}  =  (-1)  [ (s-2*b-1)(s-2*c) - 2*(s+2)(s-2*a)]  (------------------------------------------------------)
            #                 ( 3/2  c-1/2  b+1/2 )                                                  ( (2*b)(2*b+1)(2*b+2)(2*b+3)(2*c-1)(2*c)(2*c+1)(2*c+2) )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 3//2, ww.c-1//2,  ww.b+1//2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * ((s-2*ww.b-1)*(s-2*ww.c) - 2*(s+2)*(s-2*ww.a)) *   sqrt( (s-2*ww.b) * (s-2*ww.c + 1) /
                                                    (2*ww.b * (2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3) * (2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+2)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 7       ( a   b    c  )            s   (                  (s-2)(s-1)s(s+1)(s-2*a-3)(s-2*a-2)(s-2*a-1)(s-2*a) ) 1/2
            #   -----        {(             )}   =   (-1)    (---------------------------------------------------------------------)
            #                 ( 2  c-2  b-2 )                ( (2*b-3)(2*b-2)(2*b-1)(2*b)(2*b+1)(2*c-3)(2*c-2)(2*c-1)(2*c)(2*c+1)  )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-2,  ww.b-2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( (s-2) * (s-1) *s* (s-2*ww.a-3)* (s-2*ww.a-2)* (s-2*ww.a-1)* (s-2*ww.a) /
                            ((2*ww.b-3)*(2*ww.b-2)*(2*ww.b-1)*2*ww.b*(2*ww.b+1) * (2*ww.c-3)*(2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 8       ( a   b    c  )            s     (     (s-1)s(s+1)(s-2*a-2)(s-2*a-1)(s-2*a)(s-2*b)(s-2*c+1)            ) 1/2
            #   -----        {(             )}   =   (-1)    2 (---------------------------------------------------------------------)
            #                 ( 2  c-2  b-1 )                  ( (2*b-2)(2*b-1)(2*b)(2*b+1)(2*b+2)(2*c-3)(2*c-2)(2*c-1)(2*c)(2*c+1)  )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-2,  ww.b-1)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 2 + sqrt( (s-1)*s*(s+1) * (s-2*ww.a-2) * (s-2*ww.a-1) * (s-2*ww.a) * (s-2*ww.b) * (s-2*ww.c+1)/
                            ((2*ww.b-2)*(2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2) * (2*ww.c-3)*(2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 9       ( a   b    c  )            s  (       6*s*(s+1)(s-2*a-1)(s-2*b-1)(s-2*a)(s-2*b)(s-2*c+1)(s-2*c+2)  ) 1/2
            #   -----        {(             )}   =   (-1)   (--------------------------------------------------------------------)
            #                 ( 2  c-2   b  )               ( (2*b-1)(2*b)(2*b+1)(2*b+2)(2*b+3)(2*c-3)(2*c-2)(2*c-1)(2*c)(2*c+1) )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-2,  ww.b)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( 6*s * (s+1)*(s-2*ww.a-1)*(s-2*ww.b-1)*(s-2*ww.a)*(s-2*ww.b)*(s-2*ww.c+1)*(s-2*ww.c+2)/
                            ((2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3) * (2*ww.c-3)*(2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 10      ( a   b    c  )            s    (    (s+1)(s-2*a)(s-2*b-2)(s-2*b-1)(s-2*b)(s-2*c+1)(s-2*c+2)(s-2*c+3) ) 1/2
            #   -----        {(             )}   =   (-1)   2 (---------------------------------------------------------------------)
            #                 ( 2  c-2  b+1 )                 ( (2*b)(2*b+1)(2*b+2)(2*b+3)(2*b+4)(2*c-3)(2*c-2)(2*c-1)(2*c)(2*c+1)  )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-2,  ww.b+1)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 2 * sqrt( (s+1)*(s-2*ww.a)*(s-2*ww.b-2)*(s-2*ww.b-1)*(s-2*ww.b)*(s-2*ww.c+1)*(s-2*ww.c+2)*(s-2*ww.c+3) /
                                      ( 2*ww.b*(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3)*(2*ww.b+4) * (2*ww.c-3)*(2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 11      ( a   b    c  )            s   (  (s-2*b-3)(s-2*b-2)(s-2*b-1)(s-2*b)(s-2*c+1)(s-2*c+2)(s-2*c+3)(s-2*c+4) ) 1/2
            #   -----        {(             )}   =   (-1)    (-------------------------------------------------------------------------)
            #                 ( 2  c-2  b+2 )                ( (2*b+1)(2*b+2)(2*b+3)(2*b+4)(2*b+5)(2*c-3)(2*c-2)(2*c-1)(2*c)(2*c+1)    )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-2,  ww.b+2)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * sqrt( (s-2*ww.b-3)*(s-2*ww.b-2)*(s-2*ww.b-1)*(s-2*ww.b) * (s-2*ww.c+1)*(s-2*ww.c+2)*(s-2*ww.c+3)*(s-2*ww.c+4) /
                                      ( (2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3)*(2*ww.b+4)*(2*ww.b+5) * (2*ww.c-3)*(2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 12      ( a   b    c  )        s                                    (                           s*(s+1)(s-2*a-1)(s-2*a)                  ) 1/2
            #   -----        {(             )} = (-1)   4 [ (a+b)(a-b+1) - (c-1)(c-b+1)]  (--------------------------------------------------------------------)
            #                 ( 2  c-1  b-1 )                                             ( (2*b-2)(2*b-1)(2*b)(2*b+1)(2*b+2)(2*c-2)(2*c-1)(2*c)(2*c+1)(2*c+2) )
            #
            #      with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-1,  ww.b-1)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 4 * ( (ww.a+ww.b)*(ww.a-ww.b+1) - (ww.c-1)*(ww.c-ww.b+1) ) *
                                      sqrt( s*(s+1) * (s-2*ww.a-1)*(s-2*ww.a) /
                                          ( (2*ww.b-2)*(2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+2) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 13      ( a   b    c  )        s                                (                    6*(s+1)(s-2*a)(s-2*b)(s-2*c+1)                  ) 1/2
            #   -----        {(             )} = (-1)   2 [ (a+b+1)(a-b) - c*c  + 1 ] (--------------------------------------------------------------------)
            #                 ( 2  c-1   b  )                                         ( (2*b-1)(2*b)(2*b+1)(2*b+2)(2*b+3)(2*c-2)(2*c-1)(2*c)(2*c+1)(2*c+2) )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-1, ww.b)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 2 * ( (ww.a+ww.b+1)*(ww.a-ww.b) - ww.c*ww.c + 1 ) *
                                      sqrt(  6*(s+1) * (s-2*ww.a) * (s-2*ww.b) * (s-2*ww.c+1) /
                                           ( (2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+2) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 14      ( a   b    c  )        s                                     (                  (s-2*b-1)(s-2*b)(s-2*c+1)(s-2*c+2)                ) 1/2
            #   -----        {(             )} = (-1)  4 [ (a+b+2)(a-b-1) - (c-1)(b+c+2) ] (--------------------------------------------------------------------)
            #                 ( 2  c-1  b+1 )                                              ( (2*b)(2*b+1)(2*b+2)(2*b+3)(2*b+4)(2*c-2)(2*c-1)(2*c)(2*c+1)(2*c+2) )
            #
            #                 with  s = a + b + c
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c-1,  ww.b+1)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 4 * ( (ww.a+ww.b+2)*(ww.a-ww.b-1) - (ww.c-1)*(ww.b+ww.c+2) ) *
                                      sqrt( (s-2*ww.b-1)*(s-2*ww.b) * (s-2*ww.c+1)*(s-2*ww.c+2) /
                                          ( 2*ww.b*(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3)*(2*ww.b+4) * (2*ww.c-2)*(2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+2) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #   Rule: 15      ( a   b   c )        s                                     (                                                                    ) 1/2
            #   -----        {(           )} = (-1)  2 [ 3*X(X-1) - 4*b*(b+1)*c*(c+1) ]  (--------------------------------------------------------------------)
            #                 ( 2   c   b )                                              ( (2*b-1)(2*b)(2*b+1)(2*b+2)(2*b+3)(2*c-1)(2*c)(2*c+1)(2*c+2)(2*c+3) )
            #
            #                 with  s = a + b + c   and   X = b*(b+1) + c*(c+1) - a*(a+1)
            #
            specialW6j = W6j( ww.a, ww.b, ww.c, 2, ww.c,  ww.b)
            if  ww == specialW6j
                s = ww.a + ww.b + ww.c;     X = ww.b*(ww.b+1) + ww.c*(ww.c+1) - ww.a*(ww.a+1)
                wa = RacahExpression( rex.summations, integrals, rex.phase + s,
                                      rex.weight * 2 * ( 3*X*(X-1) - 4*ww.b*(ww.b+1)*ww.c*(ww.c+1) ) *
                                      sqrt( 1 / ( (2*ww.b-1)*2*ww.b*(2*ww.b+1)*(2*ww.b+2)*(2*ww.b+3) * (2*ww.c-1)*2*ww.c*(2*ww.c+1)*(2*ww.c+2)*(2*ww.c+3) ) ),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
        end

        return( (false, RacahExpression() ) )
    end

    """
    `RacahAlgebra.specialValue(w9j::RacahAlgebra.W9j)`
        ... attempts to find a special value for the Wigner 9-j symbol w9j.
            A (istrue, rex)::Tuple{Bool, RacahExpression} is returned, and where istrue determined of whether a
            special value is returned in rex. For istrue = false, rex has no meaning.
    """
    function specialValue(w9j::RacahAlgebra.W9j)
        deltas = Kronecker[];    triangles = Triangle[];   integrals = Integral[]
        w3js= W3j[];   w6js = W6j[];   w9js = W9j[];   ylms = Ylm[];    djpqs = Djpq[]

        rexList = RacahAlgebra.symmetricForms(w9j)
        for  rex in rexList
            ww = rex.w9js[1]
            #
            #  Rule:         ( ja   jb   0 )
            #  -----         (             )         D(ja,jc,je)
            #               {( jc   jd   0 )}   =  -------------- 1/2  delta(a,b)  delta(c,d)  delta(e,f)
            #                (             )       [ ja, jc, je ]
            #                ( je   jf   0 )
            #
            specialW9j = W9j( ww.a, ww.b, 0, ww.c, ww.d, 0, ww.e, ww.f, 0)
            if  ww == specialW9j
                push!( deltas, Kronecker(ww.a, ww.b) )
                push!( deltas, Kronecker(ww.c, ww.d) )
                push!( deltas, Kronecker(ww.e, ww.f) )
                push!( triangles, Kronecker(ww.a, ww.c, ww.e) )
                wa = RacahExpression( rex.summations, integrals, rex.phase, rex.weight / sqrt( (2*ww.a+1)*(2*ww.c+1)*(2*ww.c+1)),
                                      deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
            #
            #  Rule:         ( ja   jb   jc )           jb+jc+je+jf
            #  -----         (              )       (-1)                   ( ja   jb   jc )
            #               {( jd   je   jf )}   =  --------------- 1/2   {(              )}
            #                (              )          [ jc, jg ]          ( je   jd   jg )
            #                ( jg   jh   0  )
            #
            specialW9j = W9j( ww.a, ww.b, ww.c, ww.d, ww.e, ww.f, ww.g, ww.h, 0)
            if  ww == specialW9j
                push!( w6js, W6j(ww.a, ww.b, ww.c, ww.e, ww.d, ww.g) )
                wa = RacahExpression( rex.summations, integrals, rex.phase + ww.b + ww.c + ww.e + ww.f,
                                      rex.weight / sqrt( (2*ww.c+1)*(2*ww.g+1)), deltas, triangles, w3js, w6js, w9js, ylms, djpqs )
                return( (true, wa) )
            end
        end

        return( (false, RacahExpression() ) )
    end


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

"""
    update_auxvars_ss(point, model)

This function is for internal use. Do not call directly.
"""
function update_auxvars_ss(point::Vector{Float64}, model::Model)
    nvars = length(model.variables)
    nshks = length(model.shocks)
    nauxs = length(model.auxvars)

    nall = nvars + nshks + nauxs

    if nauxs == 0
        return point
    end

    trange = -model.maxlag:model.maxlead
    ntimes = length(trange)
    time0 = 1 + model.maxlag

    pt0 = [ones(ntimes) 0.0 .+ trange] * reshape(point, 2, :)
    # pt0 = [pt0[:,1:nvars] zeros(ntimes, nshks)]
    pt0 = ModelBaseEcon.update_auxvars(pt0, model)

    result = Vector{Float64}(undef, 2 * nall)
    result[1:length(point)] = point
    for i = 1:nauxs
        result[2 * (nvars + nshks + i) - 1] = pt0[time0, nvars + nshks + i]
    end

    if ! model.flags.ssZeroSlope
        shift = model.options.shift
        pt1 = [ones(ntimes) shift .+ trange] * reshape(point, 2, :)
        # pt1 = [pt1[:,1:nvars] zeros(ntimes, nshks)]
        pt1 = ModelBaseEcon.update_auxvars(pt1, model)
        for i = 1:nauxs
            result[2 * (nvars + nshks + i)] = (pt1[time0, nvars + nshks + i] - pt0[time0, nvars + nshks + i]) / shift
        end
    end

    return result
end
@assert precompile(update_auxvars_ss, (Vector{Float64}, Model))

"""
    _do_update_auxvars_presolve!(model)

Call `update_auxvars_ss`, then call `presolve_sstate!`.

This function is for internal use. Do not call directly.
"""
function _do_update_auxvars_presolve!(model::Model; verbose::Bool)
    ss = model.sstate
    # set shocks (level and slope) and steady slopes to 0.0 presolved
    for (i, v) in enumerate(model.allvars)
        if isshock(v)
            ss.values[2i .+ (-1:0)] .= 0.0
            ss.mask[2i .+ (-1:0)] .= true
        elseif issteady(v)
            ss.values[2i] = 0.0
            ss.mask[2i] = true
        elseif isexog(v)
            ss.mask[2i .+ (-1:0)] .= true
        end
    end
    # if all slopes are zero, make it so
    if model.flags.ssZeroSlope
        ss.values[2:2:end] .= 0.0
        ss.mask[2:2:end] .= true
    end
    # sometimes update_auxvars_ss might change the behaviour of presolve_sstate!
    # because it might set the values of aux variable differently and so
    # the presolve would be done on a different line parallel to the presolve-variable coordinate.
    aux_vals = update_auxvars_ss(ss.values, model)
    ss.values[ .! ss.mask ] = aux_vals[ .! ss.mask ]
    old_solved = -1
    solved = sum(ss.mask)
    while old_solved < solved
        presolve_sstate!(model; verbose=verbose)
        aux_vals .= update_auxvars_ss(ss.values, model)
        ss.values[ .! ss.mask ] = aux_vals[ .! ss.mask ]
        old_solved = solved
        solved = sum(ss.mask)
    end
    return nothing
end
@assert precompile(_do_update_auxvars_presolve!, (Model,))

"""
    clear_sstate!(model; lvl=0.1, slp=0.0, <options>)

Set the steady state values to the provided defaults and presolve.

### Arguments
  * `model` - the model instance
  * `lvl`, `slp` - the initial guess for the level and the slope. Each could be
    a number or a vector of length equal to the number of variable in the mode.

### Options
Standard options (default values are taken from `model.options`)
  * `verbose`
"""
function clear_sstate!(model::Model; lvl=0.1, slp=0.0, verbose=model.options.verbose)
    ss = model.sstate
    nvars = length(model.variables)
    nshks = length(model.shocks)
    fill!(ss.values, 0.0)
    ss.values[1:2:2nvars] .= lvl  # default initial guess for level
    ss.values[2:2:2nvars] .= slp  # default initial guess for slope
    ss.mask[:] .= false
    return _do_update_auxvars_presolve!(model; verbose=verbose)
end
export clear_sstate!
@assert precompile(clear_sstate!, (Model,))

"""
    initial_sstate!(model, init; <options>)

Set the steady state values from the given vector and presolve.

Call this function to specify initial guesses for the iterative steady state
solver. If the value of a steady state variable is known, it is better to use
[`@steadystate`](@ref ModelBaseEcon.@steadystate) to add that as a steady state
constraint.

### Arguments
  * `model` - the model.
  * `init` - a vector of length equal to twice the number of variables in the
    model. The level and slope values are staggered, i.e., the level and slope
    of variable j are in init[2j-1] and init[2j].

### Options
Standard options (default values are taken from `model.options`)
  * `verbose`

"""
function initial_sstate!(model::Model, init::AbstractVector{Float64}; verbose=model.options.verbose)
    ss = model.sstate
    nvars = length(model.variables)
    nshks = length(model.shocks)
    nauxs = length(model.auxvars)
    ninit = length(init)
    if ninit ∉ (2nvars, 2nvars + 2nauxs, 2nvars + 2nshks + 2nauxs)
        error("Incorrect dimension if initial guess: $(ninit). Expected $(2nvars) or $(2nvars + 2nauxs) or $(2nvars + 2nshks + 2nauxs)")
    end
    ss.values[1:ninit] = init
    ss.values[ninit + 1:end] .= 0.0
    ss.mask[:] .= false
    return _do_update_auxvars_presolve!(model; verbose=verbose)
end
export initial_sstate!
@assert precompile(initial_sstate!, (Model, Vector{Float64}))


text:
# the core loop for several CRC calculations.  it is used by CRC32,
# whose zlib implementation in C runs ~2.5 times faster than this
# code.

function generic_loop(::Type{D}, remainder::A, data,
        table::Vector{A}, word_size) where {D<:Unsigned, A<:Unsigned}
    for word::D in data
        word ⊻= remainder%D
        remainder = (remainder >>> word_size) ⊻ table[1 + word]
    end
    remainder
end

# try accessing the lowest byte directly to avoid the shift (slow)
struct Bytes
    a::UInt8
    b::UInt8
    c::UInt8
    d::UInt8
end

function struct_loop(::Type{UInt8}, remainder::UInt64, data::Vector{UInt8}, table::Vector{UInt64}, ignored)
    r::UInt64 = remainder
    b = reinterpret(Bytes, r)
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        word ⊻= b.a
        r = (r >>> 8) ⊻ table[1 + word]
        b = reinterpret(Bytes, r)
    end
    return r
end

# try accessing via aliased arrays

function array_loop(::Type{UInt8}, remainder::UInt64, data::Vector{UInt8}, table::Vector{UInt64}, ignored)
    r = Vector{UInt64}(undef, 1)
    @inbounds r[1] = remainder
    b = reinterpret(UInt8, r)
    for i in 1:length(data)
        @inbounds word::UInt8 = data[i]
        @inbounds word ⊻= b[1]
        @inbounds r[1] = (r[1] >>> 8) ⊻ table[1 + word]  # @inbounds no help here
    end
    r[1]
end

# according to my notes this was slightly (~15%) faster than
# generic_loop.  repeated here, however, it seems to be identical to
# the generic case.

function byte_loop(::Type{UInt8}, remainder::A, data::Vector{UInt8}, table::Vector{A}, ignored) where {A<:Unsigned}
    for i in 1:length(data)
        word = data[i]
        word ⊻= remainder%UInt8
        remainder = (remainder >>> 8) ⊻ table[1 + word]
    end
    remainder
end

# timing with a standard data set (the table is "meaningless")

Random.seed!(0)  # repeatable results
TABLE = UInt64[rand(UInt64) for _ in 1:256]
DATA = UInt8[rand(UInt8) for _ in 1:10_000_000]
#KNOWN = 17880998858342525629  # 1 million
KNOWN = 0x65280cfa186b00d5  # 10 million
# KNOWN = 10790464939549839568  # 100 million
using Printf
function run_std(loop)
    loop(UInt8, zero(UInt64), DATA[1:10], TABLE, 8)  # warmup
    @printf("%15s: ", loop)
    @time (result = loop(UInt8, zero(UInt64), DATA, TABLE, 8))
    @assert result == KNOWN
end

for loop in (generic_loop, byte_loop, array_loop)
    run_std(loop)
end
run_std(array_loop)

#  generic_loop: elapsed time: 0.023491408 seconds (16 bytes allocated)
#   struct_loop: elapsed time: 2.027252231 seconds (1680000168 bytes allocated)
#     byte_loop: elapsed time: 0.02304929 seconds (16 bytes allocated)
#    array_loop: elapsed time: 0.045178428 seconds (168 bytes allocated)

# benchmarks against Libz

SIZE = 300_000_000
# SIZE = 300_000

function time_libz()
    println("libz")
    data = rand(UInt8, SIZE)
    check = crc32(data)
    @time crc32(data)
    for tables in (Single, Multiple)
        ours = crc(CRC_32, tables=tables)
        @assert ours(data) == check
        println(tables)
        @time ours(data)
    end
end

function time_no_tables()
    println("no_tables")
    ours = crc(CRC_15, tables=NoTables)
    data = rand(UInt8, round(Int, SIZE//10))
    @assert ours(CHECK) == CRC_15.check
    @time ours(data)
end

function time_padded()
    println("padded")
    ours = crc(CRC_64)
    data = rand(UInt8, SIZE)
    @assert ours(CHECK) == CRC_64.check
    @time ours(data)
end

Random.seed!(0)  # repeatable results

time_libz()
time_no_tables()
time_padded()


text:
#========================================================================================#
"""
	NBodies

A system of N interacting bodies moving in dim dimensions over time T in res timesteps.

Author: Niall Palfreyman, 24/3/2022.
"""
module NBodies

# Externally callable methods of NBodies
export NBody

#-----------------------------------------------------------------------------------------
# Module types:

"""
	NBody

An NBody system capable of containing multiple (N) bodies that gravitationally interact
with each other.
"""
struct NBody
	N							# Number of bodies
	T							# Duration of simulation
	res							# Timestep resolution
	xi							# Initial positions of bodies
	pi							# Initial momenta of bodies
	m							# Masses of bodies

	"Construct a new NBody"
	function NBody( T=40, resolution=20000)
		# Initialise all fields of the decoding apparatus:
		new(
			0,					# Initially no bodies in the system
			T,					# Duration
			resolution,			# Timestep resolution
			zeros(1,N),			# Initial positions
			zeros(1,N),			# Initial momenta
			zeros(1,N),			# Masses of bodies
		)
	end
end

#-----------------------------------------------------------------------------------------
# Module methods:

"""
	addbody( nbody::NBody, xi, pi, m=1)

Add to the system a new body with initial position and momentum xi, pi, and with mass m.
"""
function addbody( nbody, xi, pi, m=1)
	nbody.xi = [nbody.xi xi]
	nbody.pi = [nbody.pi pi]
	nbody.m  = [nbody.m m]
	nbody.N  += 1
end

"""
	simulate( nbody::NBody)

Run a simulation of the given NBody system over duration divided into res timesteps.
"""
function simulate( nbody)
	dt = nbody.T / nbody.res
	t = 0:dt:nbody.T
	x = zeros(1,nbody.N,nbody.res+1)
	p = zeros(1,nbody.N,nbody.res+1)
	G = 1
	Gmm = repmat(G*nbody.m,[1 1 nbody.N])
	Gmm = Gmm .* permute(Gmm,[1 3 2])

	# Relative positions of bodies at positions r:
	relPos = @(r) repmat(permute(r,[1 3 2]),[1 nbody.N 1]) - repmat(r,[1 1 nbody.N])

	# Initialisation:
	x[1,:,1] = nbody.xi
	p[1,:,1] = nbody.pi

	# Simulation:
	for n = 1:nbody.res
		# Simulate timestep n:
		x[1,:,n+1] = x[1,:,n] + dt * p[1,:,n]./nbody.m

		diff = relPos(x(1,:,n))
		absDiff = sqrt(sum(diff.*diff,1))
		invSq = absDiff.^3 + permute(eye(nbody.N),[3 2 1])
		forces = Gmm.*diff ./ invSq
		forces = sum(forces,3)
		p[1,:,n+1] = p[1,:,n] + forces * dt
	end

	(t,x,p)
end

#-----------------------------------------------------------------------------------------
"""
	demo()

Demonstrate simulation of a simple 2-body problem.
"""
function demo()
	nb = NBody( 2.25, 10000)
	nb.addbody( -1, 0)
	nb.addbody( 1, 0)

	[t,x] = nb.simulate()
	plot( ...
		t,squeeze(x(1,1,:)),'b', ...
		t,squeeze(x(1,2,:)),'r' ...
	);
end

end		# of NBodies


text:
const DecayRates = Union{Vector, Matrix, Nothing}

"""
    timeevolution.master_h(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_h as derivative function.

Further information can be found at [`master`](@ref).
"""
function master_h(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
                rates::DecayRates=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:Operator{B,B}}
    check_master(rho0, H, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_nh(tspan, rho0, H, J; <keyword arguments>)

Integrate the master equation with dmaster_nh as derivative function.

In this case the given Hamiltonian is assumed to be the non-hermitian version:
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
Further information can be found at [`master`](@ref).
"""
function master_nh(tspan, rho0::T, Hnh::AbstractOperator{B,B}, J::Vector;
                rates::DecayRates=nothing,
                Hnhdagger::AbstractOperator=dagger(Hnh),
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:Operator{B,B}}
    check_master(rho0, Hnh, J, Jdagger, rates)
    tmp = copy(rho0)
    dmaster_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master(tspan, rho0, H, J; <keyword arguments>)

Time-evolution according to a master equation.

There are two implementations for integrating the master equation:

* [`master_h`](@ref): Usual formulation of the master equation.
* [`master_nh`](@ref): Variant with non-hermitian Hamiltonian.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `J`: Vector containing all jump operators which can be of any arbitrary
        operator type.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::T, H::AbstractOperator{B,B}, J::Vector;
                rates::DecayRates=nothing,
                Jdagger::Vector=dagger.(J),
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:Operator{B,B}}
    isreducible = check_master(rho0, H, J, Jdagger, rates)
    if !isreducible
        tmp = copy(rho0)
        dmaster_h_(t, rho::T, drho::T) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
        return integrate_master(tspan, dmaster_h_, rho0, fout; kwargs...)
    else
        Hnh = copy(H)
        if isa(rates, Matrix)
            for i=1:length(J), j=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im*rates[i,j])*Jdagger[i]*J[j]
            end
        elseif isa(rates, Vector)
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im*rates[i])*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= complex(float(eltype(H)))(0.5im)*Jdagger[i]*J[i]
            end
        end
        Hnhdagger = dagger(Hnh)
        tmp = copy(rho0)
        dmaster_nh_(t, rho::T, drho::T) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
        return integrate_master(tspan, dmaster_nh_, rho0, fout; kwargs...)
    end
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic non-hermitian Hamiltonian and J.

In this case the given Hamiltonian is assumed to be the non-hermitian version.
```math
H_{nh} = H - \\frac{i}{2} \\sum_k J^†_k J_k
```
The given function can either be of the form `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger)`
or `f(t, rho) -> (Hnh, Hnhdagger, J, Jdagger, rates)` For further information look
at [`master_dynamic`](@ref).
"""
function master_nh_dynamic(tspan, rho0::T, f::Function;
                rates::DecayRates=nothing,
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:Operator{B,B}}
    tmp = copy(rho0)
    dmaster_(t, rho::T, drho::T) = dmaster_nh_dynamic(t, rho, f, rates, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

"""
    timeevolution.master_dynamic(tspan, rho0, f; <keyword arguments>)

Time-evolution according to a master equation with a dynamic Hamiltonian and J.

There are two implementations for integrating the master equation with dynamic
operators:

* [`master_dynamic`](@ref): Usual formulation of the master equation.
* [`master_nh_dynamic`](@ref): Variant with non-hermitian Hamiltonian.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `f`: Function `f(t, rho) -> (H, J, Jdagger)` or `f(t, rho) -> (H, J, Jdagger, rates)`
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan, rho0::T, f::Function;
                rates::DecayRates=nothing,
                fout::Union{Function,Nothing}=nothing,
                kwargs...) where {B<:Basis,T<:Operator{B,B}}
    tmp = copy(rho0)
    dmaster_(t, rho::T, drho::T) = dmaster_h_dynamic(t, rho, f, rates, drho, tmp)
    integrate_master(tspan, dmaster_, rho0, fout; kwargs...)
end

# Automatically convert Ket states to density operators
master(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master(tspan, dm(psi0), H, J; kwargs...)
master_h(tspan, psi0::Ket{B}, H::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_h(tspan, dm(psi0), H, J; kwargs...)
master_nh(tspan, psi0::Ket{B}, Hnh::AbstractOperator{B,B}, J::Vector; kwargs...) where B<:Basis = master_nh(tspan, dm(psi0), Hnh, J; kwargs...)
master_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_dynamic(tspan, dm(psi0), f; kwargs...)
master_nh_dynamic(tspan, psi0::Ket{B}, f::Function; kwargs...) where B<:Basis = master_nh_dynamic(tspan, dm(psi0), f; kwargs...)

# Recasting needed for the ODE solver is just providing the underlying data
function recast!(x::T, rho::Operator{B,B,T}) where {B<:Basis,T}
    rho.data = x
end
recast!(rho::Operator{B,B,T}, x::T) where {B<:Basis,T} = nothing

function integrate_master(tspan, df::Function, rho0::T,
                        fout::Union{Nothing, Function}; kwargs...) where {B<:Basis,T<:Operator{B,B}}
    tspan_ = convert(Vector{float(eltype(tspan))}, tspan)
    x0 = rho0.data
    state = T(rho0.basis_l, rho0.basis_r, rho0.data)
    dstate = T(rho0.basis_l, rho0.basis_r, rho0.data)
    integrate(tspan_, df, x0, state, dstate, fout; kwargs...)
end

# Time derivative functions
#   * dmaster_h
#   * dmaster_nh
#   * dmaster_h_dynamic -> callback(t, rho) -> dmaster_h
#   * dmaster_nh_dynamic -> callback(t, rho) -> dmaster_nh
# dmaster_h and dmaster_nh provide specialized implementations depending on
# the type of the given decay rate object which can either be nothing, a vector
# or a matrix.

function dmaster_h(rho::T, H::AbstractOperator{B,B},
                    rates::Nothing, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho)
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[i],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[i],tmp,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(tmp,rho,Jdagger[i],true,false)
        QuantumOpticsBase.mul!(drho,tmp,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

function dmaster_h(rho::T, H::AbstractOperator{B,B},
                    rates::Vector, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho,eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[i],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[i],tmp,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(tmp,rho,Jdagger[i],eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

function dmaster_h(rho::T, H::AbstractOperator{B,B},
                    rates::Matrix, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,H,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,H,eltype(rho)(im),one(eltype(rho)))
    for j=1:length(J), i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho,eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[j],true,true)

        QuantumOpticsBase.mul!(drho,Jdagger[j],tmp,eltype(rho)(-0.5),one(eltype(rho)))

        QuantumOpticsBase.mul!(tmp,rho,Jdagger[j],eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,J[i],eltype(rho)(-0.5),one(eltype(rho)))
    end
    return drho
end

function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
                    rates::Nothing, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho)
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[i],true,true)
    end
    return drho
end

function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
                    rates::Vector, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho,eltype(rho)(rates[i]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[i],true,true)
    end
    return drho
end

function dmaster_nh(rho::T, Hnh::AbstractOperator{B,B}, Hnh_dagger::AbstractOperator{B,B},
                    rates::Matrix, J::Vector, Jdagger::Vector,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    QuantumOpticsBase.mul!(drho,Hnh,rho,-eltype(rho)(im),zero(eltype(rho)))
    QuantumOpticsBase.mul!(drho,rho,Hnh_dagger,eltype(rho)(im),one(eltype(rho)))
    for j=1:length(J), i=1:length(J)
        QuantumOpticsBase.mul!(tmp,J[i],rho,eltype(rho)(rates[i,j]),zero(eltype(rho)))
        QuantumOpticsBase.mul!(drho,tmp,Jdagger[j],true,true)
    end
    return drho
end

function dmaster_h_dynamic(t, rho::T, f::Function,
                    rates::DecayRates,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    result = f(t, rho)
    QO_CHECKS[] && @assert 3 <= length(result) <= 4
    if length(result) == 3
        H, J, Jdagger = result
        rates_ = rates
    else
        H, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_master(rho, H, J, Jdagger, rates_)
    dmaster_h(rho, H, rates_, J, Jdagger, drho, tmp)
end

function dmaster_nh_dynamic(t, rho::T, f::Function,
                    rates::DecayRates,
                    drho::T, tmp::T) where {B<:Basis,T<:Operator{B,B}}
    result = f(t, rho)
    QO_CHECKS[] && @assert 4 <= length(result) <= 5
    if length(result) == 4
        Hnh, Hnh_dagger, J, Jdagger = result
        rates_ = rates
    else
        Hnh, Hnh_dagger, J, Jdagger, rates_ = result
    end
    QO_CHECKS[] && check_master(rho, Hnh, J, Jdagger, rates_)
    dmaster_nh(rho, Hnh, Hnh_dagger, rates_, J, Jdagger, drho, tmp)
end

function check_master(rho0::Operator{B,B}, H::AbstractOperator{B,B}, J::Vector, Jdagger::Vector, rates::DecayRates) where B<:Basis
    # TODO: clean up type checks by dispatch; make type of J known
    isreducible = true # test if all operators are sparse or dense
    if !(isa(H, DenseOpType) || isa(H, SparseOpType))
        isreducible = false
    end
    for j=J
        @assert isa(j, AbstractOperator{B,B})
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    for j=Jdagger
        @assert isa(j, AbstractOperator{B,B})
        if !(isa(j, DenseOpType) || isa(j, SparseOpType))
            isreducible = false
        end
        check_samebases(rho0, j)
    end
    @assert length(J)==length(Jdagger)
    if isa(rates,Matrix)
        @assert size(rates, 1) == size(rates, 2) == length(J)
    elseif isa(rates,Vector)
        @assert length(rates) == length(J)
    end
    isreducible
end

get_type(rho, H, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype.(J)...,eltype.(Jdagger)...)
get_type(rho, H, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(H),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)
get_type(rho, Hnh, Hnhdagger, rates::Nothing, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger), eltype.(J)...,eltype.(Jdagger)...)
get_type(rho, Hnh, Hnhdagger, rates::Union{Vector,Matrix}, J, Jdagger) = promote_type(eltype(rho),eltype(Hnh),eltype(Hnhdagger),eltype(rates),eltype.(J)...,eltype.(Jdagger)...)


text:
# ---
# title: 1009. Complement of Base 10 Integer
# id: problem1009
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/complement-of-base-10-integer/description/>
# hidden: true
# ---
#
# Every non-negative integer `N` has a binary representation.  For example, `5`
# can be represented as `"101"` in binary, `11` as `"1011"` in binary, and so
# on.  Note that except for `N = 0`, there are no leading zeroes in any binary
# representation.
#
# The _complement_  of a binary representation is the number in binary you get
# when changing every `1` to a `0` and `0` to a `1`.  For example, the
# complement of `"101"` in binary is `"010"` in binary.
#
# For a given number `N` in base-10, return the complement of it's binary
# representation as a base-10 integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: 5
#     Output: 2
#     Explanation: 5 is "101" in binary, with complement "010" in binary, which is 2 in base-10.
#
#
# **Example 2:**
#
#
#
#     Input: 7
#     Output: 0
#     Explanation: 7 is "111" in binary, with complement "000" in binary, which is 0 in base-10.
#
#
# **Example 3:**
#
#
#
#     Input: 10
#     Output: 5
#     Explanation: 10 is "1010" in binary, with complement "0101" in binary, which is 5 in base-10.
#
#
#
#
# **Note:**
#
#   1. `0 <= N < 10^9`
#   2. This question is the same as 476: <https://leetcode.comhttps://leetcode.com/problems/number-complement/>
#
#
## @lc code=start
using LeetCode

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


text:
using RandomMatrices
using Random: seed!
using Test

@testset "RandomMatrices" begin
	seed!(1)
	include("GaussianEnsembles.jl")
	include("FormalPowerSeries.jl")
	include("Haar.jl")
	include("StochasticProcess.jl")

	@testset "densities" begin
		include("densities/Semicircle.jl")
		include("densities/TracyWidom.jl")
	end
end


text:
optimal = DiscretizeMODL_Optimal()
greedy  = DiscretizeMODL_Greedy()
post_greedy = DiscretizeMODL_PostGreedy()

# Case1
continuous = [1.0]; class_value = [1];
@test binedges(optimal,continuous,class_value)     == [1.0,1.0]
@test binedges(greedy,continuous,class_value)      == [1.0,1.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,1.0]

# Case2
continuous = collect(1.0:1.0:4.0); class_value = [1,1,2,2];
@test binedges(optimal,continuous,class_value)     == [1.0,4.0]
@test binedges(greedy,continuous,class_value)      == [1.0,4.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,4.0]

# Case3
continuous = collect(1.0:1.0:10.0); class_value = [1,1,1,2,1,2,1,2,2,2];
@test binedges(optimal,continuous,class_value)     == [1.0,10.0]
@test binedges(greedy,continuous,class_value)      == [1.0,10.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,10.0]

# Case4
continuous = collect(1.0:1.0:20.0)
class_value = [1,1,1,2,2,2,2,2,2,2,1,1,1,1,2,2,2,2,2,2];
@test binedges(optimal,continuous,class_value)     == [1.0,20.0]
@test binedges(greedy,continuous,class_value)      == [1.0,3.5,10.5,14.5,20.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,20.0]

# Case5
continuous = collect(1.0:1.0:20.0)
class_value = [1,1,1,2,2,2,2,2,2,2,3,3,3,3,3,2,2,2,2,2];
@test binedges(optimal,continuous,class_value)     == [1.0,3.5,20.0]
@test binedges(greedy,continuous,class_value)      == [1.0,3.5,10.5,15.5,20.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,3.5,20.0]

#Case6
continuous = collect(1.0:1.0:100.0)
class_value = Array{Int64}(100)
for i in 1 : 100
    class_value[i] = div((i-1),20)
end
@test binedges(optimal,continuous,class_value)     == [1.0,20.5,40.5,60.5,80.5,100.0]
@test binedges(greedy,continuous,class_value)      == [1.0,20.5,40.5,60.5,80.5,100.0]
@test binedges(post_greedy,continuous,class_value) == [1.0,20.5,40.5,60.5,80.5,100.0]
post_greedy.max_bin_count = 2
@test binedges(post_greedy,continuous,class_value) == [1.0,40.5,100.0]
post_greedy.max_bin_count = 3
@test binedges(post_greedy,continuous,class_value) == [1.0,40.5,80.5,100.0]
post_greedy.max_bin_count = 4
@test binedges(post_greedy,continuous,class_value) == [1.0,20.5,40.5,60.5,100.0]

#Case7
continuous = collect(1.0:1.0:100.0)
class_value = [1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
               1,1,2,3,2,1,2,2,1,1,2,2,2,2,2,2,2,2,3,3,
               3,3,3,3,3,3,3,2,2,3,2,3,2,3,3,3,3,2,2,2,
               1,1,1,1,2,2,2,2,1,1,1,1,3,3,3,3,4,4,4,4,
               4,4,4,4,4,4,4,4,1,1,4,1,4,1,4,4,4,4,4,1];

@test binedges(optimal,continuous,class_value)     == [1.0,12.5,38.5,57.5,72.5,76.5,100.0]
@test binedges(greedy,continuous,class_value)      == [1.0,12.5,38.5,60.5,72.5,76.5,100.0]
post_greedy.max_bin_count = 0
@test binedges(post_greedy,continuous,class_value) == [1.0,12.5,38.5,57.5,72.5,76.5,100.0]
post_greedy.max_bin_count = 2
@test binedges(post_greedy,continuous,class_value) == [1.0,76.5,100.0]
post_greedy.max_bin_count = 3
@test binedges(post_greedy,continuous,class_value) == [1.0,12.5,76.5,100.0]
post_greedy.max_bin_count = 4
@test binedges(post_greedy,continuous,class_value) == [1.0,12.5,38.5,76.5,100.0]
post_greedy.max_bin_count = 5
@test binedges(post_greedy,continuous,class_value) == [1.0,12.5,38.5,57.5,76.5,100.0]


text:
# EP simulations with i.i.d. and correlated Gaussian patterns and noisy labels

using EPDilutedClassifier

N=128; ρ=0.25; α=6.0;
epsconv=1e-4; maxiter=10000; η=0.95;
damp=0.99;
λ0=1e4
w=generate_weights(N,ρ);
report=Reporter(w)

# iid patterns
println("Test with i.i.d. Gaussian patterns")
run_instance_gaussian_iid(w,ρ,α,η,callback=report,ρ0=ρ,λ0=λ0,δρ=0.0,η0=η,δη=0.0,
        damp=damp,epsconv=epsconv,maxiter=maxiter)

# gaussian correlated patterns
println("Test with correlated Gaussian patterns")
run_instance_gaussian_cor(w,ρ,α,η,callback=report,ρ0=ρ,λ0=λ0,δρ=0.0,η0=η,δη=0.0,
        damp=damp,epsconv=epsconv,maxiter=maxiter)

max_stars_repo_name: mpieropan/EPDilutedClassifier

max_issues_repo_name: mpieropan/EPDilutedClassifier

max_forks_repo_name: mpieropan/EPDilutedClassifier


text:
using FastInverseSqrt
using Test

@testset "FastInverseSqrt.jl" begin
    # Write your tests here.
end