text:
using AxisIndices: OffsetArray, OffsetVector

@testset "OffsetAxis" begin
    inds = Base.OneTo(3)
    ks = 1:3
    offset = 0
    @test @inferred(OffsetAxis(ks)) == OffsetAxis(1:3, Base.OneTo(3))
    @test @inferred(OffsetAxis(ks, inds)) == OffsetAxis(1:3, Base.OneTo(3))
    @test @inferred(OffsetAxis(offset, inds)) == OffsetAxis(1:3, Base.OneTo(3))
    @test @inferred(OffsetAxis(OffsetAxis(ks))) == OffsetAxis(1:3, Base.OneTo(3))

    @test @inferred(OffsetAxis{Int16}(ks)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(3))
    @test @inferred(OffsetAxis{Int16}(ks, inds)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(3))
    @test @inferred(OffsetAxis{Int16}(offset, inds)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(3))
    @test @inferred(OffsetAxis{Int16}(OffsetAxis(ks))) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(3))

    @test @inferred(OffsetAxis{Int16}(ks)) == OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))
    @test @inferred(OffsetAxis{Int16}(ks, inds)) == OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))
    @test @inferred(OffsetAxis{Int16}(offset, inds)) == OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))
    @test @inferred(OffsetAxis{Int16}(OffsetAxis(ks))) == OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))

    @test @inferred(OffsetAxis{Int16,SimpleAxis{Int16,Base.OneTo{Int16}}}(ks)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))
    @test @inferred(OffsetAxis{Int16,SimpleAxis{Int16,Base.OneTo{Int16}}}(ks, inds)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))
    @test @inferred(OffsetAxis{Int16,SimpleAxis{Int16,Base.OneTo{Int16}}}(offset, inds)) ==
        OffsetAxis(Int16(1):Int16(3), Base.OneTo(Int16(3)))

    @test OffsetArray(OffsetArray(ones(2, 2))) isa OffsetArray
    #=
    function same_value(r1, r2)
        length(r1) == length(r2) || return false
        for (v1, v2) in zip(r1, r2)
            v1 == v2 || return false
        end
        return true
    end

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

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

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

    # @test_throws ArgumentError oftype(ro, 3:4)
    #
    @testset "values -> keys" begin
        # firstindex(keys(axis)) == 3, firstindex(axis) == 2
        @test @inferred(AxisIndices.v2k(Axis(OffsetAxis(2, 1:10), OffsetAxis(1, 1:10)), 2)) == 3
        # firstindex(keys(axis)) == 3, firstindex(axis) == 1
        @test @inferred(AxisIndices.v2k(Axis(OffsetAxis(2, 1:10), 1:10), 1)) == 3
        # firstindex(keys(axis)) == 1, firstindex(axis) == 2
        @test @inferred(AxisIndices.v2k(Axis(1:10, OffsetAxis(1, 1:10)), 2)) == 1
        # firstindex(keys(axis)) == 1, firstindex(axis) == 1
        @test @inferred(AxisIndices.v2k(Axis(1:10, 1:10), 1)) == 1

        # firstindex(keys(axis)) == 3, firstindex(axis) == 2
        @test @inferred(AxisIndices.v2k(Axis(OffsetAxis(2, 1:10), OffsetAxis(1, 1:10)), 2:3)) == 3:4
        # firstindex(keys(axis)) == 3, firstindex(axis) == 1
        @test @inferred(AxisIndices.v2k(Axis(OffsetAxis(2, 1:10), 1:10), 1:2)) == 3:4
        # firstindex(keys(axis)) == 1, firstindex(axis) == 2
        @test @inferred(AxisIndices.v2k(Axis(1:10, OffsetAxis(1, 1:10)), 2:3)) == 1:2
        # firstindex(keys(axis)) == 1, firstindex(axis) == 1
        @test @inferred(AxisIndices.v2k(Axis(1:10, 1:10), 1:2)) == 1:2
    end

    @testset "keys -> values" begin
        # firstindex(keys(axis)) == 3, firstindex(axis) == 2
        @test @inferred(AxisIndices.k2v(Axis(OffsetAxis(2, 1:10), OffsetAxis(1, 1:10)), 3)) == 2
        # firstindex(keys(axis)) == 3, firstindex(axis) == 1
        @test @inferred(AxisIndices.k2v(Axis(OffsetAxis(2, 1:10), 1:10), 3)) == 1
        # firstindex(keys(axis)) == 1, firstindex(axis) == 2
        @test @inferred(AxisIndices.k2v(Axis(1:10, OffsetAxis(1, 1:10)), 1)) == 2
        # firstindex(keys(axis)) == 1, firstindex(axis) == 1
        @test @inferred(AxisIndices.k2v(Axis(1:10, 1:10), 1)) == 1

        # firstindex(keys(axis)) == 3, firstindex(axis) == 2
        @test @inferred(AxisIndices.k2v(Axis(OffsetAxis(2, 1:10), OffsetAxis(1, 1:10)), 3:4)) == 2:3
        # firstindex(keys(axis)) == 3, firstindex(axis) == 1
        @test @inferred(AxisIndices.k2v(Axis(OffsetAxis(2, 1:10), 1:10), 3:4)) == 1:2
        # firstindex(keys(axis)) == 1, firstindex(axis) == 2
        @test @inferred(AxisIndices.k2v(Axis(1:10, OffsetAxis(1, 1:10)), 1:2)) == 2:3
        # firstindex(keys(axis)) == 1, firstindex(axis) == 1
        @test @inferred(AxisIndices.k2v(Axis(1:10, 1:10), 1:2)) == 1:2
    end
    =#

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

using Base: IdentityUnitRange
# differences from OffsetArrays
# * indexing preserves the offset so things like:
#!!! note this is different than what OffsetArryas does because offset arrays produced 0:1
# which doesn't make sense given the axis doesn't have a length of 2
#    @test S[0] == 3
#    @test axes(S) == (0:0,)

#=@testset "offset" begin
    @test @inferred(offset(OneTo(10))) == 0
    @test @inferred(offset(2:3)) == -1
end=#

getoffset(axis::AbstractAxis) = first(keys(axis)) - first(parent(axis))

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

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

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

    y = OffsetArray{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)
    @test eltype(y) === Float64

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

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

@testset "OffsetArray of OffsetArray construction" begin
    # guarantee no unnecessary nesting of `OffsetArray`s
    r = -2:5
    d = collect(r)
    y = OffsetArray(d, r)
    @test axes(y, 1) == -2:5

    # range constructor
    y0 = OffsetArray(y, 0:7)
    @test axes(y0, 1) == 0:7
    @test y0[0] == r[1]
    @test typeof(parent(y0)) <: Array

    # FIXME
    #= offset constructor
    y1 = OffsetArray(y, +2)
    @test y1[0] == r[1]
    @test typeof(parent(y1)) <: Array
    =#
end

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

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

    @test @inferred(A[0,3]) ==
          @inferred(A[0,3,1]) ==
          @inferred(A[1]) ==
          @inferred(S[0,3]) ==
          @inferred(S[0,3,1]) ==
          @inferred(S[1]) == 1
    @test A[1,3] == A[1,3,1] == A[2] == S[1,3] == S[1,3,1] == S[2] == 2
    @test A[0,4] == A[0,4,1] == A[3] == S[0,4] == S[0,4,1] == S[3] == 3
    @test A[1,4] == A[1,4,1] == A[4] == S[1,4] == S[1,4,1] == S[4] == 4
    @test @inbounds(A[0,3]) == @inbounds(A[0,3,1]) == @inbounds(A[1]) == @inbounds(S[0,3]) == @inbounds(S[0,3,1]) == @inbounds(S[1]) == 1
    @test @inbounds(A[1,3]) == @inbounds(A[1,3,1]) == @inbounds(A[2]) == @inbounds(S[1,3]) == @inbounds(S[1,3,1]) == @inbounds(S[2]) == 2
    @test @inbounds(A[0,4]) == @inbounds(A[0,4,1]) == @inbounds(A[3]) == @inbounds(S[0,4]) == @inbounds(S[0,4,1]) == @inbounds(S[3]) == 3
    @test @inbounds(A[1,4]) == @inbounds(A[1,4,1]) == @inbounds(A[4]) == @inbounds(S[1,4]) == @inbounds(S[1,4,1]) == @inbounds(S[4]) == 4
    @test_throws BoundsError A[1,1]
    @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

    y = OffsetArray{Float64}(undef, -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 "Vector indexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

    @test A[:, 3] == S[:, 3] == OffsetArray([1,2], (getoffset(axes(A, 1)),))
    @test A[:, 4] == S[:, 4] == OffsetArray([3,4], (getoffset(axes(A, 1)),))
    @test_throws BoundsError A[:, 1]
    @test_throws BoundsError S[:, 1]
    @test A[0, :] == S[0, :] == OffsetArray([1,3], (getoffset(axes(A, 2)),))
    @test A[1, :] == S[1, :] == OffsetArray([2,4], (getoffset(axes(A, 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 "CartesianIndexing" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))
    S = OffsetArray(view(A0, 1:2, 1:2), (-1,2))

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

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

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

    S = view(A, :, 3)
    @test S == OffsetArray([1,2], (getoffset(axes(A, 1)),))
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[2]
    @test axes(S) == (IdentityUnitRange(0:1),)
    S = view(A, 0, :)
    @test S == OffsetArray([1,3], (getoffset(axes(A, 2)),))
    @test S[3] == 1
    @test S[4] == 3
    @test_throws BoundsError S[1]
    @test axes(S) == (IdentityUnitRange(3:4),)

    # NOTE: These differ from OffsetArrays in that offset axes are preserved
    @testset "preserve offset axes" begin
        S = view(A, 0:0, 4)
        @test S == [3]
        @test S[0] == 3
        @test_throws BoundsError S[1]
        @test axes(S) == (OffsetAxis(0:0, Base.OneTo(1)),)
        S = view(A, 1, 3:4)
        @test S == [2,4]
        @test S[3] == 2
        @test S[4] == 4
        @test_throws BoundsError S[1]
        @test axes(S) == (OffsetAxis(3:4, Base.OneTo(2)),)
    end

    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) == (OffsetAxis(0:1, Base.OneTo(2)), OffsetAxis(3:4, Base.OneTo(2)))  # OffsetArrays uses == IdentityUnitRange.((0:1, 3:4))
    S = view(A, axes(A)...)
    @test S == A
    @test S[0,3] == S[1] == 1
    @test S[1,3] == S[2] == 2
    @test S[0,4] == S[3] == 3
    @test S[1,4] == S[4] == 4
    @test_throws BoundsError S[1,1]
    @test axes(S) == (OffsetAxis(0:1, Base.OneTo(2)), OffsetAxis(3:4, Base.OneTo(2)))
    # issue 100
    S = view(A, axes(A, 1), 3)
    @test S == A[:, 3]
    @test S[0] == 1
    @test S[1] == 2
    @test_throws BoundsError S[length(S)]
    @test axes(S) == (OffsetAxis(0:1, Base.OneTo(2)),)
    # issue 100
    S = view(A, 1, axes(A, 2))
    @test S == A[1, :]
    @test S[3] == 2
    @test S[4] == 4
    @test_throws BoundsError S[1]
    @test axes(S) == (OffsetAxis(3:4, Base.OneTo(2)),)

    # this bit also differs from OffsetArrays
    A0 = collect(reshape(1:24, 2, 3, 4))
    A = OffsetArray(A0, (-1,2,1))
    S = view(A, axes(A, 1), 3:4, axes(A, 3))
    @test S == A[:, 3:4, :]
    @test S[0, 3, 2] == A[0, 3, 2]
    @test S[0, 4, 2] == A[0, 4, 2]
    @test S[1, 3, 2] == A[1, 3, 2]
    @test axes(S) == (OffsetAxis(0:1, Base.OneTo(2)), OffsetAxis(3:4, Base.OneTo(2)), OffsetAxis(2:5, Base.OneTo(4)))
end

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

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

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

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

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

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

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

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

    @test A.+1 == OffsetArray(parent(A).+1, map(getoffset, axes(A)))
    @test 2*A == OffsetArray(2*parent(A), map(getoffset, axes(A)))
    @test A+A == OffsetArray(parent(A)+parent(A), map(getoffset, axes(A)))
    @test A.*A == OffsetArray(parent(A).*parent(A), map(getoffset, axes(A)))
end

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

    @test rotl90(A) == OffsetArray(rotl90(parent(A)), reverse(map(getoffset, axes(A))))
    @test rotr90(A) == OffsetArray(rotr90(parent(A)), (getoffset(axes(A, 2)), getoffset(axes(A, 1))))
    @test reverse(A, dims = 1) == OffsetArray(reverse(parent(A), dims = 1), map(getoffset, axes(A)))
    @test reverse(A, dims = 2) == OffsetArray(reverse(parent(A), dims = 2), map(getoffset, axes(A)))
end

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

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

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

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

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

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

    b = 1:2    # copy between AbstractArrays
    bo = OffsetArray(1:2, (-3,))
    # why should this error
    #@test_throws BoundsError copyto!(a, b)
    fill!(a, -1)
    copyto!(a, bo)
    @test a[-3] == 1
    @test a[-2] == 2
    @test a[-1] == -1
    # the behavior of copyto! is corrected as the documentation says "first n element"
    # https://github.com/JuliaLang/julia/pull/34049
    fill!(a, -1)
    copyto!(a, bo)
    @test a[-3] == 1
    @test a[-2] == 2
    @test a[-1] == -1

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

# TODO it seems odd to me that indexing with an offset array would change the axes
# of the resulting returned array
#=
@testset "Vector indexing with offset ranges" begin
    r = OffsetArray(8:10, -1:1)
    r1 = r[0:1]
    @test r1 == 9:10  # FIXME for some reason this doesn't work with `===` like OffsetArrays.jl
    r1 = (8:10)[OffsetArray(1:2, -5:-4)]
    @test axes(r1) == (OffsetAxis(-5:-4),)
    @test parent(r1) === 8:9
    r1 = OffsetArray(8:10, -1:1)[OffsetArray(0:1, -5:-4)]
    @test axes(r1) == (IdentityUnitRange(-5:-4),)
    @test parent(r1) === 9:10
end
=#

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

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

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

    @test mapslices(v->sort(v), A, dims = 1) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 1), map(getoffset, axes(A)))
    @test mapslices(v->sort(v), A, dims = 2) == OffsetArray(mapslices(v->sort(v), parent(A), dims = 2), map(getoffset, axes(A)))
end

# TODO move this test to an appropriate file
@testset "fill" begin
    B = fill(5, OffsetAxis(1:3), OffsetAxis(-1:1))
    @test axes(B) == (1:3,-1:1)
    @test all(B.==5)
end

# TODO I have a different approach to this. The offset is preserved
#=
@testset "similar" begin
    A0 = [1 3; 2 4]
    A = OffsetArray(A0, (-1,2))

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

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

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

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

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

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

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

=#


text:
using OrdinaryDiffEq, DiffEqBase

function MLplant(n,x0::Vector,t::Vector,U::Matrix,t0::Float64,tf::Float64)
  spU = linearSpline(t,U[:,1])   # create spline
  f = (dx,x,p,t) -> begin        # diff eqs.
    dx[1] = x[2]
    dx[2] = spU[t] - 1.625
  end
  tspan = (t0,tf)
  prob = ODEProblem(f,x0,tspan)
  sol = DiffEqBase.solve(prob,Tsit5())
  U = [spU]
  return sol, U
end

```
integrationConfig=:lgrExplicit
interpolationOn=true
linearInterpolation=false
numInterpPts=250
AlltpolyPts=false
predictX0=false
maxSim=100
lastOptimal=false
tolOn=true
```
function MoonLander(integrationConfig;interpolationOn::Bool=true,linearInterpolation::Bool=false,numInterpPts::Int64=250,AlltpolyPts::Bool=false,predictX0::Bool=false,maxSim::Int64=100,lastOptimal::Bool=true,tolOn::Bool=false)
  n=define(numStates=2,numControls=1,X0=[10.,-2],XF=[0.,0.],CL=[0.],CU=[3.])
  n.s.ocp.numInterpPts = numInterpPts
  n.s.ocp.linearInterpolation = linearInterpolation
  n.s.ocp.interpolationOn = interpolationOn
  if tolOn
    X0_tol = [0.1,0.1]
    XF_tol = [0.1,0.1]
    defineTolerances!(n;X0_tol=X0_tol,XF_tol=XF_tol)
  end
  states!(n,[:h,:v];descriptions=["h(t)","v(t)"])
  controls!(n,[:T];descriptions=["T(t)"])
  dx=[:(v[j]),:(T[j]-1.625)]
  dynamics!(n,dx)
  configure!(n;(:integrationScheme=>integrationConfig),(:finalTimeDV=>true))
  obj=integrate!(n,:(T[j]))
  @NLobjective(n.ocp.mdl, Min, obj)
  initOpt!(n)
  defineMPC!(n;predictX0=predictX0,lastOptimal=lastOptimal,tex=0.1,printLevel=0,maxSim=maxSim)
  defineIP!(n,MLplant)
  simMPC!(n)

  if AlltpolyPts
    return (isequal(length(n.r.ocp.AlltpolyPts),length(n.r.ocp.dfs)) && n.f.mpc.goalReached)
  else
    return n.f.mpc.goalReached
  end
end

@testset "MoonLander MPC Tests" begin
  # Test #1
  @test MoonLander(:lgrExplicit)

  # Test #2
  @test MoonLander(:lgrExplicit;lastOptimal=false,tolOn=true)

  # Test #3
  # a) :ps methods with several infeasibe solutions, but gets to the goal.
      # make sure that data lines up for plots even with infeasibe solutions
  # b) If using :ps methods and there is an infeasibe solution and the lastOptimal==true
      # if you try to use the interpolateLagrange!() function an error like this will occur:
        # WARNING: Not solved to optimality, status: Infeasible
        # ERROR: LoadError: BoundsError: attempt to access 39-element Array{Float64,1} at index [31:41]
      # To mitigate this error, if it was infeasibe and lastOptimal==true then a linearInterpolation is used.
  @test MoonLander(:lgrExplicit;AlltpolyPts=true,numInterpPts=25,predictX0=true,maxSim=67)

  # Test #4
  @test MoonLander(:lgrExplicit;linearInterpolation=true,predictX0=true,maxSim=72)

  # Test #5
  @test MoonLander(:lgrExplicit;interpolationOn=false,predictX0=true)

  # Test #6
  @test MoonLander(:bkwEuler;interpolationOn=false,predictX0=true)

  # Test #7
  @test MoonLander(:trapezoidal;interpolationOn=false,predictX0=true)

  # Test #8
  # :tm methods with an infeasibe solution and a restoration error
  @test MoonLander(:bkwEuler;predictX0=false)
end

ext: jl


text:
"""
# Description

Determine if the input expression is solely the identifier of a column in a
table. Downstream systems can use this to avoid evaluating an anonymous
function by making a copy of the column instead.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `_::Bool`: `true` if the input object is a symbol; `false` otherwise.

# Examples

```
julia> is_column(:x)
true

julia> is_column(1)
false

julia> is_column(:(a + b))
false
```
"""
is_column(@nospecialize(e::Any)) = isa(e, Symbol)

"""
# Description

Determine if the input expression is solely a fixed constant. Downstream
systems can use this to avoid evaluating an anonymous function by making a new
column that contains many copies of that value instead.

# Arguments

1. `e::Any`: An expression.

# Return Values

1. `_::Bool`: `true` if the input object is a constant; `false` otherwise.

# Examples

```
julia> is_constant(:x)
false

julia> is_constant(1)
true

julia> is_constant(:(a + b))
false

julia> is_constant(Expr(:\$, :z))
true
```
"""
function is_constant(@nospecialize(e::Any))
    if isa(e, Expr) && e.head == :$ && length(e.args) == 1 && isa(e.args[1], Symbol)
        true
    elseif isa(e, Expr)
        false
    elseif isa(e, Symbol)
        false
    else
        true
    end
end


text:
"""
    isstandard(kr::KRatio)::Boolean

Does this k-ratio have all the necessary basic properties required for use as a standard
(:TakeOffAngle, and :BeamEnergy for both `stdProps` and `unkProps` and :Composition for
`unkProps`.)
"""
function isstandard(std::KRatio)::Bool
    b1 = haskey(std.unkProps, :Composition)
    b1 || @warn "The :Composition property needs to be specified for the `unknown` in $std."
    b2 = haskey(std.unkProps, :BeamEnergy)
    b2 || @warn "The :BeamEnergy property needs to be specified for the `unknown` in $std."
    b3 = haskey(std.unkProps, :TakeOffAngle)
    b3 || @warn "The :TakeOffAngle property needs to be specified for the `unknown` in $std."
    b4 = haskey(std.stdProps, :BeamEnergy)
    b4 || @warn "The :BeamEnergy property needs to be specified for the `standard` in $std."
    b5 = haskey(std.stdProps, :TakeOffAngle)
    b5 || @warn "The :TakeOffAngle property needs to be specified for the `standard` in $std."
    return b1 && b2 && b3 && b4 && b5
end

"""
    matches(kr::Union{KRatio, KRatios}, std::Standard)::Bool

Is `std` a match for `kr`? (Same element, same standard, same lines, same :BeamEnergy & :TakeOffAngle )

stdProps[:Composition] items match if 1) are Materials with almost the same mass-fractions; 2) they are equivalent
AbstractString values.  The later is for
"""
function matches(unk::Union{KRatio, KRatios}, std::KRatio)::Bool
    up, sp = unk.stdProps, std.stdProps
    return element(unk)==element(std) && #  Same element
        isapprox(unk.standard, std.standard, atol=1.0e-5) && # Relative to same material
        issetequal(unk.xrays, std.xrays) && # Same lines
        isapprox(sp[:BeamEnergy], up[:BeamEnergy], rtol=0.001) && # Same E0
        isapprox(sp[:TakeOffAngle], up[:TakeOffAngle], atol=deg2rad(0.1)) # Same take-off angle
end

"""
    standardize(kr::KRatio, std::KRatio)::KRatio
    standardize(kr::KRatios, std::KRatio)::KRatios
    standardize(kratios::Union{AbstractVector{KRatio},AbstractVector{<:KRatios}}, stds::AbstractVector{KRatio})

If the `std::KRatio` is a suitable match for `kr` then `kr` is restandardized using `std`.  Otherwise, the original
`KRatio` or `KRatios` is returned.
"""
function standardize(kr::KRatio, std::KRatio)::KRatio
    div(n::AbstractFloat, d::AbstractFloat) = value(n)/value(d)
    div(n::UncertainValue, d::UncertainValue) = divide(n,d,0.0)
    @assert isstandard(std) "$std does not contain all the necessary properties (see details in the above warnings)."
    return if matches(kr, std)
        newkr = div(kr.kratio, std.kratio)
        KRatio(kr.xrays, kr.unkProps, std.unkProps, std.unkProps[:Composition], newkr)
    else
        kr
    end
end
function standardize(krs::KRatios, std::KRatio)::KRatios
    @assert isstandard(std) "$std does not contain all the necessary properties (see details in the above warnings)."
    div(n::AbstractFloat, d::AbstractFloat) = value(n)/value(d)
    div(n::UncertainValue, d::UncertainValue) = divide(n, d, 0.0)
    return if matches(krs, std)
        newkrs = map(krs.kratios) do kr
            div(kr, std.kratio)
        end
        KRatios(krs.xrays, krs.unkProps, std.unkProps, std.unkProps[:Composition], newkrs)
    else
        krs
    end
end
function standardize(kratios::Union{AbstractVector{KRatio},AbstractVector{<:KRatios}}, stds::AbstractVector{KRatio})
    return map(kratios) do kr
        i = findfirst(std->matches(kr,std), stds)
        isnothing(i) ? kr : standardize(kr, stds[i])
    end
end


text:
using Test
prec = 15
@testset "slowstorage" begin
    for i in 1:100000
        Ks = rand(0.0001:0.001:0.1)[1]
        Ratio_Rip = rand(0.001:0.01:0.8)[1]
        Storage = rand(0.1:0.01:20)[1]
        GWflow = rand(0:0.1:10)[1]
        AreaRiparian = rand(0.01:0.01:1)[1]
        #Preferential = rand(0:0.1:10)[1]
        # test that output 0, if input 0
        @test slowstorage(0., 0.,AreaRiparian, Ks, Ratio_Rip) == (0, 0, 0)
        # test that storage decreases if no input
        Riparian_Discharge, Slow_Discharge, Slowstorage = slowstorage(0., Storage, AreaRiparian, Ks, 0.)
        @test -eps(Float64)*10 <= Storage - Slowstorage - Slow_Discharge <= eps(Float64)*10
        Riparian_Discharge, Slow_Discharge, Slowstorage = slowstorage(0., Storage, AreaRiparian, 0., Ratio_Rip)
        @test round(Storage - Slowstorage, digits=prec) == round(Riparian_Discharge, digits=prec)
        # test that discharge should be the Ks-ratio of sum of storage and overlandflow
        Riparian_Discharge, Slow_Discharge, Slowstorage = slowstorage(GWflow, Storage, AreaRiparian, Ks, Ratio_Rip)
        @test round((Storage + GWflow) * Ks, digits=prec)  == round(Slow_Discharge + Riparian_Discharge / AreaRiparian, digits = prec)
        @test round(Storage + GWflow - Slow_Discharge - Riparian_Discharge, digits = prec) == round(Slowstorage, digits = prec)
        # test right behavior if storage = 0
        Riparian_Discharge, Slow_Discharge, Slowstorage = slowstorage(GWflow, 0, AreaRiparian, Ks, Ratio_Rip)
        @test GWflow * Ks * (1 - Ratio_Rip) == Slow_Discharge
        @test round(GWflow * Ks, digits=prec) == round(Slow_Discharge + Riparian_Discharge / AreaRiparian, digits=prec)
        # if storage empty before, input = output+ new storage
        @test -eps(Float64)*10 <= Slow_Discharge + Riparian_Discharge + Slowstorage - GWflow <= eps(Float64)*10
    end
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
function Base.isapprox(x::DoubleFloat{T}, y::T; atol::Real=0, rtol::Real=atol>0 ? 0 : eps(T)^(5/8), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    return isapprox(x, DoubleFloat{T}(y), atol=atol, rtol=rtol, nans=nans, norm=norm)
end
function Base.isapprox(x::T, y::DoubleFloat{T}; atol::Real=0, rtol::Real=atol>0 ? 0 : eps(T)^(5/8), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    return isapprox(DoubleFloat{T}(x), y, atol=atol, rtol=rtol, nans=nans, norm=norm)
end
function Base.isapprox(x::DoubleFloat{T}, y::DoubleFloat{T}; atol::Real=0, rtol::Real=atol>0 ? 0 : eps(DoubleFloat{T})^(5/8), nans::Bool=false, norm::Function=norm) where {T<:IEEEFloat}
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= max(atol, rtol*max(abs(x), abs(y)))) || (nans && isnan(x) && isnan(y))
end

function Base.lerpi(j::Integer, d::Integer, a::DoubleFloat{T}, b::DoubleFloat{T}) where {T}
    t = DoubleFloat{T}(j)/d
    a = fma(-t, a, a)
    return fma(t, b, a)
end

# for compatibility with old or unrevised outside linalg functions
function Base.:(+)(v::Vector{DoubleFloat{T}}, x::T) where {T}
    return v .+ x
end
function Base.:(-)(v::Vector{DoubleFloat{T}}, x::T) where {T}
    return v .- x
end
function Base.:(+)(m::Matrix{DoubleFloat{T}}, x::T) where {T}
    return m .+ x
end
function Base.:(-)(m::Matrix{DoubleFloat{T}}, x::T) where {T}
    return m .- x
end


text:
module Alchemy

using Makie
using GLMakie
using Molecules
using Colors

include("Constants.jl")
include("Molecules.jl")
include("AtomicOrbitals.jl")
include("MolecularOrbitals.jl")
include("Interfaces/Fermi.jl")

end # module

max_stars_repo_name: gustavojra/Alchemy

max_issues_repo_name: gustavojra/Alchemy

max_forks_repo_name: gustavojra/Alchemy


text:
using ModelingToolkit
using OrdinaryDiffEq

"""
solves STG calcium model (Liu 1998)

The log(3000.0/Ca) in the algorithm doesn't play nicely with Zygote. Don't know why. I've maxed it with zero to make it work.

All neurons have following channels (with some conductances at 0:)
NaV
CaT
CaS
H
Kd
Leak
"""

function  CalciumNeuron(input)

    @parameters t
    D = Differential(t)

    @parameters eNa eh eK eleak tauCa gNabar gCaSbar gCaTbar gKabar gKCabar gKdrbar ghbar gleak

    paramvars = [eNa, eh, eK, eleak, tauCa, gNabar, gCaSbar, gCaTbar, gKabar, gKCabar, gKdrbar, ghbar, gleak]

    @variables V(t) Ca(t) mNa(t) hNa(t) mCaS(t) hCaS(t) mCaT(t) hCaT(t) mIh(t) mKa(t) hKa(t) mKCa(t) mKdr(t)
    statevars = [V, Ca, mNa, hNa, mCaS, hCaS, mCaT, hCaT, mIh, mKa, hKa, mKCa, mKdr]

    # calcium reversal potential
    # eCa = (500.0)*(8.6174e-5)*(283.15)*(log(max((3000.0/Ca),1e-5)))
    eCa = (500.0)*(8.6174e-5)*(283.15)*(log((3000.0/Ca)))

    # currents
    INa=  gNabar*mNa^3*hNa*(eNa - V)
    ICaS= gCaSbar*mCaS^3*hCaS*(eCa - V)
    ICaT= gCaTbar*mCaT^3*hCaT*(eCa - V)
    Ih= ghbar*mIh*(eh - V)
    IKa= gKabar*mKa^3*hKa*(eK - V)
    IKCa=  gKCabar*mKCa^4*(eK - V)
    IKdr=  gKdrbar*mKdr^4*(eK - V)
    Ileak= gleak*(eleak - V)

    # calcium
    Ca_inf = 0.05 + 0.94*(ICaS + ICaT);

    # gating
    Na_m_inf=  1.0./(1.0+exp((V+25.5)./-5.29));
    Na_h_inf=  1.0./(1.0+exp((V+48.9)./5.18));
    Na_tau_m=  1.32 - 1.26./(1+exp((V+120.0)./-25.0));
    Na_tau_h=  (0.67./(1.0+exp((V+62.9)./-10.0))).*(1.5+1.0./(1.0+exp((V+34.9)./3.6)));

    CaS_m_inf=  1.0./(1.0+exp((V+33.0)./-8.1));
    CaS_h_inf=  1.0./(1.0+exp((V+60.0)./6.2));
    CaS_tau_m=  1.4 + 7.0./(exp((V+27.0)./10.0) + exp((V+70.0)./-13.0));
    CaS_tau_h=  60.0 + 150.0./(exp((V+55.0)./9.0) + exp((V+65.0)./-16.0));

    CaT_m_inf=  1.0./(1.0 + exp((V+27.1)./-7.2));
    CaT_h_inf=  1.0./(1.0 + exp((V+32.1)./5.5));
    CaT_tau_m=  21.7 - 21.3./(1.0 + exp((V+68.1)./-20.5));
    CaT_tau_h=  105.0 - 89.8./(1.0 + exp((V+55.0)./-16.9));

    Ih_m_inf=  1.0./(1.0+exp((V+70.0)./6.0));
    Ih_tau_m=  (272.0 + 1499.0./(1.0+exp((V+42.2)./-8.73)));

    Ka_m_inf=  1.0./(1.0+exp((V+27.2)./-8.7));
    Ka_h_inf=  1.0./(1.0+exp((V+56.9)./4.9));
    Ka_tau_m=  11.6 - 10.4./(1.0+exp((V+32.9)./-15.2));
    Ka_tau_h=  38.6 - 29.2./(1.0+exp((V+38.9)./-26.5));

    KCa_m_inf=  (Ca./(Ca+3.0))./(1.0+exp((V+28.3)./-12.6));
    KCa_tau_m=  90.3 - 75.1./(1.0+exp((V+46.0)./-22.7));

    Kdr_m_inf=  1.0./(1.0+exp((V+12.3)./-11.8));
    Kdr_tau_m=  7.2 - 6.4./(1.0+exp((V+28.3)./-19.2));

    eqs = [
    D(V) ~              INa + ICaS + ICaT + Ih + IKa + IKCa + IKdr + Ileak,
    D(Ca) ~                   (1/tauCa)*(Ca_inf - Ca),
    D(mNa) ~                   (1/Na_tau_m)*(Na_m_inf - mNa),
    D(hNa) ~                   (1/Na_tau_h)*(Na_h_inf - hNa),
    D(mCaS) ~                   (1/CaS_tau_m)*(CaS_m_inf - mCaS),
    D(hCaS) ~                   (1/CaS_tau_h)*(CaS_h_inf - hCaS),
    D(mCaT) ~                   (1/CaT_tau_m)*(CaT_m_inf - mCaT),
    D(hCaT) ~                   (1/CaT_tau_h)*(CaT_h_inf - hCaT),
    D(mIh) ~                   (1/Ih_tau_m)*(Ih_m_inf - mIh),
    D(mKa) ~                  (1/Ka_tau_m)*(Ka_m_inf - mKa),
    D(hKa) ~                   (1/Ka_tau_h)*(Ka_h_inf - hKa),
    D(mKCa) ~                  (1/KCa_tau_m)*(KCa_m_inf - mKCa),
    D(mKdr) ~                   (1/Kdr_tau_m)*(Kdr_m_inf - mKdr)
    ]

    ps = paramvars .=> [50,-20,-80,-50,20,100,3,1.3,5,10,20,0.5,0.01]
    ics = statevars .=> [-60,0.05,0,0,0,0,0,0,0,0,0,0,0]

    od = ODESystem(eqs, t, statevars, paramvars)
    tspan = (0.,500.)
    return od, ics, tspan, ps
end

# function predict(p)
#     return Array(solve(prob, Tsit5(), p=p, saveat=tsteps, sensitivity=BacksolveAdjoint(;autojacvec=true)))
# end

# function loss(p)
#     prediction = predict(p)
#     return sum(abs2, prediction - data)
# end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

SubWorker{T,A,PT} = RemoteChannel{Channel{Vector{SubProblem{T,A,PT}}}}
ScenarioProblemChannel{S} = RemoteChannel{Channel{ScenarioProblems{S}}}
Work = RemoteChannel{Channel{Int}}
Progress{T <: AbstractFloat} = Tuple{Int,Int,SubproblemSolution{T}}
ProgressQueue{T <: AbstractFloat} = RemoteChannel{Channel{Progress{T}}}

function initialize_subproblems!(ph::AbstractProgressiveHedging,
                                 subworkers::Vector{SubWorker{T,A,PT}},
                                 scenarioproblems::DistributedScenarioProblems,
                                 penaltyterm::AbstractPenaltyTerm) where {T <: AbstractFloat,
                                                                          A <: AbstractVector,
                                                                          PT <: AbstractPenaltyTerm}
    # Create subproblems on worker processes
    @sync begin
        for w in workers()
            subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{T,A,PT}}}(1), w)
            prev = map(2:(w-1)) do p
                scenarioproblems.scenario_distribution[p-1]
            end
            start_id = isempty(prev) ? 0 : sum(prev)
            @async remotecall_fetch(initialize_subworker!,
                                    w,
                                    subworkers[w-1],
                                    scenarioproblems[w-1],
                                    penaltyterm,
                                    start_id)
        end
    end
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedging,
                          subworkers::Vector{<:SubWorker{T}}) where T <: AbstractFloat
    # Update δ₂
    partial_δs = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_δs[i] = remotecall_fetch(
                w,
                subworkers[w-1],
                ph.ξ) do sw, ξ
                    subproblems = fetch(sw)
                    return mapreduce(+, subproblems, init = zero(T)) do subproblem
                        π = subproblem.probability
                        x = subproblem.x
                        π * norm(x - ph.ξ, 2) ^ 2
                    end
                end
        end
    end
    ph.data.δ₂ = sum(partial_δs)
    return nothing
end

function initialize_subworker!(subworker::SubWorker{T,A,PT},
                               scenarioproblems::ScenarioProblemChannel,
                               penaltyterm::AbstractPenaltyTerm,
                               start_id::Integer) where {T <: AbstractFloat,
                                                         A <: AbstractArray,
                                                         PT <: AbstractPenaltyTerm}
    sp = fetch(scenarioproblems)
    subproblems = Vector{SubProblem{T,A,PT}}(undef, num_subproblems(sp))
    for i = 1:num_subproblems(sp)
        subproblems[i] = SubProblem(
            subproblem(sp, i),
            start_id + i,
            T(probability(sp, i)),
            copy(penaltyterm))
    end
    put!(subworker, subproblems)
    return nothing
end

function restore_subproblems!(subworkers::Vector{<:SubWorker})
    @sync begin
        for w in workers()
            @async remotecall_fetch(w, subworkers[w-1]) do sw
                for subproblem in fetch(sw)
                    restore_subproblem!(subproblem)
                end
            end
        end
    end
    return nothing
end

function resolve_subproblems!(subworker::SubWorker{T,A,PT},
                              ξ::AbstractVector,
                              r::AbstractFloat) where {T <: AbstractFloat,
                                                       A <: AbstractArray,
                                                       PT <: AbstractPenaltyTerm}
    subproblems::Vector{SubProblem{T,A,PT}} = fetch(subworker)
    Qs = Vector{SubproblemSolution{T}}(undef, length(subproblems))
    # Reformulate and solve sub problems
    for (i,subproblem) in enumerate(subproblems)
        reformulate_subproblem!(subproblem, ξ, r)
        Qs[i] = subproblem(ξ)
    end
    # Return current objective value
    return sum(Qs)
end

function collect_primals(subworker::SubWorker{T,A,PT}, n::Integer) where {T <: AbstractFloat,
                                                                          A <: AbstractArray,
                                                                          PT <: AbstractPenaltyTerm}
    subproblems::Vector{SubProblem{T,A,PT}} = fetch(subworker)
    return mapreduce(+, subproblems, init = zeros(T, n)) do subproblem
        π = subproblem.probability
        x = subproblem.x
        π * x
    end
end

function calculate_objective_value(subworkers::Vector{<:SubWorker{T}}) where T <: AbstractFloat
    partial_objectives = Vector{Float64}(undef, nworkers())
    @sync begin
        for (i, w) in enumerate(workers())
            @async partial_objectives[i] = remotecall_fetch(w, subworkers[w-1]) do sw
                return mapreduce(+, fetch(sw), init = zero(T)) do subproblem
                    _objective_value(subproblem)
                end
            end
        end
    end
    return sum(partial_objectives)
end

function work_on_subproblems!(subworker::SubWorker{T,A,PT},
                              work::Work,
                              finalize::Work,
                              progress::ProgressQueue{T},
                              x̄::RemoteRunningAverage{A},
                              δ::RemoteRunningAverage{T},
                              iterates::RemoteIterates{A},
                              r::IteratedValue{T}) where {T <: AbstractFloat,
                                                          A <: AbstractArray,
                                                          PT <: AbstractPenaltyTerm}
    subproblems::Vector{SubProblem{T,A,PT}} = fetch(subworker)
    if isempty(subproblems)
       # Workers has nothing do to, return.
       return nothing
    end
    x̄ = fetch(x̄)
    δ = fetch(δ)
    quit = false
    while true
        t::Int = try
            if isready(finalize)
                quit = true
                take!(finalize)
            else
                wait(work)
                take!(work)
            end
        catch err
            if err isa InvalidStateException
                # Master closed the work/finalize channel. Worker finished
                return nothing
            end
        end
        t == -1 && continue
        ξ::A = fetch(iterates, t)
        if t > 1
            update_subproblems!(subproblems, ξ, fetch(r,t-1))
        end
        for (i,subproblem) in enumerate(subproblems)
            !quit && subtract!(δ, i)
            !quit && subtract!(x̄, i)
            x = subproblem.x
            π = subproblem.probability
            !quit && add!(δ, i, norm(x - ξ, 2) ^ 2, π)
            reformulate_subproblem!(subproblem, ξ, fetch(r, t))
            Q::SubproblemSolution{T} = subproblem(ξ)
            !quit && add!(x̄, i, π)
            !quit && put!(progress, (t, subproblem.id, Q))
        end
        if quit
            # Worker finished
            return nothing
        end
    end
end


text:
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using ShapeReconstructionPaLS.Utils
using ParamLevelSet
using ShapeReconstructionPaLS.ShapeFromSilhouette;
using Statistics
using Distributed
using LinearAlgebra
using Test

plotting = false;
if plotting
	using PyPlot
	using jInvVisPyPlot
end

println("Testing getData for all the methods:")
n = [64,64,64];
Mesh = getRegularMesh([0.0;3.0;0.0;3.0;0.0;3.0],n);
ScreenMeshX2X3 = getRegularMesh([0.0 3.0 0.0 3.0],n[2:3]);
XScreenLoc = 5.0;
CamLoc = [-10.0 ; 0.0 ;0.0];

alpha = [1.5;2.5;2.0;-1.0];
beta = [2.0;4.0;2.6;2.5];
Xs = [1.5 1.5 1.5 ; 2.0 2.0 2.0; 1.2 2.3 1.5; 2.2 1.5 2.0];
# alpha = [1.0;];
# beta = [1.0;];
# Xs = [1.5 1.5 1.5 ;];

m = wrapTheta(alpha,beta,Xs);
sigmaH = getDefaultHeaviside();
u, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0,sigma = sigmaH);
theta_phi = deg2rad.([37.0 24.0 ; 27.0 14.0;11.0 11.0]);
b = [-8.0 5.2 5.5; 6.0 3.0 5.0 ;4.5 -2.3 -5.6]*mean(Mesh.h);

methods = [MATFree,RBFBasedSimple2,RBFBased];
D = Array{Array{Float32,2}}(undef,length(methods));

for k=1
	println("test GetData MATFree");
	pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],nworkers());
	Dremote, = getData(u[:],pFor);
	D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);
end

for k=2
	println("test GetData RBFBasedSimple2")
	pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],nworkers());
	Dremote, = getData(m[:],pFor);
	D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);
end

k=3
println("test GetData RBFBased")
m_wraped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b)
pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],nworkers());
Dremote, = getData(m_wraped[:],pFor);
D[k] = arrangeRemoteCallDataIntoLocalData(Dremote);
println("Done!")

if plotting
	close("all");
	jj = 3;
	figure();
	plotModel(reshape(u,n[1],n[2],n[3]));
	figure();
	imshow(reshape(D[1][:,jj],tuple(ScreenMeshX2X3.n...))); colorbar();
	figure();
	imshow(reshape(D[2][:,jj],tuple(ScreenMeshX2X3.n...))); colorbar();
	figure();
	imshow(reshape(D[3][:,jj],tuple(ScreenMeshX2X3.n...))); colorbar();
	figure();
	imshow(reshape(D[1][:,jj] .- D[3][:,jj],tuple(ScreenMeshX2X3.n...))); colorbar();
	println(norm(D[1][:,jj] .- D[3][:,jj]))
	println(norm(D[1][:,jj]))
end

for k=2:length(D)
	r = norm(D[k] - D[1])./norm(D[1]);
	# figure(); imshow(reshape(D[k][:,1],tuple(ScreenMeshX2X3.n...)))
	println("Comparing 1 and ",k,". Relative diff: ",r);
	@test r <= 0.25
end

println("Testing Sensitivities: MATFree");
du = 0.01*randn(size(u));
for k=1
	pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
	Dremote,pFor = getData(u[:],pFor);
	pFor = fetch(pFor[1]);
	Dk = arrangeRemoteCallDataIntoLocalData(Dremote);
	for ii = 1:8
		dm = (0.5^ii)*du;
		pFor_t = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
		Dremote, = getData(u[:] + dm[:],pFor_t);
		Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
		println("norm(dt-d0): ",norm(Dk_t[:]-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] .- Dk[:] .- getSensMatVec(dm[:],u[:],pFor)[:]));
	end
end
println("Testing Sensitivities: RBFBasedSimple2");
dmm = 0.01*randn(size(m));
for k = 2
	pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
	Dremote, = getData(m[:],pFor);
	pFor = fetch(pFor[1]);
	Dk = arrangeRemoteCallDataIntoLocalData(Dremote);

	for ii = 1:8
		dm = (0.5^ii)*dmm;
		pFor_t = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
		Dremote, = getData(m[:] + dm[:],pFor_t);
		Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
		println("norm(dt-d0): ",norm(Dk_t[:].-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] .- Dk[:] .- getSensMatVec(dm[:],m[:],pFor)[:]));
	end
end
println("Testing Sensitivities: RBFBased");
k=3
m_wraped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b)
dm_wraped = wrapRBFparamAndRotationsTranslations(dmm,0.05*deg2rad.(randn(size(theta_phi))),0.05*mean(Mesh.h)*randn(size(b)));
pFor = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
Dremote, = getData(m_wraped[:],pFor);
Dk = arrangeRemoteCallDataIntoLocalData(Dremote);
V = 0;

for ii = 1:8
	dm = (0.5^ii)*dm_wraped;
	pFor_t = getSfSParam(Mesh,theta_phi,b,ScreenMeshX2X3,XScreenLoc,CamLoc,methods[k],1);
	Dremote, = getData(m_wraped[:] + dm[:],pFor_t);
	Dk_t = arrangeRemoteCallDataIntoLocalData(Dremote);
	V = getSensMatVec(dm[:],m_wraped[:],fetch(pFor[1]))[:];
	println("norm(dt-d0): ",norm(Dk_t[:]-Dk[:]),", norm(dt - d0 - J0*dm): ",norm(Dk_t[:] - Dk[:] - V));
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
import ForwardDiff: Dual, partials, value

######################
# Extraction methods #
######################

# Extractions are supposed to unpack the value and the partials
# The partials should be put into a tensor of higher order.
# The extraction methods need to know the input type to the function
# that generated the result. The reason for this is that there is no
# difference in the output type (except the number of partials) for
# norm(v) and det(T) where v is a vector and T is a second order tensor.
# For dim = 1 the type is exactly the same.

# Scalar
@inline function _extract_value{N, T}(v::Dual{N, T}, ::Vec{N})
    return value(v)
end
@inline function _extract_gradient{N, T}(v::Dual{N, T}, ::Vec{N})
    return Vec{N, T}(partials(v).values)
end

@inline function _extract_value{N, T, dim, T2}(v::Dual{N, T}, ::SymmetricTensor{2, dim, T2, N})
    return value(v)
end
@inline function _extract_gradient{N, T, dim, T2}(v::Dual{N, T}, ::SymmetricTensor{2, dim, T2, N})
    return SymmetricTensor{2, dim, T}(partials(v).values)
end

@inline function _extract_value{N, T, dim, T2}(v::Dual{N, T}, ::Tensor{2, dim, T2, N})
    return value(v)
end
@inline function _extract_gradient{N, T, dim, T2}(v::Dual{N, T}, ::Tensor{2, dim, T2, N})
    return Tensor{2, dim, T}(partials(v).values)
end

# Vec
@inline function _extract_value{D <: Dual}(v::Vec{1, D}, ::Any)
    @inbounds begin
        v1 = value(v[1])
        f = Vec{1}((v1,))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Vec{1, D}, ::Any)
    @inbounds begin
        p1 = partials(v[1])
        ∇f = Tensor{2, 1}((p1[1],))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::Vec{2, D}, ::Any)
    @inbounds begin
        v1, v2 = value(v[1]), value(v[2])
        f = Vec{2}((v1, v2))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Vec{2, D}, ::Any)
    @inbounds begin
        p1, p2 = partials(v[1]), partials(v[2])
        ∇f = Tensor{2, 2}((p1[1], p2[1], p1[2], p2[2]))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::Vec{3, D}, ::Any)
    @inbounds begin
        v1, v2, v3 = value(v[1]), value(v[2]), value(v[3])
        f = Vec{3}((v1, v2, v3))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Vec{3, D}, ::Any)
    @inbounds begin
        p1, p2, p3 = partials(v[1]), partials(v[2]), partials(v[3])
        ∇f = Tensor{2, 3}((p1[1], p2[1], p3[1], p1[2], p2[2], p3[2], p1[3], p2[3], p3[3]))
    end
    return ∇f
end

# Second order tensor
@inline function _extract_value{D <: Dual}(v::Tensor{2, 1, D}, ::Any)
    @inbounds begin
        v1 = value(v[1,1])
        f = Tensor{2, 1}((v1,))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 1, D}, ::Any)
    @inbounds begin
        p1 = partials(v[1,1])
        ∇f = Tensor{4, 1}((p1[1],))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 1, D}, ::Any)
    @inbounds begin
        v1 = value(v[1,1])
        f = SymmetricTensor{2, 1}((v1,))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 1, D}, ::Any)
    @inbounds begin
        p1 = partials(v[1,1])
        ∇f = SymmetricTensor{4, 1}((p1[1],))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::Tensor{2, 2, D}, ::Any)
    @inbounds begin
        v1, v2, v3, v4 = value(v[1,1]), value(v[2,1]), value(v[1,2]), value(v[2,2])
        f = Tensor{2, 2}((v1, v2, v3, v4))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 2, D}, ::Any)
    @inbounds begin
        p1, p2, p3, p4 = partials(v[1,1]), partials(v[2,1]), partials(v[1,2]), partials(v[2,2])
        ∇f = Tensor{4, 2}((p1[1], p2[1], p3[1], p4[1],
                           p1[2], p2[2], p3[2], p4[2],
                           p1[3], p2[3], p3[3], p4[3],
                           p1[4], p2[4], p3[4], p4[4]))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 2, D}, ::Any)
    @inbounds begin
        v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[2,2])
        f = SymmetricTensor{2, 2}((v1, v2, v3))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 2, D}, ::Any)
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[2,2])
        ∇f = SymmetricTensor{4, 2}((p1[1], p2[1], p3[1],
                                    p1[2], p2[2], p3[2],
                                    p1[3], p2[3], p3[3]))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::Tensor{2, 3, D}, ::Any)
    @inbounds begin
        v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[3,1])
        v4, v5, v6 = value(v[1,2]), value(v[2,2]), value(v[3,2])
        v7, v8, v9 = value(v[1,3]), value(v[2,3]), value(v[3,3])
        f = Tensor{2,3}((v1, v2, v3, v4, v5, v6, v7, v8, v9))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::Tensor{2, 3, D}, ::Any)
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1])
        p4, p5, p6 = partials(v[1,2]), partials(v[2,2]), partials(v[3,2])
        p7, p8, p9 = partials(v[1,3]), partials(v[2,3]), partials(v[3,3])
        ∇f = Tensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1], p7[1], p8[1], p9[1],
                           p1[2], p2[2], p3[2], p4[2], p5[2], p6[2], p7[2], p8[2], p9[2],  #    ###  #
                           p1[3], p2[3], p3[3], p4[3], p5[3], p6[3], p7[3], p8[3], p9[3],  #    # #  #
                           p1[4], p2[4], p3[4], p4[4], p5[4], p6[4], p7[4], p8[4], p9[4],  ###  ###  ###
                           p1[5], p2[5], p3[5], p4[5], p5[5], p6[5], p7[5], p8[5], p9[5],
                           p1[6], p2[6], p3[6], p4[6], p5[6], p6[6], p7[6], p8[6], p9[6],
                           p1[7], p2[7], p3[7], p4[7], p5[7], p6[7], p7[7], p8[7], p9[7],
                           p1[8], p2[8], p3[8], p4[8], p5[8], p6[8], p7[8], p8[8], p9[8],
                           p1[9], p2[9], p3[9], p4[9], p5[9], p6[9], p7[9], p8[9], p9[9]))
    end
    return ∇f
end

@inline function _extract_value{D <: Dual}(v::SymmetricTensor{2, 3, D}, ::Any)
    @inbounds begin
        v1, v2, v3 = value(v[1,1]), value(v[2,1]), value(v[3,1])
        v4, v5, v6 = value(v[2,2]), value(v[3,2]), value(v[3,3])
        f = SymmetricTensor{2, 3}((v1, v2, v3, v4, v5, v6))
    end
    return f
end
@inline function _extract_gradient{D <: Dual}(v::SymmetricTensor{2, 3, D}, ::Any)
    @inbounds begin
        p1, p2, p3 = partials(v[1,1]), partials(v[2,1]), partials(v[3,1])
        p4, p5, p6 = partials(v[2,2]), partials(v[3,2]), partials(v[3,3])
        ∇f = SymmetricTensor{4, 3}((p1[1], p2[1], p3[1], p4[1], p5[1], p6[1],
                                    p1[2], p2[2], p3[2], p4[2], p5[2], p6[2],
                                    p1[3], p2[3], p3[3], p4[3], p5[3], p6[3],
                                    p1[4], p2[4], p3[4], p4[4], p5[4], p6[4],
                                    p1[5], p2[5], p3[5], p4[5], p5[5], p6[5],
                                    p1[6], p2[6], p3[6], p4[6], p5[6], p6[6]))
    end
    return ∇f
end

##################
# Load functions #
##################

# Loaders are supposed to take a tensor of real values and convert it
# into a tensor of dual values where the seeds are correctly defined.

@inline function _load{T}(v::Vec{1, T})
    @inbounds v_dual = Vec{1}((Dual(v[1], one(T)),))
    return v_dual
end

@inline function _load{T}(v::Vec{2, T})
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Vec{2}((Dual(v[1], o, z),
                               Dual(v[2], z, o)))
    return v_dual
end

@inline function _load{T}(v::Vec{3, T})
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Vec{3}((Dual(v[1], o, z, z),
                               Dual(v[2], z, o, z),
                               Dual(v[3], z, z, o)))
    return v_dual
end

# Second order tensors
@inline function _load{T}(v::Tensor{2, 1, T})
    @inbounds v_dual = Tensor{2, 1}((Dual(v.data.data[1], one(T)),))
    return v_dual
end

@inline function _load{T}(v::SymmetricTensor{2, 1, T})
    @inbounds v_dual = SymmetricTensor{2, 1}((Dual(v.data.data[1], one(T)),))
    return v_dual
end

@inline function _load{T}(v::Tensor{2, 2, T})
    data = v.data.data
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Tensor{2, 2}((Dual(data[1], o, z, z, z),
                                     Dual(data[2], z, o, z, z),
                                     Dual(data[3], z, z, o, z),
                                     Dual(data[4], z, z, z, o)))
    return v_dual
end

@inline function _load{T}(v::SymmetricTensor{2, 2, T})
    data = v.data.data
    o = one(T)
    o2 = T(1/2)
    z = zero(T)
    @inbounds v_dual = SymmetricTensor{2, 2}((Dual(data[1], o, z, z),
                                              Dual(data[2], z, o2, z),
                                              Dual(data[3], z, z, o)))
    return v_dual
end

@inline function _load{T}(v::Tensor{2, 3, T})
    data = v.data.data
    o = one(T)
    z = zero(T)
    @inbounds v_dual = Tensor{2, 3}((Dual(data[1], o, z, z, z, z, z, z, z, z),
                                     Dual(data[2], z, o, z, z, z, z, z, z, z),
                                     Dual(data[3], z, z, o, z, z, z, z, z, z),
                                     Dual(data[4], z, z, z, o, z, z, z, z, z),
                                     Dual(data[5], z, z, z, z, o, z, z, z, z),
                                     Dual(data[6], z, z, z, z, z, o, z, z, z),
                                     Dual(data[7], z, z, z, z, z, z, o, z, z),
                                     Dual(data[8], z, z, z, z, z, z, z, o, z),
                                     Dual(data[9], z, z, z, z, z, z, z, z, o)))
    return v_dual
end

@inline function _load{T}(v::SymmetricTensor{2, 3, T})
    data = v.data.data
    o = one(T)
    o2 = T(1/2)
    z = zero(T)
    @inbounds v_dual = SymmetricTensor{2, 3}((Dual(data[1], o, z, z, z, z, z),
                                              Dual(data[2], z, o2, z, z, z, z),
                                              Dual(data[3], z, z, o2, z, z, z),
                                              Dual(data[4], z, z, z, o, z, z),
                                              Dual(data[5], z, z, z, z, o2, z),
                                              Dual(data[6], z, z, z, z, z, o)))
    return v_dual
end

"""
```julia
gradient(f::Function, v::Union{SecondOrderTensor, Vec})
gradient(f::Function, v::Union{SecondOrderTensor, Vec}, :all)
```
Computes the gradient of the input function. If the (pseudo)-keyword `all`
is given, the value of the function is also returned as a second output argument.

**Example:**

```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> ∇f = gradient(norm, A)
2×2 ContMechTensors.SymmetricTensor{2,2,Float64,3}:
 0.434906  0.56442
 0.56442   0.416793

julia> ∇f, f = gradient(norm, A, :all);
```
"""
function Base.gradient{F}(f::F, v::Union{SecondOrderTensor, Vec})
    v_dual = _load(v)
    res = f(v_dual)
    return _extract_gradient(res, v)
end
function Base.gradient{F}(f::F, v::Union{SecondOrderTensor, Vec}, ::Symbol)
    v_dual = _load(v)
    res = f(v_dual)
    return _extract_gradient(res, v), _extract_value(res, v)
end

"""
```julia
hessian(f::Function, v::Union{SecondOrderTensor, Vec})
hessian(f::Function, v::Union{SecondOrderTensor, Vec}, :all)
```
Computes the hessian of the input function. If the (pseudo)-keyword `all`
is given, the lower order results (gradient and value) of the function is
also returned as a second and third output argument.

**Example:**

```jldoctest
julia> A = rand(SymmetricTensor{2, 2});

julia> ∇∇f = hessian(norm, A)
2×2×2×2 ContMechTensors.SymmetricTensor{4,2,Float64,9}:
[:, :, 1, 1] =
  0.596851  -0.180684
 -0.180684  -0.133425

[:, :, 2, 1] =
 -0.180684   0.133546
  0.133546  -0.173159

[:, :, 1, 2] =
 -0.180684   0.133546
  0.133546  -0.173159

[:, :, 2, 2] =
 -0.133425  -0.173159
 -0.173159   0.608207

julia> ∇∇f, ∇f, f = hessian(norm, A, :all);
```
"""
function hessian{F}(f::F, v::Union{SecondOrderTensor, Vec})
    gradf = y -> gradient(f, y)
    return gradient(gradf, v)
end

function hessian{F}(f::F, v::Union{SecondOrderTensor, Vec}, ::Symbol)
    gradf = y -> gradient(f, y)
    return gradient(gradf, v), gradient(f, v, :all)...
end


text:
# 32-bit fixed point; parameter `f` is the number of fraction bits
immutable Fixed{T <: Signed,f} <: FixedPoint{T,  f}
    i::T

    # constructor for manipulating the representation;
    # selected by passing an extra dummy argument
    Fixed(i::Integer,_) = new(i % T)

    Fixed(x) = convert(Fixed{T,f}, x)
end

typealias Fixed16 Fixed{Int32, 16}

  rawtype{T,f}(::Type{Fixed{T,f}}) = T
nbitsfrac{T,f}(::Type{Fixed{T,f}}) = f

# basic operators
-{T,f}(x::Fixed{T,f}) = Fixed{T,f}(-x.i,0)
abs{T,f}(x::Fixed{T,f}) = Fixed{T,f}(abs(x.i),0)

+{T,f}(x::Fixed{T,f}, y::Fixed{T,f}) = Fixed{T,f}(x.i+y.i,0)
-{T,f}(x::Fixed{T,f}, y::Fixed{T,f}) = Fixed{T,f}(x.i-y.i,0)

# with truncation:
#*{f}(x::Fixed32{f}, y::Fixed32{f}) = Fixed32{f}(Base.widemul(x.i,y.i)>>f,0)
# with rounding up:
*{T,f}(x::Fixed{T,f}, y::Fixed{T,f}) = Fixed{T,f}((Base.widemul(x.i,y.i) + (convert(widen(T), 1) << (f-1) ))>>f,0)

/{T,f}(x::Fixed{T,f}, y::Fixed{T,f}) = Fixed{T,f}(div(convert(widen(T), x.i) << f, y.i), 0)

# # conversions and promotions
convert{T,f}(::Type{Fixed{T,f}}, x::Integer) = Fixed{T,f}(round(T, convert(widen1(T),x)<<f),0)
convert{T,f}(::Type{Fixed{T,f}}, x::AbstractFloat) = Fixed{T,f}(round(T, trunc(widen1(T),x)<<f + rem(x,1)*(1<<f)),0)
convert{T,f}(::Type{Fixed{T,f}}, x::Rational) = Fixed{T,f}(x.num)/Fixed{T,f}(x.den)

rem{T,f}(x::Integer, ::Type{Fixed{T,f}}) = Fixed{T,f}(rem(x,T)<<f,0)
rem{T,f}(x::Real,    ::Type{Fixed{T,f}}) = Fixed{T,f}(rem(Integer(trunc(x)),T)<<f + rem(Integer(round(rem(x,1)*(1<<f))),T),0)

convert{T,f}(::Type{BigFloat}, x::Fixed{T,f}) =
    convert(BigFloat,x.i>>f) + convert(BigFloat,x.i&(1<<f - 1))/convert(BigFloat,1<<f)
convert{TF<:AbstractFloat,T,f}(::Type{TF}, x::Fixed{T,f}) =
    convert(TF,x.i>>f) + convert(TF,x.i&(1<<f - 1))/convert(TF,1<<f)

convert{T,f}(::Type{Bool}, x::Fixed{T,f}) = x.i!=0
function convert{TI<:Integer, T,f}(::Type{TI}, x::Fixed{T,f})
    isinteger(x) || throw(InexactError())
    convert(TI, x.i>>f)
end

convert{TR<:Rational,T,f}(::Type{TR}, x::Fixed{T,f}) =
    convert(TR, x.i>>f + (x.i&(1<<f-1))//(1<<f))

promote_rule{T,f,TI<:Integer}(ft::Type{Fixed{T,f}}, ::Type{TI}) = Fixed{T,f}
promote_rule{T,f,TF<:AbstractFloat}(::Type{Fixed{T,f}}, ::Type{TF}) = TF
promote_rule{T,f,TR}(::Type{Fixed{T,f}}, ::Type{Rational{TR}}) = Rational{TR}

# TODO: Document and check that it still does the right thing.
decompose{T,f}(x::Fixed{T,f}) = x.i, -f, 1


text:
# In this test we classify integers accoring to their divisor
function divisortest()
    @testset "divisortest" begin
        ar=[1:100...]
        r=Repo(ar,4)
        for i in 1:100
            a=mod(r[i],2)
            b=mod(r[i],3)
            if a==0 && b!=0 setclass!(r,i,2); end    #dividble by 2 but not 3
            if a!=0 && b==0 setclass!(r,i,3); end    #dividble by 3 but not 2
            if a==0 && b==0 setclass!(r,i,4); end    #dividble by 6
        end
        divble=class(r,4)
        @test length(divble)==16    # 16 numbers<=100 muliples of 6
        N=0
        for x in divble
            N+=1
            @test mod(x,6)==0
        end
        @test N==16
    end
    return true
end


text:
import Base

include("grammar.jl")

HAIR_SPACE = '\u200a'

"""
A *derivation*.
"""
struct Derivation
    "The grammar to which the derivations refers to."
    G::Grammar
    "The derivation steps"
    steps::Array
    "The sentential form of the derivations"
    sf::Array
    "The string representation of the derivation"
    repr::AbstractString
end

### Constuctors ###

"""
    Derivation(G::Grammar)::Derivation
Builds a [`Derivation`](@ref) from a given [`Grammar`](@ref).
"""
Derivation(G::Grammar)::Derivation = Derivation(G, [], [G.S], G.S)

### Operators ###

"""
    next(d::Derivation, prod::Int, pos::Int)::Derivation
Applies the specified production to the given position in the sentential form.
"""
function next(d::Derivation, prod::Int, pos::Int)::Derivation
    sf = d.sf
    prod = ensure_production_index(d, prod)
    P = astype0(d.G.P[prod])
    if sf[pos:pos+length(P.left)-1] != P.left throw(ArgumentError("Cannot apply " * string(P) * " at position " * string(pos) * " of " * string(sf))) end
    sf = [c for c ∈ [sf[begin:pos-1]; P.right; sf[pos + length(P.left):end]] if c ≠ "ε"]
    steps = [d.steps; [(prod, pos)]]
    repr = string(d.repr," -> ", join(sf,HAIR_SPACE))
    clone = Derivation(d.G, steps, sf, repr)
    return clone
end

"""
    next(d::Derivation, prod::Production, pos::Int)::Derivation
Applies the specified [`Production`](@ref) to the given position in the sentential form.
"""
function next(d::Derivation, prod::Production, pos::Int)::Derivation
    p = ensure_production_index(d, prod)
    return next(d, p, pos)
end

"""
    next(d::Derivation, prod::Int, pos::Int)::Derivation
Applies the specified productions to the given position in the sentential form, one by one
"""
function next(d::Derivation, prod::AbstractArray{Tuple{Int, Int}})::Derivation
    res = d
    for (nprod, pos) ∈ prod
        res = next(res, nprod, pos)
    end
    return res
end

"""
    leftmost(d::Derivation, prod::Int)::Derivation
Performs a *leftmost* derivation step.
Applies the specified production to the current leftmost nonterminal in the sentential form.
"""
function leftmost(d::Derivation, prod::Int)::Derivation
    if ~d.G.iscontextfree
        throw(ArgumentError("Cannot perform a leftmost derivation on a non context-free grammar"))
    end
    if length(d.sf) == 0
        throw(ArgumentError("Cannot apply " * string(d.G.P[prod]) * ": there are non terminals in " * string(d.sf)))
    end
    for (pos, symbol) ∈ enumerate(d.sf)
        if symbol ∈ d.G.N
            if d.G.P[prod].left == symbol
                return next(d, prod, pos)
            else
                throw(ArgumentError("Cannot apply " * string(d.G.P[prod]) * ": the leftmost nonterminal of " * string(d.sf) * " is " * string(symbol)))
            end
        end
    end
    throw(ArgumentError("The specified production could not be found"))
end

"""
    leftmost(d::Derivation, prod::Production)::Derivation
Performs a *leftmost* derivation step.
Applies the specified [`Production`](@ref) to the current leftmost nonterminal in the sentential form.
"""
function leftmost(d::Derivation, prod::Production)::Derivation
    p = ensure_production_index(d,prod)
    return leftmost(d, p)
end

"""
    leftmost(d::Derivation, prod::AbstractArray{Production})::Derivation
Performs a *leftmost* derivation step.
Applies the specified productions to the current leftmost nonterminal in the sentential form, one by one.
"""
function leftmost(d::Derivation, prod::AbstractArray{Production})::Derivation
    res = d
    for i ∈ prod
        res = leftmost(res, i)
    end
    return res
end

"""
    leftmost(d::Derivation, prod::AbstractArray{Int})::Derivation
Performs a *leftmost* derivation step.
Applies the specified productions to the current leftmost nonterminal in the sentential form, one by one.
"""
function leftmost(d::Derivation, prod::AbstractArray{Int})::Derivation
    res = d
    for i ∈ prod
        res = leftmost(res, i)
    end
    return res
end

"""
    rightmost(d::Derivation, prod::Int)::Derivation
Performs a *rightmost* derivation step.
Applies the specified production(s) to the current rightmost nonterminal in the sentential form.
"""
function rightmost(d::Derivation, prod::Int)::Derivation
    if ~d.G.iscontextfree
        throw(ArgumentError("Cannot perform a rightmost derivation on a non context-free grammar"))
    end
    if length(d.sf) == 0
        throw(ArgumentError("Cannot apply " * string(d.G.P[prod]) * ": there are non terminals in " * string(d.sf)))
    end
    for (pos, symbol) ∈ (enumerate(d.sf) |> collect |> reverse)
        if symbol ∈ d.G.N
            if d.G.P[prod].left == symbol
                return next(d, prod, pos)
            else
                throw(ArgumentError("Cannot apply " * string(d.G.P[prod]) * ": the rightmost nonterminal of " * string(d.sf) * " is " * string(symbol)))
            end
        end
    end
    throw(ArgumentError("The specified production could not be found"))
end

"""
    rightmost(d::Derivation, prod::Production)::Derivation
Performs a *rightmost* derivation step.
Applies the specified [`Production`](@ref) to the current rightmost nonterminal in the sentential form.
"""
function rightmost(d::Derivation, prod::Production)::Derivation
    p = ensure_production_index(d,prod)
    return rightmost(d, p)
end

"""
    rightmost(d::Derivation, prod::Int)::Derivation
Performs a *rightmost* derivation step.
Applies the specified productions to the current rightmost nonterminal in the sentential form, one by one.
"""
function rightmost(d::Derivation, prod::AbstractArray{Int})::Derivation
    res = d
    for i ∈ prod
        res = rightmost(res, i)
    end
    return res
end

"""
    rightmost(d::Derivation, prod::AbstractArray{Production})::Derivation
Performs a *rightmost* derivation step.
Applies the specified productions to the current rightmost nonterminal in the sentential form, one by one.
"""
function rightmost(d::Derivation, prod::AbstractArray{Production})::Derivation
    res = d
    for i ∈ prod
        res = rightmost(res, i)
    end
    return res
end

"""
    possiblesteps(d::Derivation; prod::Union{Int, Nothing} = nothing, pos::Union{Int, Nothing} = nothing)
Returns all the possible steps that can be performed given the grammar and current *sentential form*.

Determines all the position of the *sentential form* that correspond to the left-hand side of
one of the production in the grammar, returning the position and production number. If a
production is specified, it yields only the pairs referring to it; similarly, if a position
is specified, it yields only the pairs referring to it.
"""
function possiblesteps(d::Derivation; prod::Union{Int, Nothing} = nothing, pos::Union{Int, Nothing} = nothing)
    res = []
    type0prods = map(x -> astype0(x), d.G.P)
    for (n, P) ∈ (prod === nothing ? enumerate(type0prods) : [(prod, type0prods[prod])])
        for p ∈ (pos === nothing ? range(1, length=(length(d.sf) - length(P.left) + 1)) : [pos])
            if [d.sf[p:p + length(P.left) - 1]] == P.left || P.left ∈ [d.sf[p:p + length(P.left) - 1]]
                push!(res, (n, p))
            end
        end
    end
    return res
end

"""
    sententialform(d::Derivation)
Returns the sentential form of the [`Derivation`](@ref).
"""
sententialform(d::Derivation) = d.sf

"""
    steps(d::Derivation)
Returns the steps performed by the [`Derivation`](@ref)
"""
steps(d::Derivation) = d.steps

ensure_production_index(d::Derivation, prod::Int)::Int = if 1 <= prod <= length(d.G.P) return prod else throw(ArgumentError("There is no production of index " * string(prod) * " in G")) end

ensure_production_index(d::Derivation, prod::Production)::Int = if prod ∈ d.G.P return findfirst(x -> x == prod, d.G.P)[1] else throw(ArgumentError("Production " * string(prod) * " does not belong to G")) end

### Operators ###

Base.show(io::IO, d::Derivation) = Base.show(io, d.repr)

Base.:(==)(x::Derivation, y::Derivation) = (x.G, x.steps) == (y.G, y.steps)

Base.hash(x::Derivation) = Base.hash((x.G, x.steps))


text:
using SearchablePDFs
using SearchablePDFs: searchable, require_extension, require_no_file
using Test
using Poppler_jll
using Aqua

TEST_PDF_PATH = joinpath(@__DIR__, "test.pdf")
TEST_PDF_RASTERIZED_PATH = joinpath(@__DIR__, "test_rasterized.pdf")

@testset "SearchablePDFs.jl" begin
    @test SearchablePDFs.num_pages(TEST_PDF_PATH) == 3

    for verbose in (false, true),
        apply_unpaper in (false, true),
        f in (searchable, ocr),
        opt in (true, false)

        kwargs = f === searchable ?
                 (; logfile=joinpath(@__DIR__, "test_logs.csv"), quiet=!verbose,
                  keep_intermediates=opt) :
                 (; verbose=verbose, max_files_per_unite=opt ? 2 : 100)

        result = f(TEST_PDF_RASTERIZED_PATH, joinpath(@__DIR__, "out.pdf"); apply_unpaper,
                   kwargs...)
        # make sure we delete the generated files eventually, even if the tests throw
        atexit(() -> rm(result.output_path; force=true))
        atexit(() -> rm(result.tmp; recursive=true, force=true))
        atexit(() -> rm(joinpath(@__DIR__, "test_logs.csv"); force=true))

        @test isfile(result.output_path)

        text = pdftotext() do exe
            return read(`$exe $(result.output_path) -`, String)
        end

        @test occursin("Chapter 9", text)
        @test occursin("evaluate expressions written in a source file", text)

        if f === searchable
            rm(joinpath(@__DIR__, "test_logs.csv"))
        end
        rm(result.output_path)
    end

    @testset "Errors" begin
        @test require_no_file("DOES_NOT_EXIST.pdf"; exception=true) === nothing
        @test_throws ArgumentError require_no_file("runtests.jl"; exception=true)

        @test require_extension("DOES_NOT_EXIST.pdf", ".pdf"; exception=true) === nothing
        @test_throws ArgumentError require_extension("runtests.jl", ".pdf"; exception=true)
    end
end

@testset "Aqua tests" begin
    Aqua.test_all(SearchablePDFs; ambiguities=false)
end


text:
matchesprefix(c::AbstractString, pre::AbstractString) = isempty(pre) || lowercase(c[1]) == lowercase(pre[1])
matchesprefix(c::Dict, pre::AbstractString) = matchesprefix(c[:text], pre)
matchesprefix(c, ::Nothing) = true

handle("completions") do data
  @destruct [path || nothing, mod || "Main", line, force] = data
  withpath(path) do
    pre = CodeTools.prefix(line)
    pre = isempty(pre) ? nothing : pre[end]
    m = getthing(mod)
    m = isa(m, Module) ? m : Main
    cs = CodeTools.completions(line, m, default = false)
    cs == nothing && pre == nothing && !force && (cs = [])
    d(:completions => cs,
      :prefix      => pre,
      :mod         => mod)
  end
end

handle("cacheCompletions") do mod
  m = getthing(mod)
  m = isa(m, Module) ? m : Main
  CodeTools.completions(m)
end


text:
module Io

export save_csv

"""
Save the real and imaginary parts of level np1
"""
function save_csv(;
      tc::Int64,
      mv::Int64,
      Rv::Vector{Float64},
      Yv::Vector{Float64},
      outdir::String,
      f)

   nx, ny = f.nx, f.ny

   open("$(outdir)/$(f.name)_$(mv)_re_$tc.csv","a") do out
      write(out,"R,Y,val\n")
      saved_once = false
      for i=1:nx
         for j=1:ny
            if abs(f.np1[i,j].re)>1e-16
               write(out,"$(Rv[i]),$(Yv[j]),$(f.np1[i,j].re)\n")
               saved_once=true
            end
         end
      end
      if !saved_once
         write(out,"$(Rv[1]),$(Yv[1]),0.0\n")
      end
   end

   open("$(outdir)/$(f.name)_$(mv)_im_$tc.csv","a") do out
      write(out,"R,Y,val\n")
      saved_once = false
      for i=1:nx
         for j=1:ny
            if abs(f.np1[i,j].im)>1e-16
               write(out,"$(Rv[i]),$(Yv[j]),$(f.np1[i,j].im)\n")
               saved_once=true
            end
         end
      end
      if !saved_once
         write(out,"$(Rv[1]),$(Yv[1]),0.0\n")
      end
   end

   return nothing
end

end


text:
"""
    generate([pkg::AbstractString]) -> Template

Shortcut for `Template(; interactive=true)(pkg)`.
If no package name is supplied, you will be prompted for one.
"""
function generate(pkg::AbstractString=prompt(Template, String, :pkg))
    t = Template(; interactive=true)
    t(pkg)
    return t
end

"""
    interactive(T::Type{<:Plugin}) -> T

Interactively create a plugin of type `T`. Implement this method and ignore other
related functions only if you want completely custom behaviour.
"""
function interactive(T::Type)
    pairs = Vector{Pair{Symbol, Type}}(interactive_pairs(T))

    # There must be at least 2 MultiSelectMenu options.
    # If there are none, return immediately.
    # If there's just one, add a "dummy" option.
    isempty(pairs) && return T()
    just_one = length(pairs) == 1
    just_one && push!(pairs, :None => Nothing)

    menu = MultiSelectMenu(
        collect(map(pair -> string(first(pair)), pairs));
        pagesize=length(pairs),
    )
    println("$(nameof(T)) keywords to customize:")
    customize = sort!(collect(request(menu)))

    # If the "None" option was selected, don't customize anything.
    just_one && lastindex(pairs) in customize && return T()

    kwargs = Dict{Symbol, Any}()
    foreach(pairs[customize]) do (name, F)
        kwargs[name] = prompt(T, F, name)
    end
    return T(; kwargs...)
end

struct NotCustomizable end

"""
    customizable(::Type{<:Plugin}) -> Vector{Pair{Symbol, DataType}}

Return a list of keyword arguments that the given plugin type accepts,
which are not fields of the type, and should be customizable in interactive mode.
For example, for a constructor `Foo(; x::Bool)`, provide `[x => Bool]`.
If `T` has fields which should not be customizable, use `NotCustomizable` as the type.
"""
customizable(::Type) = ()

function pretty_message(s::AbstractString)
    replacements = [
        r"Array{(.*?),1}" => s"Vector{\1}",
        r"Union{Nothing, (.*?)}" => s"Union{\1, Nothing}",
    ]
    return reduce((s, p) -> replace(s, p), replacements; init=s)
end

"""
    input_tips(::Type{T}) -> Vector{String}

Provide some extra tips to users on how to structure their input for the type `T`,
for example if multiple delimited values are expected.
"""
input_tips(::Type{Vector{T}}) where T = ["comma-delimited", input_tips(T)...]
input_tips(::Type{Nothing}) = String[]
input_tips(::Type{Union{T, Nothing}}) where T = ["'nothing' for nothing", input_tips(T)...]
input_tips(::Type{Secret}) = ["name only"]
input_tips(::Type) = String[]

"""
    convert_input(::Type{P}, ::Type{T}, s::AbstractString) -> T

Convert the user input `s` into an instance of `T` for plugin of type `P`.
A default implementation of `T(s)` exists.
"""
convert_input(::Type, T::Type{<:Real}, s::AbstractString) = parse(T, s)
convert_input(::Type, T::Type, s::AbstractString) = T(s)

function convert_input(P::Type, ::Type{Union{T, Nothing}}, s::AbstractString) where T
    # This is kind of sketchy because technically, there might be some other input
    # whose value we want to instantiate with the string "nothing",
    # but I think that would be a pretty rare occurrence.
    # If that really happens, they can just override this method.
    return s == "nothing" ? nothing : convert_input(P, T, s)
end

function convert_input(::Type, ::Type{Bool}, s::AbstractString)
    s = lowercase(s)
    return if startswith(s, 't') || startswith(s, 'y')
        true
    elseif startswith(s, 'f') || startswith(s, 'n')
        false
    else
        throw(ArgumentError("Unrecognized boolean response"))
    end
end

function convert_input(P::Type, T::Type{<:Vector}, s::AbstractString)
    startswith(s, '[') && endswith(s, ']') && (s = s[2:end-1])
    xs = map(strip, split(s, ","))
    return map(x -> convert_input(P, eltype(T), x), xs)
end

"""
    prompt(::Type{P}, ::Type{T}, ::Val{name::Symbol}) -> Any

Prompts for an input of type `T` for field `name` of plugin type `P`.
Implement this method to customize particular fields of particular types.
"""
prompt(P::Type, T::Type, name::Symbol) = prompt(P, T, Val(name))

# The trailing `nothing` is a hack for `fallback_prompt` to use, ignore it.
function prompt(P::Type, ::Type{T}, ::Val{name}, ::Nothing=nothing) where {T, name}
    tips = join([T; input_tips(T); "default=$(repr(defaultkw(P, name)))"], ", ")
    default = defaultkw(P, name)
    input = Base.prompt(pretty_message("Enter value for '$name' ($tips)"))
    input === nothing && throw(InterruptException())
    return if isempty(input)
        default
    else
        try
            # Working around what appears to be a bug in Julia 1.0:
            # #145#issuecomment-623049535
            if VERSION < v"1.1" && T isa Union && Nothing <: T
                if input == "nothing"
                    nothing
                else
                    convert_input(P, T.a === Nothing ? T.b : T.a, input)
                end
            else
                convert_input(P, T, input)
            end
        catch ex
            ex isa InterruptException && rethrow()
            @warn "Invalid input" ex
            prompt(P, T, name)
        end
    end
end

# Compute all the concrete subtypes of T.
concretes_rec(T::Type) = isabstracttype(T) ? vcat(map(concretes_rec, subtypes(T))...) : Any[T]
concretes(T::Type) = sort!(concretes_rec(T); by=nameof)

# Compute name => type pairs for T's interactive options.
function interactive_pairs(T::Type)
    pairs = collect(map(name -> name => fieldtype(T, name), fieldnames(T)))

    # Use prepend! here so that users can override field types if they wish.
    prepend!(pairs, reverse(customizable(T)))
    uniqueby!(first, pairs)
    filter!(p -> last(p) !== NotCustomizable, pairs)
    sort!(pairs; by=first)

    return pairs
end

# unique!(f, xs) added here: https://github.com/JuliaLang/julia/pull/30141
if VERSION >= v"1.1"
    const uniqueby! = unique!
else
    function uniqueby!(f, xs)
        seen = Set()
        todelete = Int[]
        foreach(enumerate(map(f, xs))) do (i, out)
            out in seen && push!(todelete, i)
            push!(seen, out)
        end
        return deleteat!(xs, todelete)
    end
end


text:
using Test, LinearAlgebra, OrderedCollections
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Single element" begin
    pde  = op = Helmholtz(dim=2,k=1)
    el   = LagrangeLine((1,0),(1.1,0))
    q    = GaussLegendre(20)
    shand = Kress()
    qs  = SingularQuadratureRule(q,shand)
    x̂,ŵ = q()
    yi  = map(el,x̂)
    νi  = map(x->normal(el,x),x̂)
    G   = SingleLayerKernel(op)
    ##
    ϕ  = y->cos(y[1])  # regular density
    s  = 0.51 # location of singularity in parameter space
    x  = el(s)
    ν  = normal(el,s)
    k  = let el = el, G=G
        v -> G(x,el(v))*measure(el,v)
    end
    w  = Nystrom.singular_weights(k,q,qs,s)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia  = transpose(w)*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
    y₀  = el(s)
    x   = x + 1e-2*ν # outside
    w  = Nystrom.singular_weights(k,q,qs,s)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia = transpose(w)*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
    x   = x - 1e-2*ν # inside
    w  = Nystrom.singular_weights(k,q,qs,s)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia = transpose(w)*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
end

@testset "Single obstacle" begin
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    # geo  = Circle()
    Geometry.clear!()
    Ω,M   = WaveProp.IO.gmsh_disk()
    Γ     = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    mesh  = NystromMesh(view(M,Γ);order=5)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    Sa    = Nystrom.assemble(S)
    Da    = Nystrom.assemble(D)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Sa,Da,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
    # K     = AdjointDoubleLayerOperator(pde,mesh)
    # H     = HyperSingularOperator(pde,mesh)
    # e0    = WaveProp.Nystrom.error_interior_derivative_green_identity(K,H,γ₀u,γ₁u)
    # Kldim  = Nystrom.assemble_ldim(K)
    # Hldim  = Nystrom.assemble_ldim(H)
    # e1    = WaveProp.Nystrom.error_interior_derivative_green_identity(Kldim,Hldim,γ₀u,γ₁u)
    # @test 10*norm(e1,Inf) < norm(e0,Inf)
    # @test norm(e1,Inf) < 1e-4
end

@testset "Single obstacle" begin
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    Geometry.clear!()
    geo  = Kite()
    Ω    = Domain(geo)
    M    = meshgen(Ω,n=100)
    Γ    = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    mesh    = NystromMesh(view(M,Γ);order=5)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    Sa    = Nystrom.assemble(S)
    Da    = Nystrom.assemble(D)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Sa,Da,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp