text:
function estimate_noise_level(observations::Observations, 𝛉::AbstractVector)
    @unpack data = observations
    𝛉 = 𝛉 / norm(𝛉)
    ℳ = data[1]
    N = length(ℳ)
    aml = 0.0
    for n = 1:N
        𝐦 = ℳ[n]
        𝐮ₙ = SA_F64[𝐦[1]^2, 𝐦[1]*𝐦[2], 𝐦[2]^2, 𝐦[1], 𝐦[2], 1]
        ∂𝐮ₙ = SA_F64[2*𝐦[1]  𝐦[2]  0  1  0  0; 0 𝐦[1] 2*𝐦[2] 0 1 0]'
        𝐀ₙ = 𝐮ₙ * 𝐮ₙ'
        𝚲ₙ =  SA_F64[1 0 ; 0 1]
        𝐁ₙ = ∂𝐮ₙ * 𝚲ₙ * ∂𝐮ₙ'
        aml = aml + abs((𝛉' * 𝐀ₙ * 𝛉)/(𝛉' * 𝐁ₙ * 𝛉))
    end
    σ² = aml / (N-5)
    return σ²
end

function determine_algebraic_covariance(observations::Observations, 𝛉₀::AbstractVector, estimator::GuaranteedEllipseFit)
    @unpack data = observations
    # We derive the covariance matrix of the estimate under the unit-norm Gauge constraint.
    𝛉₀ = 𝛉₀ / norm(𝛉₀)
    N = length(first(data))
    σ² = estimate_noise_level(observations, 𝛉₀)
    Λ = [SA_F64[σ² 0.0; 0.0 σ²] for n = 1:N]
    uncertain_observations = UncertainObservations(data, tuple(Λ))
    return determine_algebraic_covariance(uncertain_observations, 𝛉₀, estimator)
end

function determine_algebraic_covariance(observations::UncertainObservations, 𝛉₀::AbstractVector, estimator::GuaranteedEllipseFit)
     # Convert observations and covariance matrices to a data-driven normalised coordinate system.
    normalise = NormaliseDataContext(observations, IsotropicScalingTranslation())
    𝒯 = matrices(normalise)
    𝐓 = 𝒯[1]
    𝛉₁ = normalise(ToNormalisedSpace(), 𝛉₀)
    @unpack data, covariance_matrices = normalise(observations)
    ℳ = data[1]
    Λ = first(covariance_matrices)
    𝐌 = zeros(6,6)
    N = length(first(data))
    for n = 1:N
        𝐦 = ℳ[n]
        𝚲ₙ = Λ[n]
        𝐮ₙ = SA_F64[𝐦[1]^2, 𝐦[1]*𝐦[2], 𝐦[2]^2, 𝐦[1], 𝐦[2], 1]
        ∂𝐮ₙ = SA_F64[2*𝐦[1]  𝐦[2]  0  1  0  0; 0 𝐦[1] 2*𝐦[2] 0 1 0]'
        𝐀ₙ = 𝐮ₙ * 𝐮ₙ'
        𝐁ₙ = ∂𝐮ₙ * 𝚲ₙ * ∂𝐮ₙ'
        𝐌 = 𝐌 + 𝐀ₙ / (𝛉₁' * 𝐁ₙ * 𝛉₁)
    end
    𝐏ₜ = UniformScaling(1) - (𝛉₁*𝛉₁') / norm(𝛉₁)^2
    # Compute rank-5 constrained pseudo-inverse of 𝐌.
    F = svd(𝐌)
    𝐒 = vcat(SVector([1 / F.S[i] for i = 1:5]...), 0)
    𝐌⁻¹ = F.U * Diagonal(𝐒) * F.V'
    𝚺₁ = 𝐏ₜ * 𝐌⁻¹ * 𝐏ₜ

    # Matrices used to transform from normalised to unnormalised (original) coordinate system.
    𝐄 = Diagonal(SVector(1, 2^-1, 1, 2^-1, 2^-1, 1))
    # Permutation matrix for interchanging the 3rd and 4th entries of a length-6 vector.
    𝐏₃₄ = Diagonal(SVector(0,1,0)) ⊗ SMatrix{2,2,Float64}(0,1,1,0) + Diagonal(SVector(1,0,1)) ⊗ SMatrix{2,2,Float64}(1,0,0,1)
    # 9 x 6 duplication matrix
    𝐃₃ = [1 0 0 0 0 0;
          0 1 0 0 0 0;
          0 0 1 0 0 0;
          0 1 0 0 0 0;
          0 0 0 1 0 0;
          0 0 0 0 1 0;
          0 0 1 0 0 0;
          0 0 0 0 1 0;
          0 0 0 0 0 1]
    𝐅 = 𝐄 \ 𝐏₃₄ * pinv(𝐃₃) * kron(𝐓, 𝐓)' * 𝐃₃ * 𝐏₃₄ * 𝐄
    # 𝛉₀ is not unit-normalised
    𝛉₀ = 𝐅 * 𝛉₁
    𝐏ₜ = UniformScaling(1) - (𝛉₀*𝛉₀') / norm(𝛉₀)^2
    𝚺₀ = norm(𝛉₀)^-2 * 𝐏ₜ * 𝐅 * 𝚺₁ * 𝐅' * 𝐏ₜ
    return 𝚺₀
end

function determine_geometric_covariance(observations::Union{Observations,UncertainObservations}, 𝛏::AbstractVector, estimator::GuaranteedEllipseFit)
    geo_to_alg = GeometricToAlgebraic()
    𝛉 =  geo_to_alg(𝛏)
    𝛉 = 𝛉 / norm(𝛉)
    a, b, c, d, e ,f = 𝛉

    Δ = b^2 - 4*a*c
    λ₊ = 0.5*(a + c - (b^2 + (a - c)^2)^0.5)
    λ₋ = 0.5*(a + c + (b^2 + (a - c)^2)^0.5)

    ψ = b*d*e - a*e^2 - b^2*f + c*(4*a*f - d^2)
    V₊ = (ψ/(λ₊*Δ))^0.5
    V₋ = (ψ/(λ₋*Δ))^0.5

    A = max(V₊, V₋)
    B = min(V₊, V₋)

    ∂A = A == V₊ ? ∂V₊(𝛉, ψ, λ₊, Δ) : ∂V₋(𝛉, ψ, λ₋, Δ)
    ∂B = A == V₊ ? ∂V₋(𝛉, ψ, λ₋, Δ) : ∂V₊(𝛉, ψ, λ₊, Δ)

    ∂𝛏 = vcat(∂A', ∂B', ∂H(𝛉, Δ)', ∂K(𝛉, Δ)',  ∂τ(𝛉)')

    𝚺₀ = determine_algebraic_covariance(observations, 𝛉, estimator)
    return ∂𝛏*𝚺₀*∂𝛏'
end

function ∂H(𝛉::AbstractVector, Δ::Number)
    a, b, c, d, e, f = 𝛉
    ∂H_a = (4*c*(2*c*d-b*e))/Δ^2
    ∂H_b = (b^2*e+4*a*c*e-4*b*c*d)/Δ^2
    ∂H_c = (2*b*(b*d-2*a*e))/Δ^2
    ∂H_d = (2*c)/Δ
    ∂H_e = -b/Δ
    ∂H_f =  0
    return SVector(∂H_a, ∂H_b, ∂H_c, ∂H_d, ∂H_e, ∂H_f)
end

function ∂K(𝛉::AbstractVector, Δ::Number)
    a, b, c, d, e, f = 𝛉
    ∂K_a = (2*b*(b*e-2*c*d))/Δ^2
    ∂K_b = (b^2*d+4*a*c*d-4*a*b*e)/Δ^2
    ∂K_c = (4*a*(2*a*e-b*d))/Δ^2;
    ∂K_d = -b/Δ
    ∂K_e =  2*a/Δ
    ∂K_f =  0
    return SVector(∂K_a, ∂K_b, ∂K_c, ∂K_d, ∂K_e, ∂K_f)
end

function ∂τ(𝛉::AbstractVector)
    a, b, c, d, e, f = 𝛉
    ∂τ_a = -b/(2*(b^2+(a-c)^2))
    ∂τ_b = (a-c)/(2*(b^2+(a-c)^2))
    ∂τ_c = b/(2*(b^2+(a-c)^2))
    ∂τ_d = 0
    ∂τ_e = 0
    ∂τ_f = 0
    return SVector(∂τ_a, ∂τ_b, ∂τ_c, ∂τ_d, ∂τ_e, ∂τ_f)
end

function ∂V₊(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    return SVector(∂V₊_a(𝛉, ψ, λ₊, Δ),
                   ∂V₊_b(𝛉, ψ, λ₊, Δ),
                   ∂V₊_c(𝛉, ψ, λ₊, Δ),
                   ∂V₊_d(𝛉, ψ, λ₊, Δ),
                   ∂V₊_e(𝛉, ψ, λ₊, Δ),
                   ∂V₊_f(𝛉, ψ, λ₊, Δ))
end

function ∂V₊_a(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₊*Δ);
    part2 = (ψ/(λ₊*Δ))^(-0.5);
    part3 = 4*c*f - e^2 + 4*Δ^(-1)*c*ψ;
    part4 = (ψ/(2*λ₊))*(1 + ((c - a)/((a - c)^2 + b^2)^0.5));
    return part1*part2*(part3 - part4)
end

function ∂V₊_b(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₊*Δ)
    part2 = (ψ/(λ₊*Δ))^(-0.5)
    part3 = d*e - 2*b*f - 2*Δ^(-1)*b*ψ
    part4 = (b*ψ)/((2*λ₊)*((a - c)^2 + b^2)^0.5)
    return  part1*part2*(part3 + part4)
end

function ∂V₊_c(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₊*Δ)
    part2 = (ψ/(λ₊*Δ))^(-0.5)
    part3 = 4*a*f - d^2 + 4*Δ^(-1)*a*ψ
    part4 = (ψ/(2*λ₊))*(1 + ((a - c)/((a - c)^2 + b^2)^0.5))
    return part1*part2*(part3 - part4)
end

function ∂V₊_d(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = (ψ/(λ₊*Δ))^(0.5)
    part2 = (b*e - 2*c*d)/(2*ψ)
    return part1*part2;
end

function ∂V₊_e(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = (ψ/(λ₊*Δ))^(0.5)
    part2 = (b*d - 2*a*e)/(2*ψ)
    return part1*part2;
end

function ∂V₊_f(𝛉::AbstractVector, ψ::Number, λ₊::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = -(2*λ₊)^-1
    part2 = (ψ/(λ₊*Δ))^(-0.5)
    return part1*part2
end

function ∂V₋(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    return SVector(∂V₋_a(𝛉, ψ, λ₋, Δ),
                   ∂V₋_b(𝛉, ψ, λ₋, Δ),
                   ∂V₋_c(𝛉, ψ, λ₋, Δ),
                   ∂V₋_d(𝛉, ψ, λ₋, Δ),
                   ∂V₋_e(𝛉, ψ, λ₋, Δ),
                   ∂V₋_f(𝛉, ψ, λ₋, Δ))
end

function ∂V₋_a(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₋*Δ)
    part2 = (ψ/(λ₋*Δ))^(-0.5)
    part3 = 4* c*f - e^2 + 4*Δ^(-1)*c*ψ
    part4 = (ψ/(2*λ₋))*(1 - ((c - a)/((a - c)^2 + b^2)^0.5))
    return part1*part2*(part3 - part4)
end

function ∂V₋_b(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₋*Δ)
    part2 = (ψ/(λ₋*Δ))^(-0.5)
    part3 = d*e - 2*b*f - 2*Δ^(-1)*b*ψ
    part4 = (b*ψ)/((2*λ₋)*((a - c)^2 + b^2)^0.5)
    return part1*part2*(part3 - part4)
end

function ∂V₋_c(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = 1/(2*λ₋*Δ)
    part2 = (ψ/(λ₋*Δ))^(-0.5)
    part3 = 4*a*f - d^2 + 4*Δ^(-1)*a*ψ
    part4 = (ψ/(2*λ₋))*(1 - ((a - c)/((a - c)^2 + b^2)^0.5))
    return  part1*part2*(part3 - part4)
end

function ∂V₋_d(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = (ψ/(λ₋*Δ))^(0.5)
    part2 = (b*e - 2*c*d)/(2*ψ)
    return part1*part2
end

function ∂V₋_e(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = (ψ/(λ₋*Δ))^(0.5)
    part2 = (b*d - 2*a*e)/(2*ψ)
    return part1*part2
end

function ∂V₋_f(𝛉::AbstractVector, ψ::Number, λ₋::Number, Δ::Number)
    a, b, c, d, e, f = 𝛉
    part1 = -(2*λ₋)^-1
    part2 = (ψ/(λ₋*Δ))^(-0.5)
    return part1*part2
end


text:
dummy = cosd.(latvec(grd)) + sqrt.(depthvec(grd)) / 30
xs = (dummy, dummy*u"mol/m^3")

nokwplots = (surfacemap, plotverticalintegral, plotverticalaverage, plotzonalaverage, plotzonalintegral, plotmeridionalaverage, plotmeridionalintegral, plothorizontalintegral, plothorizontalaverage)
@testset "plots without kwargs" begin
    @testset "$f(x, grd)" for f in nokwplots, x in xs
        @test f(x, grd) isa Plots.Plot
    end
end

@testset "Plots with kwargs" begin
    lons = (330, 330°)
    lats = (30, 30°)
    depths = (500, 500m, 5km)
    lonlats = ((-30,30), (30,30), (-30°,30°))
    @testset "plothorizontalslice depth=$depth" for x in xs, depth in depths
        @test plothorizontalslice(x, grd, depth=depth) isa Plots.Plot
    end
    @testset "plotzonalslice lat=$lat" for x in xs, lat in lats
        @test plotzonalslice(x, grd, lat=lat) isa Plots.Plot
    end
    @testset "plotmeridionalslice lon=$lon" for x in xs, lon in lons
        @test plotmeridionalslice(x, grd, lon=lon) isa Plots.Plot
    end
    @testset "plotdepthprofile lonlat=$lonlat" for x in xs, lonlat in lonlats
        @test plotdepthprofile(x, grd, lonlat=lonlat) isa Plots.Plot
    end
end

ext: jl


text:
################## implement AbstractTrees interface #######################
# refer to https://github.com/JuliaCollections/AbstractTrees.jl for more details
function AbstractTrees.children(ver4::Ver4)
    return ver4.bubble
end

function AbstractTrees.children(bubble::Bubble)
    return (bubble.Lver, bubble.Rver)
end

function iterate(ver4::Ver4{W}) where {W}
    if length(ver4.bubble) == 0
        return nothing
    else
        return (ver4.bubble[1], 1)
    end
end

function iterate(bub::Bubble)
    return (bub.Lver, false)
end

function iterate(ver4::Ver4{W}, state) where {W}
    if state >= length(ver4.bubble) || length(ver4.bubble) == 0
        return nothing
    else
        return (ver4.bubble[state+1], state + 1)
    end
end

function iterate(bub::Bubble, state::Bool)
    state && return nothing
    return (bub.Rver, true)
end

Base.IteratorSize(::Type{Ver4{W}}) where {W} = Base.SizeUnknown()
Base.eltype(::Type{Ver4{W}}) where {W} = Ver4{W}

Base.IteratorSize(::Type{Bubble{Ver4{W}}}) where {W} = Base.SizeUnknown()
Base.eltype(::Type{Bubble{Ver4{W}}}) where {W} = Bubble{Ver4{W}}

AbstractTrees.printnode(io::IO, ver4::Ver4) = print(io, tpair(ver4))
AbstractTrees.printnode(io::IO, bub::Bubble) = print(io,
    "\u001b[32m$(bub.id): $(bub.chan) $(bub.Lver.para.innerLoopNum)Ⓧ $(bub.Rver.para.innerLoopNum)\u001b[0m")

function tpair(ver4, MaxT = 18)
    s = "\u001b[31m$(ver4.id):\u001b[0m"
    if ver4.para.innerLoopNum > 0
        s *= "$(ver4.para.innerLoopNum)lp, T$(length(ver4.Tpair))⨁ "
    else
        s *= "⨁ "
    end
    # if ver4.loopNum <= 1
    for (ti, T) in enumerate(ver4.Tpair)
        if ti <= MaxT
            s *= "($(T[1]),$(T[2]),$(T[3]),$(T[4]))"
        else
            s *= "..."
            break
        end
    end
    # end
    return s
end

##### pretty print of Bubble and Ver4  ##########################
Base.show(io::IO, bub::Bubble) = AbstractTrees.printnode(io::IO, bub)
Base.show(io::IO, ver4::Ver4) = AbstractTrees.printnode(io::IO, ver4)

"""
convert Ver4 tree struct to a string in the newick format
"""
function newick(ver4::Ver4)
    return "$(newickVer4(ver4));"
end

"""
convert Bubble tree struct to a string in the newick format
"""
function newick(bub::Bubble)
    return "$(newickBuble(bub));"
end

function newickBubble(bub::Bubble)
    # Recursive version
    # Practically a postorder tree traversal
    left = newickVer4(bub.Lver)
    right = newickVer4(bub.Rver)
    return "($left,$right)$(bub.id)_$(ChanName[bub.chan])_$(bub.Lver.para.innerLoopNum)Ⓧ$(bub.Rver.para.innerLoopNum)"
end

function newickVer4(ver4::Ver4)
    # Recursive version
    # Practically a postorder tree traversal

    function tpairNewick(ver4)
        if ver4.para.innerLoopNum > 0
            s = "$(ver4.id):lp$(ver4.para.innerLoopNum)_T$(length(ver4.Tpair))⨁"
        else
            s = "$(Ver4.id):⨁"
        end
        # if ver4.loopNum <= 1
        for (ti, T) in enumerate(ver4.Tpair)
            if ti <= 5
                s *= "⟨$(T[1])-$(T[2])-$(T[3])-$(T[4])⟩"
            else
                s *= "…"
                break
            end
        end
        # end
        return s
    end

    if ver4.para.innerLoopNum == 0
        return tpairNewick(ver4)
    else
        s = "("
        for (bi, bub) in enumerate(ver4.bubble)
            s *= newickBubble(bub)
            if bi != length(ver4.bubble)
                s *= ","
            else
                s *= ")"
            end
        end
        return s * tpairNewick(ver4)
    end
end

"""
    showTree(ver4, para::Para; verbose=0, depth=999)

    Visualize the diagram tree using ete3 python package

#Arguments
- `ver4`: the 4-vertex diagram tree to visualize
- `para`: parameters
- `verbose=0`: the amount of information to show
- `depth=999`: deepest level of the diagram tree to show
"""
function showTree(ver4; verbose = 0, depth = 999)

    # pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__) #comment this line if no need to load local python module
    ete = pyimport("ete3")

    function tpairETE(ver4, depth)
        s = "$(ver4.id):"
        if ver4.para.innerLoopNum > 0
            s *= "$(ver4.para.innerLoopNum)lp, T$(length(ver4.Tpair))⨁ "
        else
            s *= "⨁ "
        end
        if ver4.para.innerLoopNum == 0 || ver4.level > depth
            MaxT = Inf
        else
            MaxT = 1
        end

        # if ver4.loopNum <= 1
        for (ti, T) in enumerate(ver4.Tpair)
            if ti <= MaxT
                s *= "($(T[1]),$(T[2]),$(T[3]),$(T[4]))"
            else
                s *= "..."
                break
            end
        end
        # end
        return s
    end

    function treeview(ver4, t = nothing)
        if isnothing(t)
            t = ete.Tree(name = " ")
        end

        if ver4.para.innerLoopNum == 0 || ver4.level > depth
            nt = t.add_child(name = tpairETE(ver4, depth))
            return t
        else
            # prefix = "$(ver4.id): $(ver4.loopNum) lp, $(length(ver4.Tpair)) elem"
            # nt = t.add_child(name=prefix * ", ⨁")
            nt = t.add_child(name = tpairETE(ver4, depth))
            name_face = ete.TextFace(nt.name, fgcolor = "black", fsize = 10)
            nt.add_face(name_face, column = 0, position = "branch-top")
        end

        for bub in ver4.bubble
            nnt = nt.add_child(name = "$(bub.id): $(bub.chan) $(bub.Lver.para.innerLoopNum)Ⓧ$(bub.Rver.para.innerLoopNum)")

            name_face = ete.TextFace(nnt.name, fgcolor = "black", fsize = 10)
            nnt.add_face(name_face, column = 0, position = "branch-top")

            treeview(bub.Lver, nnt)
            treeview(bub.Rver, nnt)
        end

        return t
    end

    t = treeview(ver4)
    # style = ete.NodeStyle()
    # style["bgcolor"] = "Khaki"
    # t.set_style(style)

    ts = ete.TreeStyle()
    ts.show_leaf_name = true
    # ts.show_leaf_name = True
    # ts.layout_fn = my_layout
    ####### show tree vertically ############
    # ts.rotation = 90 #show tree vertically

    ####### show tree in an arc  #############
    # ts.mode = "c"
    # ts.arc_start = -180
    # ts.arc_span = 180
    # t.write(outfile="/home/kun/test.txt", format=8)
    t.show(tree_style = ts)
end


text:
using OffsetArrays
using OffsetArrays: IdentityUnitRange, no_offset_view
using OffsetArrays: IdOffsetRange
using Test, Aqua, Documenter
using LinearAlgebra
using DelimitedFiles
using CatIndices: BidirectionalVector
using EllipsisNotation
using Adapt
using StaticArrays

DocMeta.setdocmeta!(OffsetArrays, :DocTestSetup, :(using OffsetArrays); recursive=true)

# https://github.com/JuliaLang/julia/pull/29440
if VERSION < v"1.1.0-DEV.389"
    Base.:(:)(I::CartesianIndex{N}, J::CartesianIndex{N}) where N =
        CartesianIndices(map((i,j) -> i:j, Tuple(I), Tuple(J)))
end

# Custom index types
struct ZeroBasedIndexing end
struct NewColon end
struct TupleOfRanges{N}
    x ::NTuple{N, UnitRange{Int}}
end

function same_value(r1, r2)
    length(r1) == length(r2) || return false
    for (v1, v2) in zip(r1, r2)
        v1 == v2 || return false
    end
    return true
end

@testset "Project meta quality checks" begin
    # Not checking compat section for test-only dependencies
    Aqua.test_all(OffsetArrays; project_extras=true, deps_compat=true, stale_deps=true, project_toml_formatting=true)
    if VERSION >= v"1.2"
        doctest(OffsetArrays, manual = false)
    end
end

@testset "IdOffsetRange" begin

    function check_indexed_by(r, rindx)
        for i in rindx
            r[i]
        end
        @test_throws BoundsError r[minimum(rindx)-1]
        @test_throws BoundsError r[maximum(rindx)+1]
        return nothing
    end

    ro = OffsetArrays.IdOffsetRange(Base.OneTo(3))
    rs = OffsetArrays.IdOffsetRange(3:5, -2)
    @test typeof(ro) !== typeof(rs)
    @test same_value(ro, 1:3)
    check_indexed_by(ro, 1:3)
    @test same_value(rs, 1:3)
    check_indexed_by(rs, -1:1)
    @test firstindex(ro) == 1
    @test lastindex(ro) == 3
    @test firstindex(rs) == -1
    @test lastindex(rs) == 1
    @test @inferred(typeof(ro)(ro)) === ro
    @test @inferred(OffsetArrays.IdOffsetRange{Int}(ro))   === ro
    @test @inferred(OffsetArrays.IdOffsetRange{Int16}(ro)) === OffsetArrays.IdOffsetRange(Base.OneTo(Int16(3)))
    @test @inferred(OffsetArrays.IdOffsetRange(ro))        === ro
    @test parent(ro) === ro.parent
    @test parent(rs) === rs.parent
    # construction/coercion preserves the values, altering the axes if needed
    r2 = @inferred(typeof(rs)(ro))
    @test typeof(r2) === typeof(rs)
    @test same_value(ro, 1:3)
    check_indexed_by(ro, 1:3)
    r2 = @inferred(typeof(ro)(rs))
    @test typeof(r2) === typeof(ro)
    @test same_value(r2, 1:3)
    check_indexed_by(r2, 1:3)
    # check the example in the comments
    r = OffsetArrays.IdOffsetRange{Int,UnitRange{Int}}(3:4)
    @test same_value(r, 3:4)
    check_indexed_by(r, 1:2)
    r = OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}(3:4)
    @test same_value(r, 3:4)
    check_indexed_by(r, 3:4)
    r = OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}(3:4, -2)
    @test same_value(r, 1:2)
    check_indexed_by(r, 1:2)

    r = OffsetArrays.IdOffsetRange{Int32, Base.OneTo{Int32}}(Base.OneTo(Int64(2)), 3)
    @test same_value(r, 4:5)
    check_indexed_by(r, 4:5)

    r = IdOffsetRange{Int, UnitRange{Int}}(IdOffsetRange(3:5, 2), 2)
    @test typeof(r) == IdOffsetRange{Int, UnitRange{Int}}
    @test same_value(r, 7:9)
    check_indexed_by(r, 5:7)

    r = IdOffsetRange{Int, Base.OneTo{Int}}(IdOffsetRange(Base.OneTo(3), 1), 1)
    @test typeof(r) == IdOffsetRange{Int,Base.OneTo{Int}}
    @test same_value(r, 3:5)
    check_indexed_by(r, 3:5)

    rp = Base.OneTo(3)
    r = IdOffsetRange(rp)
    r2 = IdOffsetRange{Int,typeof(r)}(r, 1)
    @test same_value(r2, 2:4)
    check_indexed_by(r2, 2:4)

    r2 = IdOffsetRange{Int32,IdOffsetRange{Int32,Base.OneTo{Int32}}}(r, 1)
    @test typeof(r2) == IdOffsetRange{Int32,IdOffsetRange{Int32,Base.OneTo{Int32}}}
    @test same_value(r2, 2:4)
    check_indexed_by(r2, 2:4)

    # conversion preserves both the values and the axes, throwing an error if this is not possible
    @test @inferred(oftype(ro, ro)) === ro
    @test @inferred(convert(OffsetArrays.IdOffsetRange{Int}, ro)) === ro
    @test @inferred(convert(OffsetArrays.IdOffsetRange{Int}, rs)) === rs
    @test @inferred(convert(OffsetArrays.IdOffsetRange{Int16}, ro)) === OffsetArrays.IdOffsetRange(Base.OneTo(Int16(3)))
    r2 = @inferred(oftype(rs, ro))
    @test typeof(r2) === typeof(rs)
    @test same_value(r2, 1:3)
    check_indexed_by(r2, 1:3)
    # These two broken tests can be fixed by uncommenting the `convert` definitions
    # in axes.jl, but unfortunately Julia may not quite be ready for this. (E.g. `reinterpretarray.jl`)
    @test_broken try oftype(ro, rs); false catch err true end  # replace with line below
    # @test_throws ArgumentError oftype(ro, rs)
    @test @inferred(oftype(ro, Base.OneTo(2))) === OffsetArrays.IdOffsetRange(Base.OneTo(2))
    @test @inferred(oftype(ro, 1:2)) === OffsetArrays.IdOffsetRange(Base.OneTo(2))
    @test_broken try oftype(ro, 3:4); false catch err true end
    # @test_throws ArgumentError oftype(ro, 3:4)

    # broadcasting behavior with scalars (issue #104)
    r3 = (1 .+ OffsetArrays.IdOffsetRange(3:5, -1) .+ 1) .- 1
    @test same_value(r3, 3:5)
    check_indexed_by(r3, 0:2)

    @testset "Idempotent indexing" begin
        @testset "Indexing into an IdOffsetRange" begin
            r = OffsetArrays.IdOffsetRange(3:5, -1)
            # Indexing with IdentityUnitRange
            s = IdentityUnitRange(0:2)
            @test axes(r[s]) == axes(s)
            for i in eachindex(s)
                @test r[s[i]] == r[s][i]
            end

            # Indexing with IdOffsetRange
            s = OffsetArrays.IdOffsetRange(-4:-2, 4)
            @test axes(r[s]) == axes(s)
            for i in eachindex(s)
                @test r[s[i]] == r[s][i]
            end

            # Indexing with UnitRange
            s = 0:2
            @test axes(r[s]) == axes(s)
            for i in eachindex(s)
                @test r[s[i]] == r[s][i]
            end
        end
        @testset "Indexing using an IdOffsetRange" begin
            r = OffsetArrays.IdOffsetRange(3:5, -1)
            # Indexing into an IdentityUnitRange
            s = IdentityUnitRange(-1:5)
            @test axes(s[r]) == axes(r)
            for i in eachindex(r)
                @test s[r[i]] == s[r][i]
            end

            # Indexing into an UnitRange
            s = -3:6
            @test axes(s[r]) == axes(r)
            for i in eachindex(r)
                @test s[r[i]] == s[r][i]
            end
        end
    end

    # Test reduced index
    rred = Base.reduced_index(r)
    @test typeof(rred) == typeof(r)
    @test length(rred) == 1
    @test first(rred) == first(r)
end

# used in testing the constructor
struct WeirdInteger{T} <: Integer
    x :: T
end
# assume that it doesn't behave as expected
Base.convert(::Type{Int}, a::WeirdInteger) = a

@testset "Constructors" begin
    @testset "Single-entry arrays in dims 0:5" begin
        for n = 0:5
            for z in (OffsetArray(ones(Int,ntuple(d->1,n)), ntuple(x->x-1,n)),
                    fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)), 1),
                    fill!(OffsetArray{Float64}(undef, ntuple(x->x:x, n)...), 1),
                    fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)), 1),
                    fill!(OffsetArray{Float64,n}(undef, ntuple(x->x:x, n)...), 1))
                @test length(LinearIndices(z)) == 1
                @test axes(z) == ntuple(x->x:x, n)
                @test z[1] == 1
            end
        end
        a0 = reshape([3])
        a = OffsetArray(a0)
        @test axes(a) == ()
        @test ndims(a) == 0
        @test a[] == 3
        @test a === OffsetArray(a, ())
        @test_throws ArgumentError OffsetArray(a, 0)
        @test_throws ArgumentError OffsetArray(a0, 0)
    end

    @testset "OffsetVector" begin
        # initialization
        one_based_axes = Any[
            (Base.OneTo(4), ),
            (1:4, ),
            (big(1):big(4), ),
            (CartesianIndex(1):CartesianIndex(4), ),
            (IdentityUnitRange(1:4), ),
            (IdOffsetRange(1:4),),
            (IdOffsetRange(3:6, -2),)
        ]

        offset_axes = Any[
            (-1:2, ),
            (big(-1):big(2), ),
            (CartesianIndex(-1):CartesianIndex(2), ),
            (IdentityUnitRange(-1:2), ),
            (IdOffsetRange(-1:2),),
            (IdOffsetRange(3:6, -4),)
        ]

        offsets = size.(one_based_axes[1], 1)
        offsets_big = map(big, offsets)

        for inds in Any[offsets, offsets_big, one_based_axes...]
            # test indices API
            a = OffsetVector{Float64}(undef, inds)
            @test eltype(a) === Float64
            @test axes(a) === axes(OffsetVector{Float64}(undef, inds...)) === axes(OffsetArray{Float64, 1}(undef, inds)) === axes(OffsetArray{Float64}(undef, inds))
            @test axes(a) === (IdOffsetRange(Base.OneTo(4), 0), )
            @test a.offsets === (0, )
            @test axes(a.parent) == (Base.OneTo(4), )

            a = OffsetVector{Nothing}(nothing, inds)
            @test eltype(a) === Nothing
            @test axes(a) === axes(OffsetVector{Nothing}(nothing, inds...)) === axes(OffsetArray{Nothing, 1}(nothing, inds))
            @test axes(a) === (IdOffsetRange(Base.OneTo(4), 0), )

            a = OffsetVector{Missing}(missing, inds)
            @test eltype(a) === Missing
            @test axes(a) === axes(OffsetVector{Missing}(missing, inds...)) === axes(OffsetArray{Missing, 1}(missing, inds))
            @test axes(a) === (IdOffsetRange(Base.OneTo(4), 0), )
        end

        # offset indexing
        for inds in offset_axes
            # test offsets
            a = OffsetVector{Float64}(undef, inds)
            ax = (IdOffsetRange(Base.OneTo(4), -2), )
            @test a.offsets === (-2, )
            @test axes(a.parent) == (Base.OneTo(4), )
            @test axes(a) === ax
            a = OffsetVector{Nothing}(nothing, inds)
            @test axes(a) === ax
            a = OffsetVector{Missing}(missing, inds)
            @test axes(a) === ax

            for (T, t) in Any[(Nothing, nothing), (Missing, missing)]
                a = OffsetVector{Union{T, Vector{Int}}}(undef, inds)
                @test !isassigned(a, -1)
                @test eltype(a) === Union{T, Vector{Int}}
                @test axes(a) === ax

                a = OffsetVector{Union{T, Vector{Int}}}(t, inds)
                @test a[-1] === t
            end
        end
        @test_throws Union{ArgumentError, ErrorException} OffsetVector{Float64}(undef, -2) # only positive number works

        # convenient constructors
        a = rand(4)
        for inds in offset_axes
            oa1 = OffsetVector(a, inds...)
            oa2 = OffsetVector(a, inds)
            oa3 = OffsetArray(a, inds...)
            oa4 = OffsetArray(a, inds)
            @test oa1 === oa2 === oa3 === oa4
            @test axes(oa1) === (IdOffsetRange(Base.OneTo(4), -2), )
            @test parent(oa1) === a
            @test oa1.offsets === (-2, )
        end

        oa = OffsetArray(a, :)
        @test oa === OffsetArray(a, (:, )) === OffsetArray(a, axes(a)) === OffsetVector(a, :) === OffsetVector(a, axes(a))
        @test oa == a
        @test axes(oa) == axes(a)
        @test axes(oa) !== axes(a)

        # nested offset array
        a = rand(4)
        oa = OffsetArray(a, -1)
        for inds in Any[.-oa.offsets, one_based_axes...]
            ooa = OffsetArray(oa, inds)
            @test typeof(parent(ooa)) <: Vector
            @test ooa === OffsetArray(oa, inds...) === OffsetVector(oa, inds) === OffsetVector(oa, inds...)
            @test ooa == a
            @test axes(ooa) == axes(a)
            @test axes(ooa) !== axes(a)
        end

        # overflow bounds check
        v = rand(5)
        @test axes(OffsetVector(v, typemax(Int)-length(v))) == (IdOffsetRange(axes(v)[1], typemax(Int)-length(v)), )
        @test_throws OverflowError OffsetVector(v, typemax(Int)-length(v)+1)
        ao = OffsetArray(v, typemin(Int))
        @test_nowarn OffsetArray{Float64, 1, typeof(ao)}(ao, (-1, ))
        @test_throws OverflowError OffsetArray{Float64, 1, typeof(ao)}(ao, (-2, )) # inner Constructor
        @test_throws OverflowError OffsetArray(ao, (-2, )) # convinient constructor accumulate offsets

        # disallow OffsetVector(::Array{<:Any, N}, offsets) where N != 1
        @test_throws ArgumentError OffsetVector(zeros(2,2), (2, 2))
        @test_throws ArgumentError OffsetVector(zeros(2,2), 2, 2)
        @test_throws ArgumentError OffsetVector(zeros(2,2), (1:2, 1:2))
        @test_throws ArgumentError OffsetVector(zeros(2,2), 1:2, 1:2)
        @test_throws ArgumentError OffsetVector(zeros(), ())
        @test_throws ArgumentError OffsetVector(zeros())
        @test_throws ArgumentError OffsetVector(zeros(2,2), ())
        @test_throws ArgumentError OffsetVector(zeros(2,2))
        @test_throws ArgumentError OffsetVector(zeros(2,2), 2)
        @test_throws ArgumentError OffsetVector(zeros(2,2), (2,))
        @test_throws ArgumentError OffsetVector(zeros(2:3,2:3), 2, 3)
        @test_throws ArgumentError OffsetVector(zeros(2:3,2:3), (2, 4))
        @test_throws ArgumentError OffsetVector(zeros(2:3,2:3), ())
        @test_throws ArgumentError OffsetVector(zeros(2:3,2:3))

        # eltype of an OffsetArray should match that of the parent (issue #162)
        @test_throws TypeError OffsetVector{Float64,Vector{ComplexF64}}
        # ndim of an OffsetArray should match that of the parent
        @test_throws TypeError OffsetVector{Float64,Matrix{Float64}}
    end

    @testset "OffsetMatrix" begin
        # initialization

        one_based_axes = Any[
                (Base.OneTo(4), Base.OneTo(3)),
                (1:4, 1:3),
                (big(1):big(4), big(1):big(3)),
                (CartesianIndex(1, 1):CartesianIndex(4, 3), ),
                (CartesianIndex(1):CartesianIndex(4), CartesianIndex(1):CartesianIndex(3)),
                (CartesianIndex(1):CartesianIndex(4), 1:3),
                (IdentityUnitRange(1:4), IdentityUnitRange(1:3)),
                (IdOffsetRange(1:4), IdOffsetRange(1:3)),
                (IdOffsetRange(3:6, -2), IdOffsetRange(3:5, -2)),
                (IdOffsetRange(3:6, -2), IdentityUnitRange(1:3)),
                (IdOffsetRange(3:6, -2), 1:3),
        ]

        offset_axes = Any[
                (-1:2, 0:2),
                (big(-1):big(2), big(0):big(2)),
                (CartesianIndex(-1, 0):CartesianIndex(2, 2), ),
                (-1:2, CartesianIndex(0):CartesianIndex(2)),
                (CartesianIndex(-1):CartesianIndex(2), CartesianIndex(0):CartesianIndex(2)),
                (CartesianIndex(-1):CartesianIndex(2), 0:2),
                (IdentityUnitRange(-1:2), 0:2),
                (IdOffsetRange(-1:2), IdOffsetRange(0:2)),
                (IdOffsetRange(3:6, -4), IdOffsetRange(2:4, -2)),
                (IdOffsetRange(3:6, -4), IdentityUnitRange(0:2)),
                (IdOffsetRange(-1:2), 0:2),
        ]

        offsets = size.(one_based_axes[1], 1)
        offsets_big = map(big, offsets)

        for inds in Any[offsets, offsets_big, one_based_axes...]
            # test API
            a = OffsetMatrix{Float64}(undef, inds)
            ax = (IdOffsetRange(Base.OneTo(4), 0), IdOffsetRange(Base.OneTo(3), 0))
            @test eltype(a) === Float64
            @test axes(a) === axes(OffsetMatrix{Float64}(undef, inds...)) === axes(OffsetArray{Float64, 2}(undef, inds)) === axes(OffsetArray{Float64, 2}(undef, inds...)) === axes(OffsetArray{Float64}(undef, inds))
            @test axes(a) === ax
            @test a.offsets === (0, 0)
            @test axes(a.parent) == (Base.OneTo(4), Base.OneTo(3))

            a = OffsetMatrix{Nothing}(nothing, inds)
            @test eltype(a) === Nothing
            @test axes(a) === axes(OffsetMatrix{Nothing}(nothing, inds...)) === axes(OffsetArray{Nothing, 2}(nothing, inds)) === axes(OffsetArray{Nothing, 2}(nothing, inds...))
            @test axes(a) === ax

            a = OffsetMatrix{Missing}(missing, inds)
            @test eltype(a) === Missing
            @test axes(a) === axes(OffsetMatrix{Missing}(missing, inds...)) === axes(OffsetArray{Missing, 2}(missing, inds)) === axes(OffsetArray{Missing, 2}(missing, inds...))
            @test axes(a) === ax
        end
        @test_throws Union{ArgumentError, ErrorException} OffsetMatrix{Float64}(undef, 2, -2) # only positive numbers works

        for inds in offset_axes
            # test offsets
            a = OffsetMatrix{Float64}(undef, inds)
            ax = (IdOffsetRange(Base.OneTo(4), -2), IdOffsetRange(Base.OneTo(3), -1))
            @test a.offsets === (-2, -1)
            @test axes(a.parent) == (Base.OneTo(4), Base.OneTo(3))
            @test axes(a) === ax
            a = OffsetMatrix{Nothing}(nothing, inds)
            @test axes(a) === ax
            a = OffsetMatrix{Missing}(missing, inds)
            @test axes(a) === ax

            for (T, t) in Any[(Nothing, nothing), (Missing, missing)]
                a = OffsetMatrix{Union{T, Vector{Int}}}(undef, inds)
                @test !isassigned(a, -1, 0)
                @test eltype(a) === Union{T, Vector{Int}}
                @test axes(a) === ax

                a = OffsetMatrix{Union{T, Vector{Int}}}(t, inds)
                @test a[-1, 0] === t
            end
        end

        # convenient constructors
        a = rand(4, 3)
        for inds in offset_axes
            ax = (IdOffsetRange(Base.OneTo(4), -2), IdOffsetRange(Base.OneTo(3), -1))
            oa1 = OffsetMatrix(a, inds...)
            oa2 = OffsetMatrix(a, inds)
            oa3 = OffsetArray(a, inds...)
            oa4 = OffsetArray(a, inds)
            @test oa1 === oa2 === oa3 === oa4
            @test axes(oa1) === ax
            @test parent(oa1) === a
            @test oa1.offsets === (-2, -1)
        end
        oa = OffsetArray(a, :, axes(a, 2))
        @test oa === OffsetArray(a, (axes(oa, 1), :)) === OffsetArray(a, axes(a)) === OffsetMatrix(a, (axes(oa, 1), :)) === OffsetMatrix(a, axes(a))
        @test oa == a
        @test axes(oa) == axes(a)
        @test axes(oa) !== axes(a)

        oa = OffsetMatrix(a, :, 2:4)
        @test oa === OffsetMatrix(a, axes(a, 1), 2:4) === OffsetMatrix(a, (axes(oa, 1), 2:4))

        # nested offset array
        a = rand(4, 3)
        oa = OffsetArray(a, -1, -2)
        for inds in Any[.-oa.offsets, one_based_axes...]
            ooa = OffsetArray(oa, inds)
            @test ooa === OffsetArray(oa, inds...) === OffsetMatrix(oa, inds) === OffsetMatrix(oa, inds...)
            @test typeof(parent(ooa)) <: Matrix
            @test ooa == a
            @test axes(ooa) == axes(a)
            @test axes(ooa) !== axes(a)
        end

        # overflow bounds check
        a = rand(4, 3)
        @test axes(OffsetMatrix(a, typemax(Int)-size(a, 1), 0)) == (IdOffsetRange(axes(a)[1], typemax(Int)-size(a, 1)), axes(a, 2))
        @test_throws OverflowError OffsetMatrix(a, typemax(Int)-size(a,1)+1, 0)
        @test_throws OverflowError OffsetMatrix(a, 0, typemax(Int)-size(a, 2)+1)

        # disallow OffsetMatrix(::Array{<:Any, N}, offsets) where N != 2
        @test_throws ArgumentError OffsetMatrix(zeros(2), (2,))
        @test_throws ArgumentError OffsetMatrix(zeros(2), 2)
        @test_throws ArgumentError OffsetMatrix(zeros(2), (1:2,))
        @test_throws ArgumentError OffsetMatrix(zeros(2), 1:2)
        @test_throws ArgumentError OffsetMatrix(zeros(), ())
        @test_throws ArgumentError OffsetMatrix(zeros())
        @test_throws ArgumentError OffsetMatrix(zeros(2), ())
        @test_throws ArgumentError OffsetMatrix(zeros(2))
        @test_throws ArgumentError OffsetMatrix(zeros(2), (1, 2))
        @test_throws ArgumentError OffsetMatrix(zeros(2), 1, 2)
        @test_throws ArgumentError OffsetMatrix(zeros(2:3), (2,))
        @test_throws ArgumentError OffsetMatrix(zeros(2:3), 2)
        @test_throws ArgumentError OffsetMatrix(zeros(2:3, 1:2, 1:2), (2,0,0))
        @test_throws ArgumentError OffsetMatrix(zeros(2:3, 1:2, 1:2), 2,0,0)
        @test_throws ArgumentError OffsetMatrix(zeros(2:3, 1:2, 1:2), ())
        @test_throws ArgumentError OffsetMatrix(zeros(2:3, 1:2, 1:2))

        # eltype of an OffsetArray should match that of the parent (issue #162)
        @test_throws TypeError OffsetMatrix{Float64,Matrix{ComplexF64}}
        # ndim of an OffsetArray should match that of the parent
        @test_throws TypeError OffsetMatrix{Float64,Vector{Float64}}
    end

    # no need to duplicate the 2D case here,
    # only add some special test cases
    @testset "OffsetArray" begin
        a = rand(2, 2, 2)
        oa = OffsetArray(a, 0:1, 3:4, 2:3)
        @test OffsetArray(a, CartesianIndices(axes(oa))) == oa
        @test axes(OffsetArray(a, :, CartesianIndices((3:4, 2:3)))) == (1:2, 3:4, 2:3)
        @test axes(OffsetArray(a, 10:11, CartesianIndices((3:4, 2:3)) )) == (10:11, 3:4, 2:3)
        @test axes(OffsetArray(a, CartesianIndices((3:4, 2:3)), :)) == (3:4, 2:3, 1:2)
        @test axes(OffsetArray(a, CartesianIndices((3:4, 2:3)), 10:11)) == (3:4, 2:3, 10:11)
        @test axes(OffsetArray(a, :, :, CartesianIndices((3:4,)) )) == (1:2, 1:2, 3:4)
        @test axes(OffsetArray(a, 10:11, :, CartesianIndices((3:4,)) )) == (10:11, 1:2, 3:4)
        @test axes(OffsetArray(a, 10:11, 2:3, CartesianIndices((3:4,)) )) == (10:11, 2:3, 3:4)

        # ignore empty CartesianIndices
        @test OffsetArray(a, CartesianIndices(()), 0:1, :, 2:3) == OffsetArray(a, 0:1, :, 2:3)
        @test OffsetArray(a, 0:1, CartesianIndices(()), :, 2:3) == OffsetArray(a, 0:1, :, 2:3)
        @test OffsetArray(a, 0:1, :,  CartesianIndices(()), 2:3) == OffsetArray(a, 0:1, :, 2:3)
        @test OffsetArray(a, 0:1, :, 2:3, CartesianIndices(())) == OffsetArray(a, 0:1, :, 2:3)

        indices = (-1:1, -7:7, -1:2, -5:5, -1:1, -3:3, -2:2, -1:1)
        y = OffsetArray{Float64}(undef, indices...);
        @test axes(y) === axes(OffsetArray{Float64}(undef, indices))
        @test axes(y) === axes(OffsetArray{Float64, length(indices)}(undef, indices...))
        @test axes(y) == (-1:1, -7:7, -1:2, -5:5, -1:1, -3:3, -2:2, -1:1)
        @test eltype(y) === Float64

        @test_throws ArgumentError OffsetArray{Float64, 2}(undef, indices)
        @test_throws ArgumentError OffsetArray(y, indices[1:2])

        @test ndims(OffsetArray(zeros(), ())) == 0
        @test Base.axes1(OffsetArray(zeros(), ())) === OffsetArrays.IdOffsetRange(Base.OneTo(1))

        @testset "convenience constructors" begin
            ax = (2:3, 4:5)

            for f in (zeros, ones)
                a = f(Float64, ax)
                @test axes(a) == ax
                @test eltype(a) == Float64
            end

            for f in (trues, falses)
                a = f(ax)
                @test axes(a) == ax
                @test eltype(a) == Bool
            end
        end

        # eltype of an OffsetArray should match that of the parent (issue #162)
        @test_throws TypeError OffsetArray{Float64,2,Matrix{ComplexF64}}
        # ndim of an OffsetArray should match that of the parent
        @test_throws TypeError OffsetArray{Float64,3,Matrix{Float64}}

        # should throw a TypeError if the offsets can not be converted to Ints
        @test_throws TypeError OffsetVector{Int,Vector{Int}}(zeros(Int,2), (WeirdInteger(1),))
    end

    @testset "custom range types" begin
        @testset "EllipsisNotation" begin
            @testset "Vector" begin
                v = rand(5)
                @test axes(OffsetArray(v, ..)) == axes(v)
                @test OffsetArray(v, ..) == OffsetArray(v, :)
                @test axes(OffsetVector(v, ..)) == axes(v)
                @test OffsetVector(v, ..) == OffsetVector(v, :)

                @test axes(OffsetArray(v, .., 2:6)) == (2:6, )
                @test OffsetArray(v, .., 2:6) == OffsetArray(v, 2:6)
                @test axes(OffsetVector(v, .., 2:6)) == (2:6, )
                @test OffsetVector(v, .., 2:6) == OffsetVector(v, 2:6)
            end
            @testset "Matrix" begin
                m = rand(2, 2)
                @test axes(OffsetArray(m, ..)) == axes(m)
                @test OffsetArray(m, ..) == OffsetArray(m, :, :)
                @test axes(OffsetMatrix(m, ..)) == axes(m)
                @test OffsetMatrix(m, ..) == OffsetMatrix(m, :, :)

                @test axes(OffsetArray(m, .., 2:3)) == (axes(m, 1), 2:3)
                @test OffsetArray(m, .., 2:3) == OffsetArray(m, :, 2:3)
                @test axes(OffsetMatrix(m, .., 2:3)) == (axes(m, 1), 2:3)
                @test OffsetMatrix(m, .., 2:3) == OffsetMatrix(m, :, 2:3)

                @test axes(OffsetArray(m, .., 2:3, 3:4)) == (2:3, 3:4)
                @test OffsetArray(m, .., 2:3, 3:4) == OffsetArray(m, 2:3, 3:4)
                @test axes(OffsetMatrix(m, .., 2:3, 3:4)) == (2:3, 3:4)
                @test OffsetMatrix(m, .., 2:3, 3:4) == OffsetMatrix(m, 2:3, 3:4)
            end
            @testset "3D Array" begin
                a = rand(2, 2, 2)
                @test axes(OffsetArray(a, ..)) == axes(a)
                @test OffsetArray(a, ..) == OffsetArray(a, :, :, :)

                @test axes(OffsetArray(a, .., 2:3)) == (axes(a)[1:2]..., 2:3)
                @test OffsetArray(a, .., 2:3) == OffsetArray(a, :, :, 2:3)

                @test axes(OffsetArray(a, .., 2:3, 3:4)) == (axes(a, 1), 2:3, 3:4)
                @test OffsetArray(a, .., 2:3, 3:4) == OffsetArray(a, :, 2:3, 3:4)

                @test axes(OffsetArray(a, 2:3, .., 3:4)) == (2:3, axes(a, 2), 3:4)
                @test OffsetArray(a, 2:3, .., 3:4) == OffsetArray(a, 2:3, :, 3:4)

                @test axes(OffsetArray(a, .., 4:5, 2:3, 3:4)) == (4:5, 2:3, 3:4)
                @test OffsetArray(a, .., 4:5, 2:3, 3:4) == OffsetArray(a, 4:5, 2:3, 3:4)
            end
        end
        @testset "ZeroBasedIndexing" begin
            Base.to_indices(A, inds, ::Tuple{ZeroBasedIndexing}) = map(x -> 0:length(x) - 1, inds)

            a = zeros(3,3)
            oa = OffsetArray(a, ZeroBasedIndexing())
            @test axes(oa) == (0:2, 0:2)
        end
        @testset "TupleOfRanges" begin
            Base.to_indices(A, inds, t::Tuple{TupleOfRanges{N}}) where {N} = t
            OffsetArrays.AxisConversionStyle(::Type{TupleOfRanges{N}}) where {N} =
                OffsetArrays.TupleOfRanges()

            Base.convert(::Type{Tuple{Vararg{AbstractUnitRange{Int}}}}, t::TupleOfRanges) = t.x

            a = zeros(3,3)
            inds = TupleOfRanges((3:5, 2:4))
            oa = OffsetArray(a, inds)
            @test axes(oa) == inds.x
        end
        @testset "NewColon" begin
            Base.to_indices(A, inds, t::Tuple{NewColon,Vararg{Any}}) =
                (_uncolon(inds, t), to_indices(A, Base.tail(inds), Base.tail(t))...)

            _uncolon(inds::Tuple{}, I::Tuple{NewColon, Vararg{Any}}) = OneTo(1)
            _uncolon(inds::Tuple, I::Tuple{NewColon, Vararg{Any}}) = inds[1]

            a = zeros(3, 3)
            oa = OffsetArray(a, (NewColon(), 2:4))
            @test axes(oa) == (axes(a,1), 2:4)
        end
    end

    @testset "Offset range construction" begin
        r = -2:5
        for AT in Any[OffsetArray, OffsetVector]
            y = AT(r, r)
            @test axes(y) == (r,)
            @test step(y) == step(r)
            y = AT(r, (r,))
            @test axes(y) == (r,)
            y = AT(r, CartesianIndices((r, )))
            @test axes(y) == (r, )
        end
    end
end

@testset "Axes supplied to constructor correspond to final result" begin
    # Ref https://github.com/JuliaArrays/OffsetArrays.jl/pull/65#issuecomment-457181268
    B = BidirectionalVector([1, 2, 3], -2)
    A = OffsetArray(B, -1:1)
    @test axes(A) == (IdentityUnitRange(-1:1),)
end

@testset "Traits" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))                   # IndexLinear
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))   # IndexCartesian
    @test axes(A) === axes(S)
    @test axes(A) == axes(S) == (0:1, 3:4)
    @test axes(A, 1) === OffsetArrays.IdOffsetRange(Base.OneTo(2), -1)
    @test size(A) == size(A0)
    @test size(A, 1) == size(A0, 1)
    @test length(A) == length(A0)
    @test A == OffsetArray(A0, 0:1, 3:4)
    @test_throws DimensionMismatch OffsetArray(A0, 0:2, 3:4)
    @test_throws DimensionMismatch OffsetArray(A0, 0:1, 2:4)
    @test eachindex(IndexLinear(), A) == eachindex(IndexLinear(), parent(A))
    @test eachindex(IndexCartesian(), A) == CartesianIndices(A)
end

@testset "Scalar indexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    @test @inferred(A[0,3]) == @inferred(A[0,3,1]) == @inferred(A[1]) == @inferred(S[0,3]) == @inferred(S[0,3,1]) == @inferred(S[1]) == 1
    @test A[1,3] == A[1,3,1] == A[2] == S[1,3] == S[1,3,1] == S[2] == 2
    @test A[0,4] == A[0,4,1] == A[3] == S[0,4] == S[0,4,1] == S[3] == 3
    @test A[1,4] == A[1,4,1] == A[4] == S[1,4] == S[1,4,1] == S[4] == 4
    @test @inbounds(A[0,3]) == @inbounds(A[0,3,1]) == @inbounds(A[1]) == @inbounds(S[0,3]) == @inbounds(S[0,3,1]) == @inbounds(S[1]) == 1
    @test @inbounds(A[1,3]) == @inbounds(A[1,3,1]) == @inbounds(A[2]) == @inbounds(S[1,3]) == @inbounds(S[1,3,1]) == @inbounds(S[2]) == 2
    @test @inbounds(A[0,4]) == @inbounds(A[0,4,1]) == @inbounds(A[3]) == @inbounds(S[0,4]) == @inbounds(S[0,4,1]) == @inbounds(S[3]) == 3
    @test @inbounds(A[1,4]) == @inbounds(A[1,4,1]) == @inbounds(A[4]) == @inbounds(S[1,4]) == @inbounds(S[1,4,1]) == @inbounds(S[4]) == 4
    @test_throws BoundsError(A, (1,1)) A[1,1]
    @test_throws BoundsError(A, (1,1)) A[1,1] = 4
    @test_throws BoundsError(S, (1,1)) S[1,1]
    @test_throws BoundsError(S, (1,1)) S[1,1] = 4
    @test_throws BoundsError(A, (0,3,2)) A[0,3,2]
    @test_throws BoundsError(A, (0,3,2)) A[0,3,2] = 4
    @test_throws BoundsError(A, (0,3,0)) A[0,3,0]
    @test_throws BoundsError(A, (0,3,0)) A[0,3,0] = 4
    Ac = copy(A)
    Ac[0,3] = 10
    @test Ac[0,3] == 10
    Ac[0,3,1] = 11
    @test Ac[0,3] == 11
    @inbounds Ac[0,3,1] = 12
    @test Ac[0,3] == 12

    y = OffsetArray{Float64}(undef, -1:1, -7:7, -3:-1, -5:5, -1:1, -3:3, -2:2, -1:1)
    y[-1,-7,-3,-5,-1,-3,-2,-1] = 14
    y[-1,-7,-3,-5,-1,-3,-2,-1] += 5
    @test y[-1,-7,-3,-5,-1,-3,-2,-1] == 19

    @testset "setindex!" begin
        A = OffsetArray(ones(2,2), 1:2, 1:2)
        @test setindex!(A, 2, 1, 1) === A
        @test A[1,1] == 2
        @test setindex!(A, 2, 1) === A
        @test A[1] == 2

        v = OffsetArray(ones(3), 4:6)
        @test setindex!(A, 2, 4) === A
        @test A[4] == 2
    end
end

@testset "Vector indexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    @test A[:, 3] == S[:, 3] == OffsetArray([1,2], (A.offsets[1],))
    @test A[:, 4] == S[:, 4] == OffsetArray([3,4], (A.offsets[1],))
    @test_throws BoundsError A[:, 1]
    @test_throws BoundsError S[:, 1]
    @test A[0, :] == S[0, :] == OffsetArray([1,3], (A.offsets[2],))
    @test A[1, :] == S[1, :] == OffsetArray([2,4], (A.offsets[2],))
    @test_throws BoundsError A[2, :]
    @test_throws BoundsError S[2, :]
    @test A[0:1, 3] == S[0:1, 3] == [1,2]
    @test A[[1,0], 3] == S[[1,0], 3] == [2,1]
    @test A[0, 3:4] == S[0, 3:4] == [1,3]
    @test A[1, [4,3]] == S[1, [4,3]] == [4,2]
    @test A[:, :] == S[:, :] == A
end

@testset "Vector indexing with offset ranges" begin
    r = OffsetArray(8:10, -1:1)
    r1 = r[0:1]
    @test r1 === 9:10
    r1 = (8:10)[OffsetArray(1:2, -5:-4)]
    @test axes(r1) == (IdentityUnitRange(-5:-4),)
    @test parent(r1) === 8:9
    r1 = OffsetArray(8:10, -1:1)[OffsetArray(0:1, -5:-4)]
    @test axes(r1) == (IdentityUnitRange(-5:-4),)
    @test parent(r1) === 9:10

    a = OffsetVector(3:4, 10:11)
    ax = OffsetArrays.IdOffsetRange(5:6, 5)
    @test axes(a[ax]) == axes(ax)
    for i in axes(ax,1)
        @test a[ax[i]] == a[ax][i]
    end

    ax = IdentityUnitRange(10:11)
    @test axes(a[ax]) == axes(ax)
    for i in axes(ax,1)
        @test a[ax[i]] == a[ax][i]
    end
end

@testset "CartesianIndexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    @test A[CartesianIndex((0,3))] == S[CartesianIndex((0,3))] == 1
    @test A[CartesianIndex((0,3)),1] == S[CartesianIndex((0,3)),1] == 1
    @test @inbounds(A[CartesianIndex((0,3))]) == @inbounds(S[CartesianIndex((0,3))]) == 1
    @test @inbounds(A[CartesianIndex((0,3)),1]) == @inbounds(S[CartesianIndex((0,3)),1]) == 1
    @test_throws BoundsError A[CartesianIndex(1,1)]
    @test_throws BoundsError A[CartesianIndex(1,1),0]
    @test_throws BoundsError A[CartesianIndex(1,1),2]
    @test_throws BoundsError S[CartesianIndex(1,1)]
    @test_throws BoundsError S[CartesianIndex(1,1),0]
    @test_throws BoundsError S[CartesianIndex(1,1),2]
    @test eachindex(A) == 1:4
    @test eachindex(S) == CartesianIndices(IdentityUnitRange.((0:1,3:4)))
end

@testset "IdentityUnitRange indexing" begin
    # 155
    a = OffsetVector(3:4, 2:3)
    ax = IdentityUnitRange(2:3)
    @test a[ax[2]] == a[ax][2]

    s = -2:2:4
    r = 5:8
    y = OffsetArray(s, r)
    @test axes(y) == (r,)
    @test step(y) == step(s)

    a = OffsetVector(3:4, 10:11)
    ax = OffsetArrays.IdOffsetRange(5:6, 5)
    @test axes(a[ax]) == axes(ax)
    for i in axes(ax,1)
        @test a[ax[i]] == a[ax][i]
    end

    ax = IdentityUnitRange(10:11)
    @test axes(a[ax]) == axes(ax)
    for i in axes(ax,1)
        @test a[ax[i]] == a[ax][i]
    end
end

@testset "view" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    S = view(A, :, 3)
    @test S == OffsetArray([1,2], (A.offsets[1],))
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[2]
    @test axes(S) == (IdentityUnitRange(0:1),)
    S = view(A, 0, :)
    @test S == OffsetArray([1,3], (A.offsets[2],))
    @test S[3] == 1
    @test S[4] == 3
    @test_throws BoundsError S[1]
    @test axes(S) == (IdentityUnitRange(3:4),)
    S = view(A, 0:0, 4)
    @test S == [3]
    @test S[1] == 3
    @test_throws BoundsError S[0]
    @test axes(S) === (Base.OneTo(1),)
    S = view(A, 1, 3:4)
    @test S == [2,4]
    @test S[1] == 2
    @test S[2] == 4
    @test_throws BoundsError S[3]
    @test axes(S) === (Base.OneTo(2),)
    S = view(A, :, :)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) == IdentityUnitRange.((0:1, 3:4))
    S = view(A, axes(A)...)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) == OffsetArrays.IdOffsetRange.((0:1, 3:4))
    # issue 100
    S = view(A, axes(A, 1), 3)
    @test S == A[:, 3]
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[length(S)]
    @test axes(S) == (OffsetArrays.IdOffsetRange(0:1), )
    # issue 100
    S = view(A, 1, axes(A, 2))
    @test S == A[1, :]
    @test S[3] == 2
    @test S[4] == 4
    @test_throws BoundsError S[1]
    @test axes(S) == (OffsetArrays.IdOffsetRange(3:4), )

    # issue 133
    r = OffsetArrays.IdOffsetRange(1:2, -1)
    v1 = view(A, r, 3)
    @test v1[0] == 1
    @test v1[1] == 2
    @test axes(v1, 1) == axes(r, 1)
    v2 = view(A, UnitRange(r), 3)
    for (indflat, indoffset) in enumerate(r)
        @test v1[indoffset] == v2[indflat]
    end

    # issue 133
    r = OffsetArrays.IdOffsetRange(1:2, 2)
    v1 = view(A, 1, r)
    @test v1[3] == 2
    @test v1[4] == 4
    @test axes(v1, 1) == axes(r, 1)
    v2 = view(A, 1, UnitRange(r))
    for (indflat, indoffset) in enumerate(r)
        @test v1[indoffset] == v2[indflat]
    end

    # issue 133
    a12 = zeros(3:8, 3:4)
    r = OffsetArrays.IdOffsetRange(Base.OneTo(3), 5)
    a12[r, 4] .= 3
    @test all(a12[r, 4] .== 3)
    @test all(a12[UnitRange(r), 4] .== 3)

    A0 = collect(reshape(1:24, 2, 3, 4))
    A = OffsetArray(A0, (-1,2,1))
    S = view(A, axes(A, 1), 3:4, axes(A, 3))
    @test S == A[:, 3:4, :]
    @test S[0, 1, 2] == A[0, 3, 2]
    @test S[0, 2, 2] == A[0, 4, 2]
    @test S[1, 1, 2] == A[1, 3, 2]
    @test axes(S) == (OffsetArrays.IdOffsetRange(0:1), Base.OneTo(2), OffsetArrays.IdOffsetRange(2:5))

    # issue #186
    a = reshape(1:12, 3, 4)
    r = OffsetArrays.IdOffsetRange(3:4)
    av = view(a, :, r)
    @test av == a[:, 3:4]
    @test axes(av) == (axes(a,1), axes(r,1))
    r = OffsetArrays.IdOffsetRange(1:2,2)
    av = view(a, :, r)
    @test no_offset_view(av) == a[:, 3:4]
    @test axes(av) == (axes(a,1), axes(r,1))
    r = OffsetArrays.IdOffsetRange(2:3)
    av1d = view(a, r, 3)
    @test av1d == a[2:3, 3]
    @test axes(av1d) == (axes(r,1),)
    r = OffsetArrays.IdOffsetRange(Base.OneTo(2), 1)
    av1d = view(a, r, 3)
    @test no_offset_view(av1d) == a[2:3, 3]
    @test axes(av1d) == (axes(r,1),)

    # fix IdOffsetRange(::IdOffsetRange, offset) nesting from #178
    b = 1:20
    bov = OffsetArray(view(b, 3:4), 3:4)
    c = @view b[bov]
    @test same_value(c, 3:4)
    @test axes(c,1) == 3:4
    d = OffsetArray(c, 1:2)
    @test same_value(d, c)
    @test axes(d,1) == 1:2
end

@testset "iteration" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    let a
        for (a,d) in zip(A, A0)
            @test a == d
        end
    end
end

@testset "show/summary" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    @test sprint(show, A) == "[1 3; 2 4]"
    @test sprint(show, S) == "[1 3; 2 4]"
    strs = split(strip(sprint(show, MIME("text/plain"), A)), '\n')
    @test strs[2] == " 1  3"
    @test strs[3] == " 2  4"
    v = OffsetArray(rand(3), (-2,))
    @test sprint(show, v) == sprint(show, parent(v))
    io = IOBuffer()
    function cmp_showf(printfunc, io, A)
        ioc = IOContext(io, :limit=>true, :compact=>true)
        printfunc(ioc, A)
        str1 = String(take!(io))
        printfunc(ioc, parent(A))
        str2 = String(take!(io))
        @test str1 == str2
    end
    cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,5), (10,-9)))       # rows&cols fit
    cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,5), (10,-9)))    # columns fit
    cmp_showf(Base.print_matrix, io, OffsetArray(rand(5,10^3), (10,-9)))    # rows fit
    cmp_showf(Base.print_matrix, io, OffsetArray(rand(10^3,10^3), (10,-9))) # neither fits

    a = OffsetArray([1 2; 3 4], -1:0, 5:6)
    shownsz = VERSION >= v"1.2.0-DEV.229" ? Base.dims2string(size(a))*' ' : ""
    @test summary(a) == "$(shownsz)OffsetArray(::$(typeof(parent(a))), -1:0, 5:6) with eltype $(Int) with indices -1:0×5:6"
    shownsz = VERSION >= v"1.2.0-DEV.229" ? Base.dims2string(size(view(a, :, 5)))*' ' : ""
    @test summary(view(a, :, 5)) == "$(shownsz)view(OffsetArray(::$(typeof(parent(a))), -1:0, 5:6), :, 5) with eltype $(Int) with indices -1:0"
    a = OffsetArray(reshape([1]))
    @test summary(a) == "0-dimensional OffsetArray(::$(typeof(parent(a)))) with eltype $(Int)"

    a = OffsetArray([1 2; 3 4], -1:0, 5:6)
    io = IOBuffer()
    show(io, axes(a, 1))
    @test String(take!(io)) == "OffsetArrays.IdOffsetRange(-1:0)"
    show(io, axes(a, 2))
    @test String(take!(io)) == "OffsetArrays.IdOffsetRange(5:6)"

    @test Base.inds2string(axes(a)) == Base.inds2string(map(UnitRange, axes(a)))

    show(io, OffsetArray(3:5, 0:2))
    @test String(take!(io)) == "3:5 with indices 0:2"

    show(io, MIME"text/plain"(), OffsetArray(3:5, 0:2))
    @test String(take!(io)) == "3:5 with indices 0:2"

    d = Diagonal([1,2,3])
    Base.print_array(io, d)
    s1 = String(take!(io))
    od = OffsetArray(d, -1:1, 3:5)
    Base.print_array(io, od)
    s2 = String(take!(io))
    @test s1 == s2

    @test Base.replace_in_print_matrix(od, -1, 3, " ") == Base.replace_in_print_matrix(d, 1, 1, " ")
    @test Base.replace_in_print_matrix(od, -1, 4, " ") == Base.replace_in_print_matrix(d, 1, 2, " ")

    v = rand(3)
    ov = OffsetArray(v, (-2,))
    @test Base.replace_in_print_matrix(ov, -1, 1, " ") == Base.replace_in_print_matrix(v, 1, 1, " ")

    # Avoid returning the value of toplevel if it is false
    # showarg should only print values, it shouldn't return anything
    @test Base.showarg(io, a, false) === nothing
    # check the other case too for good measure
    @test Base.showarg(io, a, true) === nothing
end

@testset "readdlm/writedlm" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    io = IOBuffer()
    writedlm(io, A)
    seek(io, 0)
    @test readdlm(io, eltype(A)) == parent(A)
end

@testset "similar" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    B = similar(A, Float32)
    @test isa(B, OffsetArray{Float32,2})
    @test axes(B) === axes(A)
    B = similar(A, (3,4))
    @test isa(B, Array{Int,2})
    @test size(B) == (3,4)
    @test axes(B) === (Base.OneTo(3), Base.OneTo(4))
    B = similar(A, (-3:3,1:4))
    @test isa(B, OffsetArray{Int,2})
    @test axes(B) == IdentityUnitRange.((-3:3, 1:4))
    B = similar(parent(A), (-3:3,1:4))
    @test isa(B, OffsetArray{Int,2})
    @test axes(B) == IdentityUnitRange.((-3:3, 1:4))
    @test isa([x for x in [1,2,3]], Vector{Int})
    @test similar(Array{Int}, (0:0, 0:0)) isa OffsetArray{Int, 2}
    @test similar(Array{Int}, (1, 1)) isa Matrix{Int}
    @test similar(Array{Int}, (Base.OneTo(1), Base.OneTo(1))) isa Matrix{Int}
    B = similar(Array{Int}, (0:0, 3))
    @test isa(B, OffsetArray{Int, 2})
    @test axes(B) == (0:0, 1:3)

    @test_throws MethodError similar(A, (:,))
    @test_throws MethodError similar(A, (: ,:))
    @test_throws MethodError similar(A, (: ,2))
    @test_throws MethodError similar(A, Float64, (: ,:))
    @test_throws MethodError similar(A, Float64, (: ,2))

    function testsimilar(args...)
        try
           similar(args...)
        catch e
            @test e isa MethodError
            io = IOBuffer()
            showerror(io, e)
            s = split(String(take!(io)),'\n')[1]
            @test occursin(repr(similar), s)
        end
    end

    testsimilar(typeof(A), (:, :))
    testsimilar(typeof(A), (:, 2))
    testsimilar(typeof(A), (:, 1:3))
end

@testset "reshape" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    B = reshape(A0, -10:-9, 9:10)
    @test isa(B, OffsetArray{Int,2})
    @test parent(B) === A0
    @test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
    B = reshape(A, -10:-9, 9:10)
    @test isa(B, OffsetArray{Int,2})
    @test pointer(parent(B)) === pointer(A0)
    @test axes(B) == IdentityUnitRange.((-10:-9, 9:10))
    b = reshape(A, -7:-4)
    @test axes(b) == (IdentityUnitRange(-7:-4),)
    @test isa(parent(b), Vector{Int})
    @test pointer(parent(b)) === pointer(parent(A))
    @test parent(b) == A0[:]
    a = OffsetArray(rand(3,3,3), -1:1, 0:2, 3:5)
    # Offset axes are required for reshape(::OffsetArray, ::Val) support
    b = reshape(a, Val(2))
    @test isa(b, OffsetArray{Float64,2})
    @test pointer(parent(b)) === pointer(parent(a))
    @test axes(b) == IdentityUnitRange.((-1:1, 1:9))
    b = reshape(a, Val(4))
    @test isa(b, OffsetArray{Float64,4})
    @test pointer(parent(b)) === pointer(parent(a))
    @test axes(b) == (axes(a)..., IdentityUnitRange(1:1))

    @test reshape(OffsetArray(-1:0, -1:0), :, 1) == reshape(-1:0, 2, 1)
    @test reshape(OffsetArray(-1:2, -1:2), -2:-1, :) == reshape(-1:2, -2:-1, 2)

    @test reshape(OffsetArray(-1:0, -1:0), :) == OffsetArray(-1:0, -1:0)
    @test reshape(A, :) == reshape(A0, :)

    # pop the parent
    B = reshape(A, size(A))
    @test B == A0
    @test parent(B) === A0
    B = reshape(A, (Base.OneTo(2), 2))
    @test B == A0
    @test parent(B) === A0
    B = reshape(A, (2,:))
    @test B == A0
    @test parent(B) === A0

    # julialang/julia #33614
    A = OffsetArray(-1:0, (-2,))
    @test reshape(A, :) === A
    Arsc = reshape(A, :, 1)
    Arss = reshape(A, 2, 1)
    @test Arsc[1,1] == Arss[1,1] == -1
    @test Arsc[2,1] == Arss[2,1] == 0
    @test_throws BoundsError Arsc[0,1]
    @test_throws BoundsError Arss[0,1]
    A = OffsetArray([-1,0], (-2,))
    Arsc = reshape(A, :, 1)
    Arsc[1,1] = 5
    @test first(A) == 5
end

@testset "Indexing with OffsetArray axes" begin
    A0 = [1 3; 2 4]

    i1 = OffsetArray([2,1], (-5,))
    i1 = OffsetArray([2,1], -5)
    b = A0[i1, 1]
    @test axes(b) == (IdentityUnitRange(-4:-3),)
    @test b[-4] == 2
    @test b[-3] == 1
    b = A0[1,i1]
    @test axes(b) == (IdentityUnitRange(-4:-3),)
    @test b[-4] == 3
    @test b[-3] == 1
    v = view(A0, i1, 1)
    @test axes(v) == (IdentityUnitRange(-4:-3),)
    v = view(A0, 1:1, i1)
    @test axes(v) == (Base.OneTo(1), IdentityUnitRange(-4:-3))

    for r in (1:10, 1:1:10, StepRangeLen(1, 1, 10), LinRange(1, 10, 10), 0.1:0.2:0.9)
        for s in (IdentityUnitRange(2:3), OffsetArray(2:3, 2:3))
            @test axes(r[s]) == axes(s)
        end
    end
end

@testset "logical indexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

    @test A[A .> 2] == [3,4]
end

@testset "copyto!" begin
    a = OffsetArray{Int}(undef, (-3:-1,))
    fill!(a, -1)
    copyto!(a, (1,2))   # non-array iterables
    @test a[-3] == 1
    @test a[-2] == 2
    @test a[-1] == -1
    fill!(a, -1)
    copyto!(a, -2, (1,2))
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    @test_throws BoundsError copyto!(a, 1, (1,2))
    fill!(a, -1)
    copyto!(a, -2, (1,2,3), 2)
    @test a[-3] == -1
    @test a[-2] == 2
    @test a[-1] == 3
    @test_throws BoundsError copyto!(a, -2, (1,2,3), 1)
    fill!(a, -1)
    copyto!(a, -2, (1,2,3), 1, 2)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2

    b = 1:2    # copy between AbstractArrays
    bo = OffsetArray(1:2, (-3,))
    if VERSION < v"1.5-"
        @test_throws BoundsError copyto!(a, b)
        fill!(a, -1)
        copyto!(a, bo)
        @test a[-3] == -1
        @test a[-2] == 1
        @test a[-1] == 2
    else
        # the behavior of copyto! is corrected as the documentation says "first n element"
        # https://github.com/JuliaLang/julia/pull/34049
        fill!(a, -1)
        copyto!(a, bo)
        @test a[-3] == 1
        @test a[-2] == 2
        @test a[-1] == -1
    end
    fill!(a, -1)
    copyto!(a, -2, bo)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    @test_throws BoundsError copyto!(a, -4, bo)
    @test_throws BoundsError copyto!(a, -1, bo)
    fill!(a, -1)
    copyto!(a, -3, b, 2)
    @test a[-3] == 2
    @test a[-2] == a[-1] == -1
    @test_throws BoundsError copyto!(a, -3, b, 1, 4)
    am = OffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
    fill!(am, -1)
    copyto!(am, b)
    @test am[1] == 1
    @test am[2] == 2
    @test am[3] == -1
    @test am[1,7] == 1
    @test am[1,8] == 2
    @test am[1,9] == -1
end

@testset "map" begin
    am = OffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
    fill!(am, -1)
    copyto!(am, 1:2)

    dest = similar(am)
    map!(+, dest, am, am)
    @test dest[1,7] == 2
    @test dest[1,8] == 4
    @test dest[1,9] == -2
end

@testset "reductions" begin
    A = OffsetArray(rand(4,4), (-3,5))
    @test maximum(A) == maximum(parent(A))
    @test minimum(A) == minimum(parent(A))
    @test extrema(A) == extrema(parent(A))
    @test sum(A) == sum(parent(A))
    @test sum(A, dims=1) == OffsetArray(sum(parent(A), dims=1), A.offsets)
    @test sum(A, dims=2) == OffsetArray(sum(parent(A), dims=2), A.offsets)
    @test sum(A, dims=(1,2)) == OffsetArray(sum(parent(A), dims=(1,2)), A.offsets)
    @test sum(view(OffsetArray(reshape(1:27, 3, 3, 3), 0, 0, 0), :, :, 1:2), dims=(2,3)) == reshape([51,57,63], 3, 1, 1)
    C = similar(A)
    cumsum!(C, A, dims = 1)
    @test parent(C) == cumsum(parent(A), dims = 1)
    @test parent(cumsum(A, dims = 1)) == cumsum(parent(A), dims = 1)
    cumsum!(C, A, dims = 2)
    @test parent(C) == cumsum(parent(A), dims = 2)
    R = similar(A, (1:1, 6:9))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 1)
    R = similar(A, (-2:1, 1:1))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 2)
    amin, iamin = findmin(A)
    pmin, ipmin = findmin(parent(A))
    @test amin == pmin
    @test A[iamin] == amin
    @test amin == parent(A)[ipmin]
    amax, iamax = findmax(A)
    pmax, ipmax = findmax(parent(A))
    @test amax == pmax
    @test A[iamax] == amax
    @test amax == parent(A)[ipmax]

    amin, amax = extrema(parent(A))
    @test clamp.(A, (amax+amin)/2, amax) == OffsetArray(clamp.(parent(A), (amax+amin)/2, amax), axes(A))
end

# v  = OffsetArray([1,1e100,1,-1e100], (-3,))*1000
# v2 = OffsetArray([1,-1e100,1,1e100], (5,))*1000
# @test isa(v, OffsetArray)
# cv  = OffsetArray([1,1e100,1e100,2], (-3,))*1000
# cv2 = OffsetArray([1,-1e100,-1e100,2], (5,))*1000
# @test isequal(cumsum_kbn(v), cv)
# @test isequal(cumsum_kbn(v2), cv2)
# @test isequal(sum_kbn(v), sum_kbn(parent(v)))

@testset "Collections" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test unique(A, dims=1) == OffsetArray(parent(A), 0, first(axes(A, 2)) - 1)
    @test unique(A, dims=2) == OffsetArray(parent(A), first(axes(A, 1)) - 1, 0)
    v = OffsetArray(rand(8), (-2,))
    @test sort(v) == OffsetArray(sort(parent(v)), v.offsets)
    @test sortslices(A; dims=1) == OffsetArray(sortslices(parent(A); dims=1), A.offsets)
    @test sortslices(A; dims=2) == OffsetArray(sortslices(parent(A); dims=2), A.offsets)
    @test sort(A, dims = 1) == OffsetArray(sort(parent(A), dims = 1), A.offsets)
    @test sort(A, dims = 2) == OffsetArray(sort(parent(A), dims = 2), A.offsets)

    @test mapslices(sort, A, dims = 1) == OffsetArray(mapslices(sort, parent(A), dims = 1), A.offsets)
    @test mapslices(sort, A, dims = 2) == OffsetArray(mapslices(sort, parent(A), dims = 2), A.offsets)
end

@testset "rot/reverse" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test rotl90(A) == OffsetArray(rotl90(parent(A)), A.offsets[[2,1]])
    @test rotr90(A) == OffsetArray(rotr90(parent(A)), A.offsets[[2,1]])
    @test reverse(A, dims = 1) == OffsetArray(reverse(parent(A), dims = 1), A.offsets)
    @test reverse(A, dims = 2) == OffsetArray(reverse(parent(A), dims = 2), A.offsets)
end

@testset "fill" begin
    B = fill(5, 1:3, -1:1)
    @test axes(B) == (1:3,-1:1)
    @test all(B.==5)

    B = fill(5, (1:3, -1:1))
    @test axes(B) == (1:3,-1:1)
    @test all(B.==5)

    B = fill(5, 3, -1:1)
    @test axes(B) == (1:3,-1:1)
    @test all(B.==5)
end

@testset "broadcasting" begin
    A = OffsetArray(rand(4,4), (-3,5))

    @test A.+1 == OffsetArray(parent(A).+1, A.offsets)
    @test 2*A == OffsetArray(2*parent(A), A.offsets)
    @test A+A == OffsetArray(parent(A)+parent(A), A.offsets)
    @test A.*A == OffsetArray(parent(A).*parent(A), A.offsets)
end

@testset "@inbounds" begin
    a = OffsetArray(zeros(7), -3:3)
    unsafe_fill!(x) = @inbounds(for i in axes(x,1); x[i] = i; end)
    function unsafe_sum(x)
        s = zero(eltype(x))
        @inbounds for i in axes(x,1)
            s += x[i]
        end
        s
    end
    unsafe_fill!(a)
    for i = -3:3
        @test a[i] == i
    end
    @test unsafe_sum(a) == 0
end

@testset "Resizing OffsetVectors" begin
    local a = OffsetVector(rand(5),-3)
    axes(a,1) == -2:2
    length(a) == 5
    resize!(a,3)
    length(a) == 3
    axes(a,1) == -2:0
    @test_throws ArgumentError resize!(a,-3)
end

####
#### type defined for testing no_offset_view
####

struct NegativeArray{T,N,S <: AbstractArray{T,N}} <: AbstractArray{T,N}
    parent::S
end

# Note: this defines the axes-of-the-axes to be OneTo.
# In general this isn't recommended, because
#    positionof(A, i, j, ...) == map(getindex, axes(A), (i, j, ...))
# is quite desirable, and this requires that the axes be "identity" ranges, i.e.,
# `r[i] == i`.
# Nevertheless it's useful to test this on a "broken" implementation
# to make sure we still get the right answer.
Base.axes(A::NegativeArray) = map(n -> (-n):(-1), size(A.parent))

Base.size(A::NegativeArray) = size(A.parent)

function Base.getindex(A::NegativeArray{T,N}, I::Vararg{Int,N}) where {T,N}
    getindex(A.parent, (I .+ size(A.parent) .+ 1)...)
end

struct PointlessWrapper{T,N, A <: AbstractArray{T,N}} <: AbstractArray{T,N}
    parent :: A
end
Base.parent(x::PointlessWrapper) = x.parent
Base.size(x::PointlessWrapper) = size(parent(x))
Base.axes(x::PointlessWrapper) = axes(parent(x))

@testset "no offset view" begin
    # OffsetArray fallback
    A = randn(3, 3)
    @inferred no_offset_view(A)
    O1 = OffsetArray(A, -1:1, 0:2)
    O2 = OffsetArray(O1, -2:0, -3:(-1))
    @test no_offset_view(O2) ≡ A
    @inferred no_offset_view(O1)
    @inferred no_offset_view(O2)

    P = PointlessWrapper(A)
    @test no_offset_view(P) ≡ P

    # generic fallback
    A = collect(reshape(1:12, 3, 4))
    N = NegativeArray(A)
    @test N[-3, -4] == 1
    V = no_offset_view(N)
    @test collect(V) == A

    # bidirectional
    B = BidirectionalVector([1, 2, 3])
    pushfirst!(B, 0)
    OB = OffsetArrays.no_offset_view(B)
    @test axes(OB, 1) == 1:4
    @test collect(OB) == 0:3
end

@testset "no nesting" begin
    A = randn(2, 3)
    x = A[2, 2]
    O1 = OffsetArray(A, -1:0, -1:1)
    O2 = OffsetArray(O1, 0:1, 0:2)
    @test parent(O1) ≡ parent(O2)
    @test eltype(O1) ≡ eltype(O2)
    O2[1, 1] = x + 1            # just a sanity check
    @test A[2, 2] == x + 1
end

@testset "mutating functions for OffsetVector" begin
    # push!
    o = OffsetVector(Int[], -1)
    @test push!(o) === o
    @test axes(o, 1) == 0:-1
    @test push!(o, 1) === o
    @test axes(o, 1) == 0:0
    @test o[end] == 1
    @test push!(o, 2, 3) === o
    @test axes(o, 1) == 0:2
    @test o[end-1:end] == [2, 3]
    # pop!
    o = OffsetVector([1, 2, 3], -1)
    @test pop!(o) == 3
    @test axes(o, 1) == 0:1
    # append!
    o = OffsetVector([1, 2, 3], -1)
    append!(o, [4, 5])
    @test axes(o, 1) == 0:4
    # empty!
    o = OffsetVector([1, 2, 3], -1)
    @test empty!(o) === o
    @test axes(o, 1) == 0:-1
end

@testset "searchsorted (#85)" begin
    o = OffsetVector([1,3,4,5],-2)
    @test searchsortedfirst(o,-2) == -1
    @test searchsortedfirst(o, 1) == -1
    @test searchsortedfirst(o, 2) ==  0
    @test searchsortedfirst(o, 5) ==  2
    @test searchsortedfirst(o, 6) ==  3
    @test searchsortedlast(o, -2) == -2
    @test searchsortedlast(o,  1) == -1
    @test searchsortedlast(o,  2) == -1
    @test searchsortedlast(o,  5) ==  2
    @test searchsortedlast(o,  6) ==  2
    @test searchsorted(o, -2) == -1:-2
    @test searchsorted(o,  1) == -1:-1
    @test searchsorted(o,  2) ==  0:-1
    @test searchsorted(o,  5) ==  2:2
    @test searchsorted(o,  6) ==  3:2
end

@testset "Adapt" begin
    # We need another storage type, CUDA.jl defines one but we can't use that for CI
    # let's define an appropriate method for SArrays
    Adapt.adapt_storage(::Type{SA}, xs::AbstractArray) where SA<:SArray = convert(SA, xs)
    arr = OffsetArray(rand(3, 3), -1:1, -1:1)
    s_arr = adapt(SMatrix{3,3}, arr)
    @test parent(s_arr) isa SArray
    @test arr == adapt(Array, s_arr)
end

@testset "Pointer" begin
    a = OffsetVector(collect(10:20), 9);
    @test 12 == a[12] == unsafe_load(pointer(a), 12 + (1 - firstindex(a))) == unsafe_load(pointer(a, 12))

    A = OffsetArray(reshape(collect(10:130), (11,11)), 9, 9);
    @test 21 == A[12] == unsafe_load(pointer(A), 12) == unsafe_load(pointer(A, 12))
    @test 61 == A[52] == unsafe_load(pointer(A), 52) == unsafe_load(pointer(A, 52))

    @test pointer(a) === pointer(parent(a))
    @test pointer(A) === pointer(parent(A))
    @test pointer(a, 12) === pointer(parent(a), 12 + (1 - firstindex(a)))
    @test pointer(A, 12) === pointer(parent(A), 12)
    @test pointer(a) === pointer(a, firstindex(a))
    @test pointer(A) === pointer(A, firstindex(A))
    if VERSION ≥ v"1.5"
        @test pointer(a') === pointer(parent(a))
        @test pointer(A') === pointer(parent(A))
        @test pointer(a', 5) === pointer(parent(a), 5)
        @test pointer(A', 15) === pointer(parent(A)', 15)
    end
end

include("origin.jl")


text:
include("forecast_types.jl")
include("helper_functions.jl")

function count_cases_in_windows!(
    forecast::Forecast,
    case_counts,
    window_lengths,
    D_total,
    sim
)
    """
    Count the cases in each observation period and store in case_counts.
    """

    t_start = 1
    for (i,t) in enumerate(window_lengths)
        case_counts[i] = sum(@view D_total[t_start:t,sim])
        t_start = t
    end

    return nothing
end

function check_sim!(
    forecast::Forecast,
    forecast_start_date,
    omicron_dominant_date,
    case_counts,
    sim,
    local_cases,
    window_lengths,
    min_cases,
    max_cases,
    reinitialise_allowed;
    day = 0,
)
    """
    Check for consistency of the simulations against the data. This will also
    check for instances of superspreading events and will insert cases if conditions
    are met.
    """

    print_status = false

    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U

    # days forecast observed for
    T_observed = forecast.sim_features.T_observed
    max_forecast_cases = forecast.sim_features.max_forecast_cases

    consistency_multiplier = forecast.sim_constants.consistency_multiplier

    # initialise the bad sim
    bad_sim = false
    injected_cases = false

    days_delta = (Dates.Date(omicron_dominant_date) - Dates.Date(forecast_start_date)).value

    cases_pre_backcast = sum(@view D[1:days_delta,1:2,sim])
    cases_backcast = sum(@view D[days_delta+1:T_observed,1:2,sim])
    cases_pre_nowcast = sum(@view D[21:T_observed,1:2,sim])
    cases_nowcast = sum(@view D[10:T_observed,1:2,sim])

    case_counts[1] = cases_pre_backcast
    case_counts[2] = cases_backcast
    case_counts[3] = cases_pre_nowcast
    case_counts[4] = cases_nowcast

    # this is just the total cases over the forecast horizon
    D_forecast = sum(@view D[T_observed+1:end,1:2,sim])

    # if we've exceeded the max cases over a given period
    if D_forecast > max_forecast_cases
        bad_sim = true
        print_status && println("too many forecast cases")
    end

    if !bad_sim
        for (i, (c,m)) in enumerate(zip(case_counts,max_cases))
            if c > m
                bad_sim = true
                print_status && println("too many cases", c, " ", m, " day: ", day)
                break
            end
        end

        if day == 0
            for (i, (c,m)) in enumerate(zip(case_counts,min_cases))
                if c < m
                    bad_sim = true
                    print_status && println("too few cases", c, " ", m, " day: ", day)
                    break
                end
            end
        end
    end

    # see whether we need to inject any cases
    injected_cases = false

    if !bad_sim && reinitialise_allowed
        # calculate the number of detections missed over the 5 day period
        actual_3_day_cases = sum(@view local_cases[day-2:day])
        sim_3_day_cases = sum(@view D[day-2:day,:,sim])
        sim_3_day_threshold = consistency_multiplier*max(1, sim_3_day_cases)
        # calculate total number of cases over week
        sim_week_cases = sum(@view D[day-6:day,:,sim])
        missing_detections = 0

        if actual_3_day_cases > sim_3_day_threshold || (sim_week_cases == 0 && actual_3_day_cases > 0)
            print_status && println("actual: ", actual_3_day_cases, " 3day: ", sim_3_day_cases, " 3day thresh: ", sim_3_day_threshold, " week: ", sim_week_cases, " day added: ", day)

            if sim_week_cases == 0
                missing_detections = rand(1:actual_3_day_cases)
            else
                missing_detections = rand(1:(actual_3_day_cases - sim_3_day_cases))
            end

            injected_cases = true
            inject_cases!(
                forecast::Forecast,
                missing_detections,
                day,
                sim,
            )
        end
    end

    return (bad_sim, injected_cases)
end

function inject_cases!(
    forecast::Forecast,
    missing_detections,
    day,
    sim,
)
    """
    A model for injecting cases into the simulation following long periods of low
    incidence.
    """

    Z = forecast.sim_realisations.Z
    D = forecast.sim_realisations.D
    U = forecast.sim_realisations.U
    individual_type_map = forecast.individual_type_map

    total_injected = missing_detections

    p_symp = forecast.sim_constants.p_symp
    p_detect_given_symp = forecast.sim_constants.p_detect_given_symp
    p_symp_given_detect = forecast.sim_constants.p_symp_given_detect

    omicron_dominant_day = forecast.sim_features.omicron_dominant_day

    # sample number detected
    num_symptomatic_detected = sample_binomial_limit(missing_detections, p_symp_given_detect)
    num_asymptomatic_detected = missing_detections - num_symptomatic_detected
    # infer some undetected symptomatic
    num_symptomatic_undetected = sample_negative_binomial_limit(missing_detections, p_detect_given_symp)

    # infer some undetected asumptomatic
    num_symptomatic_total = num_symptomatic_detected + num_symptomatic_undetected

    if num_symptomatic_total == 0
        num_asymptomatic_undetected = sample_negative_binomial_limit(1, p_symp)
    else
        num_asymptomatic_undetected = sample_negative_binomial_limit(num_symptomatic_total, p_symp)
    end

    # detection times are determined triangularly
    num_symptomatic_each_day = (
        num_symptomatic_detected ÷ 2,
        num_symptomatic_detected ÷ 2 + num_symptomatic_detected ÷ 3,
        num_symptomatic_detected ÷ 2 + num_symptomatic_detected ÷ 3 + num_symptomatic_detected ÷ 6
    )
    counter = 0
    onset_time = day

    # infection time is then 1 generation interval earlier
    for _ in 1:num_symptomatic_detected
        total_injected += 1
        if counter >= num_symptomatic_each_day[2]
            onset_time = day-1
        elseif counter >= num_symptomatic_each_day[3]
            onset_time = day-2
        end
        infection_time = onset_time - sample_onset_time(omicron = onset_time >= omicron_dominant_day)
        Z[infection_time+36,individual_type_map.S,sim] += 1
        D[onset_time,individual_type_map.S,sim] += 1
        counter += 1
    end

    # detection times are determined triangularly
    num_asymptomatic_each_day = (
        num_asymptomatic_detected ÷ 2,
        num_asymptomatic_detected ÷ 2 + num_asymptomatic_detected ÷ 3,
        num_asymptomatic_detected ÷ 2 + num_asymptomatic_detected ÷ 3 + num_asymptomatic_detected ÷ 6
    )
    counter = 0
    onset_time = day

    # infection time is then 1 generation interval earlier
    for _ in 1:num_asymptomatic_detected
        total_injected += 1
        if counter >= num_asymptomatic_each_day[2]
            onset_time = day-1
        elseif counter >= num_asymptomatic_each_day[3]
            onset_time = day-2
        end
        infection_time = onset_time - sample_onset_time(omicron = onset_time >= omicron_dominant_day)
        Z[infection_time+36,individual_type_map.A,sim] += 1
        D[onset_time,individual_type_map.A,sim] += 1
        counter += 1
    end

    # now add in the inferred undetected cases in the same way. Note that the
    # onset times are kinda irrelevant here, but just help with consistent
    # simulation recording.
    num_symptomatic_each_day = (
        num_symptomatic_undetected ÷ 2,
        num_symptomatic_undetected ÷ 2 + num_symptomatic_undetected ÷ 3,
        num_symptomatic_undetected ÷ 2 + num_symptomatic_undetected ÷ 3 + num_symptomatic_undetected ÷ 6
    )
    counter = 0
    onset_time = day

    for _ in 1:num_symptomatic_undetected
        total_injected += 1
        if counter >= num_symptomatic_each_day[2]
            onset_time = day-1
        elseif counter >= num_symptomatic_each_day[3]
            onset_time = day-2
        end
        infection_time = onset_time - sample_onset_time(omicron = onset_time >= omicron_dominant_day)
        Z[infection_time+36,individual_type_map.S,sim] += 1
        U[onset_time,individual_type_map.S,sim] += 1
        counter += 1
    end

    # detection times are determined triangularly
    num_asymptomatic_each_day = (
        num_asymptomatic_undetected ÷ 2,
        num_asymptomatic_undetected ÷ 2 + num_asymptomatic_undetected ÷ 3,
        num_asymptomatic_undetected ÷ 2 + num_asymptomatic_undetected ÷ 3 + num_asymptomatic_undetected ÷ 6
    )
    counter = 0
    onset_time = day

    # now add in the inferred undetected cases
    for _ in 1:num_asymptomatic_undetected
        total_injected += 1
        if counter >= num_asymptomatic_each_day[2]
            onset_time = day-1
        elseif counter >= num_asymptomatic_each_day[3]
            onset_time = day-2
        end
        infection_time = onset_time - sample_onset_time(omicron = onset_time >= omicron_dominant_day)
        Z[infection_time+36,individual_type_map.A,sim] += 1
        U[onset_time,individual_type_map.A,sim] += 1
        counter += 1
    end

    # println("total injected is: ", total_injected, " at time: ", day)

    return nothing
end


text:
# This program computes Feigenbaum's alpha constant using
# Feigenbaum's functional equation derived via renormalization
# group methods.  This version uses a Taylor's expansion for
# the function g(z).  The expansion is evaluated on a set of
# Chebyshev nodes.
#
# Started: Nov 2017 -- Stuart Brorson, sdb@cloud9.net

using ForwardDiff
using JLD

#++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Helper fcns

#-------------------------------------------------------
# Now create the Chebyshev points
function generate_cheb_points(N)
  #println("---->  Entering generate_cheb_points\n")
  local k = big(1.0):big(N)
  local bigpi = big(pi);
  zi = cos.(bigpi*(2*k-1)./(4*N));
  #println(zi)
  #println("<----  Leaving generate_cheb_points\n")
  return zi;
end

#-------------------------------------------------------
# Define universal fcn g(z, a).  It is expressed as a
# Taylor's series, coefficients a[i].  Only even powers
# are needed.  Use Horner's rule to perform sum.
function g{T<:Real}(z, a::Vector{T})
  #println("  ---->  Entering g, z = $z\n")
  local s = zero(z);
  for i=length(a):-1:1
    s = z*z*(s + a[i]);
  end
  #println("  <----  Leaving g, s = $s\n")
  return (one(z)+s)
end

#---------------------------------------------------------
# Define relation obeyed by g(z).  Name it f(z, a)
function f{T<:Real}(z, a::Vector{T})
  #println("---->  Entering f, a = $a\n")
  local myone = one(z);
  local alpha = myone/g(myone, a);
  r = g(z, a) - alpha*g(g(z/alpha, a), a);
  #println("<----  Leaving f\n")
  return r
end

#--------------------------------------------------------
# Define fcn which takes gradient of f w.r.t. a.
# Returns row vector which is gradient of f w.r.t. a,
# evaluated at z, a.
function gradf(z, a::Vector{BigFloat})
  # println("---> Entered gradf\n")
  # Specialize to position z.  Result is fcn of coeffs a only.
  fbeta(a) = f(z, a);

  # Compute gradient w.r.t. beta and return it
  gf1 = ForwardDiff.gradient(fbeta, a)
  # println("<--- Leaving gradf\n")
  return gf1
end

# ------------------------------------------------------------
# The main event -- a computation of alpha using Newton's method.
function compute_alpha(Numdigs, Npts, betain)

  # First set up variables, including BigFloat settings
  # Setting precision to 5x Numdigs seems to keep Newton's method
  # from wandering away.  Smaller values don't.  More fidding might
  # improve things here.
  setprecision(3*Numdigs);

  # Set stopping tol for Newton's method so I get the number of digits
  # I want.
  const tol = BigFloat(10.0)^-(Npts);

  zi = generate_cheb_points(Npts);      # Grid of sample points zi

  fn = zeros(BigFloat, Npts);             # Function f vector
  Jn = zeros(BigFloat, Npts, Npts);        # Jacobian

  # Need expansion coeffs beta.  Start with the one passed in, betain.
  # N should always be larger than length(betain), so we are adding to the
  # betas we computed last time.
  betan = zeros(BigFloat, Npts);
  for i = 1:length(betain)
    betan[i] = betain[i];
  end

  # println("-----------------------------\n")
  sttim = now();  # Keep track of loop timing.

  # Now enter Newton loop.
  for cnt = 1:10

    # Compute new f and Jacobian upon each iteration
    #println("Computing fn and Jn...")
    for i = 1:Npts
      fn[i] = f(zi[i], betan);
      Jn[i,:] = transpose(gradf(zi[i], betan));
    end

    #println(fn)
    #println(Jn)

    # Now compute step to take, then take it.
    #println("Computing step sn...")
    sn = Jn\fn;
    #println("Step = $sn\n")
    betanp1 = betan - sn;
    #println("betanp1 = $betanp1\n")

    nowtim = now();  # Keep track of loop timing
    println("After Newton iteration = $cnt, time elapsed = $(nowtim-sttim)\n")

    # Write out betanp1 and Npts values in case we need to restart the
    # calculation
    save("./betan.jld", "betan", betanp1, "Npts", Npts)

    # Check for convergence
    if (norm(sn) < tol)
      println("\nDone!  Converged after $cnt iterations.\n")
      # compute alpha
      myalpha = big(1.0)/g(big(1.0), betanp1);
      return myalpha, betanp1;
    end

    # Move values back
    betan = betanp1

  end
  println("!!!!!!!  Didn't converge!  !!!!!!!!\n")

end

#--------------------------------------------------
function iterate_alpha()
  # This is a runner -- call this fcn.  This fcn calls
  # compute_alpha asking for different number of digits
  # at different precision/tolerance levels.
  # It calls with the beta coefficients computed in the
  # last iteration as the starting point to help Newton's
  # method from wandering away.

  # Seed betan, either with known start values, or with
  # the values computed at the last Newton iteration.
  println("Use betan value from file [y/N]? ==>")
  resp = readline(STDIN)
  if (resp == "")
    betan = zeros(BigFloat, 3);
    betan[1] = big(-1.5);
    betan[2] = big(1e-1);
    betan[3] = big(2e-2);
    N0 = 4;
  else
    betan = load("./betan.jld", "betan")
    N0 = load("./betan.jld", "Npts")
  end

  # Step up number of pts to compute on.  Change this to get
  # more digits
  for Npts = N0:2*N0:2500
    println("=================================\n")
    println("Npts = $Npts\n")
    Numdigs = Int(floor(2*Npts));

    # Call the function which computes alpha at this level
    # of precision.
    myalpha, betan = compute_alpha(Numdigs, Npts, betan)

    # Print out myalpha
    a = string(myalpha);
    println("myalpha = $(a)\n")

    f = open("alpha.dat", "a");
    write(f, "==========\nNpts = $Npts\n$(a)\n")
    close(f)

    # Compare against the version from
    # http://www.plouffe.fr/simon/constants/feigenbaum.txt
    oeisalpha = oeis_alpha(3*Numdigs);
    diff = Float64(abs(myalpha) - abs(oeisalpha));
    println("diff = $diff\n")
  end

end

#------------------------------------------------
function oeis_alpha(Numdigs)
  # Return up to 1018 digits of alpha from
  # http://www.plouffe.fr/simon/constants/feigenbaum.txt

  if (Numdigs <= 1018)
    setprecision(Numdigs);
  else
    setprecision(1018);
  end

  str = "-2."
  str = str*"5029078750958928222839028732182157863812713767271499773361920567"
  str = str*"7923546317959020670329964974643383412959523186999585472394218237"
  str = str*"7785445179272863314993372578112163594879503744781260997380598671"
  str = str*"2397117373289276654044010306698313834600094139322364490657889951"
  str = str*"2205843172507873377463087853424285351988587500042358246918740820"
  str = str*"4281700901714823051821621619413199856066129382742649709844084470"
  str = str*"1008054549677936760888126446406885181552709324007542506497157047"
  str = str*"0475419932831783645332562415378693957125097066387979492654623137"
  str = str*"6745918909813116752434221110130913127837160951158341230841503716"
  str = str*"4997020224681219644081216686527458043026245782561067150138521821"
  str = str*"6449532543349873487413352795815351016583605455763513276501810781"
  str = str*"1948369459574850237398235452625632779475397269902012891516645793"
  str = str*"9420198920248803394051699686551494477396533876979741232354061781"
  str = str*"9896112494095990353128997733611849847377946108428833293833903950"
  str = str*"9008914086351525626803381414669279913310743349705143545201344643"
  str = str*"4264752001621384610729922641994332772918977769053802596851"

  alpha = parse(BigFloat, str)

end

max_stars_repo_name: brorson/FeigenbaumConstants

max_issues_repo_name: brorson/FeigenbaumConstants

max_forks_repo_name: brorson/FeigenbaumConstants


text:
using Revise
using Flipping
using CSV
##
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Pharmacology/PharmacologyFlipping"
pokes, bouts, streaks, DataIndex  = Flipping.create_exp_dataframes(Directory_path)
bouts
##
groupA = ["WJ0"*string(x) for x in 1:8]
groupB = ["WJ09"]
append!(groupB,["WJ1"*string(x) for x in 0:6])
groupB
groupdic = Dict()
for x in groupA
    groupdic[x]="A"
end
for x in groupB
    groupdic[x]="B"
end
drugDic = Dict(
    "A1a" => "None",
    "A1b" => "None",
    "A2a" => "Veh",
    "A2b" => "Veh",
    "A3a" => "PreVeh",
    "A3b" => "Altanserin_0.5",
    "A4a" => "PostVeh",
    "A5a" => "None",
    "A5b" => "None",
    "A6a" => "Veh",
    "A6b" => "PreVeh",
    "A7a" => "Altanserin_0.5",
    "A7b" => "PostVeh",
    "A8a" => "Veh",
    "B1a" => "None",
    "B1b" => "None",
    "B2a" => "Veh",
    "B2b" => "PreVeh",
    "B3a" => "Altanserin_0.5",
    "B3b" => "PostVeh",
    "B4a" => "Veh",
    "B5a" => "None",
    "B5b" => "None",
    "B6a" => "Veh",
    "B6b" => "Veh",
    "B7a" => "PreVeh",
    "B7b" => "Altanserin_0.5",
    "B8a" => "PostVeh",
    "A9a" => "None",
    "A9b" => "None",
    "A10a" => "Veh",
    "A10b" => "Veh",
    "A11a" => "PreVeh",
    "A11b" => "Altanserin_0.25",
    "A12a" => "PostVeh",
    "A13a" => "None",
    "A13b" => "None",
    "A14a" => "Veh",
    "A14b" => "PreVeh",
    "A15a" => "Altanserin_0.25",
    "A15b" => "PostVeh",
    "A16a" => "Veh",
    "B9a" => "None",
    "B9b" => "None",
    "B10a" => "Veh",
    "B10b" => "PreVeh",
    "B11a" => "Altanserin_0.25",
    "B11b" => "PostVeh",
    "B12a" => "Veh",
    "B13a" => "None",
    "B13b" => "None",
    "B14a" => "Veh",
    "B14b" => "Veh",
    "B15a" => "PreVeh",
    "B15b" => "Altanserin_0.25",
    "B16a" => "PostVeh",)
##
streaks.Group = [get(groupdic,x,"missing") for x in streaks.MouseID]
streaks.ExpSession = string.(streaks.Exp_Day).*streaks.Daily_Session
streaks.GroupSession = streaks.Group .* streaks.ExpSession
println(union(streaks.GroupSession))
streaks.Drug = [get(drugDic,x,missing) for x in streaks.GroupSession]
findall(ismissing,streaks.Drug)
filetosave = joinpath(Directory_path,"Results","Full_info_streaks.csv")
CSV.write(filetosave,streaks)
##
bouts.Group = [get(groupdic,x,"missing") for x in bouts.MouseID]
bouts.ExpSession = string.(bouts.Exp_Day).*bouts.Daily_Session
bouts.GroupSession = bouts.Group .* bouts.ExpSession
println(union(bouts.GroupSession))
bouts.Drug = [get(drugDic,x,missing) for x in bouts.GroupSession]
findall(ismissing,bouts.Drug)
filetosave = joinpath(Directory_path,"Results","Full_info_bouts.csv")
CSV.write(filetosave,bouts)
##
using FlippingModel


text:
using Tables

sym(ptr) = ccall(:jl_symbol, Ref{Symbol}, (Ptr{UInt8},), ptr)

struct Query{NT}
    stmt::Stmt
    status::Base.RefValue{Cint}
end

Tables.istable(::Type{<:Query}) = true
Tables.rowaccess(::Type{<:Query}) = true
Tables.rows(q::Query) = q
Tables.schema(q::Query{NamedTuple{names, types}}) where {names, types} = Tables.Schema(names, types)

Base.IteratorSize(::Type{<:Query}) = Base.SizeUnknown()
Base.eltype(q::Query{NT}) where {NT} = NT

function reset!(q::Query)
    sqlite3_reset(q.stmt.handle)
    q.status[] = execute!(q.stmt)
    return
end

function done(q::Query)
    st = q.status[]
    st == SQLITE_DONE && return true
    st == SQLITE_ROW || sqliteerror(q.stmt.db)
    return false
end

function getvalue(q::Query, col::Int, ::Type{T}) where {T}
    handle = q.stmt.handle
    t = sqlite3_column_type(handle, col)
    if t == SQLITE_NULL
        return missing
    else
        TT = juliatype(t) # native SQLite Int, Float, and Text types
        return sqlitevalue(ifelse(TT === Any && !isbitstype(T), T, TT), handle, col)
    end
end

function generate_namedtuple(::Type{NamedTuple{names, types}}, q) where {names, types}
    if @generated
        vals = Tuple(:(getvalue(q, $i, $(fieldtype(types, i)))) for i = 1:fieldcount(types))
        return :(NamedTuple{names, types}(($(vals...),)))
    else
        return NamedTuple{names, types}(Tuple(getvalue(q, i, fieldtype(types, i)) for i = 1:fieldcount(types)))
    end
end

function Base.iterate(q::Query{NT}) where {NT}
    done(q) && return nothing
    nt = generate_namedtuple(NT, q)
    return nt, nothing
end

function Base.iterate(q::Query{NT}, ::Nothing) where {NT}
    q.status[] = sqlite3_step(q.stmt.handle)
    done(q) && return nothing
    nt = generate_namedtuple(NT, q)
    return nt, nothing
end

"""
`SQLite.Query(db, sql::String; values=[]; stricttypes::Bool=true, nullable::Bool=true)`

Constructs a `SQLite.Query` object by executing the SQL query `sql` against the sqlite database `db` and querying
the columns names and types of the result set, if any.

Will bind `values` to any parameters in `sql`.
`stricttypes=false` will remove strict column typing in the result set, making each column effectively `Vector{Any}`; in sqlite, individual
column values are only loosely associated with declared column types, and instead each carry their own type information. This can lead to
type errors when trying to query columns when a single type is expected.
`nullable` controls whether `NULL` (`missing` in Julia) values are expected in a column.

An `SQLite.Query` object will iterate NamedTuple rows by default, and also supports the Tables.jl interface for integrating with
any other Tables.jl implementation. Due note however that iterating an sqlite result set is a forward-once-only operation. If you need
to iterate over an `SQLite.Query` multiple times, but can't store the iterated NamedTuples, call `SQLite.reset!(q::SQLite.Query)` to
re-execute the query and position the iterator back at the begining of the result set.
"""
function Query(db::DB, sql::AbstractString; values=[], stricttypes::Bool=true, nullable::Bool=true)
    stmt = Stmt(db, sql)
    bind!(stmt, values)
    status = execute!(stmt)
    cols = sqlite3_column_count(stmt.handle)
    header = Vector{Symbol}(undef, cols)
    types = Vector{Type}(undef, cols)
    for i = 1:cols
        header[i] = sym(sqlite3_column_name(stmt.handle, i))
        if nullable
            types[i] = stricttypes ? Union{juliatype(stmt.handle, i), Missing} : Any
        else
            types[i] = stricttypes ? juliatype(stmt.handle, i) : Any
        end
    end
    return Query{NamedTuple{Tuple(header), Tuple{types...}}}(stmt, Ref(status))
end

# as a sink
function createtable!(db::DB, nm::AbstractString, ::Tables.Schema{names, types}; temp::Bool=false, ifnotexists::Bool=true) where {names, types}
    temp = temp ? "TEMP" : ""
    ifnotexists = ifnotexists ? "IF NOT EXISTS" : ""
    typs = [types === nothing ? "BLOB" : sqlitetype(fieldtype(types, i)) for i = 1:length(names)]
    columns = [string(esc_id(String(names[i])), ' ', typs[i]) for i = 1:length(names)]
    return execute!(db, "CREATE $temp TABLE $ifnotexists $nm ($(join(columns, ',')))")
end

"""
    source |> SQLite.load!(db::SQLite.DB, tablename::String; temp::Bool=false, ifnotexists::Bool=false)
    SQLite.load!(source, db, tablename; temp=false, ifnotexists=false)

Load a Tables.jl input `source` into an SQLite table that will be named `tablename` (will be auto-generated if not specified).

`temp=true` will create a temporary SQLite table that will be destroyed automatically when the database is closed
`ifnotexists=false` will throw an error if `tablename` already exists in `db`
"""
function load! end

load!(db::DB, table::AbstractString="sqlitejl_"*Random.randstring(5); kwargs...) = x->load!(x, db, table; kwargs...)

function load!(itr, db::DB, name::AbstractString="sqlitejl_"*Random.randstring(5); kwargs...)
    # check if table exists
    nm = esc_id(name)
    checkstmt = Stmt(db, "pragma table_info($nm)")
    status = execute!(checkstmt)
    rows = Tables.rows(itr)
    sch = Tables.schema(rows)
    return load!(sch, rows, db, nm, name, status == SQLITE_DONE; kwargs...)
end

function load!(sch::Tables.Schema, rows, db::DB, nm::AbstractString, name, shouldcreate; temp::Bool=false, ifnotexists::Bool=false)
    # create table if needed
    shouldcreate && createtable!(db, nm, sch; temp=temp, ifnotexists=ifnotexists)
    # build insert statement
    params = chop(repeat("?,", length(sch.names)))
    stmt = Stmt(db, "INSERT INTO $nm VALUES ($params)")
    # start a transaction for inserting rows
    transaction(db) do
        for row in rows
            Tables.eachcolumn(sch, row) do val, col, _
                bind!(stmt, col, val)
            end
            sqlite3_step(stmt.handle)
            sqlite3_reset(stmt.handle)
        end
    end
    execute!(db, "ANALYZE $nm")
    return name
end

# unknown schema case
function load!(::Nothing, rows, db::DB, nm::AbstractString, name, shouldcreate; temp::Bool=false, ifnotexists::Bool=false)
    state = iterate(rows)
    state === nothing && return nm
    row, st = state
    names = propertynames(row)
    sch = Tables.Schema(names, nothing)
    # create table if needed
    shouldcreate && createtable!(db, nm, sch; temp=temp, ifnotexists=ifnotexists)
    # build insert statement
    params = chop(repeat("?,", length(names)))
    stmt = Stmt(db, "INSERT INTO $nm VALUES ($params)")
    # start a transaction for inserting rows
    transaction(db) do
        while true
            Tables.eachcolumn(sch, row) do val, col, _
                bind!(stmt, col, val)
            end
            sqlite3_step(stmt.handle)
            sqlite3_reset(stmt.handle)
            state = iterate(rows)
            state === nothing && break
            row, st = state
        end
    end
    execute!(db, "ANALYZE $nm")
    return name
end