text:
module TestOffsetArrays

using Test, DelimitedFiles
using ModifiedArrays
using ModifiedArrays:
    ModifiedOffsetArray,
    IdentityUnitRange, no_offset_view
# using CatIndices: BidirectionalVector

# These tests were copied verbatim from the OffsetArrays package,
# https://github.com/JuliaArrays/OffsetArrays.jl
# See the package for the various authors of the tests.
#
# The tests were adapted to the ModifiedArrays package in the following way:
# - occurrences of the OffsetArray constructor are modified to use ModifiedOffsetArray
# - explicit references to the offsets fields of an OffsetArray were changed to
#   refer to the offsets field of the OffsetMod in this package
# - various "isa" checks that check for the type OffsetArray now check for
#   the type SimpleOffsetArray

# @test isempty(detect_ambiguities(ModifiedArrays, Core))
# println(detect_ambiguities(ModifiedArrays, Base, Core))
@test isempty(detect_ambiguities(ModifiedArrays, Base, Core))

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

@testset "OffsetVector constructors" begin
    local v = rand(5)
    @test ModifiedOffsetVector(v, -2) == ModifiedOffsetArray(v, -2)
    @test ModifiedOffsetVector(v, -2:2) == ModifiedOffsetArray(v, -2:2)
    @test typeof(ModifiedOffsetVector{Float64}(undef, -2:2)) == typeof(ModifiedOffsetArray{Float64}(undef, -2:2))
end

@testset "undef, missing, and nothing constructors" begin
    y = ModifiedOffsetArray{Float32}(undef, (IdentityUnitRange(-1:1),))
    @test axes(y) === (IdentityUnitRange(-1:1),)

    for (T, t) in ((Missing, missing), (Nothing, nothing))
        @test !isassigned(ModifiedOffsetArray{Union{T,Vector{Int}}}(undef, -1:1, -1:1), -1, -1)
        @test ModifiedOffsetArray{Union{T,Vector{Int}}}(t, -1:1, -1:1)[-1, -1] === t
        @test !isassigned(ModifiedOffsetVector{Union{T,Vector{Int}}}(undef, -1:1), -1)
        @test ModifiedOffsetVector{Union{T,Vector{Int}}}(t, -1:1)[-1] === t
    end
end

@testset "high dimensionality" begin
    y = ModifiedOffsetArray{Float64}(undef, -1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
    @test axes(y) == (-1:1, -7:7, -128:512, -5:5, -1:1, -3:3, -2:2, -1:1)
    y[-1,-7,-128,-5,-1,-3,-2,-1] = 14
    y[-1,-7,-128,-5,-1,-3,-2,-1] += 5
    @test y[-1,-7,-128,-5,-1,-3,-2,-1] == 19
end

@testset "Offset range construction" begin
    r = -2:5
    y = ModifiedOffsetArray(r, r)
    @test axes(y) == (r,)
    y = ModifiedOffsetArray(r, (r,))
    @test axes(y) == (r,)
end

@testset "Traits" begin
    A0 = [1 3; 2 4]
    A = ModifiedOffsetArray(A0, (-1,2))                   # IndexLinear
    S = ModifiedOffsetArray(view(A0, 1:2, 1:2), (-1,2))   # IndexCartesian
    @test axes(A) == axes(S) == (0:1, 3:4)
    @test size(A) == size(A0)
    @test size(A, 1) == size(A0, 1)
    @test length(A) == length(A0)
    @test A == ModifiedOffsetArray(A0, 0:1, 3:4)
    @test_throws DimensionMismatch ModifiedOffsetArray(A0, 0:2, 3:4)
    @test_throws DimensionMismatch ModifiedOffsetArray(A0, 0:1, 2:4)
end

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

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

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

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

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

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

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

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

    S = view(A, :, 3)
    @test S == ModifiedOffsetArray([1,2], (A.modifier.offsets[1],))
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[2]
    @test axes(S) === (IdentityUnitRange(0:1),)
    S = view(A, 0, :)
    @test S == ModifiedOffsetArray([1,3], (A.modifier.offsets[2],))
    @test S[3] == 1
    @test S[4] == 3
    @test_throws BoundsError S[1]
    @test axes(S) === (IdentityUnitRange(3:4),)
    S = view(A, 0:0, 4)
    @test S == [3]
    @test S[1] == 3
    @test_throws BoundsError S[0]
    @test axes(S) === (Base.OneTo(1),)
    S = view(A, 1, 3:4)
    @test S == [2,4]
    @test S[1] == 2
    @test S[2] == 4
    @test_throws BoundsError S[3]
    @test axes(S) === (Base.OneTo(2),)
    S = view(A, :, :)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) === IdentityUnitRange.((0:1, 3:4))
end

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    b = 1:2    # copy between AbstractArrays
    bo = ModifiedOffsetArray(1:2, (-3,))
    @test_throws BoundsError copyto!(a, b)
    fill!(a, -1)
    copyto!(a, bo)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    fill!(a, -1)
    copyto!(a, -2, bo)
    @test a[-3] == -1
    @test a[-2] == 1
    @test a[-1] == 2
    @test_throws BoundsError copyto!(a, -4, bo)
    @test_throws BoundsError copyto!(a, -1, bo)
    fill!(a, -1)
    copyto!(a, -3, b, 2)
    @test a[-3] == 2
    @test a[-2] == a[-1] == -1
    @test_throws BoundsError copyto!(a, -3, b, 1, 4)
    am = ModifiedOffsetArray{Int}(undef, (1:1, 7:9))  # for testing linear indexing
    fill!(am, -1)
    copyto!(am, b)
    @test am[1] == 1
    @test am[2] == 2
    @test am[3] == -1
    @test am[1,7] == 1
    @test am[1,8] == 2
    @test am[1,9] == -1
end

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

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

@testset "reductions" begin
    A = ModifiedOffsetArray(rand(4,4), (-3,5))
    @test maximum(A) == maximum(parent(A))
    @test minimum(A) == minimum(parent(A))
    @test extrema(A) == extrema(parent(A))
    C = similar(A)
    cumsum!(C, A, dims = 1)
    @test parent(C) == cumsum(parent(A), dims = 1)
    @test parent(cumsum(A, dims = 1)) == cumsum(parent(A), dims = 1)
    cumsum!(C, A, dims = 2)
    @test parent(C) == cumsum(parent(A), dims = 2)
    R = similar(A, (1:1, 6:9))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 1)
    R = similar(A, (-2:1, 1:1))
    maximum!(R, A)
    @test parent(R) == maximum(parent(A), dims = 2)
    amin, iamin = findmin(A)
    pmin, ipmin = findmin(parent(A))
    @test amin == pmin
    @test A[iamin] == amin
    @test amin == parent(A)[ipmin]
    amax, iamax = findmax(A)
    pmax, ipmax = findmax(parent(A))
    @test amax == pmax
    @test A[iamax] == amax
    @test amax == parent(A)[ipmax]

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

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

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

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

    @test mapslices(v->sort(v), A, dims = 1) == ModifiedOffsetArray(mapslices(v->sort(v), parent(A), dims = 1), A.modifier.offsets)
    @test mapslices(v->sort(v), A, dims = 2) == ModifiedOffsetArray(mapslices(v->sort(v), parent(A), dims = 2), A.modifier.offsets)
end

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

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

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

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

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

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

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

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

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

Base.axes(A::NegativeArray) = map(n -> (-n):(-1), size(A.parent))

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

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

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

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

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

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

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

end #module


text:
module WaterLily

include("util.jl")
export L₂,BC!,@inside,inside,δ,apply,apply!

include("Poisson.jl")
export AbstractPoisson,Poisson,solver!,mult

include("MultiLevelPoisson.jl")
export MultiLevelPoisson,solver!,mult

include("Flow.jl")
export Flow,mom_step!

include("Body.jl")
export AbstractBody,measure!

include("AutoBody.jl")
export AutoBody,measure

include("Metrics.jl")
using LinearAlgebra: norm2

"""
    Simulation(dims::Tuple, u_BC::Vector, L::Number;
               U=norm2(u_BC), Δt=0.25, ν=0., ϵ = 1,
               uλ::Function=(i,x)->u_BC[i],
               body::AbstractBody=NoBody())

Constructor for a WaterLily.jl simulation:
    - `dims`: Simulation domain dimensions.
    - `u_BC`: Simulation domain velocity boundary conditions, `u_BC[i]=uᵢ, i=1,2...`.
    - `L`: Simulation length scale.
    - `U`: Simulation velocity scale.
    - `ϵ`: BDIM kernel width.
    - `Δt`: Initial time step.
    - `ν`: Scaled viscosity (`Re=UL/ν`)
    - `uλ`: Function to generate the initial velocity field.
    - `body`: Embedded geometry

See files in `examples` folder for examples.
"""
struct Simulation
    U :: Number # velocity scale
    L :: Number # length scale
    ϵ :: Number # kernel width
    flow :: Flow
    body :: AbstractBody
    pois :: AbstractPoisson
    function Simulation(dims::Tuple, u_BC::Vector, L::Number;
                        Δt=0.25, ν=0., U=norm2(u_BC), ϵ = 1,
                        uλ::Function=(i,x)->u_BC[i],
                        body::AbstractBody=NoBody(),T=Float64)
        flow = Flow(dims,u_BC;uλ,Δt,ν,T)
        measure!(flow,body;ϵ)
        new(U,L,ϵ,flow,body,MultiLevelPoisson(flow.μ₀))
    end
end

time(sim::Simulation) = sum(sim.flow.Δt[1:end-1])
"""
    sim_time(sim::Simulation)

Return the current dimensionless time of the simulation `tU/L`
where `t=sum(Δt)`, and `U`,`L` are the simulation velocity and length
scales.
"""
sim_time(sim::Simulation) = time(sim)*sim.U/sim.L

"""
    sim_step!(sim::Simulation,t_end;verbose=false)

Integrate the simulation `sim` up to dimensionless time `t_end`.
If `verbose=true` the time `tU/L` and adaptive time step `Δt` are
printed every time step.
"""
function sim_step!(sim::Simulation,t_end;verbose=false,remeasure=false)
    t = time(sim)
    while t < t_end*sim.L/sim.U
        remeasure && measure!(sim,t)
        mom_step!(sim.flow,sim.pois) # evolve Flow
        t += sim.flow.Δt[end]
        verbose && println("tU/L=",round(t*sim.U/sim.L,digits=4),
            ", Δt=",round(sim.flow.Δt[end],digits=3))
    end
end

"""
    measure!(sim::Simulation,t=time(sim))

Measure a dynamic `body` to update the `flow` and `pois` coefficients.
"""
function measure!(sim::Simulation,t=time(sim))
    measure!(sim.flow,sim.body;t,ϵ=sim.ϵ)
    update!(sim.pois,sim.flow.μ₀)
end

export Simulation,sim_step!,sim_time,measure!
end # module

ext: jl


text:
################################################################################
#
#  Steinitz number
#
################################################################################

function steinitz_number(a::T) where T <: Union{gfp_poly, gfp_fmpz_poly}
  p =  characteristic(parent(a))
  pq = fmpz(1)
  st_a = fmpz(0)
  for i = 0:degree(a)
    st_a += lift(coeff(a, i))*pq
    pq *= p
  end
  return st_a
end

function steinitz_number(a::fq_nmod)
  q = characteristic(parent(a))
  pq = fmpz(1)
  st_a = fmpz(0)
  for i = 0:degree(parent(a))
    st_a += coeff(a, i)*pq
    pq *= q
  end
  return st_a
end

function steinitz_number(a::fq)
  q = characteristic(parent(a))
  pq = fmpz(1)
  st_a = fmpz(0)
  for i = 0:degree(parent(a))
    st_a += steinitz_number(coeff(a, i))*pq
    pq *= q
  end
  return st_a
end

function steinitz_number(a::T) where T <: RelFinFieldElem
  q = order(base_field(parent(a)))
  pq = fmpz(1)
  st_a = fmpz(0)
  for i = 0:degree(parent(a))
    st_a += steinitz_number(coeff(a, i))*pq
    pq *= q
  end
  return st_a
end

function steinitz_number(a::T) where T <: Union{gfp_elem, gfp_fmpz_elem}
  return lift(a)
end

################################################################################
#
#  Presentation for extensions of prime power degree
#
################################################################################
@doc Markdown.doc"""
    presentation(p, r, n) -> FinField

Computes a presentation for the finite field of order p^(r^n) as defined
by Plesken.
"""
function presentation(p::T, r::T, n::Int) where T <: Union{fmpz, Int}
  F = GF(p, cached = false)
  if r == p
    return _presentation_artin_schreier(F, n)
  elseif iszero(mod(p - 1, r))
    return _presentation_kummer(F, r, n)
  else
    return _presentation_generic(F, r, n)
  end
end

function _presentation_artin_schreier(F, n)
  p = characteristic(F)
  Fx, x = PolynomialRing(F, "x", cached = false)
  F1, a = FiniteField(x^p-x-1, "a", cached = false, check = false)
  F1y, y = PolynomialRing(F1, "y1", cached = false)
  el = a
  for i = 2:n
    pol = y^p-y-el^(p-1)
    Frel = FiniteField(pol, "a$i", cached = false, check = false)[1]
    abs_def_pol = norm(pol)
    F1, gF1 = FiniteField(abs_def_pol, "a", check = false, cached = false)
    mp = hom(F1, Frel, gen(Frel))
    el = mp\(gen(Frel)*el)
    F1y, y = PolynomialRing(F1, "y1", cached = false)
  end
  return F1
end

function smallest_pkth_root(F, r)
  #First, I find a rk-th root of unity
  #by taking random elements
  p = order(F)
  nit, cop = ppio(p-1, fmpz(r))
  exp_order = divexact(nit, r)
  x = rand(F)^cop
  while iszero(x) || isone(x^exp_order)
    x = rand(F)^cop
  end
  k = valuation(p-1, r)
  #Now, iteratively, I search for
  #the smallest r-th root
  pow_ind = 1
  for i = 1:k
    root = x^(pow_ind*fmpz(r)^(k-i))
    roots = powers(root, Int(r-1))
    new_ind = 2
    st = steinitz_number(roots[2])
    for j = 3:length(roots)
      stj = steinitz_number(roots[j])
      if stj < st
        st = stj
        new_ind = j
      end
    end
    pow_ind *= new_ind-1
  end
  return x^pow_ind
end

function _presentation_kummer(F, r::T, n::Int) where T <: Union{fmpz, Int}
  pr_root = smallest_pkth_root(F, r)
  Fx, gFx = PolynomialRing(F, cached = false)

  def_pol1 = Fx()
  setcoeff!(def_pol1, 0, -pr_root)
  setcoeff!(def_pol1, r^n, one(F))
  F1, gF1 = FiniteField(def_pol1, "a1", cached = false, check = false)
  return F1
end

function _presentation_generic(F, r::T, n::Int) where T <: Union{fmpz, Int}
  #First, I add the right roots of unity to the field.
  xZx = ZZ["x"][2]
  phi = cyclotomic(r, xZx)
  lf = factor(phi, F)
  lF = collect(keys(lf.fac))
  ind = 1
  st = steinitz_number(lF[1])
  for i = 2:length(lF)
    st_i = steinitz_number(lF[i])
    if st_i < st
      st = st_i
      ind = i
    end
  end
  F0, gF0 = FiniteField(lF[ind], "a0", cached = false)
  f = degree(F0)
  Fn = _presentation_kummer(F0, r, n)
  #Now, I need to take the trace.
  p = characteristic(F)
  e = _find_exponent(f, p, fmpz(r), n)
  t = gen(Fn)
  g = gen(Fn)
  for i = 1:f-1
    g = g^e
    t = add!(t, t, g)
  end
  return FiniteField(_minpoly(t, r^n), "a", cached = false, check = false)[1]
end

function _find_exponent(f::Int, p::fmpz, r::fmpz, n::Int)
  xZx = ZZ["x"][2]
  phi = cyclotomic(f, xZx)
  R = ResidueRing(FlintZZ, r^(n+1), cached = false)
  Rx = PolynomialRing(R, "x", cached = false)[1]
  phiR = Rx(phi)
  phiR1 = derivative(phiR)
  a = p
  it = max(clog(fmpz(n), 2)+1, 2)
  for i = 1:it
    ev1 = phiR(a)
    ev2 = phiR1(a)
    a = a - divexact(ev1, ev2)
  end
  return lift(a)
end


text:
##### Beginning of file

module CommandLine # Begin submodule MirrorUpdater.CommandLine

__precompile__(true)

import ArgParse
import Dates
import HTTP
import Pkg
import TimeZones

import ..Types
import ..Utils
import ..Common
import ..Run

import ..process_delayed_error_list

function run_mirror_updater_command_line!!(
        ;
        registry_list::Vector{Types.Registry},
        arguments::Vector{String} =
            String[],
        git_hosting_providers::AbstractVector =
            Any[],
        additional_repos::Vector{Types.SrcDestPair} =
            Types.SrcDestPair[],
        do_not_push_to_these_destinations::Vector{String} =
            String[],
        do_not_try_url_list::Vector{String} =
            String[],
        try_but_allow_failures_url_list::Vector{String} =
            String[],
        time_zone::TimeZones.TimeZone =
            TimeZones.TimeZone("America/New_York"),
        )::Nothing
    @info(
        "Running MirrorUpdater.CommandLine.run_mirror_updater_command_line!!"
        )
    @info("parsing command line arguments...")
    parsed_arguments::Dict = _parse_arguments(
        arguments
        )
    @info("processing parsed command line arguments...")
    processed_arguments::Dict = _process_parsed_arguments(
        parsed_arguments
        )
    task::String = processed_arguments[:task]
    has_gist_description::Bool =  processed_arguments[:has_gist_description]
    gist_description::String = processed_arguments[:gist_description]
    is_dry_run::Bool = processed_arguments[:is_dry_run]
    delete_gists_older_than_minutes::Int =
        processed_arguments[:delete_gists_older_than_minutes]
    Run.run_mirror_updater!!(
        ;
        git_hosting_providers = git_hosting_providers,
        task = task,
        gist_description = gist_description,
        is_dry_run = is_dry_run,
        registry_list = registry_list,
        additional_repos = additional_repos,
        time_zone = time_zone,
        do_not_push_to_these_destinations =
            do_not_push_to_these_destinations,
        do_not_try_url_list =
            do_not_try_url_list,
        try_but_allow_failures_url_list =
            try_but_allow_failures_url_list,
        delete_gists_older_than_minutes =
            delete_gists_older_than_minutes,
        )
    process_delayed_error_list()
    return nothing
end

function _parse_arguments(arguments::Vector{String})::Dict
    s = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table s begin
        "--task"
            help = "which task to run"
            arg_type = String
            default = ""
        "--gist-description"
            help = "description for the temporary gist"
            arg_type = String
            default = ""
        "--dry-run"
            help = "do everything except actually pushing the repos"
            action = :store_true
        "--delete-gists-older-than-minutes"
            help = "delete all gists older than N minutes"
            arg_type = Int
            default = 0
    end
    result::Dict = ArgParse.parse_args(arguments, s)
    return result
end

function _process_parsed_arguments(parsed_arguments::Dict)::Dict{Symbol, Any}
    task_argument::String = strip(
        convert(String, parsed_arguments["task"])
        )
    if length(task_argument) > 0
        task = task_argument
    else
        task = "all"
    end

    gist_description::String = strip(
        convert(String, parsed_arguments["gist-description"])
        )
    if length(gist_description) > 0
        has_gist_description = true
    else
        has_gist_description = false
    end

    is_dry_run::Bool = parsed_arguments["dry-run"]
    delete_gists_older_than_minutes::Int =
        parsed_arguments["delete-gists-older-than-minutes"]

    result::Dict{Symbol, Any} = Dict{Symbol, Any}()
    result[:task] = task
    result[:has_gist_description]=has_gist_description
    result[:gist_description]=gist_description
    result[:is_dry_run]=is_dry_run
    result[:delete_gists_older_than_minutes]=delete_gists_older_than_minutes

    return result
end

end # End submodule MirrorUpdater.CommandLine

##### End of file


text:
cube, angles, _psf = BetaPictoris[:cube, :pa, :psf]
cube .-= minimum(cube) # rescale so minimum is >0
psf = _psf ./ maximum(_psf)

@testset "throughput - $(nameof(typeof(alg)))" for alg in [Classic(), PCA(10), GreeDS(3)]
    tt, meta = throughput(alg, cube, angles, psf; fwhm=4)

    @test length(tt) == 10
    @test length(meta.distance) == length(meta.noise) == 10
    @test size(meta.fake_comps) == size(cube)[2:3]
    @test keys(meta) == (:distance, :fake_comps, :noise)

    t1 = throughput(alg, cube, angles, psf, (51, 61); fwhm=4)
    t2 = throughput(alg, cube, angles, psf, Polar(10, deg2rad(90)); fwhm=4)
    @test t1 ≈ t2
end

@testset "contrast - $(nameof(typeof(alg)))" for alg in [Classic(), PCA(10), GreeDS(3)]
    cc = contrast_curve(alg, cube, angles, psf; fwhm=4)

    @test keys(cc) == (:distance, :throughput, :contrast, :contrast_corr, :noise)
    ls = map(length, values(cc))
    @test all(ls .== 38)

    reduced_empty = alg(cube, angles)
    sp = Metrics.estimate_starphot(cube, 4)
    cc_raw = contrast_curve(alg, cube, angles, psf; fwhm=4, subsample=false, starphot=sp)
    cc_sub = Metrics.subsample_contrast(reduced_empty, cc_raw.distance, cc_raw.throughput; fwhm=4, starphot=sp)

    for (v1, v2) in zip(values(cc), values(cc_sub))
        v1_ = filter(isfinite, v1)
        v2_ = filter(isfinite, v2)
        @test v1_ ≈ v2_
    end
end

@testset "throughput - AnnulusView - $(nameof(typeof(alg)))" for alg in [Classic(), PCA(10), GreeDS(3)]
    av = AnnulusView(cube; inner=8)
    tt, meta = throughput(alg, av, angles, psf; fwhm=4)

    @test length(tt) == 9
    @test length(meta.distance) == length(meta.noise) == 9
    @test size(meta.fake_comps) == size(cube)[2:3]
    @test keys(meta) == (:distance, :fake_comps, :noise)

    t1 = throughput(alg, cube, angles, psf, (51, 61); fwhm=4)
    t2 = throughput(alg, cube, angles, psf, Polar(10, deg2rad(90)); fwhm=4)
    @test t1 ≈ t2
end

@testset "throughput - MultiAnnulusView - $(nameof(typeof(alg)))" for alg in [Classic(), PCA(10)]
    av = MultiAnnulusView(cube, 4; inner=8)
    tt, meta = throughput(alg, av, angles, psf; fwhm=4)

    @test length(tt) == 9
    @test length(meta.distance) == length(meta.noise) == 9
    @test size(meta.fake_comps) == size(cube)[2:3]
    @test keys(meta) == (:distance, :fake_comps, :noise)

    t1 = throughput(alg, cube, angles, psf, (51, 61); fwhm=4)
    t2 = throughput(alg, cube, angles, psf, Polar(10, deg2rad(90)); fwhm=4)
    @test t1 ≈ t2
end


text:
function julia_to_gamess_key(julia_key)
  julia_to_gamess_mapping = Dict( 7 => 1,
                                  3 => 2,
                                  1 => 3,
                                  9 => 4,
                                  5 => 5,
                                  8 => 6,
                                  4 => 7,
                                  2 => 8,
                                  6 => 9
                                )

  ijkl_match = eachmatch(r"[0-9]+", julia_key)
  ijkl = collect(ijkl_match)
  i = parse(Int,ijkl[1].match)
  j = parse(Int,ijkl[2].match)
  k = parse(Int,ijkl[3].match)
  l = parse(Int,ijkl[4].match)

  i_new = julia_to_gamess_mapping[i]
  j_new = julia_to_gamess_mapping[j]
  k_new = julia_to_gamess_mapping[k]
  l_new = julia_to_gamess_mapping[l]

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  if k_new <= l_new
    tmp = k_new
    k_new = l_new
    l_new = tmp
  end

  ij = i_new > j_new ? i_new*(i_new-1)/2 + j_new : j_new*(j_new-1)/2 + i_new
  kl = k_new > l_new ? k_new*(k_new-1)/2 + l_new : l_new*(l_new-1)/2 + k_new

  if ij <= kl
    tmp = i_new
    i_new = k_new
    k_new = tmp

    tmp = j_new
    j_new = l_new
    l_new = tmp
  end

  return "$i_new,$j_new,$k_new,$l_new"
end

function gamess_to_julia_key(gamess_key)
  gamess_to_julia_mapping = Dict( 1 => 7,
                                  2 => 3,
                                  3 => 1,
                                  4 => 9,
                                  5 => 5,
                                  6 => 8,
                                  7 => 4,
                                  8 => 2,
                                  9 => 6
                                )

  ijkl_match = eachmatch(r"[0-9]+", gamess_key)
  ijkl = collect(ijkl_match)
  i = parse(Int,ijkl[1].match)
  j = parse(Int,ijkl[2].match)
  k = parse(Int,ijkl[3].match)
  l = parse(Int,ijkl[4].match)

  i_new = gamess_to_julia_mapping[i]
  j_new = gamess_to_julia_mapping[j]
  k_new = gamess_to_julia_mapping[k]
  l_new = gamess_to_julia_mapping[l]

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  if k_new <= l_new
    tmp = k_new
    k_new = l_new
    l_new = tmp
  end

  ij = i_new > j_new ? i_new*(i_new-1)/2 + j_new : j_new*(j_new-1)/2 + i_new
  kl = k_new > l_new ? k_new*(k_new-1)/2 + l_new : l_new*(l_new-1)/2 + k_new

  if ij <= kl
    tmp = i_new
    i_new = k_new
    k_new = tmp

    tmp = j_new
    j_new = l_new
    l_new = tmp
  end

  return "$i_new,$j_new,$k_new,$l_new"
end

function get_d_quartet_integrals(input_file_string::String)
  #== read in input file ==#
  input_file::Array{String,1} = []
  open(input_file_string) do file
    input_file = readlines(file)
  end

  #== extract integrals from input file ==#
  starting::Int64 = 1
  ending::Int64 = 1

  integrals = Vector{String}([])
  for line in input_file
    if occursin("QUARTET(",line)
      push!(integrals,line)
    end
  end

  #== set up integral dictionary ==#
  quartet_dict = Dict([])
  for integral in integrals
    #==get ijkl numbers ==#
    ijkl_match = eachmatch(r"[0-9]+", integral)
    ijkl = collect(ijkl_match)
    i = ijkl[1].match
    j = ijkl[2].match
    k = ijkl[3].match
    l = ijkl[4].match

    #== get actual integral ==#
    integral_match = match(r"0","0")
    if typeof(match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",integral)) != Nothing
      integral_match = match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",integral)
    elseif typeof(match(r"(-\.[\d]{1,})", integral)) != Nothing
      integral_match = match(r"(-\.[\d]{1,})",integral)
    else
      integral_match = match(r"(-?[\d]\.[\d]{1,})",integral)
    end

    integral_value = parse(Float64,integral_match.match)
    if !haskey(quartet_dict, "$i,$j,$k,$l")
      quartet_dict["$i,$j,$k,$l"] = [ integral_value ]
    else
      push!(quartet_dict["$i,$j,$k,$l"], integral_value)
    end
  end

  return quartet_dict
end

function check_lists(julia_dict, gamess_dict)
  julia_key_amount = length(keys(julia_dict))
  gamess_key_amount = length(keys(gamess_dict))

  same_size = true
  if julia_key_amount != gamess_key_amount
    println("The lists are of different size!: $julia_key_amount, $gamess_key_amount")
    same_size = false
  else
    println("Both lists have $julia_key_amount quartets. You are good.")
  end

  return same_size
end

function check_gamess_keys_integrals(julia_dict, gamess_dict)
  for gamess_quartet in keys(gamess_dict)
    julia_quartet = gamess_to_julia_key(gamess_quartet)
    if !haskey(julia_dict, julia_quartet)
      println("Missing Julia key! $gamess_quartet => $julia_quartet")
      continue
    end

    julia_value = julia_dict[julia_quartet]
    gamess_value = gamess_dict[gamess_quartet]

    julia_array_sorted = sort(julia_value)
    gamess_array_sorted = sort(gamess_value)

    @assert length(julia_value) == length(gamess_value)
    same = true
    for idx in 1:length(julia_value)
      same = same && isapprox(julia_array_sorted[idx],
        gamess_array_sorted[idx]; atol=1e-6)
    end
    if !same
      println("Different quartet! $julia_quartet and $gamess_quartet")
    end
  end
end

function check_julia_keys_integrals(julia_dict, gamess_dict)
  for julia_quartet in keys(julia_dict)
    gamess_quartet = julia_to_gamess_key(julia_quartet)
    if !haskey(gamess_dict, gamess_quartet)
      println("Missing GAMESS key! $julia_quartet => $gamess_quartet")
      continue
    end

    julia_value = julia_dict[julia_quartet]
    gamess_value = gamess_dict[gamess_quartet]

    julia_array_sorted = sort(julia_value)
    gamess_array_sorted = sort(gamess_value)

    @assert length(julia_value) == length(gamess_value)
    same = true
    for idx in 1:length(julia_value)
      same = same && isapprox(julia_array_sorted[idx],
        gamess_array_sorted[idx]; atol=1e-6)
    end
    if !same
      println("Different quartet! $julia_quartet and $gamess_quartet")
    end
  end
end

julia_dict = get_d_quartet_integrals(ARGS[1])
gamess_dict = get_d_quartet_integrals(ARGS[2])

@assert check_lists(julia_dict, gamess_dict)
check_gamess_keys_integrals(julia_dict, gamess_dict)
check_julia_keys_integrals(julia_dict, gamess_dict)


text:
# The @compat macro is used to implement compatibility rules that require
# syntax rewriting rather than simply new function/constant/module definitions.

using Base.Meta
export @compat

if !isdefined(Base, :UndefKeywordError)
    struct UndefKeywordError <: Exception
        kw
    end
    Base.showerror(io::IO, e::UndefKeywordError) = print(io, "UndefKeywordError: keyword argument $(e.kw) not assigned")
    export UndefKeywordError
end

"Convert a functions symbol argument to the corresponding required keyword argument."
function symbol2kw(sym::Symbol)
    Expr(:kw, sym, Expr(:call, throw, UndefKeywordError(sym)))
end
symbol2kw(arg) = arg

function _compat(ex::Expr)
    if ex.head === :call
        f = ex.args[1]
        if !isdefined(Base, :UndefKeywordError) && length(ex.args) > 1 && isexpr(ex.args[2], :parameters)
            params = ex.args[2]
            params.args = map(symbol2kw, params.args)
        end
    elseif ex.head === :quote && isa(ex.args[1], Symbol)
        # Passthrough
        return ex
    end
    if VERSION < v"0.7.0-DEV.880"
        if ex.head == :curly && ex.args[1] == :CartesianRange && length(ex.args) >= 2
            a = ex.args[2]
            if a != :CartesianIndex && !(isa(a, Expr) && a.head == :curly && a.args[1] == :CartesianIndex)
                return Expr(:curly, :CartesianRange, Expr(:curly, :CartesianIndex, ex.args[2]))
            end
        end
    end
    if VERSION < v"0.7.0-DEV.2562"
        if ex.head == :call && ex.args[1] == :finalizer
            ex.args[2], ex.args[3] = ex.args[3], ex.args[2]
        end
    end
    return Expr(ex.head, map(_compat, ex.args)...)
end

_compat(ex) = ex

macro compat(ex)
    esc(_compat(ex))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6b400190-a550-4e32-baa6-d54b66d7676e
begin
	let
		using Pkg
		Pkg.activate(mktempdir())
		Pkg.Registry.update()
		Pkg.add("PlutoUI")
		Pkg.add("Statistics")
		Pkg.add("StatsBase")
		Pkg.add("ImageMorphology")
		Pkg.add("ImageFiltering")
		Pkg.add("CSV")
		Pkg.add("DataFrames")
		Pkg.add("GLM")
		Pkg.add(url="https://github.com/JuliaHealth/DICOM.jl")
		Pkg.add(url="https://github.com/Dale-Black/DICOMUtils.jl")
		Pkg.add(url="https://github.com/Dale-Black/PhantomSegmentation.jl")
		Pkg.add(url="https://github.com/Dale-Black/CalciumScoring.jl")
		Pkg.add("CairoMakie")
		Pkg.add("HypothesisTests")
		Pkg.add("Colors")
	end

	using PlutoUI
	using Statistics
	using StatsBase: quantile!, rmsd
	using ImageMorphology
	using ImageFiltering
	using CSV
	using DataFrames
	using GLM
	using DICOM
	using DICOMUtils
	using PhantomSegmentation
	using CalciumScoring
	using CairoMakie
	using HypothesisTests
	using Colors
end

# ╔═╡ 38cd4959-e50a-4724-849d-9a246dd4414e
TableOfContents()

# ╔═╡ 5017031f-2e8d-44f9-925d-a14376050342
md"""
## Integrated Mass Scoring
"""

# ╔═╡ 787dda66-c0c9-4e2f-8646-c1c5ac95feee
canon_path1 = "/Users/daleblack/Google Drive/dev/MolloiLab/CAC-stanford-data/data/output/Canon_Aquilion_One_Vision";

# ╔═╡ 1a7206ea-4be5-4236-954e-a927f9cba993
df_i_canon1 = CSV.read(string(canon_path1, "/full.csv"), DataFrame);

# ╔═╡ 955d74f0-2ee9-40e9-9275-01555e13daa2
begin
	fInt1 = Figure()
	axInt1 = Axis(fInt1[1, 1])

	scatter!(df_i_canon1[!, :ground_truth_mass_large], df_i_canon1[!, :calculated_mass_large], label="Large Inserts")
	scatter!(df_i_canon1[!, :ground_truth_mass_medium], df_i_canon1[!, :calculated_mass_medium], label="Medium Inserts")
	scatter!(df_i_canon1[!, :ground_truth_mass_small], df_i_canon1[!, :calculated_mass_small], label="Small Inserts", color=:red)
	lines!([0, 100], [0, 100], label="Unity")

	axInt1.title = "IHU Mass vs Known Mass"
	axInt1.ylabel = "Calculated Mass (mg)"
	axInt1.xlabel = "Known Mass (mg)"
	axInt1.xticks = [0, 25, 50, 75, 100]
	axInt1.yticks = [0, 25, 50, 75, 100]

	xlims!(axInt1, 0, 100)
	ylims!(axInt1, 0, 100)

	fInt1[1, 2] = Legend(fInt1, axInt1, framevisible = false)
	fInt1
end

# ╔═╡ 8d1452f9-032c-41df-830e-ab6459e7c5ae
begin
	fInt1_small = Figure()
	axInt1_small = Axis(fInt1_small[1, 1])

	scatter!(df_i_canon1[!, :ground_truth_mass_small], df_i_canon1[!, :calculated_mass_small], label="Small Inserts", color=:red)
	lines!([0, 1.2], [0, 1.2], label="Unity")

	axInt1_small.title = "IHU Mass vs Known Mass (Small)"
	axInt1_small.ylabel = "Calculated Mass (mg)"
	axInt1_small.xlabel = "Known Mass (mg)"

	fInt1_small[1, 2] = Legend(fInt1_small, axInt1_small, framevisible = false)
	fInt1_small
end

# ╔═╡ 067b9f6b-6b57-40c0-ae5f-0465400f77bb
md"""
### RMSD measurements (large, medium, and small inserts)
"""

# ╔═╡ 7f4d43df-f51e-49c2-b47e-629e3ccedc91
int1_rms_large_c = rmsd(df_i_canon1[!, :calculated_mass_large], df_i_canon1[!, :ground_truth_mass_large])

# ╔═╡ 39886e68-c4ea-4b96-96fb-ba46eeb54ecd
int1_rms_medium_c = rmsd(df_i_canon1[!, :calculated_mass_medium], df_i_canon1[!, :ground_truth_mass_medium])

# ╔═╡ fdad40dc-5e3c-4768-837f-7dbf4fca07ac
int1_rms_small_c = rmsd(df_i_canon1[!, :calculated_mass_small], df_i_canon1[!, :ground_truth_mass_small])

# ╔═╡ 565fab58-ef10-43e8-8ba8-e2e792ef872f
md"""
## Agatston Mass Scoring
"""

# ╔═╡ 09bfb23a-db24-42c8-8c34-8234da9c0d56
canon_path_agat = "/Users/daleblack/Google Drive/dev/MolloiLab/CAC-stanford-data/data/output/Canon_Aquilion_One_Vision_agatston_mass_score_script";

# ╔═╡ 5c1bea23-972b-4867-a424-9119cbc5edbc
df_a_canon = CSV.read(string(canon_path_agat, "/output.csv"), DataFrame);

# ╔═╡ 65131928-ad9a-4649-9600-e663ab5a7a28
begin
	fAgat = Figure()
	axAgat = Axis(fAgat[1, 1])

	scatter!(df_a_canon[!, :ground_truth_mass_large], df_a_canon[!, :calculated_mass_large], label="Mass: Large Inserts")
	scatter!(df_a_canon[!, :ground_truth_mass_medium], df_a_canon[!, :calculated_mass_medium], label="Mass: Medium Inserts")
	scatter!(df_a_canon[!, :ground_truth_mass_small], df_a_canon[!, :calculated_mass_small], label="Mass: Small Inserts", color=:red)
	lines!([0, 100], [0, 100], label="Unity Line")

	axAgat.title = "AS Mass vs Known Mass"
	axAgat.ylabel = "Calculated Mass (mg)"
	axAgat.xlabel = "Known Mass (mg)"
	axAgat.xticks = [0, 25, 50, 75, 100]
	axAgat.yticks = [0, 25, 50, 75, 100]

	xlims!(axAgat, 0, 100)
	ylims!(axAgat, 0, 100)

	fAgat[1, 2] = Legend(fAgat, axAgat, framevisible = false)

	fAgat
end

# ╔═╡ c4bc4b99-7adc-48a2-a3b9-114d6c5a4e3a
begin
	fAgat_small = Figure()
	axAgat_small = Axis(fAgat_small[1, 1])

	scatter!(df_a_canon[!, :ground_truth_mass_small], df_a_canon[!, :calculated_mass_small], label="Mass: Small Inserts", color=:red)
	lines!([0, 1], [0, 1], label="Unity Line")

	axAgat_small.title = "AS Mass vs Known Mass (Small)"
	axAgat_small.ylabel = "Calculated Mass (mg)"
	axAgat_small.xlabel = "Known Mass (mg)"

	fAgat_small[1, 2] = Legend(fAgat_small, axAgat_small, framevisible = false)

	fAgat_small
end

# ╔═╡ 243c0c2b-e84d-47eb-aaf8-5fb79d76d04d
md"""
### RMSD measurements (large, medium, and small inserts)
"""

# ╔═╡ c6cd7f7d-d77b-4136-9485-913f962be85a
agat_rms_large_c = rmsd(df_a_canon[!, :calculated_mass_large], df_a_canon[!, :ground_truth_mass_large])

# ╔═╡ 7f8756e1-36e0-444a-a9a9-e20f763298af
agat_rms_medium_c = rmsd(df_a_canon[!, :calculated_mass_medium], df_a_canon[!, :ground_truth_mass_medium])

# ╔═╡ ec955132-42a3-40cd-a031-480ee6e95545
agat_rms_small_c = rmsd(df_a_canon[!, :calculated_mass_small], df_a_canon[!, :ground_truth_mass_small])

# ╔═╡ Cell order:
# ╠═6b400190-a550-4e32-baa6-d54b66d7676e
# ╠═38cd4959-e50a-4724-849d-9a246dd4414e
# ╟─5017031f-2e8d-44f9-925d-a14376050342
# ╠═787dda66-c0c9-4e2f-8646-c1c5ac95feee
# ╠═1a7206ea-4be5-4236-954e-a927f9cba993
# ╟─955d74f0-2ee9-40e9-9275-01555e13daa2
# ╟─8d1452f9-032c-41df-830e-ab6459e7c5ae
# ╟─067b9f6b-6b57-40c0-ae5f-0465400f77bb
# ╠═7f4d43df-f51e-49c2-b47e-629e3ccedc91
# ╠═39886e68-c4ea-4b96-96fb-ba46eeb54ecd
# ╠═fdad40dc-5e3c-4768-837f-7dbf4fca07ac
# ╟─565fab58-ef10-43e8-8ba8-e2e792ef872f
# ╠═09bfb23a-db24-42c8-8c34-8234da9c0d56
# ╠═5c1bea23-972b-4867-a424-9119cbc5edbc
# ╟─65131928-ad9a-4649-9600-e663ab5a7a28
# ╟─c4bc4b99-7adc-48a2-a3b9-114d6c5a4e3a
# ╟─243c0c2b-e84d-47eb-aaf8-5fb79d76d04d
# ╠═c6cd7f7d-d77b-4136-9485-913f962be85a
# ╠═7f8756e1-36e0-444a-a9a9-e20f763298af
# ╠═ec955132-42a3-40cd-a031-480ee6e95545

max_stars_repo_name: MolloiLab/CAC-stanford-data

max_issues_repo_name: MolloiLab/CAC-stanford-data

max_forks_repo_name: MolloiLab/CAC-stanford-data


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

using StatsBase
using DataStructures

" Load raw corpus into `Vector{String}` "
function loaddata(fname::AbstractString)
    corpus = Vector{String}()
    open(fname, "r") do fin
        for line in eachline(fin)
            length(line) > 10 && push!(corpus, strip(line))
        end
    end
    corpus
end

" Segment words in corpus "
function segment!(corpus::Vector{<:AbstractString})
    newcorp = Vector{Vector{String}}()
    for k in 1:length(corpus)
        push!(newcorp, string.(collect(corpus[k])))
    end
    newcorp
end

" Generate vocabulary from corpus "
function genvocab(corpus::Vector{Vector{String}}, mincount::Integer = 0)
    wordcount = Dict{String,Int64}()
    for k in 1:length(corpus)
        addcounts!(wordcount, corpus[k])
    end
    totalcount = (mincount <= 0) ? sum(values(wordcount)) :
                 sum(values(filter!(kv -> kv.second >= mincount, wordcount)))
    wordfreq = OrderedDict(map(kv -> (kv.first, kv.second / totalcount), collect(wordcount)))
    wordfreq
end


text:
# helper function
axes_min_length(idxs) = argmin([a isa Colon ? Inf : length(a) for a in idxs])

# parses FITS indices into standard Julian indices
function fits_indices(string::String)
    str = replace(string, r"[\[\]\s]" => "")
    tokens = split(str, ',')

    idxs = map(tokens) do token
        t = split(token, ':', keepempty=false)
        length(t) == 0 ? Colon() : parse(Int, t[1]):parse(Int, t[2])
    end
    return reverse(idxs)
end

# Convert using `round` for integers
convert_value(S::Type{<:Integer}, x) = round(S, x)
convert_value(S, x) = convert(S, x)

#-------------------------------------------------------------------------------
"""
    subtract_bias!(frame, bias_frame; [hdu = 1])

In-place version of [`subtract_bias`](@ref)

# See Also
[`subtract_bias`](@ref)
"""
function subtract_bias!(frame::AbstractArray, bias_frame::AbstractArray)
    frame .-= bias_frame
    return frame
end

"""
    subtract_bias(frame, bias_frame; [hdu = 1])

Subtract the `bias_frame` from `frame`.

If either are strings, they will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` as either an integer or a tuple corresponding to each file.

# Examples
```jldoctest
julia> frame = [1.0 2.2 3.3 4.5];

julia> bias = [0.0 0.2 0.3 0.5];

julia> subtract_bias(frame, bias)
1×4 Matrix{Float64}:
 1.0  2.0  3.0  4.0

```

# See Also
[`subtract_bias!`](@ref)
"""
subtract_bias(frame::AbstractArray, bias_frame::AbstractArray) = subtract_bias!(deepcopy(frame), bias_frame)

"""
    subtract_overscan!(frame, idxs; dims = axes_min_length(idxs))

In-place version of [`subtract_overscan`](@ref)

# See Also
[`subtract_overscan`](@ref)
"""
function subtract_overscan!(frame::AbstractArray{T}, idxs; dims = axes_min_length(idxs)) where T
    overscan_region = @view frame[idxs...]
    overscan_value = convert_value.(T, median(overscan_region, dims = dims))
    frame .-= overscan_value
    return frame
end

subtract_overscan!(frame::AbstractArray, idxs::String; kwargs...) = subtract_overscan!(frame, fits_indices(idxs); kwargs...)

"""
    subtract_overscan(frame, idxs; dims = axes_min_length(idxs), [hdu = 1])

Subtract the overscan frame from image.

`dims` is the dimension along which `overscan_frame` is combined. The default value
of `dims` is the axis with smaller length in overscan region. If `idxs` is a string it will be parsed as FITS-style indices.

If `frame` is a string, it will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` which by default is 1.

# Examples
```jldoctest
julia> frame = [4.0 2.0 3.0 1.0 1.0];

julia> subtract_overscan(frame, (:, 4:5), dims = 2)
1×5 Matrix{Float64}:
 3.0  1.0  2.0  0.0  0.0

julia> subtract_overscan(frame, "[4:5, 1:1]", dims = 2)
1×5 Matrix{Float64}:
 3.0  1.0  2.0  0.0  0.0

```

# See Also
[`subtract_overscan!`](@ref)
"""
subtract_overscan(frame, idxs; kwargs...) = subtract_overscan!(deepcopy(frame), idxs; kwargs...)

"""
    flat_correct!(frame, flat_frame; norm_value = mean(flat_frame), [hdu = 1])

In-place version of [`flat_correct`](@ref)

# See Also
[`flat_correct`](@ref)
"""
function flat_correct!(frame::AbstractArray, flat_frame::AbstractArray; norm_value = mean(flat_frame))
    norm_value <= 0 && error("norm_value must be positive")
    frame ./= (flat_frame ./ norm_value)
    return frame
end

"""
    flat_correct(frame, flat_frame; norm_value = mean(flat_frame), [hdu = 1])

Correct `frame` for non-uniformity using the calibrated `flat_frame`.

By default, the `flat_frame` is normalized by its mean, but this can be changed by providing a custom `norm_value`.

If either are strings, they will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` as either an integer or a tuple corresponding to each file.

!!! note
    This function may introduce non-finite values if `flat_frame` contains values very close to `0` due to dividing by zero.
    The default behavior will return `Inf` if the frame value is non-zero, and `Nan` if the frame value is `0`.

# Examples
```jldoctest
julia> frame = ones(3, 3);

julia> flat = fill(2.0, (3, 3));

julia> flat_correct(frame, flat, norm_value = 1.0)
3×3 Matrix{Float64}:
 0.5  0.5  0.5
 0.5  0.5  0.5
 0.5  0.5  0.5

julia> flat_correct(frame, flat)
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

```

# See Also
[`flat_correct!`](@ref)
"""
function flat_correct(frame::AbstractArray{T}, flat_frame::AbstractArray{S}; kwargs...) where {T, S}
    V = float(promote_type(T, S))
    return flat_correct!(V.(frame), flat_frame; kwargs...)
end

"""
    trim(frame, idxs; [hdu = 1])

Trims the `frame` to remove the region specified by idxs.

This function trims the array in a manner such that final array should be rectangular.
The indices follow standard Julia convention, so `(:, 45:60)` trims all columns from 45 to 60 and `(1:20, :)` trims all the rows from 1 to 20.
The function also supports FITS-style indices.

If `frame` is a string, it will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` which by default is 1.

# Examples
```jldoctest
julia> frame = ones(5, 5);

julia> trim(frame, (:, 2:5))
5×1 Matrix{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0

julia> trim(frame, "[2:5, 1:5]")
5×1 Matrix{Float64}:
 1.0
 1.0
 1.0
 1.0
 1.0

```

# See Also
[`trimview`](@ref)
"""
trim(frame, idxs) = copy(trimview(frame, idxs))

"""
    trimview(frame, idxs)

Trims the `frame` to remove the region specified by idxs.

This function is same as the [`trim`](@ref) function but returns a view of the frame.

!!! note
    This function returns a view of the frame, so any modification to output
    array will result in modification of frame.

# See Also
[`trim`](@ref)
"""
function trimview(frame::AbstractArray, idxs)
    # this adds the support for input indices of the form (1:size(frame, 1), ...) or (..., 1:size(frame, 2))
    # It converts 1:size(frame, 1) to : and then the same subroutine follows.
    processed_idxs = map(axes(frame), idxs) do a1, a2
        (a2 isa Colon || a1 == a2) ? Colon() : a2
    end

    # can switch to using `only` for Julia v1.4+
    ds = findall(x -> !isa(x, Colon), processed_idxs)
    length(ds) == 1 || error("Invalid trim indices $idxs")

    d = ds[1]
    full_idxs = axes(frame, d)
    # checking bounds error
    processed_idxs[d] ⊆ full_idxs || error("Trim indices $(idxs[d]) out of bounds for frame dimension $d $(full_idxs)")

    # finding the complement indices
    complement_idxs = setdiff(full_idxs, processed_idxs[d])

    return selectdim(frame, d, complement_idxs)
end

trimview(frame::AbstractArray, idxs::String) = trimview(frame, fits_indices(idxs))

"""
    crop(frame, shape; force_equal = true, [hdu = 1])

Crops `frame` to the size specified by `shape` anchored by the frame center.

This will remove rows/cols of the `frame` equally on each side. When there is an uneven difference in sizes (e.g. size 9 -> 6 can't be removed equally) the default is to
increase the output size (e.g. 6 -> 7) so there is equal removal on each side. To disable this, set `force_equal=false`, which will remove the extra slice from the end of the axis.

If `frame` is a string, it will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` which by default is 1.

# Examples
```jldoctest
julia> frame = reshape(1:25, (5, 5));

julia> crop(frame, (3, 3))
3×3 Matrix{Int64}:
 7  12  17
 8  13  18
 9  14  19

julia> crop(frame, (4, 3), force_equal = false)
4×3 Matrix{Int64}:
 6  11  16
 7  12  17
 8  13  18
 9  14  19

```

# See Also
[`cropview`](@ref)
"""
crop(frame, shape; kwargs...) = copy(cropview(frame, shape; kwargs...))

"""
    cropview(frame, shape; force_equal = true)

Crops `frame` to the size specified by `shape` anchored by the frame center.

This function is same as the [`crop`](@ref) function but returns a view of the frame.

!!! note
    This function returns a view of the frame, so any modification to output
    array will result in modification of frame.

# See Also
[`crop`](@ref)
"""
function cropview(frame::AbstractArray, shape; force_equal = true)
    # testing error
    ndims(frame) == length(shape) || error("Dimension mismatch between frame and shape")
    any(s -> !isa(s, Colon) && s < 1, shape) && error("crop size $shape cant't be less than 1")

    # generating idxs for cropped frame
    idxs = map(enumerate(size(frame)), shape) do (d, s1), s2
                diff = s2 isa Colon ? 0 : s1 - s2
                lower = iseven(diff) ? diff ÷ 2 : (diff - 1) ÷ 2
                upper = if isodd(diff) && force_equal
                            @warn "dimension $d changed from $s2 to $(s2 + 1)"
                            (diff - 1) ÷ 2
                        elseif isodd(diff)
                            (diff + 1) ÷ 2
                        else
                            diff ÷ 2
                        end
                1 + lower:s1 - upper
            end

    # returning the view
    return @view frame[idxs...]
end

"""
    combine(frames...; method = median, [hdu = 1], [header_hdu = 1])
    combine(frames; method = median, [hdu = 1], [header_hdu = 1])

Combine multiple frames using `method`. Multiple frames can also be passed in a vector or as generators for combining.

To pass a custom method, it must have a signature like `method(::AbstractArray; dims)`.

If `frames` are strings, they will be loaded into [`CCDData`](@ref)s first. The HDU indices can be specified with `hdu` as either an integer or a tuple corresponding to each file.

Header of output file (if applicable) is specified by `header_hdu` which by default is 1.

# Examples
```jldoctest
julia> frame = [reshape(1.0:4.0, (2, 2)) for i = 1:4];

julia> combine(frame)
2×2 Matrix{Float64}:
 1.0  3.0
 2.0  4.0

julia> combine(frame, method = sum)
2×2 Matrix{Float64}:
 4.0  12.0
 8.0  16.0

```
"""
function combine(frames::Vararg{<:AbstractArray{<:Number}}; method = median)
    firstframe = first(frames)
    dim = ndims(firstframe) + 1
    shape = size(firstframe)
    return reshape(method(LazyStack.stack(frames...), dims = dim), shape)
end

combine(frames; kwargs...) = combine(frames...; kwargs...)

"""
    subtract_dark!(frame, dark_frame; data_exposure = 1, dark_exposure = 1, [hdu = 1])

In-place version of [`subtract_dark`](@ref)

# See Also
[`subtract_dark`](@ref)
"""
function subtract_dark!(frame::AbstractArray, dark_frame::AbstractArray; data_exposure = 1, dark_exposure = 1)
    factor = data_exposure / dark_exposure
    @. frame -= (dark_frame * factor)
    return frame
end

"""
    subtract_dark(frame, dark_frame; data_exposure = 1, dark_exposure = 1, [hdu = 1])

Subtract the `dark_frame` from `frame`.

If either are strings, they will be loaded into [`CCDData`](@ref) first. The HDU loaded can be specified by `hdu` as either an integer or a tuple corresponding to each file.

# Examples
```jldoctest
julia> frame = ones(3, 3);

julia> dark_frame = ones(3, 3);

julia> subtract_dark(frame, dark_frame)
3×3 Matrix{Float64}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> subtract_dark(frame, dark_frame, data_exposure = 1, dark_exposure = 4)
3×3 Matrix{Float64}:
 0.75  0.75  0.75
 0.75  0.75  0.75
 0.75  0.75  0.75

```

# See Also
[`subtract_dark!`](@ref)
"""
function subtract_dark(frame::AbstractArray{T}, dark_frame::AbstractArray{S}; kwargs...) where {T, S}
    V = float(promote_type(T, S))
    return subtract_dark!(V.(frame), dark_frame; kwargs...)
end


text:
#Methods some Some can broadcast
getindex(s::Some) = s.value
getindex(s::Some, ::CartesianIndex{0}) = s.value

iterate(s::Some) = (s.value, nothing)
iterate( ::Some, s) = nothing

ndims(::Some) = 0
ndims(::Type{<:Some}) = 0

length(::Some) = 1
size(::Some) = ()
axes(::Some) = ()

IteratorSize(::Type{<:Some}) = HasShape{0}()
Broadcast.broadcastable(s::Some) = s

eltype(::Type{Some{T}}) where {T} = T

max_stars_repo_name: MasonProtter/julia

max_issues_repo_name: MasonProtter/julia

max_forks_repo_name: MasonProtter/julia


text:
using Revise
using ADCME
using AdFem
using SparseArrays
using PyPlot

function plot_disp(u, m, n, h)
    close("all")
    u = reshape(u, m+1, n+1)
    x = collect((0:m) * h)
    y = collect((0:n) * h)
    pcolormesh(x, y, u')
    colorbar()
    axis("scaled")
end

function plot_disp_line(u, m, n, h)
   close("all")
    u = reshape(u, m+1, n+1)
   plot(u[:, n÷2], label="middle")
   plot(u[:, n÷3], label="upper")
   plot(u[:, n÷3*2], label="lower")
   legend()
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
#i64o-constants.jl

#various superint constants and ways of generating them
Base.zero{FSS}(::Type{ArrayNum{FSS}}) = ArrayNum{FSS}(zeros(UInt64, __cell_length(FSS)))
Base.zero{FSS}(a::ArrayNum{FSS}) = ArrayNum{FSS}(zeros(UInt64, __cell_length(FSS)))

doc"""
`zero!` sets an `ArrayNum` to zero.
"""
function zero!{FSS}(a::ArrayNum{FSS})
  #unroll the loop that fills the contents of the a with zero.
  for idx = 1:__cell_length(FSS)
    @inbounds a.a[idx] = 0
  end
end

##create Unums.frac_zero!, out of zero(UInt64) and zero!(ArrayNum)
doc"""
`Unums.frac_zero!` sets the fraction of a unum to zero.
"""
@fracproc zero

function Base.one{FSS}(::Type{ArrayNum{FSS}})
  l = __cell_length(FSS)
  arr = zeros(UInt64, l)
  @inbounds arr[l] = 1
  ArrayNum{FSS}(arr)
end

function Base.one{FSS}(a::ArrayNum{FSS})
  one(ArrayNum{FSS})
end

top(::Type{UInt64}) = t64
function top{FSS}(::Type{ArrayNum{FSS}})
  l = __cell_length(FSS)
  arr = zeros(UInt64, l)
  @inbounds arr[1] = t64
  ArrayNum{FSS}(arr)
end

top(n::UInt64) = t64
function top{FSS}(a::ArrayNum{FSS})
  top(ArrayNum{FSS})
end
function top!{FSS}(a::ArrayNum{FSS})
  zero!(a)
  @inbounds a[1] = t64
end
export top

doc"""`Unums.frac_top!` sets the fraction of a unum to be zero, except with the top bit set."""
@fracproc top

doc"""`Unums.all(x::UInt64)` returns a UInt64 or an ArrayNum with all bits set."""
all(x::UInt64) = f64
all{FSS}(::Type{ArrayNum{FSS}}) = ArrayNum{FSS}([f64 for idx = 1:__cell_length(FSS)])
doc"""`Unums.all!(x::UInt64)` sets all bits on an arraynum."""
all!{FSS}(x::ArrayNum{FSS}) = for idx = 1:__cell_length(FSS); @inbounds x.a[idx] = f64; end

doc"""`Unums.frac_all!` sets the fraction of a unum to be all set."""
@fracproc all

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


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

module Types

export VersionInterval, VersionSet, Requires, Available, Fixed, merge_requires!, satisfies,
       ResolveBacktraceItem, ResolveBacktrace
import Base: show, isempty, in, intersect, union!, union, ==, hash, copy, deepcopy_internal, push!

struct VersionInterval
    lower::VersionNumber
    upper::VersionNumber
end
VersionInterval(lower::VersionNumber) = VersionInterval(lower,typemax(VersionNumber))
VersionInterval() = VersionInterval(typemin(VersionNumber))

show(io::IO, i::VersionInterval) = print(io, "[$(i.lower),$(i.upper))")
isempty(i::VersionInterval) = i.upper <= i.lower
in(v::VersionNumber, i::VersionInterval) = i.lower <= v < i.upper
intersect(a::VersionInterval, b::VersionInterval) = VersionInterval(max(a.lower,b.lower), min(a.upper,b.upper))
==(a::VersionInterval, b::VersionInterval) = a.lower == b.lower && a.upper == b.upper
hash(i::VersionInterval, h::UInt) = hash((i.lower, i.upper), h + (0x0f870a92db508386 % UInt))

function normalize!(ivals::Vector{VersionInterval})
    # VersionSet internal normalization:
    # removes empty intervals and fuses intervals without gaps
    # e.g.:
    #     [0.0.0,1.0.0) ∪ [1.0.0,1.5.0) ∪ [1.6.0,1.6.0) ∪ [2.0.0,∞)
    # becomes:
    #     [0.0.0,1.5.0) ∪ [2.0.0,∞)
    # (still assumes that lower bounds are sorted, and intervals do
    # not overlap)
    l = length(ivals)
    l == 0 && return ivals

    lo, up, k0 = ivals[1].lower, ivals[1].upper, 1
    fusing = false
    for k = 2:l
        lo1, up1 = ivals[k].lower, ivals[k].upper
        if lo1 == up
            up = up1
            fusing = true
            continue
        end
        if lo < up
            # The only purpose of the "fusing" check is to avoid
            # extra allocations
            ivals[k0] = fusing ? VersionInterval(lo, up) : ivals[k-1]
            k0 += 1
        end
        fusing = false
        lo, up = lo1, up1
    end
    if lo < up
        ivals[k0] = fusing ? VersionInterval(lo, up) : ivals[l]
        k0 += 1
    end
    resize!(ivals, k0 - 1)
    return ivals
end

struct VersionSet
    intervals::Vector{VersionInterval}
    VersionSet(intervals::Vector{VersionInterval}) = new(normalize!(intervals))
    # copy is defined inside the struct block to call `new` directly
    # without going through `normalize!`
    Base.copy(vset::VersionSet) = new(copy(vset.intervals))
end
function VersionSet(versions::Vector{VersionNumber})
    intervals = VersionInterval[]
    if isempty(versions)
        push!(intervals, VersionInterval())
    else
        isodd(length(versions)) && push!(versions, typemax(VersionNumber))
        while !isempty(versions)
            push!(intervals, VersionInterval(shift!(versions), shift!(versions)))
        end
    end
    VersionSet(intervals)
end
VersionSet(versions::VersionNumber...) = VersionSet(VersionNumber[versions...])

const empty_versionset = VersionSet(VersionInterval[])

# Windows console doesn't like Unicode
const _empty_symbol = @static is_windows() ? "empty" : "∅"
const _union_symbol = @static is_windows() ? " or " : " ∪ "
show(io::IO, s::VersionSet) = isempty(s) ? print(io, _empty_symbol) :
                                           join(io, s.intervals, _union_symbol)
isempty(s::VersionSet) = all(isempty, s.intervals)
in(v::VersionNumber, s::VersionSet) = any(i->in(v,i), s.intervals)
function intersect(A::VersionSet, B::VersionSet)
    (isempty(A) || isempty(B)) && return copy(empty_versionset)
    ivals = [intersect(a,b) for a in A.intervals for b in B.intervals]
    sort!(ivals, by=i->i.lower)
    VersionSet(ivals)
end

union(A::VersionSet, B::VersionSet) = union!(copy(A), B)
function union!(A::VersionSet, B::VersionSet)
    A == B && return A
    ivals = A.intervals
    for intB in B.intervals
        lB, uB = intB.lower, intB.upper
        k0 = findfirst(i->(i.upper > lB), ivals)
        if k0 == 0
            push!(ivals, intB)
            continue
        end
        lB = min(lB, ivals[k0].lower)
        for k1 = k0:length(ivals)
            intA = ivals[k1]
            if uB < intA.lower
                splice!(ivals, k0:(k1-1), (VersionInterval(lB, uB),))
                break
            elseif uB ∈ intA || k1 == length(ivals)
                splice!(ivals, k0:k1, (VersionInterval(lB, max(uB, intA.upper)),))
                break
            end
        end
    end
    normalize!(ivals)
    return A
end

==(A::VersionSet, B::VersionSet) = A.intervals == B.intervals
hash(s::VersionSet, h::UInt) = hash(s.intervals, h + (0x2fd2ca6efa023f44 % UInt))
deepcopy_internal(vs::VersionSet, ::ObjectIdDict) = copy(vs)

const Requires = Dict{String,VersionSet}

function merge_requires!(A::Requires, B::Requires)
    for (pkg,vers) in B
        A[pkg] = haskey(A,pkg) ? intersect(A[pkg],vers) : vers
    end
    return A
end

satisfies(pkg::AbstractString, ver::VersionNumber, reqs::Requires) =
    !haskey(reqs, pkg) || in(ver, reqs[pkg])

struct Available
    sha1::String
    requires::Requires
end

==(a::Available, b::Available) = a.sha1 == b.sha1 && a.requires == b.requires
hash(a::Available, h::UInt) = hash((a.sha1, a.requires), h + (0xbc8ae0de9d11d972 % UInt))
copy(a::Available) = Available(a.sha1, copy(a.requires))

show(io::IO, a::Available) = isempty(a.requires) ?
    print(io, "Available(", repr(a.sha1), ")") :
    print(io, "Available(", repr(a.sha1), ",", a.requires, ")")

struct Fixed
    version::VersionNumber
    requires::Requires
end
Fixed(v::VersionNumber) = Fixed(v,Requires())

==(a::Fixed, b::Fixed) = a.version == b.version && a.requires == b.requires
hash(f::Fixed, h::UInt) = hash((f.version, f.requires), h + (0x68628b809fd417ca % UInt))

show(io::IO, f::Fixed) = isempty(f.requires) ?
    print(io, "Fixed(", repr(f.version), ")") :
    print(io, "Fixed(", repr(f.version), ",", f.requires, ")")

# TODO: Available & Fixed are almost the same – merge them?
# Free could include the same information too, it just isn't
# required by anything that processes these things.

const VersionReq = Union{VersionNumber,VersionSet}
const WhyReq = Tuple{VersionReq,Any}

# This is used to keep track of dependency relations when propagating
# requirements, so as to emit useful information in case of unsatisfiable
# conditions.
# The `versionreq` field keeps track of the remaining allowed versions,
# intersecting all requirements.
# The `why` field is a Vector which keeps track of the requirements. Each
# entry is a Tuple of two elements:
# 1) the first element is the version requirement (can be a single VersionNumber
#    or a VersionSet).
# 2) the second element can be either :fixed (for requirements induced by
#    fixed packages), :required (for requirements induced by explicitly
#    required packages), or a Pair p=>backtrace_item (for requirements induced
#    indirectly, where `p` is the package name and `backtrace_item` is
#    another ResolveBacktraceItem.
mutable struct ResolveBacktraceItem
    versionreq::VersionReq
    why::Vector{WhyReq}
    ResolveBacktraceItem() = new(VersionSet(), WhyReq[])
    ResolveBacktraceItem(reason, versionreq::VersionReq) = new(versionreq, WhyReq[(versionreq,reason)])
end

function push!(ritem::ResolveBacktraceItem, reason, versionset::VersionSet)
    if isa(ritem.versionreq, VersionSet)
        ritem.versionreq = ritem.versionreq ∩ versionset
    elseif ritem.versionreq ∉ versionset
        ritem.versionreq = copy(empty_versionset)
    end
    push!(ritem.why, (versionset,reason))
end

function push!(ritem::ResolveBacktraceItem, reason, version::VersionNumber)
    if isa(ritem.versionreq, VersionSet)
        if version ∈ ritem.versionreq
            ritem.versionreq = version
        else
            ritem.versionreq = copy(empty_versionset)
        end
    elseif ritem.versionreq ≠ version
        ritem.versionreq = copy(empty_versionset)
    end
    push!(ritem.why, (version,reason))
end

show(io::IO, ritem::ResolveBacktraceItem) = _show(io, ritem, "", Set{ResolveBacktraceItem}([ritem]))

function _show(io::IO, ritem::ResolveBacktraceItem, indent::String, seen::Set{ResolveBacktraceItem})
    l = length(ritem.why)
    for (i,(vs,w)) in enumerate(ritem.why)
        print(io, indent, (i==l ? '└' : '├'), '─')
        if w ≡ :fixed
            @assert isa(vs, VersionNumber)
            println(io, "version $vs set by fixed requirement (package is checked out, dirty or pinned)")
        elseif w ≡ :required
            @assert isa(vs, VersionSet)
            println(io, "version range $vs set by an explicit requirement")
        else
            @assert isa(w, Pair{<:AbstractString,ResolveBacktraceItem})
            if isa(vs, VersionNumber)
                print(io, "version $vs ")
            else
                print(io, "version range $vs ")
            end
            print(io, "required by package $(w[1]), ")
            if isa(w[2].versionreq, VersionSet)
                println(io, "whose allowed version range is $(w[2].versionreq):")
            else
                println(io, "whose only allowed version is $(w[2].versionreq):")
            end
            if w[2] ∈ seen
                println(io, (i==l ? "  " : "│ ") * indent, "└─[see above for $(w[1]) backtrace]")
                continue
            end
            push!(seen, w[2])
            _show(io, w[2], (i==l ? "  " : "│ ") * indent, seen)
        end
    end
end

const ResolveBacktrace = Dict{AbstractString,ResolveBacktraceItem}

end # module