text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

A = rand(5,4,3)
@testset "Bounds checking" begin
    @test checkbounds(Bool, A, 1, 1, 1) == true
    @test checkbounds(Bool, A, 5, 4, 3) == true
    @test checkbounds(Bool, A, 0, 1, 1) == false
    @test checkbounds(Bool, A, 1, 0, 1) == false
    @test checkbounds(Bool, A, 1, 1, 0) == false
    @test checkbounds(Bool, A, 6, 4, 3) == false
    @test checkbounds(Bool, A, 5, 5, 3) == false
    @test checkbounds(Bool, A, 5, 4, 4) == false
    @test checkbounds(Bool, A, 1) == true           # linear indexing
    @test checkbounds(Bool, A, 60) == true
    @test checkbounds(Bool, A, 61) == false
    @test checkbounds(Bool, A, 2, 2, 2, 1) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, 2) == false
    @test checkbounds(Bool, A, 1, 1)  == true       # partial linear indexing (PLI)
    # @test checkbounds(Bool, A, 1, 12) == false     # PLI TODO: Re-enable after partial linear indexing deprecation
    # @test checkbounds(Bool, A, 5, 12) == false     # PLI TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, 1, 13) == false      # PLI
    # @test checkbounds(Bool, A, 6, 12) == false    # PLI TODO: Re-enable after partial linear indexing deprecation
end

@testset "single CartesianIndex" begin
    @test checkbounds(Bool, A, CartesianIndex((1, 1, 1))) == true
    @test checkbounds(Bool, A, CartesianIndex((5, 4, 3))) == true
    @test checkbounds(Bool, A, CartesianIndex((0, 1, 1))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 0, 1))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 1, 0))) == false
    @test checkbounds(Bool, A, CartesianIndex((6, 4, 3))) == false
    @test checkbounds(Bool, A, CartesianIndex((5, 5, 3))) == false
    @test checkbounds(Bool, A, CartesianIndex((5, 4, 4))) == false
    @test checkbounds(Bool, A, CartesianIndex((1,))) == true
    # @test checkbounds(Bool, A, CartesianIndex((60,))) == false     # TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, CartesianIndex((61,))) == false
    @test checkbounds(Bool, A, CartesianIndex((2, 2, 2, 1,))) == true
    @test checkbounds(Bool, A, CartesianIndex((2, 2, 2, 2,))) == false
    @test checkbounds(Bool, A, CartesianIndex((1, 1,)))  == true
    # @test checkbounds(Bool, A, CartesianIndex((1, 12,))) == false  # TODO: Re-enable after partial linear indexing deprecation
    # @test checkbounds(Bool, A, CartesianIndex((5, 12,))) == false  # TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, CartesianIndex((1, 13,))) == false
    # @test checkbounds(Bool, A, CartesianIndex((6, 12,))) == false # TODO: Re-enable after partial linear indexing deprecation
end

@testset "mix of CartesianIndex and Int" begin
    @test checkbounds(Bool, A, CartesianIndex((1,)), 1, CartesianIndex((1,))) == true
    @test checkbounds(Bool, A, CartesianIndex((5, 4)), 3)  == true
    @test checkbounds(Bool, A, CartesianIndex((0, 1)), 1)  == false
    @test checkbounds(Bool, A, 1, CartesianIndex((0, 1)))  == false
    @test checkbounds(Bool, A, 1, 1, CartesianIndex((0,))) == false
    @test checkbounds(Bool, A, 6, CartesianIndex((4, 3)))  == false
    @test checkbounds(Bool, A, 5, CartesianIndex((5,)), 3) == false
    @test checkbounds(Bool, A, CartesianIndex((5,)), CartesianIndex((4,)), CartesianIndex((4,)))  == false
end

@testset "vector indices" begin
    @test checkbounds(Bool, A, 1:5, 1:4, 1:3) == true
    @test checkbounds(Bool, A, 0:5, 1:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 0:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:4, 0:3) == false
    @test checkbounds(Bool, A, 1:6, 1:4, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:5, 1:3) == false
    @test checkbounds(Bool, A, 1:5, 1:4, 1:4) == false
    @test checkbounds(Bool, A, 1:60) == true
    @test checkbounds(Bool, A, 1:61) == false
    @test checkbounds(Bool, A, 2, 2, 2, 1:1) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, 1:2) == false
    @test checkbounds(Bool, A, 1:5, 1:4) == true
    # @test checkbounds(Bool, A, 1:5, 1:12) == false # TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, 1:5, 1:13) == false
    # @test checkbounds(Bool, A, 1:6, 1:12) == false # TODO: Re-enable after partial linear indexing deprecation
end

@testset "logical" begin
    @test checkbounds(Bool, A, trues(5), trues(4), trues(3)) == true
    @test checkbounds(Bool, A, trues(6), trues(4), trues(3)) == false
    @test checkbounds(Bool, A, trues(5), trues(5), trues(3)) == false
    @test checkbounds(Bool, A, trues(5), trues(4), trues(4)) == false
    @test checkbounds(Bool, A, trues(60)) == true
    @test checkbounds(Bool, A, trues(61)) == false
    @test checkbounds(Bool, A, 2, 2, 2, trues(1)) == true  # extra indices
    @test checkbounds(Bool, A, 2, 2, 2, trues(2)) == false
    # @test checkbounds(Bool, A, trues(5), trues(12)) == false  # TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, trues(5), trues(13)) == false
    # @test checkbounds(Bool, A, trues(6), trues(12)) == false # TODO: Re-enable after partial linear indexing deprecation
    @test checkbounds(Bool, A, trues(5, 4, 3)) == true
    @test checkbounds(Bool, A, trues(5, 4, 2)) == false
    @test checkbounds(Bool, A, trues(5, 12)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 4, 1), trues(1, 1, 3)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 4, 1), trues(1, 1, 2)) == false
    @test checkbounds(Bool, A, trues(1, 5), trues(1, 5, 1), trues(1, 1, 3)) == false
    @test checkbounds(Bool, A, trues(1, 5), :, 2) == false
end

@testset "array of CartesianIndex" begin
    @test checkbounds(Bool, A, [CartesianIndex((1, 1, 1))]) == true
    @test checkbounds(Bool, A, [CartesianIndex((5, 4, 3))]) == true
    @test checkbounds(Bool, A, [CartesianIndex((0, 1, 1))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 0, 1))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1, 0))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((6, 4, 3))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 5, 3))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 4, 4))]) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1))], 1) == true
    @test checkbounds(Bool, A, [CartesianIndex((5, 4))], 3) == true
    @test checkbounds(Bool, A, [CartesianIndex((0, 1))], 1) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 0))], 1) == false
    @test checkbounds(Bool, A, [CartesianIndex((1, 1))], 0) == false
    @test checkbounds(Bool, A, [CartesianIndex((6, 4))], 3) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 5))], 3) == false
    @test checkbounds(Bool, A, [CartesianIndex((5, 4))], 4) == false
end

@testset "sub2ind & ind2sub" begin
    @testset "0-dimensional" begin
        for i = 1:4
            @test sub2ind((), i) == i
        end
        @test sub2ind((), 2, 2) == 3
        @test ind2sub((), 1) == ()
        @test_throws BoundsError ind2sub((), 2)
    end

    @testset "1-dimensional" begin
        for i = 1:4
            @test sub2ind((3,), i) == i
            @test ind2sub((3,), i) == (i,)
        end
        @test sub2ind((3,), 2, 2) == 5
        @test_throws MethodError ind2sub((3,), 2, 2)
        #   ambiguity btw cartesian indexing and linear indexing in 1d when
        #   indices may be nontraditional
        @test_throws ArgumentError sub2ind((1:3,), 2)
        @test_throws ArgumentError ind2sub((1:3,), 2)
    end

    @testset "2-dimensional" begin
        k = 0
        for j = 1:3, i = 1:4
            @test sub2ind((4,3), i, j) == (k+=1)
            @test ind2sub((4,3), k) == (i,j)
            @test sub2ind((1:4,1:3), i, j) == k
            @test ind2sub((1:4,1:3), k) == (i,j)
            @test sub2ind((0:3,3:5), i-1, j+2) == k
            @test ind2sub((0:3,3:5), k) == (i-1, j+2)
        end
        @testset "Delete when partial linear indexing is deprecated (#14770)" begin
            @test sub2ind((4,3), 7) == 7
            @test sub2ind((1:4,1:3), 7) == 7
            @test sub2ind((0:3,3:5), 7) == 8
        end
    end

    @testset "3-dimensional" begin
        l = 0
        for k = 1:2, j = 1:3, i = 1:4
            @test sub2ind((4,3,2), i, j, k) == (l+=1)
            @test ind2sub((4,3,2), l) == (i,j,k)
            @test sub2ind((1:4,1:3,1:2), i, j, k) == l
            @test ind2sub((1:4,1:3,1:2), l) == (i,j,k)
            @test sub2ind((0:3,3:5,-101:-100), i-1, j+2, k-102) == l
            @test ind2sub((0:3,3:5,-101:-100), l) == (i-1, j+2, k-102)
        end

        A = reshape(collect(1:9), (3,3))
        @test ind2sub(size(A), 6) == (3,2)
        @test sub2ind(size(A), 3, 2) == 6
        @test ind2sub(A, 6) == (3,2)
        @test sub2ind(A, 3, 2) == 6

        @testset "PR #9256" begin
            function pr9256()
                m = [1 2 3; 4 5 6; 7 8 9]
                ind2sub(m, 6)
            end
            @test pr9256() == (3,2)
        end
    end
end

# token type on which to dispatch testing methods in order to avoid potential
# name conflicts elsewhere in the base test suite
mutable struct TestAbstractArray end

## Tests for the abstract array interfaces with minimally defined array types

# A custom linear fast array type with 24 elements that doesn't rely upon Array storage
mutable struct T24Linear{T,N,dims} <: AbstractArray{T,N}
    v1::T;  v2::T;  v3::T;  v4::T;  v5::T;  v6::T;  v7::T;  v8::T
    v9::T;  v10::T; v11::T; v12::T; v13::T; v14::T; v15::T; v16::T
    v17::T; v18::T; v19::T; v20::T; v21::T; v22::T; v23::T; v24::T
    T24Linear{T,N,d}() where {T,N,d} =
        (prod(d) == 24 || throw(DimensionMismatch("T24Linear must have 24 elements")); new())
    function T24Linear{T,N,d}(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,
                              v13,v14,v15,v16,v17,v18,v19,v20,v21,v22,v23,v24) where {T,N,d}
        prod(d) == 24 || throw(DimensionMismatch("T24Linear must have 24 elements"))
        new(v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16,v17,v18,v19,v20,v21,v22,v23,v24)
    end
end

T24Linear(::Type{T}, dims::Int...) where T = T24Linear(T, dims)
T24Linear(::Type{T}, dims::NTuple{N,Int}) where {T,N} = T24Linear{T,N,dims}()

Base.convert{T,N  }(::Type{T24Linear     }, X::AbstractArray{T,N}) = convert(T24Linear{T,N}, X)
Base.convert{T,N,_}(::Type{T24Linear{T  }}, X::AbstractArray{_,N}) = convert(T24Linear{T,N}, X)
Base.convert{T,N  }(::Type{T24Linear{T,N}}, X::AbstractArray     ) = T24Linear{T,N,size(X)}(X...)

Base.size{T,N,dims}(::T24Linear{T,N,dims}) = dims
import Base: IndexLinear
Base.IndexStyle{A<:T24Linear}(::Type{A}) = IndexLinear()
Base.getindex(A::T24Linear, i::Int) = getfield(A, i)
Base.setindex!{T}(A::T24Linear{T}, v, i::Int) = setfield!(A, i, convert(T, v))

# A custom linear slow sparse-like array that relies upon Dict for its storage
struct TSlow{T,N} <: AbstractArray{T,N}
    data::Dict{NTuple{N,Int}, T}
    dims::NTuple{N,Int}
end
TSlow{T}(::Type{T}, dims::Int...) = TSlow(T, dims)
TSlow{T,N}(::Type{T}, dims::NTuple{N,Int}) = TSlow{T,N}(Dict{NTuple{N,Int}, T}(), dims)

Base.convert{T,N  }(::Type{TSlow{T,N}}, X::TSlow{T,N}) = X
Base.convert{T,N  }(::Type{TSlow     }, X::AbstractArray{T,N}) = convert(TSlow{T,N}, X)
Base.convert{T,N,_}(::Type{TSlow{T  }}, X::AbstractArray{_,N}) = convert(TSlow{T,N}, X)
Base.convert{T,N  }(::Type{TSlow{T,N}}, X::AbstractArray     ) = begin
    A = TSlow(T, size(X))
    for I in CartesianRange(size(X))
        A[I.I...] = X[I.I...]
    end
    A
end

Base.size(A::TSlow) = A.dims
Base.similar{T}(A::TSlow, ::Type{T}, dims::Dims) = TSlow(T, dims)
import Base: IndexCartesian
Base.IndexStyle{A<:TSlow}(::Type{A}) = IndexCartesian()
# Until #11242 is merged, we need to define each dimension independently
Base.getindex{T}(A::TSlow{T,0}) = get(A.data, (), zero(T))
Base.getindex{T}(A::TSlow{T,1}, i1::Int) = get(A.data, (i1,), zero(T))
Base.getindex{T}(A::TSlow{T,2}, i1::Int, i2::Int) = get(A.data, (i1,i2), zero(T))
Base.getindex{T}(A::TSlow{T,3}, i1::Int, i2::Int, i3::Int) =
    get(A.data, (i1,i2,i3), zero(T))
Base.getindex{T}(A::TSlow{T,4}, i1::Int, i2::Int, i3::Int, i4::Int) =
    get(A.data, (i1,i2,i3,i4), zero(T))
Base.getindex{T}(A::TSlow{T,5}, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) =
    get(A.data, (i1,i2,i3,i4,i5), zero(T))

Base.setindex!{T}(A::TSlow{T,0}, v) = (A.data[()] = v)
Base.setindex!{T}(A::TSlow{T,1}, v, i1::Int) = (A.data[(i1,)] = v)
Base.setindex!{T}(A::TSlow{T,2}, v, i1::Int, i2::Int) = (A.data[(i1,i2)] = v)
Base.setindex!{T}(A::TSlow{T,3}, v, i1::Int, i2::Int, i3::Int) =
    (A.data[(i1,i2,i3)] = v)
Base.setindex!{T}(A::TSlow{T,4}, v, i1::Int, i2::Int, i3::Int, i4::Int) =
    (A.data[(i1,i2,i3,i4)] = v)
Base.setindex!{T}(A::TSlow{T,5}, v, i1::Int, i2::Int, i3::Int, i4::Int, i5::Int) =
    (A.data[(i1,i2,i3,i4,i5)] = v)

const can_inline = Base.JLOptions().can_inline != 0
function test_scalar_indexing{T}(::Type{T}, shape, ::Type{TestAbstractArray})
    N = prod(shape)
    A = reshape(collect(1:N), shape)
    B = T(A)
    @test A == B
    # Test indexing up to 5 dimensions
    i=0
    for i5 = 1:size(B, 5)
        for i4 = 1:size(B, 4)
            for i3 = 1:size(B, 3)
                for i2 = 1:size(B, 2)
                    for i1 = 1:size(B, 1)
                        i += 1
                        @test A[i1,i2,i3,i4,i5] == B[i1,i2,i3,i4,i5] == i
                        @test A[i1,i2,i3,i4,i5] ==
                              Base.unsafe_getindex(B, i1, i2, i3, i4, i5) == i
                    end
                end
            end
        end
    end
    # Test linear indexing and partial linear indexing
    i=0
    for i1 = 1:length(B)
        i += 1
        @test A[i1] == B[i1] == i
    end
    i=0
    for i2 = 1:size(B, 2)
        for i1 = 1:size(B, 1)
            i += 1
            @test A[i1,i2] == B[i1,i2] == i
        end
    end
    @test A == B
    i=0
    for i3 = 1:size(B, 3)
        for i2 = 1:size(B, 2)
            for i1 = 1:size(B, 1)
                i += 1
                @test A[i1,i2,i3] == B[i1,i2,i3] == i
            end
        end
    end
    # Test zero-dimensional accesses
    @test A[] == B[] == A[1] == B[1] == 1
    # Test multidimensional scalar indexed assignment
    C = T(Int, shape)
    D1 = T(Int, shape)
    D2 = T(Int, shape)
    D3 = T(Int, shape)
    i=0
    for i5 = 1:size(B, 5)
        for i4 = 1:size(B, 4)
            for i3 = 1:size(B, 3)
                for i2 = 1:size(B, 2)
                    for i1 = 1:size(B, 1)
                        i += 1
                        C[i1,i2,i3,i4,i5] = i
                        # test general unsafe_setindex!
                        Base.unsafe_setindex!(D1, i, i1,i2,i3,i4,i5)
                        # test for dropping trailing dims
                        Base.unsafe_setindex!(D2, i, i1,i2,i3,i4,i5, 1, 1, 1)
                        # test for expanding index argument to appropriate dims
                        Base.unsafe_setindex!(D3, i, i1,i2,i3,i4)
                    end
                end
            end
        end
    end
    @test D1 == D2 == C == B == A
    @test D3[:, :, :, :, 1] == D2[:, :, :, :, 1]
    # Test linear indexing and partial linear indexing
    C = T(Int, shape)
    fill!(C, 0)
    @test C != B && C != A
    i=0
    for i1 = 1:length(C)
        i += 1
        C[i1] = i
    end
    @test C == B == A
    C = T(Int, shape)
    i=0
    C2 = reshape(C, Val{2})
    for i2 = 1:size(C2, 2)
        for i1 = 1:size(C2, 1)
            i += 1
            C2[i1,i2] = i
        end
    end
    @test C == B == A
    C = T(Int, shape)
    i=0
    C3 = reshape(C, Val{3})
    for i3 = 1:size(C3, 3)
        for i2 = 1:size(C3, 2)
            for i1 = 1:size(C3, 1)
                i += 1
                C3[i1,i2,i3] = i
            end
        end
    end
    @test C == B == A
    # Test zero-dimensional setindex
    A[] = 0; B[] = 0
    @test A[] == B[] == 0
    @test A == B
end

function test_vector_indexing{T}(::Type{T}, shape, ::Type{TestAbstractArray})
    @testset "test_vector_indexing{$(T)}" begin
        N = prod(shape)
        A = reshape(collect(1:N), shape)
        B = T(A)
        idxs = rand(1:N, 3, 3, 3)
        @test B[idxs] == A[idxs] == idxs
        @test B[vec(idxs)] == A[vec(idxs)] == vec(idxs)
        @test B[:] == A[:] == collect(1:N)
        @test B[1:end] == A[1:end] == collect(1:N)
        # @test B[:,:] == A[:,:] == reshape(1:N, shape[1], prod(shape[2:end])) # TODO: Re-enable after partial linear indexing deprecation
        # @test B[1:end,1:end] == A[1:end,1:end] == reshape(1:N, shape[1], prod(shape[2:end])) # TODO: Re-enable after partial linear indexing deprecation

        @testset "Test with containers that aren't Int[]" begin
            @test B[[]] == A[[]] == []
            @test B[convert(Array{Any}, idxs)] == A[convert(Array{Any}, idxs)] == idxs
        end

        idx1 = rand(1:size(A, 1), 3)
        idx2 = rand(1:size(A, 2), 4, 5)
        @testset "Test adding dimensions with matrices" begin
            @test B[idx1, idx2] == A[idx1, idx2] == reshape(A[idx1, vec(idx2)], 3, 4, 5) == reshape(B[idx1, vec(idx2)], 3, 4, 5)
            @test B[1, idx2] == A[1, idx2] == reshape(A[1, vec(idx2)], 4, 5) == reshape(B[1, vec(idx2)], 4, 5)
        end
            # test removing dimensions with 0-d arrays
        @testset "test removing dimensions with 0-d arrays" begin
            idx0 = reshape([rand(1:size(A, 1))])
            @test B[idx0, idx2] == A[idx0, idx2] == reshape(A[idx0[], vec(idx2)], 4, 5) == reshape(B[idx0[], vec(idx2)], 4, 5)
            # @test B[reshape([end]), reshape([end])] == A[reshape([end]), reshape([end])] == reshape([A[end,end]]) == reshape([B[end,end]]) # TODO: Re-enable after partial linear indexing deprecation
        end

        mask = bitrand(shape)
        @testset "test logical indexing" begin
            @test B[mask] == A[mask] == B[find(mask)] == A[find(mask)] == find(mask)
            @test B[vec(mask)] == A[vec(mask)] == find(mask)
            mask1 = bitrand(size(A, 1))
            mask2 = bitrand(size(A, 2))
            @test B[mask1, mask2] == A[mask1, mask2] == B[find(mask1), find(mask2)]
            @test B[mask1, 1] == A[mask1, 1] == find(mask1)
        end
    end
end

function test_primitives{T}(::Type{T}, shape, ::Type{TestAbstractArray})
    N = prod(shape)
    A = reshape(collect(1:N), shape)
    B = T(A)

    # last(a)
    @test last(B) == B[length(B)]

    # strides(a::AbstractArray)
    @inferred strides(B)
    strides_B = strides(B)
    for (i, _stride) in enumerate(collect(strides_B))
        @test _stride == stride(B, i)
    end

    # isassigned(a::AbstractArray, i::Int...)
    j = rand(1:length(B))
    @test isassigned(B, j) == true
    if T == T24Linear
        @test isassigned(B, length(B) + 1) == false
    end

    # reshape(a::AbstractArray, dims::Dims)
    @test_throws DimensionMismatch reshape(B, (0, 1))

    # copy!(dest::AbstractArray, src::AbstractArray)
    @test_throws BoundsError copy!(Array{Int}(10), [1:11...])

    # convert{T, N}(::Type{Array}, A::AbstractArray{T, N})
    X = [1:10...]
    Y = [1 2; 3 4]
    @test convert(Array, X) == X
    @test convert(Array, Y) == Y

    # convert{T}(::Type{Vector}, A::AbstractVector{T})
    @test convert(Vector, X) == X
    @test convert(Vector, view(X, 2:4)) == [2,3,4]
    @test_throws MethodError convert(Vector, Y)

    # convert{T}(::Type{Matrix}, A::AbstractMatrix{T})
    @test convert(Matrix, Y) == Y
    @test convert(Matrix, view(Y, 1:2, 1:2)) == Y
    @test_throws MethodError convert(Matrix, X)
end

mutable struct TestThrowNoGetindex{T} <: AbstractVector{T} end
@testset "ErrorException if getindex is not defined" begin
    Base.length(::TestThrowNoGetindex) = 2
    Base.size(::TestThrowNoGetindex) = (2,)
    @test_throws ErrorException isassigned(TestThrowNoGetindex{Float64}(), 1)
end

function test_in_bounds(::Type{TestAbstractArray})
    n = rand(2:5)
    sz = rand(2:5, n)
    len = prod(sz)
    A = zeros(sz...)
    for i in 1:len
        @test checkbounds(Bool, A, i) == true
    end
    @test checkbounds(Bool, A, len + 1) == false
end

mutable struct UnimplementedFastArray{T, N} <: AbstractArray{T, N} end
Base.IndexStyle(::UnimplementedFastArray) = Base.IndexLinear()

mutable struct UnimplementedSlowArray{T, N} <: AbstractArray{T, N} end
Base.IndexStyle(::UnimplementedSlowArray) = Base.IndexCartesian()

mutable struct UnimplementedArray{T, N} <: AbstractArray{T, N} end

function test_getindex_internals{T}(::Type{T}, shape, ::Type{TestAbstractArray})
    N = prod(shape)
    A = reshape(collect(1:N), shape)
    B = T(A)

    @test getindex(A) == 1
    @test getindex(B) == 1
    @test Base.unsafe_getindex(A) == 1
    @test Base.unsafe_getindex(B) == 1
end

function test_getindex_internals(::Type{TestAbstractArray})
    U = UnimplementedFastArray{Int, 2}()
    V = UnimplementedSlowArray{Int, 2}()
    @test_throws ErrorException getindex(U, 1)
    @test_throws ErrorException Base.unsafe_getindex(U, 1)
    @test_throws ErrorException getindex(V, 1, 1)
    @test_throws ErrorException Base.unsafe_getindex(V, 1, 1)
end

function test_setindex!_internals{T}(::Type{T}, shape, ::Type{TestAbstractArray})
    N = prod(shape)
    A = reshape(collect(1:N), shape)
    B = T(A)

    Base.unsafe_setindex!(B, 1)
    @test B[1] == 1
end

function test_setindex!_internals(::Type{TestAbstractArray})
    U = UnimplementedFastArray{Int, 2}()
    V = UnimplementedSlowArray{Int, 2}()
    @test_throws ErrorException setindex!(U, 0, 1)
    @test_throws ErrorException Base.unsafe_setindex!(U, 0, 1)
    @test_throws ErrorException setindex!(V, 0, 1, 1)
    @test_throws ErrorException Base.unsafe_setindex!(V, 0, 1, 1)
end

function test_get(::Type{TestAbstractArray})
    A = T24Linear([1:24...])
    B = TSlow([1:24...])

    @test get(A, (), 0) == Int[]
    @test get(B, (), 0) == TSlow(Int, 0)
end

function test_cat(::Type{TestAbstractArray})
    A = T24Linear([1:24...])
    b_int = reshape([1:27...], 3, 3, 3)
    b_float = reshape(Float64[1:27...], 3, 3, 3)
    b2hcat = Array{Float64}(3, 6, 3)
    b1 = reshape([1:9...], 3, 3)
    b2 = reshape([10:18...], 3, 3)
    b3 = reshape([19:27...], 3, 3)
    b2hcat[:, :, 1] = hcat(b1, b1)
    b2hcat[:, :, 2] = hcat(b2, b2)
    b2hcat[:, :, 3] = hcat(b3, b3)
    b3hcat = Array{Float64}(3, 9, 3)
    b3hcat[:, :, 1] = hcat(b1, b1, b1)
    b3hcat[:, :, 2] = hcat(b2, b2, b2)
    b3hcat[:, :, 3] = hcat(b3, b3, b3)
    B = TSlow(b_int)
    B1 = TSlow([1:24...])
    B2 = TSlow([1:25...])
    C1 = TSlow([1 2; 3 4])
    C2 = TSlow([1 2 3; 4 5 6])
    C3 = TSlow([1 2; 3 4; 5 6])
    D = [1:24...]
    i = rand(1:10)

    @test cat(i) == Any[]
    @test vcat() == Any[]
    @test hcat() == Any[]
    @test hcat(1, 1.0, 3, 3.0) == [1.0 1.0 3.0 3.0]
    @test_throws ArgumentError hcat(B1, B2)
    @test_throws ArgumentError vcat(C1, C2)

    @test vcat(B) == B
    @test hcat(B) == B
    @test Base.typed_hcat(Float64, B) == TSlow(b_float)
    @test Base.typed_hcat(Float64, B, B) == TSlow(b2hcat)
    @test Base.typed_hcat(Float64, B, B, B) == TSlow(b3hcat)

    @test vcat(B1, B2) == TSlow(vcat([1:24...], [1:25...]))
    @test hcat(C1, C2) == TSlow([1 2 1 2 3; 3 4 4 5 6])
    @test hcat(C1, C2, C1) == TSlow([1 2 1 2 3 1 2; 3 4 4 5 6 3 4])

    # hvcat
    for nbc in (1, 2, 3, 4, 5, 6)
        @test hvcat(nbc, 1:120...) ==
              transpose(reshape([1:120...], nbc, round(Int, 120 / nbc)))
    end

    @test_throws ArgumentError hvcat(7, 1:20...)
    @test_throws ArgumentError hvcat((2), C1, C3)
    @test_throws ArgumentError hvcat((1), C1, C2)
    @test_throws ArgumentError hvcat((1), C2, C3)

    tup = tuple(rand(1:10, i)...)
    @test hvcat(tup) == []

    # check for shape mismatch
    @test_throws ArgumentError hvcat((2, 2), 1, 2, 3, 4, 5)
    @test_throws ArgumentError Base.typed_hvcat(Int, (2, 2), 1, 2, 3, 4, 5)
    # check for # of columns mismatch b/w rows
    @test_throws ArgumentError hvcat((3, 2), 1, 2, 3, 4, 5, 6)
    @test_throws ArgumentError Base.typed_hvcat(Int, (3, 2), 1, 2, 3, 4, 5, 6)

    # 18395
    @test isa(Any["a" 5; 2//3 1.0][2,1], Rational{Int})

    # 13665, 19038
    @test @inferred(hcat([1.0 2.0], 3))::Array{Float64,2} == [1.0 2.0 3.0]
    @test @inferred(vcat([1.0, 2.0], 3))::Array{Float64,1} == [1.0, 2.0, 3.0]

    @test @inferred(vcat(["a"], "b"))::Vector{String} == ["a", "b"]
    @test @inferred(vcat((1,), (2.0,)))::Vector{Tuple{Real}} == [(1,), (2.0,)]
end

function test_ind2sub(::Type{TestAbstractArray})
    n = rand(2:5)
    dims = tuple(rand(1:5, n)...)
    len = prod(dims)
    A = reshape(collect(1:len), dims...)
    I = ind2sub(dims, [1:len...])
    for i in 1:len
        idx = [ I[j][i] for j in 1:n ]
        @test A[idx...] == A[i]
    end
end

# A custom linear slow array that insists upon Cartesian indexing
mutable struct TSlowNIndexes{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
end
Base.IndexStyle{A<:TSlowNIndexes}(::Type{A}) = Base.IndexCartesian()
Base.size(A::TSlowNIndexes) = size(A.data)
Base.getindex(A::TSlowNIndexes, index::Int...) = error("Must use $(ndims(A)) indexes")
Base.getindex{T}(A::TSlowNIndexes{T,2}, i::Int, j::Int) = A.data[i,j]

mutable struct GenericIterator{N} end
Base.start{N}(::GenericIterator{N}) = 1
Base.next{N}(::GenericIterator{N}, i) = (i, i + 1)
Base.done{N}(::GenericIterator{N}, i) = i > N ? true : false
Base.iteratorsize{N}(::Type{GenericIterator{N}}) = Base.SizeUnknown()

function test_map(::Type{TestAbstractArray})
    empty_pool = WorkerPool([myid()])
    pmap_fallback = (f, c...) -> pmap(empty_pool, f, c...)

    for mapf in [map, asyncmap, pmap_fallback]
        for typ in (Float16, Float32, Float64,
                    Int8, Int16, Int32, Int64, Int128,
                    UInt8, UInt16, UInt32, UInt64, UInt128),
            arg_typ in (Integer,
                        Signed,
                        Unsigned)
            X = typ[1:10...]
            _typ = typeof(arg_typ(one(typ)))
            @test mapf(arg_typ, X) == _typ[1:10...]
        end

        # generic map
        f(x) = x + 1
        I = GenericIterator{10}()
        @test mapf(f, I) == Any[2:11...]

        # AbstractArray map for 2 arg case
        f(x, y) = x + y
        B = Float64[1:10...]
        C = Float64[1:10...]
        @test mapf(f, convert(Vector{Int},B), C) == Float64[ 2 * i for i in 1:10 ]
        @test mapf(f, Int[], Float64[]) == Union{}[]
        # map with different result types
        let m = mapf(x->x+1, Number[1, 2.0])
            # FIXME why is this different for asyncmap?
            @test mapf !== map || isa(m, Vector{Real})
            @test m == Real[2, 3.0]
        end

        # AbstractArray map for N-arg case
        A = Array{Int}(10)
        f(x, y, z) = x + y + z
        D = Float64[1:10...]

        @test map!(f, A, B, C, D) == Int[ 3 * i for i in 1:10 ]
        @test mapf(f, B, C, D) == Float64[ 3 * i for i in 1:10 ]
        @test mapf(f, Int[], Int[], Complex{Int}[]) == Union{}[]
    end

    # In-place map
    A = Float64[1:10...]
    map!(x -> x*x, A, A)
    @test A == map(x -> x*x, Float64[1:10...])
    B = Float64[1:10...]
    Base.asyncmap!(x->x*x, B, B)
    @test A == B

    # Map to destination collection
    map!((x,y,z)->x*y*z, A, Float64[1:10...], Float64[1:10...], Float64[1:10...])
    @test A == map(x->x*x*x, Float64[1:10...])
    Base.asyncmap!((x,y,z)->x*y*z, B, Float64[1:10...], Float64[1:10...], Float64[1:10...])
    @test A == B
end

@testset "issue #15689, mapping an abstract type" begin
    @test isa(map(Set, Array[[1,2],[3,4]]), Vector{Set{Int}})
end

@testset "mapping over scalars and empty arguments:" begin
    @test map(sin, 1) === sin(1)
    @test map(()->1234) === 1234
end

function test_UInt_indexing(::Type{TestAbstractArray})
    A = [1:100...]
    _A = Expr(:quote, A)
    for i in 1:100
        _i8 = convert(UInt8, i)
        _i16 = convert(UInt16, i)
        _i32 = convert(UInt32, i)
        for _i in (_i8, _i16, _i32)
            @eval begin
                @test $_A[$_i] == $i
            end
        end
    end
end

# Issue 13315
function test_13315(::Type{TestAbstractArray})
    U = UInt(1):UInt(2)
    @test [U;[U;]] == [UInt(1), UInt(2), UInt(1), UInt(2)]
end

# checksquare
function test_checksquare()
    @test LinAlg.checksquare(zeros(2,2)) == 2
    @test LinAlg.checksquare(zeros(2,2),zeros(3,3)) == [2,3]
    @test_throws DimensionMismatch LinAlg.checksquare(zeros(2,3))
end

#----- run tests -------------------------------------------------------------#

@testset for T in (T24Linear, TSlow), shape in ((24,), (2, 12), (2,3,4), (1,2,3,4), (4,3,2,1))
    test_scalar_indexing(T, shape, TestAbstractArray)
    test_vector_indexing(T, shape, TestAbstractArray)
    test_primitives(T, shape, TestAbstractArray)
    test_getindex_internals(T, shape, TestAbstractArray)
    test_setindex!_internals(T, shape, TestAbstractArray)
end
test_in_bounds(TestAbstractArray)
test_getindex_internals(TestAbstractArray)
test_setindex!_internals(TestAbstractArray)
test_get(TestAbstractArray)
test_cat(TestAbstractArray)
test_ind2sub(TestAbstractArray)
test_map(TestAbstractArray)
test_UInt_indexing(TestAbstractArray)
test_13315(TestAbstractArray)
test_checksquare()

A = TSlowNIndexes(rand(2,2))
@test_throws ErrorException A[1]
@test A[1,1] == A.data[1]
@test first(A) == A.data[1]

@testset "#16381" begin
    @inferred size(rand(3,2,1), 2, 1)
    @inferred size(rand(3,2,1), 2, 1, 3)

    @test @inferred(indices(rand(3,2)))    == (1:3,1:2)
    @test @inferred(indices(rand(3,2,1)))  == (1:3,1:2,1:1)
    @test @inferred(indices(rand(3,2), 1)) == 1:3
    @test @inferred(indices(rand(3,2), 2)) == 1:2
    @test @inferred(indices(rand(3,2), 3)) == 1:1
end

@testset "#17088" begin
    n = 10
    M = rand(n, n)
    @testset "vector of vectors" begin
        v = [[M]; [M]] # using vcat
        @test size(v) == (2,)
        @test !issparse(v)
    end
    @testset "matrix of vectors" begin
        m1 = [[M] [M]] # using hcat
        m2 = [[M] [M];] # using hvcat
        @test m1 == m2
        @test size(m1) == (1,2)
        @test !issparse(m1)
        @test !issparse(m2)
    end
end

@testset "isinteger and isreal" begin
    @test all(isinteger, Diagonal(rand(1:5,5)))
    @test isreal(Diagonal(rand(5)))
end

@testset "unary ops" begin
    let A = Diagonal(rand(1:5,5))
        @test +(A) == A
        @test *(A) == A
    end
end

@testset "flipdim on empty" begin
    @test flipdim(Diagonal([]),1) == Diagonal([])
end

@testset "ndims and friends" begin
    @test ndims(Diagonal(rand(1:5,5))) == 2
    @test ndims(Diagonal{Float64}) == 2
    @test Base.elsize(Diagonal(rand(1:5,5))) == sizeof(Int)
end

@testset "Issue #17811" begin
    A17811 = Integer[]
    I = [abs(x) for x in A17811]
    @test isa(I, Array{Any,1})
    push!(I, 1)
    @test I == Any[1]
    @test isa(map(abs, A17811), Array{Any,1})
end

@testset "copymutable for itrs" begin
    @test Base.copymutable((1,2,3)) == [1,2,3]
end

@testset "sub2ind for empty tuple" begin
    @test sub2ind(()) == 1
end

@testset "to_shape" begin
    @test Base.to_shape(()) === ()
    @test Base.to_shape(1) === 1
end

@testset "issue #19267" begin
    @test ndims((1:3)[:]) == 1
    @test ndims((1:3)[:,:]) == 2
    @test ndims((1:3)[:,[1],:]) == 3
    @test ndims((1:3)[:,[1],:,[1]]) == 4
    @test ndims((1:3)[:,[1],1:1,:]) == 4
    @test ndims((1:3)[:,:,1:1,:]) == 4
    @test ndims((1:3)[:,:,1:1]) == 3
    @test ndims((1:3)[:,:,1:1,:,:,[1]]) == 6
end

@testset "dispatch loop introduced in #19305" begin
    @test [(1:2) zeros(2,2); ones(3,3)] == [[1,2] zeros(2,2); ones(3,3)] == [reshape([1,2],2,1) zeros(2,2); ones(3,3)]
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""
    HexapolarField(sys::AxisymmetricOpticalSystem; collimated = true, samples = 8, wavelength = 0.55, sourcepos = (0.0, 0.0, 3.0), sourceangle = 0.0, sourcenum = 0)

Distributes rays over the entrance pupil of the system in a hexapolar pattern.
"""
function HexapolarField(sys::AxisymmetricOpticalSystem{T}; kwargs...) where {T<:Real}
    return HexapolarField(semidiameter(sys) * 0.95, zero(SVector{3,T}); kwargs...)
end

"""
    HexapolarField(semidiameter, pupilpos; collimated = true, samples = 8, wavelength = 0.55, sourcepos = (0.0, 0.0, 3.0), sourceangle = 0.0, sourcenum = 0)

Distributes rays over a circular pupil with half-diameter defined by `semidiameter`, centred at `pupilpos` in a hexapolar pattern.
`samples` is the number of rings in the hexapolar pattern, so the number of rays in total is `samples * (samples + 1) / 2) * 6 + 1`.
"""
function HexapolarField(semidiameter::T, pupilpos::SVector{3,T}; collimated::Bool = true, samples::Int = 8, wavelength::T = T(0.55), sourcepos::SVector{3,T} = SVector{3,T}(0.0, 0.0, 3.0), sourceangle::T = zero(T), sourcenum::Int = 0) where {T<:Real}
    if sourceangle != zero(T)
        sourcepos = norm(sourcepos) .* SVector{3,T}(0.0, tan(sourceangle), 1.0) + pupilpos
    end
    pointsonpupil = Vector{SVector{3,T}}(undef, 0)
    for n in 0:samples
        if n == 0
            push!(pointsonpupil, pupilpos)
        else
            ρ = n / samples
            ninring = n * 6
            for i in 1:ninring
                ϕ = (i / ninring) * 2π
                u = cos(ϕ) * semidiameter
                v = sin(ϕ) * semidiameter
                push!(pointsonpupil, SVector{3,T}(ρ * u, ρ * v, 0) + pupilpos)
            end
        end
    end
    N = length(pointsonpupil)
    rays = Vector{OpticalRay{T,3}}(undef, N)
    if collimated
        direction = pupilpos - sourcepos
        for i in 1:N
            rays[i] = OpticalRay(pointsonpupil[i] - direction, direction, one(T), wavelength, sourcenum = sourcenum)
        end
        raygenerator = RayListSource(rays)
    else
        for i in 1:N
            rays[i] = OpticalRay(sourcepos, pointsonpupil[i] - sourcepos, one(T), wavelength, sourcenum = sourcenum)
        end
        raygenerator = RayListSource(rays)
    end
    return raygenerator
end
export HexapolarField

"""
    GridField(sys::AxisymmetricOpticalSystem; collimated = true, samples = 20, wavelength = 0.55, sourcepos = (0.0, 0.0, 3.0), sourceangle = 0.0, sourcenum = 0)

Distributes rays over the entrance pupil of the system in a rectangular grid pattern.
"""
function GridField(sys::AxisymmetricOpticalSystem{T}; kwargs...) where {T<:Real}
    return GridField(semidiameter(sys) * 0.95, zero(SVector{3,T}); kwargs...)
end

"""
    GridField(semidiameter, pupilpos; collimated = true, samples = 20, wavelength = 0.55, sourcepos = (0.0, 0.0, 3.0), sourceangle = 0.0, sourcenum = 0)

Distributes rays over a circular pupil with half-diameter defined by `semidiameter`, centred at `pupilpos` in a rectangular grid pattern.
`samples` is the number of rays on each side of the grid, so there are `samples×samples` rays in total.
"""
function GridField(semidiameter::T, pupilpos::SVector{3,T}; collimated::Bool = true, samples::Int = 20, wavelength::T = T(0.55), sourcepos::SVector{3,T} = SVector{3,T}(0.0, 0.0, 3.0), sourceangle::T = zero(T), sourcenum::Int = 0) where {T<:Real}
    if sourceangle != zero(T)
        sourcepos = norm(sourcepos) .* SVector{3,T}(0.0, tan(sourceangle), 1.0) + pupilpos
    end
    pointsonpupil = Vector{SVector{3,T}}(undef, 0)
    for x in 0:samples
        for y in 0:samples
            u = (2.0 * (x / samples) - one(T)) * semidiameter
            v = (2.0 * (y / samples) - one(T)) * semidiameter
            if (u == zero(T) && v == zero(T)) || sqrt(u^2 + v^2) < semidiameter
                push!(pointsonpupil, SVector{3,T}(u, v, 0) + pupilpos)
            end
        end
    end

    N = length(pointsonpupil)
    rays = Vector{OpticalRay{T,3}}(undef, N)
    if collimated
        direction = pupilpos - sourcepos
        for i in 1:N
            rays[i] = OpticalRay(pointsonpupil[i] - direction, direction, one(T), wavelength, sourcenum = sourcenum)
        end
        raygenerator = RayListSource(rays)
    else
        for i in 1:N
            rays[i] = OpticalRay(sourcepos, pointsonpupil[i] - sourcepos, one(T), wavelength, sourcenum = sourcenum)
        end
        raygenerator = RayListSource(rays)
    end
    return raygenerator
end
export GridField

ext: jl


text:
using ArgParse

ArgParse.parse_item(::Type{Symbol}, s::AbstractString) = Symbol(s)

# Taken from the docstring of `printstyled()`
const VALID_COLORS = Set([
    :normal, :default, :bold, :black, :blink, :blue, :cyan, :green, :hidden, :light_black,
    :light_blue, :light_cyan, :light_green, :light_magenta, :light_red, :light_yellow,
    :magenta, :nothing, :red, :reverse, :underline, :white, :yellow
    ])

is_valid_color(color::String) = Symbol(color) in VALID_COLORS

const ARG_DEFAULTS = Dict(
    :include_only => nothing,
    :exclude => nothing,
    :include_only_tags => Symbol[],
    :exclude_tags => Symbol[],
    :max_fails => 0,
    :capture_output => false,
    :verbosity => 1,
    :dont_add_runtests_path => false,
    :test_file_postfix => ".test.jl",
    :color_pass => :green,
    :color_fail => :red,
    :color_error => :yellow,
    :color_broken => :green,
    :color_return => :blue,
    )

function get_default(options, key)
    get(options, key, ARG_DEFAULTS[key])
end

# Used to escape underscores in symbols so that Markdown doesn't turn them into italic.
escape_symbol(s::Symbol) = replace(string(s), "_" => "\\_")

"""
For every option, the corresponding command-line argument names are given in parentheses.
If supplied via the `options` keyword argument of [`runtests()`](@ref),
their type must be as given or `convert()`-able to it.

**`:include_only`**`:: Union{Nothing, Regex}` (`--include-only`, `-i`):
takes a regular expression; tests with full names that do not match it will not be executed.

**`:exclude`**`:: Union{Nothing, Regex}` (`--exclude`, `-e`):
takes a regular expression; tests with full names that match it will not be executed.

!!! note

    The full name refers to the name of the testcase and all its parent groups,
    separated by slashes (`/`).
    For example, the full name of the testcase "foo" in the group "bar", which, in turn,
    is defined in the group "baz", will be "baz/bar/foo".

**`:verbosity`**`:: Int` (`--verbosity`, `-v`):
`0`, `1` or `2`, defines the amount of output that will be shown. `1` is the default.

**`:include_only_tags`**`:: Array{Symbol, 1}` (`--include-only-tags`, `-t`):
include only tests with any of the specified tags.
You can pass several tags to this option, separated by spaces.

**`:exclude_tags`**`:: Array{Symbol, 1}` (`--exclude-tags`, `-n`):
exclude tests with any of the specified tags.
You can pass several tags to this option, separated by spaces.

!!! note

    Programmatically (using `options` keyword of [`runtests`](@ref))
    tags must be passed as symbols.

!!! note

    If `-` is one of the arguments to `include_only_tags` or `exclude_tags`,
    it is equivalent to not passing anything for these options at all
    (all other values passed to them will be ignored).
    This can be helpful when overriding the options passed programmatically to
    [`runtests()`](@ref).

**`:max_fails`**`:: Int` (`--max-fails`):
stop after the given amount of failed testcases
(a testcase is considered failed if at least one test in it failed,
or an unhandled exception was thrown).

**`:test_file_postifx`**`:: String` (`--test-file-postfix`):
postfix of the files which will be picked up by the automatic testcase discovery.

The flag arguments are processed in a special way if a non-default value is passed via
the `options` keyword argument of [`runtests()`](@ref).
In this case, the command-line option is replaced by its inverse (given second in the parentheses),
with the corresponding change in the default.

**`:capture_output`**`:: Bool` (`--capture-output`/`--dont-capture-output`):
capture all the output from testcases
and only show the output of the failed ones at the end of the test run.

**`:dont_add_runtests_path`**`:: Bool` (`--dont-add-runtests-path`/`--add-runtests-path`):
do not push the test root path into `LOAD_PATH` before including test files

**`:color_(pass/fail/error/broken/return)`**`:: Symbol`
(`--color-(pass/fail/error/broken/return)`):
custom colors to use in the report. Supported colors: $(join(escape_symbol.(VALID_COLORS), ", ")).

!!! note

    Programmatically (using `options` keyword of [`runtests`](@ref))
    colors must be passed as symbols.
"""
function build_parser(options)
    s = ArgParseSettings(; autofix_names=true)

    default(key) = get_default(options, key)

    # We have to insert the defaults passed programmatically in the argument table,
    # since there is currently no way in ArgParse to determine whether an argument
    # was actually passed by the user or not (which may change with PR#47).

    add_arg_table(
        s,

        ["--include-only", "-i"],
        Dict(
            :help => "include only tests (by path/name)",
            :metavar => "REGEX",
            :default => default(:include_only)),

        ["--exclude", "-e"],
        Dict(
            :help => "exclude tests (by path/name)",
            :metavar => "REGEX",
            :default => default(:exclude)),

        ["--include-only-tags", "-t"],
        Dict(
            :help => "include only tests (by tag)",
            :metavar => "TAGS",
            :arg_type => String,
            :nargs => '+',
            :default => String.(default(:include_only_tags))),

        ["--exclude-tags", "-n"],
        Dict(
            :help => "exclude tests (by tag)",
            :metavar => "TAGS",
            :arg_type => String,
            :nargs => '+',
            :default => String.(default(:exclude_tags))),

        "--max-fails",
        Dict(
            :help => "stop after a certain amount of failed testcases",
            :metavar => "NUM",
            :arg_type => Int,
            :default => default(:max_fails)),

        ["--verbosity", "-v"],
        Dict(
            :help => "the output verbosity (0-2)",
            :arg_type => Int,
            :default => default(:verbosity)),

        "--test-file-postfix",
        Dict(
            :help =>
                "postfix of the files which will be picked up " *
                "by the automatic testcase discovery.",
            :metavar => "STR",
            :arg_type => String,
            :default => default(:test_file_postfix)),
        )

    # The flag arguments require a special treatment.
    # If a non-default value was made default programmatically,
    # the corresponding flag is useless, and instead we need a flag for the negation of said value.

    if default(:capture_output)
        add_arg_table(
            s,
            "--dont-capture-output",
            Dict(
                :help =>
                    "don't capture testcase output and display it during test run",
                :nargs => 0))
    else
        add_arg_table(
            s,
            "--capture-output",
            Dict(
                :help =>
                    "capture testcase output and display only the output from failed testcases " *
                    "after all the testcases are finished",
                :nargs => 0))
    end

    if default(:dont_add_runtests_path)
        add_arg_table(
            s,
            "--add-runtests-path",
            Dict(
                :help =>
                    "push the test root path into `LOAD_PATH` " *
                    "before including test files.",
                :nargs => 0))
    else
        add_arg_table(
            s,
            "--dont-add-runtests-path",
            Dict(
                :help =>
                    "do not push the test root path into `LOAD_PATH` " *
                    "before including test files.",
                :nargs => 0))
    end

    add_arg_group(s, "Color scheme customization (supported colors: $(join(VALID_COLORS, ", ")))")

    add_arg_table(
        s,

        "--color-pass",
        Dict(
            :help => "The color used to mark passed tests",
            :metavar => "COLOR",
            :arg_type => String,
            :default => String(default(:color_pass)),
            :range_tester => is_valid_color),

        "--color-fail",
        Dict(
            :help => "The color used to mark failed tests",
            :metavar => "COLOR",
            :arg_type => String,
            :default => String(default(:color_fail)),
            :range_tester => is_valid_color),

        "--color-error",
        Dict(
            :help => "The color used to mark non-fail errors in tests",
            :metavar => "COLOR",
            :arg_type => String,
            :default => String(default(:color_error)),
            :range_tester => is_valid_color),

        "--color-broken",
        Dict(
            :help => "The color used to mark known broken tests",
            :metavar => "COLOR",
            :arg_type => String,
            :default => String(default(:color_broken)),
            :range_tester => is_valid_color),

        "--color-return",
        Dict(
            :help => "The color used to mark custom return values in tests",
            :metavar => "COLOR",
            :arg_type => String,
            :default => String(default(:color_return)),
            :range_tester => is_valid_color),
        )

    s
end

maybe_regex(::Nothing) = nothing
maybe_regex(s::String) = Regex(s)

struct ReportColorScheme

    color_pass :: Symbol
    color_fail :: Symbol
    color_error :: Symbol
    color_broken :: Symbol
    color_return :: Symbol

    function ReportColorScheme(run_options)
        new(
            Symbol(run_options[:color_pass]),
            Symbol(run_options[:color_fail]),
            Symbol(run_options[:color_error]),
            Symbol(run_options[:color_broken]),
            Symbol(run_options[:color_return]),
            )
    end
end

function normalize_options(run_options)

    run_options = deepcopy(run_options)

    run_options[:include_only] = maybe_regex(run_options[:include_only])
    run_options[:exclude] = maybe_regex(run_options[:exclude])

    # Collect tags as strings in order to make autogenerated ArgParse help more clear
    run_options[:include_only_tags] = Symbol.(run_options[:include_only_tags])
    run_options[:exclude_tags] = Symbol.(run_options[:exclude_tags])

    if :- in run_options[:include_only_tags]
        run_options[:include_only_tags] = Symbol[]
    end

    if :- in run_options[:exclude_tags]
        run_options[:exclude_tags] = Symbol[]
    end

    # Remove the possible negated flags

    if haskey(run_options, :dont_capture_output)
        run_options[:capture_output] = !run_options[:dont_capture_output]
        delete!(run_options, :dont_capture_output)
    end

    if haskey(run_options, :add_runtests_path)
        run_options[:dont_add_runtests_path] = !run_options[:add_runtests_path]
        delete!(run_options, :add_runtests_path)
    end

    run_options[:report_color_scheme] = ReportColorScheme(run_options)
    for name in fieldnames(ReportColorScheme)
        delete!(run_options, name)
    end

    run_options
end

function build_run_options_from_commandline(args, options)
    s = build_parser(options)
    normalize_options(parse_args(args, s; as_symbols=true))
end

function build_run_options(; args=nothing, options=nothing)
    if options === nothing
        options = Dict{Symbol, Any}()
    end

    if args === nothing
        args = String[]
    end

    for (key, val) in options
        if !haskey(ARG_DEFAULTS, key)
            error("Unknown option: $key")
        end
        options[key] = convert(typeof(ARG_DEFAULTS[key]), val)
    end

    build_run_options_from_commandline(args, options)
end


text:
# mbe_demo.jl
#
# This script contains the work towards a "reduced effort" control law for an
# underactuated vehicle. In includes the initial simulation, the Monte-Carlo
# runs, and TODO.
#
# To run this file, start Julia, navigate to the directory containing this
# file, and just include it in your main workspace, like so:
#
# ```
# $ julia
# julia> cd(Pkg.dir("Eratosthenes") * "/examples/5_mbe")
# julia> include("mbe_demo.jl")
# ```

# Include our custom models.
include("models.jl")

# We work inside a module because it's more flexible than the main workspace.
module MBEDemo

#############
# Resources #
#############

# Use/import the modules we'll need.
using Eratosthenes # The simulation engine (setup, simulate, mc)
using ..MBEModels  # The MBEModels module included above (and "above" this module)
import HDF5        # For loading the log file
using Plots        # For plotting our results
using LinearAlgebra

# Create a directory for the log file.
if !isdir("out"); mkdir("out"); end

##############
# Simulation #
##############

# Specify the name of the YAML file that defines our scenario.
yaml_file = joinpath(@__DIR__, "underactuated.yaml")

# Run the scenario defined by our YAML file. The returned scenario will
# have all of the final states.
scenario = simulate(yaml_file, @__MODULE__)

#########
# Plots #
#########

# Read some data from the logs. This do block takes care of
# error-handling, so we can just read whatever we want and return it.
t, q_BI, ω_BI_B = HDF5.h5open(scenario.sim.log_file, "r") do logs
    (read(logs, "/smallsat/body/state/time"),
     read(logs, "/smallsat/body/state/data/q_BI"),
     read(logs, "/smallsat/body/state/data/ω_BI_B"))
end

# Scalars are logged as 1-by-n. Convert to just n for the plot.
t = dropdims(t, dims=1)

# Choose a friendly plotting package.
pyplot()
# plotlyjs()

# Define each plot that we'll need.
display(plot(t, transpose(q_BI),
             label  = ["q1" "q2" "q3" "q4"],
             xlabel = "Time (s)",
             ylabel = "Attitude Quaternion"))
display(plot(t, 180/π * transpose(ω_BI_B),
             label  = ["ω1" "ω2" "ω3"],
             xlabel = "Time (s)",
             ylabel = "Rotation Rate (deg/s)"))

# The plot looks fine, so we can now try some Monte-Carlo runs.
# See mbe_demo_mc.jl.

end # module


text:
#!/usr/bin/env julia

"""
    nanremove(arr::Array{T,2}) where T<:Real

For a real 2D array with columns of e.g. time, stress and strain, or frequency, G' and G'',
find any row which contains NaN values and remove that row completely.
"""
function nanremove(arr::Array{T,2}) where T<:Real
    # get 1D array which is NaN for any corresponding row with NaNs in
    rowsums = sum(arr, dims=2)

    # get list of NaNs / not NaNs (true for not NaN, false for NaN)
    notnanlist = vec(.!isnan.(rowsums))

    if false in notnanlist
        # raise warning
        @warn "Please note that NaN data rows are not included in resultant data struct."
    end

    # extract only those rows whose sum did not evaluate to NaN
    arr[notnanlist, :]
end

function look_for_csv_col(var, cols, header_cells, default_options)
	col = nothing

	if var in keys(cols)
		if typeof(cols[var])<:String
			col = findfirst(e -> e == cols[var], header_cells)
		elseif typeof(cols[var])<:Integer
			col = cols[var]
		end

	else
		col = findfirst(e -> lowercase(e) in default_options, header_cells)
	end
	return(col)
end

"""
    importcsv(filepath::String; delimiter=',', header=false, comment = "Imported from csv file", savelog = true, kwargs...)

Load data from a CSV file (two/three columns, comma separated by default but
delimiter can be specified in the `delimiter` keyword argument).

The function can be used to construct either a RheoTimeData instance or a RheoFreqData
instance, depending on the parameters provided or column names;the function detects whether time or frequency has been included
and proceeds accordingly. For oscillatory data, all three columns (Gp, Gpp, Frequency)
must be provided. For regular viscoelastic data only time, or time-stress, or time-strain or
time-stress-strain data can be provided.

Column can be specified in the csv file by providing their values or header name (case insensitive) as keyword parameters.

`importcsv("filename.csv", time=1, strain=2, stress=3)` loads a `RheoTimeData` from a csv file with time in the first column, strain in the second, and stress in the third.

`importcsv("filename.csv", omega="Frequency", Gp="Storage", Gpp="Loss")` would detect the column numbers by reading the headers on the csv file.

If no information is provided, it would try to detect columns from header names, expecting standard names. Some of the recognised keywords are `time`, `stress`, `strain`, `frequency`, `storage modulus`, `loss modulus`.

If the csv file contains headers, but it is prefered to indicate columns by their numbers rather than header strings, the keyword `header` must be set to `true`.
"""
function importcsv(filepath::String; delimiter = ',', header = false, comment = "Imported from csv file", savelog = true, kwargs...)
# Convert parameter names to standard values.
	cols=symbol_to_unicode(kwargs.data)

	use_header = header
	if length(cols)==0
		use_header = true
	elseif any([typeof(e)<:String for e in cols])
		use_header = true
	end

	# read data from file
    if use_header
		dataraw,h = readdlm(filepath, delimiter, header=use_header)
        header_cells = [strip(string(e)) for e in [h...]]
	else
		dataraw = readdlm(filepath, delimiter)
        header_cells=[]
	end

    # remove and rows with NaN values in any of the columns
    data = nanremove(dataraw)

	# is it time data?
	t_col = look_for_csv_col(:t, cols, header_cells, ("t", "time"))

	if !isnothing(t_col)
		σ_col = look_for_csv_col(:σ, cols, header_cells, ("stress", "sigma", "σ"))
        ϵ_col = look_for_csv_col(:ϵ, cols, header_cells, ("strain", "epsilon", "ϵ"))

    #  @assert (!isnothing(t_col) && isnothing(ω_col)) || (isnothing(t_col) && !isnothing(ω_col)) "Data must contain either \"time\" or \"frequency\" "
    #  @assert endswith(lowercase(filepath), ".csv") "filepath must point to a .csv file."

        log = if savelog
                info = (comment=comment, folder=pwd(), stats=(t_min=data[1,t_col],t_max=data[end,t_col], n_sample=size(data[:,t_col])))
                RheoLogItem( (type=:source, funct=:importcsv, params=(filepath=filepath,), keywords=(t=t_col, σ=σ_col, ϵ=ϵ_col, header=header)), info )
              else
                nothing
              end

        # generate RheologyData struct and output
        if !isnothing(ϵ_col) && !isnothing(σ_col)
            return RheoTimeData(t = data[:, t_col], ϵ = data[:, ϵ_col], σ = data[:, σ_col], log = log)
        elseif !isnothing(ϵ_col) && isnothing(σ_col)
            return RheoTimeData(t = data[:, t_col], ϵ = data[:, ϵ_col], log = log)
        elseif !isnothing(σ_col) && isnothing(ϵ_col)
            return RheoTimeData(t = data[:, t_col], σ = data[:, σ_col], log = log)
        end
    end

	ω_col = look_for_csv_col(:ω, cols, header_cells, ("ω", "omega", "frequency", "freq"))

    if !isnothing(ω_col)

        # check colnames length is correct
        #@assert isnothing(σ_col) || isnothing(ϵ_col) "Stress and strain not allowed for frequency data"
        #@assert !isnothing(Gp_col) && !isnothing(Gpp_col) "\"Gp\" and \"Gpp\" are required."
        Gp_col = look_for_csv_col(:Gp, cols, header_cells, ("Gp", "G'", "storage modulus"))
        Gpp_col = look_for_csv_col(:Gpp, cols, header_cells, ("Gpp", "G''", "loss modulus"))

        # Add error message for missing loss moduli / error loading

        log = if savelog
                info = (comment=comment, folder=pwd(), stats=(ω_min=data[1,ω_col],ω_max=data[end,ω_col], n_sample=size(data[:,ω_col])))
                RheoLogItem( (type=:source, funct=:importcsv, params=(filepath=filepath,), keywords=(ω=ω_col, Gp=Gp_col, Gpp=Gpp_col, header=header)), info )
              else
                nothing
              end

        return RheoFreqData(ω = data[:,ω_col], Gp = data[:,Gp_col], Gpp = data[:,Gpp_col], log = log)

    end
end

"""
    exportcsv(self::Union{RheoTimeData, RheoFreqData}, filedir::String; delimiter=',', colorder=nothing)

Export `RheoTimeData` or `RheoFreqData` type to csv format. May be useful for plotting/analysis in other software.
By default, full time data will be exported with columns ordered as (t, σ, ϵ). Partial time data will be ordered
as either (t, σ) or (t, ϵ). Full frequency data will be ordered as (ω, Gp, Gpp). The order of columns can be customised
by passing a NamedTuple to the `colorder` arguments. For example (σ = 1, t = 3, ϵ = 2) would export the columns in the
order (σ, ϵ, t). As with `importcsv`, the delimiter can be set by keyword argument.
"""
function exportcsv(self::Union{RheoTimeData, RheoFreqData}, filedir::String; delimiter=',', colorder=nothing)

    #@assert endswith(lowercase(filedir), ".csv") "filedir must be a .csv file."

    # if ordering of columns not provided, get default ordering depending on data contained in struct
    if isnothing(colorder)
        if typeof(self)==RheoTimeData
            datacontained = rheotimedatatype(self)
            if datacontained == stress_only
                colorder = (t=1, σ=2)
            elseif datacontained == strain_only
                colorder = (t=1, ϵ=2)
            elseif datacontained == strain_and_stress
                colorder = (t=1, σ=2, ϵ=3)
            end

        elseif typeof(self)==RheoFreqData
            # invalid_freq_data=-1 freq_only=0 with_modulus=1
            datacontained = rheofreqdatatype(self)
            colorder = (ω = 1, Gp = 2, Gpp = 3)
        end
    end

    cols = collect(keys(colorder))
    sortedcols = sort(cols, by=i->getfield(colorder, i))
    dataraw = getfield.((self,), sortedcols)
    dataout = hcat(dataraw...)
    open(filedir, "w") do io
        writedlm(io, dataout, delimiter)
    end
end

# """
#     savedata(self::Union{RheoTimeData, RheoFreqData}, filedir::String)

# Save RheoTimeData or RheoFreqData object using JLD2 format for re-use in
# a later Julia session.
# """
# function savedata(self::Union{RheoTimeData, RheoFreqData}, filedir::String)

#     @assert endswith(lowercase(filedir), ".jld2") "filedir must be a .jld2 file."

#     @save filedir self

# end

# """
#     loaddata(filedir::String)

# Loads RheoTimeData or RheoFreqData object from a jld2 file.
# """
# function loaddata(filedir::String)

#     @assert endswith(lowercase(filedir), ".jld2") "filedir must be a .jld2 file."

#     @load filedir self

#     return self

# end

# """
#     savemodel(self::RheologyModel, filedir::String; ext = ".jld2")

# Save RheologyModel object using JLD2 format reuse in a later Julia session.
# """
# function savemodel(self::RheologyModel, filedir::String; ext = ".jld2")

#     fulldir = string(filedir, ext)

#     @save fulldir self

# end

# """
#     loaddata(filedir::String)

# Loads RheologyModel from a JLD2 file.
# """
# function loadmodel(filedir::String)

#     @load filedir self

#     return self

# end


text:
using SubpixelRegistration
using FFTW
using Test

## Testing the registration on a small 4d array
@testset "Testing the registration on a small 4d array" begin
    test4d = zeros(40,40,20,2)

    test4d[10:20,10:20,2:10,1] .= 1
    test4d[5:15,15:25,5:13,2] .= 1

    ## No subpixel resolution
    dftResults = stackDftReg(test4d,ufac=1)

    @test length(dftResults) == size(test4d)[4]
    @test dftResults[2]["shift"] == [5,-5,-3]

    ## Factor 2
    dftResults = stackDftReg(test4d,ufac=2)

    @test length(dftResults) == size(test4d)[4]
    @test dftResults[2]["shift"] == [5,-5,-3]

    # More refined (needs dftups)
    dftResults = stackDftReg(test4d,ufac=3)

    @test length(dftResults) == size(test4d)[4]
    @test dftResults[2]["shift"] == [5,-5,-3]

    ## Testing the translation
    back4d = alignFromDict(test4d,dftResults)
    @test all(round.(back4d[:,:,:,2]-back4d[:,:,:,1]).<eps())
end


text:
function listen_localhost()
    @async HTTP.listen(HTTP.Sockets.localhost, 8080) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do ws
                while !eof(ws)
                    data = readavailable(ws)
                    write(ws, data)
                end
            end
        end
    end
end

function do_nothing_localhost()
    @async HTTP.listen(HTTP.Sockets.localhost, 8081) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do ws
                while !eof(ws)
                    data = readavailable(ws)
                end
            end
        end
    end
end

@testset "timers" begin
    listen_localhost()
    HTTP.WebSockets.open("ws://127.0.0.1:8080") do ws
        # timeout
        ch = GraphQLClient.async_reader_with_timeout(ws, 0.1)
        @test take!(ch) == :timeout

        ch = GraphQLClient.async_reader_with_timeout(ws, 5)
        write(ws, "Data")
        @test String(take!(ch)) == "Data"

        # stopfn
        stop = Ref(false)
        stopfn() = stop[]
        ch = GraphQLClient.async_reader_with_stopfn(ws, stopfn, 0.5)
        @test !isready(ch)
        stop[] = true
        sleep(1.0)
        @test take!(ch) == :stopfn
        stop[] = false
        ch = GraphQLClient.async_reader_with_stopfn(ws, stopfn, 0.5)
        write(ws, "Data")
        @test String(take!(ch)) == "Data"

        # readfromwebsocket - no timeout or stopfn
        write(ws, "Data")
        @test String(GraphQLClient.readfromwebsocket(ws, nothing, 0)) == "Data"

        # readfromwebsocket - timeout
        @test GraphQLClient.readfromwebsocket(ws, nothing, 0.1) == :timeout

        # readfromwebsocket - stopfn
        count = Ref(0)
        function stopfn2()
            if count[] > 0
                return true
            else
                count[] += 1
                return false
            end
        end
        @test GraphQLClient.readfromwebsocket(ws, stopfn2, 0.1) == :stopfn
    end
end

function send_error_localhost(message, port)
    @async HTTP.listen(HTTP.Sockets.localhost, port) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do ws
                while !eof(ws)
                    data = readavailable(ws)
                    isempty(data) && continue
                    query = JSON3.read(data)
                    error_payload = """
                    {   "id": "$(query["id"])",
                        "type": "sub",
                        "payload": {
                            "errors" : [
                                {
                                    "message": "$message",
                                    "locations": [{
                                        "line": 1,
                                        "column": 10
                                    }]
                                }
                            ]
                        }
                    }
                    """
                    write(ws, error_payload)
                end
            end
        end
    end
end

function send_data_localhost(sub_name, port)
    @async HTTP.listen(HTTP.Sockets.localhost, port) do http
        if HTTP.WebSockets.is_upgrade(http.message)
            HTTP.WebSockets.upgrade(http) do ws
                while !eof(ws)
                    data = readavailable(ws)
                    isempty(data) && continue
                    query = JSON3.read(data)
                    data_payload = """
                    {   "id": "$(query["id"])",
                        "type": "sub",
                        "payload": {
                            "data": {
                                "$sub_name": {
                                    "field": 1
                                }
                            }
                        }
                    }
                    """
                    write(ws, data_payload)
                end
            end
        end
    end
end

@testset "subscriptions.jl" begin
    # Set up client
    client = Client("http", "ws://127.0.0.1:8081"; introspect=false)
    client.introspection_complete = true
    push!(client.subscriptions, "MySub")

    @test_throws GraphQLClient.GraphQLError open_subscription(()->(), client, "Doesn't Exist")

    do_nothing_localhost()

    # init_func
    val = Ref(1)
    init_arg = 2
    initfn(arg) = val[] = arg
    open_subscription(
        (val)->(sleep(0.5); false),
        client,
        "MySub",
        initfn=() -> initfn(init_arg),
        output_fields="field",
        subtimeout=0.1)
    @test val[] == 2

    # Error response - not throwing
    port = 8093
    send_error_localhost("This failed", port)
    client = Client("http", "ws://127.0.0.1:$port"; introspect=false)
    client.introspection_complete = true
    push!(client.subscriptions, "MySub")
    results = []
    open_subscription(
        (val)->(push!(results, val); true),
        client,
        "MySub",
        output_fields="field")
    @test length(results) == 1
    @test results[1] isa GraphQLClient.GQLResponse
    @test isnothing(results[1].data)
    @test !isnothing(results[1].errors)
    @test !isempty(results[1].errors)
    @test results[1].errors[1].message == "This failed"

    # Error response, thrown
    @test_throws GraphQLClient.GraphQLError open_subscription(
        (val)->(push!(results, val); true),
        client,
        "MySub",
        output_fields="field",
        throw_on_execution_error=true)

    # Actual data returned
    port = 8096
    send_data_localhost("MySub", port)
    client = Client("http", "ws://127.0.0.1:$port"; introspect=false)
    client.introspection_complete = true
    push!(client.subscriptions, "MySub")
    results = []
    open_subscription(
        (val)->(push!(results, val); true),
        client,
        "MySub",
        output_fields="field")
    @test results[1] isa GraphQLClient.GQLResponse
    @test isnothing(results[1].errors)
    @test !isnothing(results[1].data) # No point testing content as we've coded it in the test function

    # Test struct types
    struct Response
        field::Int
    end
    StructTypes.StructType(::Type{Response}) = StructTypes.OrderedStruct()

    results = []
    open_subscription(
        (val)->(push!(results, val); true),
        client,
        "MySub",
        Response,
        output_fields="field")
    @test results[1] isa GraphQLClient.GQLResponse{Response}
    @test isnothing(results[1].errors)
    @test !isnothing(results[1].data) # No point testing content as we've coded it in the test function
end


text:
using MPI
using CLIMA.Topologies
using CLIMA.Grids
using CLIMA.DGBalanceLawDiscretizations
using Printf
using LinearAlgebra

@inline function auxiliary_state_initialization!(aux, x, y, z, dim)
  @inbounds begin
    if dim == 2
      aux[1] = x^2 + y^3 - x*y
      aux[5] = 2*x - y
      aux[6] = 3*y^2 - x
      aux[7] = 0
    else
      aux[1] = x^2 + y^3 + z^2*y^2 - x*y*z
      aux[5] = 2*x - y*z
      aux[6] = 3*y^2 + 2*z^2*y - x*z
      aux[7] = 2*z*y^2 - x*y
    end
  end
end

using Test
function run(dim, Ne, N, DFloat)
  ArrayType = Array

  MPI.Initialized() || MPI.Init()
  Sys.iswindows() || (isinteractive() && MPI.finalize_atexit())

  mpicomm = MPI.COMM_WORLD

  brickrange = ntuple(j->range(DFloat(-1); length=Ne[j]+1, stop=1), dim)
  topl = BrickTopology(mpicomm, brickrange, periodicity=ntuple(j->true, dim))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = DFloat,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )

  spacedisc = DGBalanceLaw(grid = grid,
                           length_state_vector = 0,
                           flux! = (x...) -> (),
                           numerical_flux! = (x...) -> (),
                           auxiliary_state_length = 7,
                           auxiliary_state_initialization! = (x...) ->
                           auxiliary_state_initialization!(x..., dim))

  DGBalanceLawDiscretizations.grad_auxiliary_state!(spacedisc, 1, (2,3,4))

  @test spacedisc.auxstate.Q[:, 2, :] ≈ spacedisc.auxstate.Q[:, 5, :]
  @test spacedisc.auxstate.Q[:, 3, :] ≈ spacedisc.auxstate.Q[:, 6, :]
  @test spacedisc.auxstate.Q[:, 4, :] ≈ spacedisc.auxstate.Q[:, 7, :]
end

let
  numelem = (5, 5, 1)
  lvls = 1

  polynomialorder = 4

  for DFloat in (Float64,) #Float32)
    for dim = 2:3
      err = zeros(DFloat, lvls)
      for l = 1:lvls
        run(dim, ntuple(j->2^(l-1) * numelem[j], dim), polynomialorder, DFloat)
      end
    end
  end
end

isinteractive() || MPI.Finalize()

nothing

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(detect_bedgraph),IOStream})
    precompile(Tuple{typeof(detect_noometiff),IOStream})
    precompile(Tuple{typeof(detect_rdata),IOStream})
    precompile(Tuple{typeof(detect_rdata_single),IOStream})
    precompile(Tuple{typeof(detectwav),IOStream})
    precompile(Tuple{typeof(load),File})
    precompile(Tuple{typeof(load),String})

    precompile(Tuple{typeof(query),String})
    precompile(Tuple{typeof(query),IOStream})
    precompile(Tuple{typeof(query),IOStream,String})
    precompile(Tuple{typeof(hasfunction),Function})
    precompile(Tuple{typeof(hasmagic),Function})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{10,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{20,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{30,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{32,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{35,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{4,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{6,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{7,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},NTuple{8,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},Tuple{UInt8,UInt8,UInt8}})
    precompile(Tuple{typeof(iter_eq),Array{UInt8,1},Tuple{UInt8,UInt8}})
end


text:
__precompile__()
"""
An outline of the API for ACTRModels is presented below. In the REPL, use `?` to see documentation,
    i.e.,

````julia
? print_memory

to see documentation for `print_memory`.
````
# API
## Types
- `ACTR`
- `Declarative`
- `Goal`
- `Imaginal`
- `Chunk`
- `Visual`
- `VisualLocation`
- `Motor`
- `Procedural`
- `Rule`
## Functions
- `get_chunks`
- `get_iconic_memory`
- `get_visicon`
- `get_buffer`
- `set_buffer!`
- `update_lags!`
- `update_recent!`
- `update_chunk!`
- `modify!`
- `add_chunk!`
- `retrieval_prob`
- `retrieval_probs`
- `retrieve`
- `compute_activation!`
- `get_parm`
- `spreading_activation`
- `match`
- `compute_RT`
- `retrieval_request`
- `get_subset`
- `first_chunk`
- `posterior_predictive`
- `find_index`
- `find_indices`
- `import_printing`
- `print_chunk`
- `print_memory`
-  `run`
- `get_mean_activations`
"""
module ACTRModels
    using Reexport
    @reexport using Distributions, Parameters, Random, StatsBase, StatsFuns
    using  ConcreteStructs, PrettyTables
    import Distributions: pdf, logpdf
    import SequentialSamplingModels: LNR
    import Base: rand, match
    export AbstractACTR, ACTR, Declarative, Imaginal, Chunk, BufferState, Mod
    export Goal, Visual, Motor, VisualLocation, Procedural, Rule
    export AbstractVisualObject, VisualObject, Parms, AbstractParms
    export default_penalty, LNR, reduce_data, get_buffer, set_buffer!
    export get_chunks, update_lags!, update_recent!, update_chunk!, modify!, add_chunk!
    export retrieval_prob, retrieval_probs, retrieve, compute_activation!, get_parm
    export spreading_activation!, match, compute_RT, retrieval_request, get_subset
    export first_chunk, posterior_predictive, find_index, find_indices, get_mean_activations
    export get_visicon, get_iconic_memory, get_rules
    export get_time, set_time!, add_time!, reset_time!, rnd_time, get_name
    export blend_chunks, blended_activation

    include("Structs.jl")
    include("Utilities/Utilities.jl")
    include("MemoryFunctions.jl")
end


text:
function _acb_calc_func!(
    out::Ptr{acb_struct},
    inp::Ptr{acb_struct},
    param::Ptr{Cvoid}, # pointer to the actual function
    order::Clong,
    prec::Clong,
)
    @assert iszero(order) || isone(order)
    x = AcbRef(inp, nothing, prec = prec)
    res = AcbRef(out, nothing, prec = prec)
    f! = unsafe_pointer_to_objref(param)

    f!(res, x, analytic = isone(order), prec = prec)

    return zero(Cint)
end

_acb_calc_cfunc!() = @cfunction(
    _acb_calc_func!,
    Cint,
    (Ptr{acb_struct}, Ptr{acb_struct}, Ptr{Cvoid}, Clong, Clong)
)

function calc_integrate!(
    res::AcbLike,
    acb_calc_cfunc::Ptr{Cvoid}, # cfunction
    param,
    a::AcbLike,
    b::AcbLike,
    rel_goal::Clong,
    abs_tol::MagLike,
    options::calc_integrate_opt_struct,
    prec::Clong,
)
    return ccall(
        @libarb(acb_calc_integrate),
        Cint,
        (
            Ref{acb_struct},
            Ptr{Cvoid},
            Any,
            Ref{acb_struct},
            Ref{acb_struct},
            Clong,
            Ref{mag_struct},
            Ref{calc_integrate_opt_struct},
            Clong,
        ),
        res,
        acb_calc_cfunc,
        param,
        a,
        b,
        rel_goal,
        abs_tol,
        options,
        prec,
    )
end

function calc_integrate!(
    res::AcbLike,
    acb_calc_func,
    param,
    a::AcbLike,
    b::AcbLike,
    rel_goal::Clong,
    abs_tol::MagLike,
    options::Ptr{Cvoid},
    prec::Clong,
)
    return ccall(
        @libarb(acb_calc_integrate),
        Cint,
        (
            Ref{acb_struct},
            Ptr{Cvoid},
            Any,
            Ref{acb_struct},
            Ref{acb_struct},
            Clong,
            Ref{mag_struct},
            Ptr{Cvoid},
            Clong,
        ),
        res,
        acb_calc_func,
        param,
        a,
        b,
        rel_goal,
        abs_tol,
        options,
        prec,
    )
end

"""
    _integrate!(res, f, a, b; prec, rel_goal, abs_goal, opts)

Internal function for integration. Assumes that `f!` is of the type
`f!(out::AcbRef, inp::AcbRef; analytic::Bool, prec::Integer)` and
that all other parameters are set to valid values.
"""
function _integrate!(
    f!,
    res::AcbLike,
    a::AcbLike,
    b::AcbLike;
    prec::Integer = precision(res),
    rel_goal::Integer = prec,
    abs_tol::MagLike = set_ui_2exp!(Arblib.Mag(), one(UInt), -prec),
    opts::Union{Ptr{Cvoid},calc_integrate_opt_struct} = C_NULL,
)
    status = calc_integrate!(
        res,
        _acb_calc_cfunc!(),
        f!,
        a,
        b,
        rel_goal,
        abs_tol,
        opts, # passing C_NULL uses the default options
        prec,
    )

    # status:
    # ARB_CALC_SUCCESS = 0
    # ARB_CALC_NO_CONVERGENCE = 2
    status == 2 && @warn "integration did not converge"

    return res
end

"""
    integrate!(f!, res::Acb, a::Number, b::Number;
        check_analytic::Bool = false,
        take_prec::Bool = false,
        prec::Integer = precision(res),
        rtol = 0.0,
        atol = 2.0^-prec,
        opts::Union{acb_calc_integrate_opt_struct, Ptr{Cvoid}} = C_NULL,
    )

Like [`integrate`](@ref), but make use of in-place operations. In
particular, there are three differences from `integrate`:

1. The function `f!` should be of the form `f!(y, x) = set!(y, f(x))`.
   That is, it writes the return value of the integand `f(x)` in-place
   into its first argument `y`. (The return value of `f!` is ignored).

2. `res` is set to the result.

3. The default precision is taken from `res` instead of from `a` and
   `b`.

# Examples
```jldoctest
julia> Arblib.integrate!(Arblib.sin!, Acb(0), Acb(0), Acb(10)) # Integrate sin from 0 to 10
[1.83907152907645245225886394782406483451993016513316854683595373104879258687 +/- 5.15e-75]

julia> Arblib.integrate!(Arblib.inv!, Acb(0), Acb(1, -5), Acb(1, 5)) # Integrate 1/z from 1 - 5i to 1 + 5i
[+/- 2.02e-75] + [2.74680153389003172172254385288992229730199919179940161793956671182574846633 +/- 2.83e-75]im

julia> # Integrate √z from 1 - 5im to 1 + 5im, taking into account the branch cut at (-∞, 0]
julia> f! = (res, z; analytic = false) -> Arblib.sqrt_analytic!(res, z, analytic);
julia> Arblib.integrate!(f!, Acb(0), Acb(1, -5), Acb(1, 10), check_analytic = true, prec = 64)
[-9.0064084416559764 +/- 6.53e-17] + [23.8636067095598007 +/- 6.98e-17]im
```
"""
function integrate!(
    f!,
    res::AcbOrRef,
    a::AcbOrRef,
    b::AcbOrRef;
    check_analytic::Bool = false,
    take_prec::Bool = false,
    prec::Integer = precision(res),
    rtol = 0.0,
    atol = set_ui_2exp!(Mag(), one(UInt), -prec),
    opts::Union{Ptr{Cvoid},calc_integrate_opt_struct} = C_NULL,
)
    # rtol ≈ 2^(-rel_goal)
    rel_goal = iszero(rtol) ? prec : max(-floor(Int, log2(rtol)), 0)

    if !check_analytic && !take_prec
        g! = (res, x; analytic, prec) -> f!(res, x)
    elseif !check_analytic && take_prec
        g! = (res, x; analytic, prec) -> f!(res, x, prec = prec)
    elseif check_analytic && !take_prec
        g! = (res, x; analytic, prec) -> f!(res, x, analytic = analytic)
    else
        g! = f!
    end

    return _integrate!(
        g!,
        res,
        a,
        b,
        prec = prec,
        rel_goal = rel_goal,
        abs_tol = convert(Mag, atol),
        opts = opts,
    )
end

function integrate!(
    f!,
    res::AcbOrRef,
    a::Number,
    b::Number;
    check_analytic::Bool = false,
    take_prec::Bool = false,
    prec::Integer = precision(res),
    rtol = 0.0,
    atol = set_ui_2exp!(Mag(), one(UInt), -prec),
    opts::Union{Ptr{Cvoid},calc_integrate_opt_struct} = C_NULL,
)
    return integrate!(
        f!,
        res,
        Acb(a, prec = prec),
        Acb(b, prec = prec),
        check_analytic = check_analytic,
        take_prec = take_prec,
        prec = prec,
        rtol = rtol,
        atol = atol,
        opts = opts,
    )
end

"""
    integrate(f, a::Number, b::Number;
        check_analytic::Bool = false,
        take_prec::Bool = false,
        prec = max(precision(a), precision(b)),
        rtol=0.0,
        atol=2.0^-prec,
        opts::Union{acb_calc_integrate_opt_struct, Ptr{Cvoid}} = C_NULL)
Computes a rigorous enclosure of the integral ∫ₐᵇ f(x) dx where
`f(x::AcbRef)` is any (holomorphic) julia function. From Arb docs:

> The integral follows a straight-line path between the complex numbers `a` and
> `b`. For finite results, `a`, `b` must be finite and `f` must be bounded on
> the path of integration. To compute improper integrals, the user should
> therefore truncate the path of integration manually (or make a regularizing
> change of variables, if possible).

The error estimates used require that `f` is holomorphic on certain
ellipses around the path of integration.

> The integration algorithm combines direct interval enclosures,
> Gauss-Legendre quadrature where f is holomorphic, and adaptive
> subdivision. This strategy supports integrands with discontinuities
> while providing exponential convergence for typical piecewise
> holomorphic integrands.

In general the integration will work for any function which is
holomorpic or meromorphic on the whole complex plane. For functions
with branch cuts or other things which makes them non-holomorphic the
argument `check_analytic` has to be set to `true`. In this case `f`
will be given a keyword argument `analytic::Bool`, if `analytic` is
`false` then nothing special has to be done, but if `analytic` is
`true` then the output has to be non-finite (typically `Acb(NaN)`) if
`f` is not holomorphic on the whole input ball.

!!! Note: It's users responsibility to verify holomorphicity of `f`.

Parameters:
 * `take_prec` if true then `f` will be given the keyword argument
   `prec = prec`, useful for functions requiring an explicit precision
   to be given.
 * `rtol` relative tolerance
 * `atol` absolute tolerance
 * `opts` a `C_NULL` (using the default options), or an instance of
   `acb_calc_integrate_opt_struct` controlling the algorithmic aspects of integration.

!!! Note: `integrate` does not guarantee to satisfy provided
    tolerances. But the result is guaranteed to be contained in the
    resulting ball.

For more information please consider arblib documentation and the
paper
> Fredrik Johansson, Numerical integration in arbitrary-precision ball arithmetic,
> _Mathematical Software – ICMS 2018_
> https://doi.org/10.1007/978-3-319-96418-8
> https://arxiv.org/abs/1802.07942

# Examples
```jldoctest
julia> Arblib.integrate(sin, 0, 10) # Integrate sin from 0 to 10
[1.83907152907645245225886394782406483451993016513316854683595373104879258687 +/- 5.15e-75]

julia> Arblib.integrate(z -> 1/z, Acb(1, -5), Acb(1, 5)) # Integrate 1/z from 1 - 5i to 1 + 5i
[+/- 2.02e-75] + [2.74680153389003172172254385288992229730199919179940161793956671182574846633 +/- 2.83e-75]im

julia> # Integrate √z from 1 - 5im to 1 + 5im, taking into account the branch cut at (-∞, 0]
julia> f = (z; analytic = false) -> begin
           if analytic && Arblib.contains_nonpositive(real(z))
               return Acb(NaN, prec = precision(z))
           else
               return sqrt(z)
           end
       end;
julia> Arblib.integrate(f, Acb(1, -5), Acb(1, 10), check_analytic = true, prec = 64)
[-9.0064084416559764 +/- 7.40e-17] + [23.8636067095598007 +/- 9.03e-17]im
```
"""
function integrate(
    f,
    a::AcbOrRef,
    b::AcbOrRef;
    check_analytic::Bool = false,
    take_prec::Bool = false,
    prec = max(precision(a), precision(b)),
    rtol = 0.0,
    atol = set_ui_2exp!(Mag(), one(UInt), -prec),
    opts::Union{Ptr{Cvoid},calc_integrate_opt_struct} = C_NULL,
)
    f! = (res, x; kwargs...) -> set!(res, f(x; kwargs...))
    res = Acb(prec = prec)

    return integrate!(
        f!,
        res,
        a,
        b;
        check_analytic = check_analytic,
        take_prec = take_prec,
        rtol = rtol,
        atol = atol,
        opts = opts,
        prec = prec,
    )
end

function integrate(
    f,
    a::Number,
    b::Number;
    check_analytic::Bool = false,
    take_prec::Bool = false,
    prec::Integer = max(_precision(a), _precision(b)),
    rtol = 0.0,
    atol = set_ui_2exp!(Mag(), one(UInt), -prec),
    opts::Union{Ptr{Cvoid},calc_integrate_opt_struct} = C_NULL,
)
    return integrate(
        f,
        Acb(a, prec = prec),
        Acb(b, prec = prec),
        check_analytic = check_analytic,
        take_prec = take_prec,
        prec = prec,
        rtol = rtol,
        atol = atol,
        opts = opts,
    )
end


text:
using Compat
using Base.Test
using CUDArt
using KUnet
density = 0.4

iseq01(a::KUsparse,b::SparseMatrixCSC)=(for f in fieldnames(a); iseq01(a.(f),b.(f)) || (warn("$f mismatch:\n$(a.(f))\n$(b.(f))"); return false); end; return true)
iseq01(a::BaseArray,b::BaseArray)=(to_host(a)==to_host(b))
iseq01(a,b)=(a==b)

for A in (CudaArray, Array)
    for T in (Float64, Float32, Int64, Int32)
        m = rand(1:20)
        n = rand(1:20)
        @show (A,T,m,n)
        a = sprand(m,n,density,rand,T)
        s = convert(KUsparse{A}, copy(a))

        @test atype(s)==A
        @test iseq01(s,a)
        for fname in (:eltype, :length, :ndims, :size, :isempty)
            @test (@eval $fname($s)==$fname($a))
        end
        @test size(s,1)==size(a,1)
        @test size(s,2)==size(a,2)

        # copy
        s1 = cpucopy(s); @test iseq01(s1,a)
        s2 = gpucopy(s); @test iseq01(s2,a)
        s3 = copy(s); @test iseq01(s3,a)
        # b = sprand(m,rand(1:20),density,rand,T)
        # @show map(size, (a,b,s,s1,s2,s3))
        # copy!(s3,b); @test iseq01(s3,b)
        # copy!(s3,s); @test iseq01(s3,a)

        # cslice!
        r1 = rand(1:ccount(a))
        r2 = rand(1:ccount(a))
        r1 <= r2 || ((r1,r2)=(r2,r1))
        s3 = copy(s)
        cslice!(s3, a, r1:r2)
        @test iseq01(s3, a[:,r1:r2])

        # ccat!
        b = sprand(m,rand(1:20),density,rand,T)
        da = convert(KUsparse{A}, copy(a))
        db = convert(KUsparse{A}, copy(b))
        ccat!(da, db)
        @test iseq01(da, [a b])

        # uniq!
        for i=1:5; ccat!(da, db); end
        uniq!(da)
        @test iseq01(da, unique([a b], 2))
    end
end


text:
module InformationMeasures

include("../src/Discretization.jl")
include("../src/Estimators.jl")
include("../src/Formulae.jl")
include("../src/Measures.jl")

end # module


text:
using Pkg

Pkg.add(url="https://github.com/snotskie/EpistemicNetworkAnalysis.jl")
Pkg.add(url="https://github.com/dillondaudert/UMAP.jl")
Pkg.add("Plots")
Pkg.add("LinearAlgebra")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("Statistics")
Pkg.add("Distances")
Pkg.add("Colors")
Pkg.add("GLM")
Pkg.add("Clustering")


text:
@testset "motif.jl" begin
    n = 5

    profile=rand(n,2)

    model = Peak_Motif(Gaussian_Peak(),profile)

    @test parameter_size(model) == 2

    θ = rand(parameter_size(model))

    X = rand(10)
    Y = zeros(10)

    @test (@benchmark accumulate_y!($model,$Y,$X,$θ)).allocs == 0
end


text:
using Cxx

const BoostSharedPtr{T}  = cxxt"boost::shared_ptr<$T>"
const BoostSharedPtrR{T} = cxxt"boost::shared_ptr<$T>&"

function (::Type{BoostSharedPtr{T}})() where T
    icxx"$(BoostSharedPtr{T})();"
end

function make_shared(ptr::BoostSharedPtr{T}) where T
    icxx"$ptr = boost::make_shared<$T>();"
    ptr
end

set(ptr0::BoostSharedPtr{T}, ptr1::BoostSharedPtr{T}) where T = icxx"$ptr0 = $ptr1;"
Base.get(ptr::BoostSharedPtr{T}) where T = icxx"$ptr.get();"

use_count(ptr::BoostSharedPtr) = icxx"$ptr.use_count();"
Base.reset(ptr) = icxx"$ptr.reset();"

export BoostSharedPtr, make_shared, set, use_count

max_stars_repo_name: moyhig/CxxROS

max_issues_repo_name: moyhig/CxxROS

max_forks_repo_name: moyhig/CxxROS


text:
module TSNE

##############################################################################
##
## Dependencies and Reexports
##
##############################################################################

# None

##############################################################################
##
## Exported methods and types
##
##############################################################################

export tsne

##############################################################################
##
## Source code
##
##############################################################################

"""
Distance matrix D containing the pairwise squared Euclidean distances. Each
element is Dᵢⱼ = ∥xᵢ−xⱼ∥² with the diagonal elements given by the input `diagval`.
"""
function pairwise_sq_euc_dist(X::AbstractMatrix{F}, diagval::F) where {F<:AbstractFloat}
    n = size(X, 2)
    s = sum(abs2, X, 1)                # sᵢ = ||xᵢ||²
    D = BLAS.syrk('L', 'T', F(-2), X)  # Computes lower triangle of D = - 2 .* (X' * X)

    # Complete the full distance matrix: D = - 2 .* (X' * X) .+ s .+ s'
    for j in 1:n

        for i in 1:(j - 1)  # Upper triangle - copy from lower triangle (matrix is symmetric)
            D[i, j] = D[j, i]
        end

        # Diagonal - often 0, but can use for tricks e.g. set to -Inf to underflow to 0 under exp or inv
        D[j, j] = diagval

        sj = s[j]
        for i in (j + 1):n  # Lower triangle
            D[i, j] += (s[i] + sj)
        end
    end

    return D
end

"""
Calculates the entropy of the distribution around point j as well as the conditional
probability of each point i given j (pᵢ∣ⱼ) with distances Dⱼ and parameter βⱼ = 1 / 2σⱼ².
Pⱼ is changed in-place and Hⱼ is returned.
"""
function entropy_and_conditional_p!(Pj::AbstractVector{F}, Dj::AbstractVector{F},
                                    βj::F) where {F<:AbstractFloat}
    sumPj, dotDjPj = F(0), F(0)
    for i in eachindex(Pj)
        sumPj += (Pj[i] = exp(-βj * Dj[i]))
        dotDjPj += Dj[i] * Pj[i]
    end
    Hj = (βj * dotDjPj / sumPj + log(sumPj)) / log(F(2))
    Pj ./= sumPj
    return Hj
end

"""
Reduces each vector by it's minimum value. Has the effect of moving the range of
−βⱼDᵢⱼ (the input into the exponential function) to (−∞,0).
"""
function reducemin!(Dj)
    minDj = minimum(Dj)
    for i in eachindex(Dj)
        Dj[i] -= minDj
    end
end

"""
Performs, for each column (point j), a binary search for the value  βⱼ = 1 / 2σⱼ² and
returns the matrix of conditional probabilities pᵢ∣ⱼ. Each column sums to 1 i.e. ∑ᵢpᵢ∣ⱼ = 1.
"""
function conditional_p(D::AbstractMatrix{F}, perplexity::F, tol::F,
                       max_iter::T) where {T<:Integer, F<:AbstractFloat}

    # Number of observations - columns of X
    n = size(D, 1)

    # Initialise P matrix and β vector
    P = zeros(F, n, n)
    β = zeros(F, n)

    # Initialise column vectors
    Dj = zeros(F, n)
    Pj = zeros(F, n)

    # Binary search target
    target = log2(perplexity)

    # Binary search over each column
    for j in 1:n

        # Re-initialise column vectors
        fill!(Pj, F(0))
        copy!(Dj, view(D, :, j))

        # Subtract minimum value for better stability
        reducemin!(Dj)

        # Initial binary search values
        βj_min = F(0)
        βj_max = F(Inf)
        βj = F(1)

        # Binary search loop
        for iter in 1:max_iter

            # Calculate conditional probabilities and entropy given βj
            Hj = entropy_and_conditional_p!(Pj, Dj, βj)

            # Break from loop if within tolerence
            abs(target - Hj) < tol && break

            # Search branches
            if target < Hj
                βj_min = βj
                βj = isinf(βj_max) ? 2βj : (βj_min + βj_max) / 2
            else
                βj_max = βj
                βj = (βj_min + βj_max) / 2
            end

            # Warn if max iterations reached
            iter == max_iter && warn("Max iterations reached for column $j with perplexity = $(2^Hj)")
        end

        # Save the results
        P[:, j] = Pj
        β[j] = βj
    end

    # Return conditional probabilities
    return P
end

"""
Replaces the conditional P matrix in-place with the joint P matrix
"""
function joint_p!(P::AbstractMatrix{F}) where {F<:AbstractFloat}
    n = size(P, 2)
    for j = 1:n
        for i = (j + 1):n  # Lower triangle
            P[i, j] += P[j, i]
            P[j, i] = P[i, j]
        end
    end
    scale!(P, F(inv(2n)))
end

"""
Create the matrix P of joint probabilities in high dimension space
"""
function create_P(X::AbstractMatrix{F}, perplexity::F, tol::F, max_iter::T) where {T<:Integer, F<:AbstractFloat}
    D = pairwise_sq_euc_dist(X,  prevfloat(F(Inf)))
    P = conditional_p(D, perplexity, tol, max_iter)
    joint_p!(P)
    return P
end

"""
Create the matrix Q of joint probabilities in low dimension space
For computation reasons `joint_q!` doesn't return the full Q matrix but instead returns
the matrix of the numerators (1 + ∥yᵢ−yⱼ∥²)⁻¹ - called Qtop - and separately the
denominator ∑ₖ∑ₗ≠ₖ(1 + ∥yₖ−yₗ∥²)⁻¹ - called sumQ - which is identical for all Qᵢⱼ.
"""
function joint_q!(Qtop::AbstractMatrix{F}, s::AbstractMatrix{F},
                  Y::AbstractMatrix{F}) where {F<:AbstractFloat}
    n = size(Qtop, 2)
    sum!(abs2, s, Y)                            # sᵢ = ||Yᵢ||² (in-place)
    BLAS.syrk!('L', 'T', F(-2), Y, F(0), Qtop)  # Lower triangle of Qtop = - 2 .* (Y' * Y) (in-place)

    sumQ = zero(F)
    @inbounds for j = 1:n
        Qtop[j, j] = F(0)  # Diagonal = 0
        sj = s[j]
        @simd for i = (j + 1):n  # Lower triangle
            Qtop[i, j] += (s[i] + sj)
            sumQ += (Qtop[i, j] = inv(one(F) + Qtop[i, j]))
        end
    end
    sumQ = 2 * sumQ
end

"""
Cost function - KL Divergence
"""
function kl_div(P::AbstractMatrix{F}, Qtop::AbstractMatrix{F}, sumQ::F, exag::F) where {F<:AbstractFloat}
    n = size(P, 2)
    cost = zero(F)
    for j = 1:n
        for i = (j + 1):n  # Lower triangle
            @inbounds cost += P[i, j] * log(P[i, j] / Qtop[i, j])
        end
    end
    return (2 * cost) / exag + log(sumQ / exag)
end

"""
Gradients
"""
function grad!(ΔY::AbstractMatrix{F}, Y::AbstractMatrix{F}, P::AbstractMatrix{F},
               Qtop::AbstractMatrix{F}, sumQ::F)  where {F<:AbstractFloat}
    fill!(ΔY, zero(F))
    invsumQ = inv(sumQ)
    n = size(P, 2)
    for j in 1:n
        for i = (j+1):n  # Lower Triangle
            @inbounds tmp1 = 4 * (P[i, j] - Qtop[i, j] * invsumQ) * Qtop[i, j]
            for d in indices(ΔY, 1)
                tmp2 = tmp1 * (Y[d, j] - Y[d, i])
                ΔY[d, j] += tmp2
                ΔY[d, i] -= tmp2
            end
        end
    end
end

"""
Check the analytic gradients with a numeric approximation
"""
function check_gradients(X, d)

    # Number of observations
    n = size(X, 2)

    # Create P Matrix
    P = create_P(X, 30.0, 1e-5, 50)

    # Initialise Y matrix - sample from Normal(0, 1e-4)
    Y = randn(d, n) * 1e-4

    # Initialise re-use matrices
    s, Qtop, ΔY = zeros(1, n), zeros(n, n), zeros(d, n)

    # Calculate gradients analytically
    sumQ = joint_q!(Qtop, s, Y)
    grad!(ΔY, Y, P, Qtop, sumQ)

    # Size of perturbation
    ϵ = 1e-5

    # Initialise numeric gradient array
    numeric_grad = zeros(ΔY)

    # Estimate gradient for each parameter separately
    for i in eachindex(Y)

        # Store value of parameter
        tmp = Y[i]

        # Forward
        Y[i] += ϵ
        sumQ = joint_q!(Qtop, s, Y)
        cost1 = kl_div(P, Qtop, sumQ, 1.0)

        # Backward
        Y[i] -= 2ϵ
        sumQ = joint_q!(Qtop, s, Y)
        cost2 = kl_div(P, Qtop, sumQ, 1.0)

        # Calculate numeric gradient using central difference
        numeric_grad[i] = (cost1 - cost2) ./ 2ϵ

        # Reset the parameter
        Y[i] = tmp
    end

    # Display results
    # Display the difference between numeric and analytic gradient for each parameter
    @printf "\e[1m%s\e[0m" " Numeric Gradient    Analytic Gradient           Difference\n"
    for i=1:min(15, length(numeric_grad))
        @printf "%17e %20e %20e \n" numeric_grad[i] ΔY[i] numeric_grad[i] - ΔY[i]
    end
    @printf "%17s %20s %20s \n\n" "..." "..." "..."
    @printf "\e[1m%s\e[0m" " Largest differences:\n"
    @printf "%17e %20e \n" extrema(numeric_grad .- ΔY)...
end

"""
Perform PCA to reduce the number of dimensions before applying t-SNE
"""
function pca(X, k)
    Σ = X * X' ./ size(X, 1)         # Covariance natrix
    F = svdfact(Σ)                   # Factorise into Σ = U * diagm(S) * V'
    Xrot = F.U' * X                  # Rotate onto the basis defined by U
    pvar = sum(F.S[1:k]) / sum(F.S)  # Percentage of variance retained with top k vectors
    X̃ = Xrot[1:k, :]                 # Keep top k vectors
    return X̃, pvar
end

"""
Gradient descent update
Uses both momentum and adaptive learning
"""
function update!(Y::AbstractArray{F}, Yμ::AbstractArray{F}, ∇Y::AbstractArray{F},
                 Yg::AbstractArray{F}, μ::F, η::F, min_gain::F) where {F<:AbstractFloat}
    for i in eachindex(Y)

        # Update gains - adaptive learning params
        ∇Y[i] * Yμ[i] > 0 ? Yg[i] = max(Yg[i] * 0.8, min_gain) : Yg[i] += 0.2

        # Momentum gradient descent update with adaptive learning
        Yμ[i] = μ * Yμ[i] - η * Yg[i] * ∇Y[i]
        Y[i] += Yμ[i]
    end
end

"""
`t()` is for displaying the time during the gradient descent optimisation.
"""
t() = Dates.format(now(), "HH:MM:SS")

"""
    tsne(X, d, [perplexity = 30.0, perplexity_tol = 1e-5, perplexity_max_iter = 50,
                pca_init = true, pca_dims = 30, exag = 12.0, stop_exag = 250,
                μ_init = 0.5, μ_final = 0.8, μ_switch = 250,
                η = 100.0, min_gain = 0.01, num_iter = 1000])

### Input types
- `F::AbstractFloat`
- `T::Integer`

### Arguments
- `X::Matrix{F}`: Data matrix - where each row is a feature and each column is an observation / point
- `d::T`: The number of dimensions to reduce X to (e.g. 2 or 3)

### Keyword arguments and default values
- `perplexity::F = 30.0`: User specified perplexity for each conditional distribution
- `perplexity_tol::F = 1e-5`: Tolerence for binary search of bandwidth
- `perplexity_max_iter::T = 50`: Maximum number of iterations used in binary search for bandwidth
- `pca_init::Bool = true`: Choose whether to perform PCA before running t-SNE
- `pca_dims::T = 30`: Number of dimensions to reduce to using PCA before applying t-SNE
- `exag::T = 12`: Early exaggeration - multiply all pᵢⱼ's by this constant
- `stop_exag::T = 250`: Stop the early exaggeration after this many iterations
- `μ_init::F = 0.5`: Initial momentum parameter
- `μ_final::F = 0.8`: Final momentum parameter
- `μ_switch::T = 250`: Switch from initial to final momentum parameter at this iteration
- `η::F = 100.0`: Learning rate
- `min_gain::F = 0.01`: Minimum gain for adaptive learning
- `num_iter::T = 1000`: Number of iterations
- `show_every::T = 100`: Display progress at intervals of this number of iterations
"""
function tsne(X::AbstractMatrix{F}, d::T;
              perplexity::F = 30.0, perplexity_tol::F = 1e-5, perplexity_max_iter::T = 50,
              pca_init::Bool = true, pca_dims::T = 30,
              exag::F = 12.0, stop_exag::T = 250,
              μ_init::F = 0.5, μ_final::F = 0.8, μ_switch::T = 250,
              η::F = 100.0, min_gain::F = 0.01, num_iter::T = 1000,
              show_every::T = 100) where {T<:Integer, F<:AbstractFloat}

    # Number of observations
    n = size(X, 2)

    # Run PCA
    if pca_init
        print(t())
        @printf(" Running PCA ...")
        X, pvar = pca(X, pca_dims)
        @printf(" completed ... percentage of variance retained = %.1f%%\n", pvar*100)
    end

    # Create P Matrix
    print(t())
    @printf(" Computing high dimension joint probabilities ...")
    P = create_P(X, perplexity, perplexity_tol, perplexity_max_iter)
    @printf(" completed\n\n")
    scale!(P, exag)             # early exaggeration

    # Initialise matrices
    Y = randn(F, d, n) * 1e-4   # low dimension map - initialise by sampling from Normal(0, 1e-4)
    s = zeros(F, 1, n)          # sum(abs2, Y, 1)
    Qtop = zeros(F, n, n)       # ||Yⱼ - Yⱼ||² etc
    ∇Y = zeros(F, d, n)         # gradient of Y
    Yμ = zeros(F, d, n)         # velocity of Y - momentum gradient descent
    Yg = ones(F, d, n)          # gains of Y - adaptive learning rate params

    # Initial momentum
    μ = μ_init

    # Gradient descent updates
    println("Running gradient descent updates ...")
    for iter in 1:num_iter

        # Calculate Q probabilities
        sumQ = joint_q!(Qtop, s, Y)

        # Calculate gradients
        grad!(∇Y, Y, P, Qtop, sumQ)

        # Update Y
        update!(Y, Yμ, ∇Y, Yg, μ, η, min_gain)

        # Show progress
        if mod(iter, show_every) == 0
            # print progress
            cost = kl_div(P, Qtop, sumQ, exag)
            print(t())
            @printf(" Iteration %d: cost = %.5f, gradient norm = %.5f\n", iter, cost, norm(∇Y))
        end

        # Change momentum param after μ_switch iterations
        iter == μ_switch && (μ = μ_final)

        # Stop exaggeration
        if iter == stop_exag
            scale!(P, inv(exag))
            exag = F(1)
        end
    end

    # Completed results
    return Y
end

end # module


text:
using Aquarius
using Test

@testset "Aquarius.jl" begin
    # Write your tests here.
    Aquarius.greet()
end