text:
using DomainSets: ScalarAffineMap,
    VectorAffineMap,
    StaticAffineMap,
    GenericAffineMap,
    ScalarLinearMap,
    VectorLinearMap,
    StaticLinearMap,
    GenericLinearMap,
    ScalarTranslation,
    VectorTranslation,
    StaticTranslation,
    GenericTranslation,
    TupleProductMap

maps_to_test(T) = [
    StaticIdentityMap{T}(),
    VectorIdentityMap{T}(10),
    ConstantMap{T}(one(T)),
    ConstantMap{T}(SVector{2,T}(1,2)),
    ZeroMap{T}(),
    UnityMap{T}(),
    AffineMap(T(1.2), T(2.4)), # scalar map
    AffineMap(-T(1.2), T(2.4)), # scalar map with negative A
    AffineMap(randvec(T, 2, 2), randvec(T, 2)), # static map
    AffineMap(randvec(T, 3, 2), randvec(T, 3)), # static map, rectangular
    AffineMap(rand(T, 2, 2), rand(T, 2)), # vector map
    AffineMap(rand(T, 3, 2), rand(T, 3)), # vector map, rectangular
    AffineMap(LinearAlgebra.I, one(T)/2),  # use UniformScaling object as A
    AffineMap(LinearAlgebra.I, randvec(T,2)),  # use UniformScaling object as A
    DomainSets.GenericAffineMap(randvec(T, 2, 2), randvec(T, 2)),
    DomainSets.GenericAffineMap(T(1.2), randvec(T, 2)),
    DomainSets.GenericAffineMap(randvec(T, 3, 2), randvec(T, 3)),
    Translation(randvec(T, 3)),
    LinearMap(randvec(T, 2, 2)),
    LinearMap(randvec(T, 2)),
    LinearMap(randvec(T, 2, 2)) ∘ AffineMap(T(1.2), randvec(T, 2)),
    mapfrom_parameterdomain(T(5) * ComplexUnitCircle{T}() .+ T(2))
]

randvec(T,n) = SVector{n,T}(rand(n))
randvec(T,m,n) = SMatrix{m,n,T}(rand(m,n))

suitable_point_to_map(m::Map) = suitable_point_to_map(m, domaintype(m))

suitable_point_to_map(m::Map, ::Type{SVector{N,T}}) where {N,T} = SVector{N,T}(rand(N))
suitable_point_to_map(m::Map, ::Type{T}) where {T<:Number} = rand(T)
suitable_point_to_map(m::Map, ::Type{<:AbstractVector{T}}) where {T} = rand(T, mapsize(m,2))

suitable_point_to_map(m::DomainSets.ProductMap) =
    map(suitable_point_to_map, components(m))
suitable_point_to_map(m::DomainSets.VcatMap{N,T}) where {N,T} =
    SVector{N,T}(rand(T,N))

suitable_point_to_map(::CartToPolarMap{T}) where {T} = randvec(T,2)
suitable_point_to_map(::PolarToCartMap{T}) where {T} = randvec(T,2)

widertype(T) = widen(T)
widertype(::Type{SVector{N,T}}) where {N,T} = SVector{N,widen(T)}
widertype(::Type{Vector{T}}) where {T} = Vector{widen(T)}
widertype(::Type{Tuple{A}}) where {A} = Tuple{widen(A)}
widertype(::Type{Tuple{A,B}}) where {A,B} = Tuple{widen(A),widen(B)}
widertype(::Type{Tuple{A,B,C}}) where {A,B,C} = Tuple{widen(A),widen(B),widen(C)}
widertype(::Type{NTuple{N,T}}) where {N,T} = NTuple{N,widen(T)}

issquarematrix(A) = false
issquarematrix(A::AbstractArray) = size(A,1)==size(A,2)

function test_generic_inverse(m)
    M = mapsize(m,1)
    N = mapsize(m,2)
    x = suitable_point_to_map(m)
    y = m(x)

    if M==N
        inverse(m, y)      # trigger exception outside test if not implemented

        minv = inverse(m)
        @test minv(y) ≈ x
        @test inverse(m)(y) ≈ x
        @test inverse(m, y) ≈ x
        @test m\y ≈ x
        @test DomainSets.LazyInverse(m)(y) ≈ inverse(m, y)
    end
    if M >= N && numtype(m)!=BigFloat
        leftinverse(m, y)      # trigger exception outside test if not implemented

        mli = leftinverse(m)
        @test mli(y) ≈ x
        @test leftinverse(m, y) ≈ x
    end
    if M <= N && numtype(m)!=BigFloat
        rightinverse(m, y)      # trigger exception outside test if not implemented

        mri = rightinverse(m)
        @test m(mri(y)) ≈ y
        @test m(rightinverse(m, y)) ≈ y
    end
end

function test_generic_jacobian(m)
    x = suitable_point_to_map(m)
    # Trigger exception outside of test if jacobian(m, x) is not implemented.
    # test_generic_jacobian(m) should be called within a try/catch block
    jacobian(m, x)

    δ = sqrt(eps(prectype(m)))
    x2 = x .+ δ
    # we intentionally test jacobian(m, x) before testing jacobian(m)
    if !(m isa ProductMap)
        @test norm(m(x2) - m(x) + jacobian(m, x)*(x-x2)) < 100δ
    end
    jac = jacobian(m)
    @test jac(x) ≈ jacobian(m, x)
    j = jacobian(m, x)
    @test size(j) == mapsize(m)
    if j isa AbstractArray{Float64}
        y = similar(j)
        DomainSets.jacobian!(y, m, x)
        @test y ≈ jac(x)
    end
    if issquarematrix(jac(x))
        @test jacdet(m, x) ≈ det(jacobian(m, x))
    end
    if mapsize(m) == ()
        d = diffvolume(m, x)
        s = sqrt(det(transpose(jacobian(m,x))*jacobian(m,x)))
        @test d ≈ s || d ≈ -s
    else
        @test abs(diffvolume(m, x) - sqrt(det(transpose(jacobian(m,x))*jacobian(m,x)))) < 1e-7
    end
end

# generic map test suite
function test_generic_map(m)
    @test convert(Map{domaintype(m)}, m) == m

    x = suitable_point_to_map(m)
    @test applymap(m, x) == m(x)
    y = m(x)
    @test isreal(m) == (isreal(x) && isreal(y))

    S = domaintype(m)
    U = codomaintype(m)
    @test x isa S
    @test y isa U

    if isaffine(m) && !isconstant(m) && !(prectype(m) == BigFloat)
        test_generic_inverse(m)
    else
        try
            # The map may not support an inverse, let's try
            test_generic_inverse(m)
        catch e
            # make sure it failed because the inverse is not defined
            @test e isa MethodError || prectype(m) == BigFloat
        end
    end

    if isaffine(m)
        M = matrix(m)
        @test size(M) == mapsize(m)
        test_generic_jacobian(m)
    else
        try # jacobian may not be implemented
            test_generic_jacobian(m)
        catch e
        end
    end

    if domaintype(m) == Float64
        @test convert(Map{BigFloat}, m) isa Map{BigFloat}
        @test convert(Map{BigFloat}, m) == m
    end
    if prectype(m) == Float64
        U = widertype(domaintype(m))
        @test convert(Map{U}, m) isa Map{U}
        @test convert(Map{U}, m) == m
    end

    if isaffine(m)
        A = matrix(m)
        b = vector(m)
        @test m(x) == A*x+b
    end
end

function test_isomorphisms(T)
    m1 = DomainSets.VectorToNumber{T}()
    @test m1(SA[1.0]) == 1.0
    m1b = DomainSets.NumberToVector{T}()
    @test inverse(m1) == m1b
    @test inverse(m1, 0.4) == m1b(0.4)
    @test inverse(m1b) == m1
    @test inverse(m1b, SA[0.4]) == m1(SA[0.4])
    @test m1b(1.0) == SA[1.0]
    @test mapsize(m1) == (1,1)
    @test mapsize(m1b) == (1,)
    @test mapsize(m1 ∘ inverse(m1)) == ()
    @test mapsize(inverse(m1) ∘ m1) == (1,1)
    @test jacobian(m1 ∘ inverse(m1), 1) == one(T)
    @test jacobian(inverse(m1) ∘ m1, [1]) == ones(T,1,1)
    @test jacobian(m1) isa ConstantMap
    @test jacobian(m1b) isa ConstantMap

    m2 = DomainSets.VectorToComplex{T}()
    @test m2(SA[one(T), one(T)]) == 1 + im
    m2b = DomainSets.ComplexToVector{T}()
    @test inverse(m2) == m2b
    @test inverse(m2, zero(T)+0im) == m2b(zero(T)+0im)
    @test inverse(m2b) == m2
    @test inverse(m2b, SA[one(T),one(T)]) == one(T)+one(T)*im
    @test m2b(one(T)+one(T)*im) == SA[one(T),one(T)]
    @test mapsize(m2) == (1,2)
    @test mapsize(m2b) == (2,)
    @test mapsize(m2 ∘ m2b) == ()
    @test mapsize(m2b ∘ m2) == (2,2)

    m3 = DomainSets.VectorToTuple{2,T}()
    m3b = DomainSets.TupleToVector{2,T}()
    @test m3(SA[one(T), one(T)]) == (one(T),one(T))
    @test inverse(m3b, SA[one(T), one(T)]) == (one(T),one(T))
    @test inverse(m3) == m3b
    @test inverse(m3b) == m3
    @test m3b( (one(T),one(T)) ) == SA[one(T),one(T)]
    @test inverse(m3, (one(T),one(T)) ) == SA[one(T),one(T)]

    m4 = DomainSets.NestedToFlat{3,T,Tuple{Tuple{T,T},T},(2,1)}()
    m4b = DomainSets.FlatToNested{3,T,Tuple{Tuple{T,T},T},(2,1)}()
    x4 = ([T(1),T(2)],T(3))
    y4 = T[1,2,3]
    @test m4(x4) == y4
    @test m4b(y4) == x4
    @test inverse(m4) == m4b
    @test inverse(m4b) == m4
    @test inverse(m4, y4) == x4
    @test inverse(m4b, x4) == y4
end

function test_maps(T)
    generic_tests(T)

    # Test special maps
    test_identity_map(T)
    test_basic_maps(T)
    test_affine_maps(T)
    test_composite_map(T)
    test_product_map(T)
    test_wrapped_maps(T)
    test_scaling_maps(T)
    test_isomorphisms(T)
end

function generic_tests(T)
    for map in maps_to_test(T)
        @test prectype(map) == T
        test_generic_map(map)
    end
    # Test an additional composite map
    m1 = LinearMap(SMatrix{2,2}(1,2,3,4.0))
    m2 = CartToPolarMap()
    cmap = m1 ∘ m2 ∘ m1
    test_generic_map(cmap)
end

function test_affine_maps(T)
    A = rand(T,2,2)
    @test DomainSets.to_matrix(Vector{T}, A) == A
    @test DomainSets.to_matrix(T, 2) == 2
    @test DomainSets.to_matrix(SVector{2,T}, 2) == SMatrix{2,2}(2,0,0,2)
    @test DomainSets.to_matrix(SVector{2,T}, LinearAlgebra.I) == SMatrix{2,2}(1,0,0,1)
    @test DomainSets.to_matrix(Vector{T}, 2) == UniformScaling(2)
    @test DomainSets.to_matrix(Vector{T}, LinearAlgebra.I) == LinearAlgebra.I
    # test fallback with nonsensical call
    @test DomainSets.to_matrix(Tuple{Int}, 2) == 2

    @test DomainSets.to_matrix(T, A, 2) == A
    @test DomainSets.to_matrix(T, 2, 3) == 2
    @test DomainSets.to_matrix(T, UniformScaling(2), 3) == 2
    @test DomainSets.to_matrix(SVector{2,T}, 2, SVector(1,1)) == SMatrix{2,2}(2,0,0,2)
    @test DomainSets.to_matrix(Vector{T}, 2, [1,2]) == [2 0 ; 0 2]

    @test DomainSets.to_vector(T, 2) == 0
    @test DomainSets.to_vector(SVector{2,T}, 2) == SVector(0,0)
    @test DomainSets.to_vector(Vector{T}, A) == [0,0]
    @test DomainSets.to_vector(T, 2, 3) == 3

    if T != BigFloat    # BigFloat's make pinv fail for StaticArrays
        @test DomainSets.matrix_pinv(SMatrix{2,2}(rand(T),rand(T),rand(T),rand(T))) isa SMatrix{2,2}
        @test DomainSets.matrix_pinv(SVector(rand(T),rand(T))) isa Transpose{T,SVector{2,T}}
    end

    test_linearmap(T)
    test_translation(T)
    test_affinemap(T)
end

function test_linearmap(T)
    m1 = LinearMap(2one(T))
    @test m1 isa LinearMap{T}
    @test domaintype(m1) == T
    @test islinear(m1)
    @test isaffine(m1)
    @test m1(3) == 6
    @test m1(3one(T)) == 6
    @test m1(3one(widen(T))) isa widen(T)
    @test m1(SVector(1,2)) == SVector(2,4)
    @test m1([1.0,2.0]) == [2.0,4.0]
    mw1 = convert(Map{widen(T)}, m1)
    @test mw1 isa Map{widen(T)}
    @test mw1.A isa widen(T)
    @test jacobian(m1) isa ConstantMap{T}
    @test jacobian(m1, 1) == 2
    @test LinearMap(one(T)) == StaticIdentityMap{T}()

    m2 = LinearMap(2)
    @test domaintype(m2) == Int
    @test m2 === ScalarLinearMap(2)
    @test m2(one(T)) isa T
    @test m2 == convert(Map, 2)
    @test jacobian(m2, 1) == 2
    @test jacobian(m2) isa ConstantMap{Int}
    @test jacobian(m2, 1) == 2
    @test LinearMap(1) == StaticIdentityMap{T}()

    m3 = LinearMap(SMatrix{2,2}(one(T), 2one(T), 3one(T), 4one(T)))
    @test m3 isa LinearMap{SVector{2,T}}
    @test m3 === StaticLinearMap(m3.A)
    @test m3(SVector(1,2)) == SVector(7, 10)
    @test m3(SVector{2,T}(1,2)) == SVector{2,T}(7, 10)
    @test m3 ∘ m3 isa LinearMap
    @test LinearMap(SA[1 0; 0 1]) == StaticIdentityMap{domaintype(m3)}()

    A = rand(T,2,2)
    m4 = LinearMap(A)
    @test m4 isa LinearMap{Vector{T}}
    @test m4([1,2]) ==  A * [1,2]
    y = zeros(T,2)
    @test (DomainSets.applymap!(y, m4, [1,2]); y == A * [1,2])
    @test jacobian(m4, [1,2]) == A

    m5 = LinearMap{Vector{T}}(UniformScaling(2*one(T)))
    @test m5 isa LinearMap{Vector{T}}
    @test m5([1,2]) ==  2 * [1,2]
    @test jacobian(m5, [1,2]) == UniformScaling(2)

    m6 = LinearMap{SVector{2,T}}(2)
    @test m6 isa GenericLinearMap{SVector{2,T},T}
    @test m6(SA[one(T),one(T)]) == [2,2]
    @test leftinverse(m6) isa LinearMap{SVector{2,T}}
    @test leftinverse(m6)([2,2]) == [1,1]

    # Test construction and conversion
    @test LinearMap{T}(1) isa ScalarLinearMap{T}
    @test LinearMap{Vector{T}}(rand(Int,5,5)) isa VectorLinearMap{T}
    @test LinearMap{SVector{2,T}}(SMatrix{2,2}(1, 2, 3, 4)) isa StaticLinearMap{T,2,2,4}

    @test GenericLinearMap(one(T)) isa GenericLinearMap{T,T}
    @test GenericLinearMap(SMatrix{2,2,T}(1,1,1,1)) isa GenericLinearMap{SVector{2,T},SMatrix{2,2,T,4}}

    @test VectorLinearMap(rand(T,3,3)) isa VectorLinearMap{T}
    @test StaticLinearMap(SMatrix{2,2,T}(1,2,3,4)) isa StaticLinearMap{T,2,2,4}

    @test convert(Map{SVector{2,T}}, LinearMap(rand(T,2,2))) isa StaticLinearMap{T,2,2,4}
    @test convert(Map{T}, 1) isa LinearMap{T}
end

function test_translation(T)
    v = randvec(T,3)
    m = Translation(v)
    @test !islinear(m)
    @test isaffine(m)
    @test inverse(inverse(m)) == m
    @test jacobian(m) isa ConstantMap
    @test vector(jacobian(m)) == [1 0 0; 0 1 0; 0 0 1]
    @test jacdet(m) isa UnityMap{SVector{3,T},T}

    @test mapsize(Translation(one(T))) == ()

    # Test construction and conversion
    @test Translation(one(T)) isa ScalarTranslation{T}
    @test Translation{T}(1) isa ScalarTranslation{T}
    @test Translation(SVector{2,T}(1,2)) isa StaticTranslation{T,2}
    @test Translation{SVector{2,T}}(rand(T,2)) isa StaticTranslation{T,2}
    @test Translation(rand(Int,2)) isa VectorTranslation{Int}
    @test Translation{Vector{T}}(rand(Int,2)) isa VectorTranslation{T}
    @test Translation(1:10) isa GenericTranslation{Vector{Int}}
    @test Translation{Vector{T}}(1:10) isa GenericTranslation{Vector{T}}

    @test StaticTranslation(MVector(one(T),2)) isa StaticTranslation{T,2}
    @test StaticTranslation{T}(SVector(1,2)) isa StaticTranslation{T,2}

    @test VectorTranslation(1:10) isa VectorTranslation{Int}
    @test VectorTranslation{T}(1:10) isa VectorTranslation{T}

    @test GenericTranslation(1:10) isa GenericTranslation{Vector{Int}}
    @test GenericTranslation{Vector{T}}(1:10) isa GenericTranslation{Vector{T}}

    @test GenericTranslation(one(T)) isa GenericTranslation{T,T}
    @test GenericTranslation{T}(1) isa GenericTranslation{T,T}
    @test_throws InexactError GenericTranslation{Int}(one(T)/2)

    @test convert(Map{SVector{2,T}}, Translation(rand(T,2))) isa StaticTranslation{T,2}
end

function test_affinemap(T)
    m1 = AffineMap(T(2), T(3))
    @test m1 isa ScalarAffineMap{T}
    @test !islinear(m1)
    @test isaffine(m1)
    @test m1(2) == 7

    @test convert(ScalarAffineMap{BigFloat}, m1) == m1

    @test m1 ∘ m1 isa AffineMap
    @test (m1 ∘ m1)(2) == 2*(2*2+3)+3

    m2 = AffineMap(T(2), SVector{2,T}(1,2))
    @test m2 isa GenericAffineMap{SVector{2,T}}
    @test m2(SVector(1,2)) == SVector(3,6)
    @test mapsize(m2) == (2,2)

    m3 = AffineMap(UniformScaling(2*one(T)), [one(T),2*one(T)])
    @test m3 isa AffineMap{Vector{T}}
    @test mapsize(m3) == (2,2)
    @test m3([1,2]) ==  2 * [1,2] + [1,2]
    y = zeros(T,2)
    @test (DomainSets.applymap!(y, m3, [1,2]); y == m3([1,2]))
    @test jacobian(m3, [1,2]) == [2 0; 0 2]
    @test jacdet(m3, [1,2]) == 4

    @test AffineMap(LinearAlgebra.I, 2.0) isa ScalarAffineMap
    m4 = GenericAffineMap(LinearAlgebra.I, 2.0)
    @test m4.A isa UniformScaling{Float64}
    @test jacdet(m4, 3.0) == 1.0

    @test mapsize(AffineMap(rand(3),rand(3))) == (3,)
    @test mapsize(AffineMap(LinearAlgebra.I,2)) == ()

    # Test construction and conversion
    @test AffineMap(1, 2*one(T)) isa ScalarAffineMap{T}
    @test AffineMap{BigFloat}(1, 2.0) isa ScalarAffineMap{BigFloat}
    @test AffineMap(rand(T,2),rand(T,2)) isa GenericAffineMap{T}
    @test AffineMap(rand(T,2,2),rand(T,2)) isa VectorAffineMap{T}
    @test AffineMap{Vector{T}}(rand(Int,2,2),rand(Int,2)) isa VectorAffineMap{T}
    @test AffineMap{SVector{2,T}}(rand(T,2,2),rand(T,2)) isa StaticAffineMap{T,2,2,4}
    @test AffineMap{SVector{2,T}}(rand(Int,2,2),rand(T,2)) isa StaticAffineMap{T,2,2,4}
    @test AffineMap{SVector{2,T}}(rand(Int,2,2),rand(Int,2)) isa StaticAffineMap{T,2,2,4}
    @test AffineMap{SVector{2,T}}(1,rand(Int,2)) isa GenericAffineMap{SVector{2,T},T,SVector{2,T}}
    @test AffineMap{SVector{2,T}}(SMatrix{3,2}(3,2,1,4,5,6),SVector(1,2,3)) isa StaticAffineMap{T,2,3,6}

    @test GenericAffineMap(SMatrix{2,2,T}(1,2,3,4), [1,2]) isa GenericAffineMap{SVector{2,T},SMatrix{2,2,T,4}}
    @test GenericAffineMap(SMatrix{1,1,Int}(1),SVector{1,T}(1)) isa GenericAffineMap{SVector{1,T}}
    @test GenericAffineMap(rand(T,2,2),rand(Int,2)) isa GenericAffineMap{Vector{T},Matrix{T},Vector{T}}
    @test GenericAffineMap{Vector{T}}(1, [1,2]) isa GenericAffineMap{Vector{T},T,Vector{Int}}
    @test GenericAffineMap{Vector{T}}(rand(Int,2,2), rand(Int,2)) isa GenericAffineMap{Vector{T},Matrix{T},Vector{T}}
    @test convert(GenericAffineMap{Vector{T}}, GenericAffineMap(rand(Int,2,2), rand(Int,2))) isa GenericAffineMap{Vector{T}}
    @test GenericAffineMap(rand(T,2),rand(T,2)) isa GenericAffineMap{T}
    @test GenericAffineMap{T}(rand(Int,2),rand(Int,2)) isa GenericAffineMap{T}

    @test VectorAffineMap(rand(T,2,2),rand(T,2)) isa VectorAffineMap{T}

    @test StaticAffineMap(SMatrix{1,1,T}(1),SVector{1,T}(1)) isa StaticAffineMap{T}
    @test StaticAffineMap(SMatrix{1,1,T}(1),SVector{1,Int}(1)) isa StaticAffineMap{T}

    @test StaticAffineMap(rand(T,2,2),SVector{2,T}(1,1)) isa StaticAffineMap{T}
    @test_throws AssertionError StaticAffineMap(rand(T,3,2),SVector{2,T}(1,1))
    @test StaticAffineMap(SMatrix{2,2,T}(1,2,3,4),[1,1]) isa StaticAffineMap{T,2}
    @test StaticAffineMap(SMatrix{3,2,T}(1,2,3,4,5,6),SVector{3,T}(1,2,3)) isa StaticAffineMap{T,2,3,6}

    @test convert(Map{SVector{2,T}}, AffineMap(rand(2,2),rand(2))) isa StaticAffineMap{T,2,2,4}
end

function test_scaling_maps(T)
    test_generic_map(LinearMap(T(2)))
    test_generic_map(LinearMap(T(2), T(3)))
    test_generic_map(LinearMap(T(2), T(3), T(4)))
    test_generic_map(LinearMap(T(2), T(3), T(4), T(5)))
end

function test_identity_map(T)
    i1 = StaticIdentityMap{T}()
    i2 = StaticIdentityMap{SVector{2,T}}()
    test_generic_map(i1)
    test_generic_map(i2)
    @test i1 == i2
    @test hash(i1) == hash(i2)
    @test islinear(i1)
    @test isaffine(i1)
    @test convert(StaticIdentityMap{SVector{2,T}}, i1) === i2
    @test jacobian(i1) isa ConstantMap
    @test jacobian(i1, 1) == 1
    @test jacdet(i1, 1) == 1
    m1 = convert(DomainSets.AbstractAffineMap{T}, i1)
    @test m1 isa LinearMap{T}
    @test jacdet(m1, 1) == 1
    @test convert(DomainSets.AbstractAffineMap, i1) isa LinearMap{T}
    m2 = convert(DomainSets.LinearMap{T}, i1)
    @test m2 isa LinearMap{T}
    @test jacdet(m2, 1) == 1

    i3 = VectorIdentityMap{T}(10)
    test_generic_map(i3)
    r = rand(T, 10)
    @test i3(r) ≈ r
    @test hash(i3) == hash(VectorIdentityMap{Int}(10))

    @test IdentityMap() ∘ LinearMap(2) == LinearMap(2.0)
end

function test_basic_maps(T)
    @test UnityMap{SVector{2,Float64}}() == UnityMap{SVector{2,Float64},Float64}()
    @test hash(ConstantMap(2)) == hash(ConstantMap(2.0))
    @test DomainSets.absmap(ConstantMap(-2)) == ConstantMap(2)
end

function test_composite_map(T)
    a = T(0)
    b = T(1)
    c = T(2)
    d = T(3)
    ma = StaticIdentityMap{T}()
    mb = DomainSets.interval_map(a, b, c, d)

    r = suitable_point_to_map(ma)
    m1 = ma∘mb
    test_generic_map(m1)
    @test m1(r) ≈ ma(mb(r))
    m2 = m1∘mb
    test_generic_map(m2)
    @test m2(r) ≈ m1(mb(r))
    m3 = mb∘m2
    test_generic_map(m3)
    @test m3(r) ≈ mb(m2(r))
    m = m2∘m3
    test_generic_map(m)
    @test m(r) ≈ m2(m3(r))

    m5 = ComposedMap(LinearMap(rand(T,2,2)), AffineMap(rand(T,2,2),rand(T,2)))
    test_generic_map(m5)
    @test jacobian(m5) isa ConstantMap
    @test m5[Component(1)] isa LinearMap
    @test m5[Component(2)] isa AffineMap
    @test ComposedMap(m5[Component(1:2)]...) == m5
    @test_throws BoundsError m5[Component(3)]

    m6 = DomainSets.multiply_map(ma,ma)
    @test m6(one(T)/2) == one(T)/4
    @test jacobian(m6) isa DomainSets.SumMap
    @test jacobian(m6)(one(T)) == 2
    @test jacobian(m6, one(T)) == 2

    m7 = DomainSets.sum_map(ma,ma)
    @test m7(one(T)) == 2
    @test jacobian(m7) isa ConstantMap
    @test jacobian(m7, one(T)) == 2

    @test mapsize(ComposedMap(LinearMap(2),LinearMap(rand(T,2)),LinearMap(rand(T,2,2)))) == (2,)
    @test mapsize(ComposedMap(LinearMap(rand(T,2)),LinearMap(rand(T,2)'),LinearMap(rand(T,2)))) == (2,)
    @test mapsize(ComposedMap(LinearMap(rand(T,2,2)),LinearMap(rand(T,2)'),LinearMap(2))) == (1,2)
    @test mapsize(ComposedMap(LinearMap(one(T)),LinearMap(one(T)))) == ()

    @test DomainSets.composedmap() == ()
    @test DomainSets.composedmap(ma) == ma
    @test DomainSets.composedmap(ma,ma) == ma
    @test DomainSets.composedmap(ma,ma,ma) == ma

    @test DomainSets.composite_jacobian(ma) == jacobian(ma)

    @test DomainSets.multiply_map() == ()
    @test DomainSets.multiply_map(ma) == ma
    @test DomainSets.multiply_map(ma,ma)(2*one(T)) == 4
    @test DomainSets.multiply_map(ma,ma,ma)(2*one(T)) == 8

    @test DomainSets.sum_jacobian() == ()
    @test DomainSets.sum_jacobian(ma) == jacobian(ma)
end

function test_product_map(T)
    ma = StaticIdentityMap{T}()
    mb = DomainSets.interval_map(T(0), T(1), T(2), T(3))

    r1 = suitable_point_to_map(ma)
    r2 = suitable_point_to_map(ma)
    r3 = suitable_point_to_map(ma)
    r4 = suitable_point_to_map(ma)
    r5 = suitable_point_to_map(ma)

    m1 = productmap(ma,mb)
    test_generic_map(m1)
    @test m1(SVector(r1,r2)) ≈ SVector(ma(r1),mb(r2))
    m2 = productmap(m1,mb)
    test_generic_map(m2)
    @test m2(SVector(r1,r2,r3)) ≈ SVector(ma(r1),mb(r2),mb(r3))
    m3 = productmap(mb,m2)
    test_generic_map(m3)
    @test m3(SVector(r1,r2,r3,r4)) ≈ SVector(mb(r1),ma(r2),mb(r3),mb(r4))
    m4 = productmap(m1,m2)
    test_generic_map(m4)
    @test m4(SVector(r1,r2,r3,r4,r5)) ≈ SVector(m1(SVector(r1,r2))...,m2(SVector(r3,r4,r5))...)

    @test ProductMap(ma, mb) == ProductMap([ma,mb])
    @test hash(ProductMap(ma, mb)) == hash(ProductMap([ma,mb]))

    m5 = productmap(AffineMap(SMatrix{2,2,T}(1.0,2,3,4), SVector{2,T}(1,3)), LinearMap{T}(2.0))
    @test domaintype(component(m5,1)) == SVector{2,T}
    @test domaintype(component(m5,2)) == T
    test_generic_map(m5)
    x = SVector{3,T}(rand(T,3))
    @test m5(x) ≈ SVector(component(m5,1)(SVector(x[1],x[2]))...,component(m5,2)(x[3]))

    m6 = ProductMap([ma,mb])
    @test m6 isa DomainSets.VectorProductMap
    @test mapsize(m6) == (2,2)
    @test convert(Map{SVector{2,T}}, m6) isa DomainSets.VcatMap
    test_generic_map(m6)

    d1 = ProductDomain(T(1)..T(2), T(1)..T(2))
    d2 = ProductDomain(T(2)..T(4), T(2)..T(4))
    @test mapto(d1, d2) isa DomainSets.VcatMap
    @test jacobian(mapto(d1,d2)) isa ConstantMap
    @test jacdet(mapto(d1,d2)) == ConstantMap{SVector{2,T}}(4)
    d1v = ProductDomain([T(1)..T(2), T(1)..T(2)])
    d2v = ProductDomain([T(2)..T(4), T(2)..T(4)])
    @test mapto(d1v, d2v) isa DomainSets.VectorProductMap
    @test jacobian(mapto(d1v,d2v)) isa ConstantMap
    @test jacdet(mapto(d1v,d2v)) == ConstantMap{Vector{T}}(4)

    @test components(productmap(LinearMap(1), LinearMap(one(T)))) isa Tuple{<:Map{T},<:Map{T}}
end

using DomainSets: WrappedMap
function test_wrapped_maps(T)
    m1 = WrappedMap{T}(cos)
    m2 = WrappedMap{T}(sin)
    @test m1(one(T)) ≈ cos(one(T))
    @test m2(one(T)) ≈ sin(one(T))
    @test m1 == cos
    @test sin == m2
    m3 = m1 ∘ m2
    @test m3(one(T)) ≈ cos(sin(one(T)))

    @test WrappedMap(cos) isa WrappedMap{Float64}
    @test convert(Map, cos) isa WrappedMap
    @test convert(Map{BigFloat}, m1) isa WrappedMap{BigFloat}

    @test convert(Map{T}, cos) isa DomainSets.WrappedMap{T,typeof(cos)}

    @test convert(Map, LinearAlgebra.I) isa GenericLinearMap{Vector{Any}}
    @test convert(Map{Vector{T}}, LinearAlgebra.I) isa GenericLinearMap{Vector{T}}
end

function test_mixed_maps()
    m1 = composedmap(cos, sin)
    @test domaintype(m1) == Any
    @test m1(0.4) == sin(cos(0.4))

    m2 = composedmap(AffineMap(2.0, 3.0), cos)
    @test domaintype(m2) == Float64
    @test m2(0.4) ≈ cos(2*0.4+3)
    @inferred m2(0.4)
    @test repr(m2) == "cos ∘ (x -> 2.0 * x + 3.0)"

    m3 = composedmap(sin, AffineMap(2.0, 3.0), cos)
    @test m3(0.5) ≈ cos(2*sin(0.5)+3)
    @test repr(m3) == "cos ∘ (x -> 2.0 * x + 3.0) ∘ sin"
    @test domaintype(m3) == Any
    @inferred m3(0.5)

    m4 = productmap(sin, cos)
    @test m4 isa TupleProductMap
    @test domaintype(m4) == Tuple{Any,Any}
    @test m4(0.3,0.5) == (sin(0.3), cos(0.5))
end

@testset "maps" begin
    test_maps(Float64)
    test_maps(BigFloat)
    test_mixed_maps()
end


text:
#print
#

function printterm(tm)
 if isa(tm, Symbol)
   print(tm)
 elseif isa(tm, Number)
   print(tm)
 else # term
   print(tm.args[1])
   print_termlist(tm.args[2:end])
 end
end

function print_termlist(tl)
 if length(tl)==0; print("()")
 elseif length(tl)==1; print("(");printterm(tl[1]); print(")")
 else
  print("(")
  printterm(tl[1])
  for tm in tl[2:end]
   print(",")
   printterm(tm)
  end
  print(")")
 end
end

#

function printliteral(lit)
 print(lit)
end

function printlid(lid, core)
  print(" $lid.")
  printliteral(literalof(lid, core).body)
end

function printlids(lids, core)
 if isempty(lids)
  print("□")
 else
  for i in 1:length(lids)
    i!=1 && print(",")
#   println()
   printlid(lids[i], core)
  end
 end
end

function printbody(cls)
 if !isempty(cls)
  for i in 1:length(cls)
   println()
   print(" ")
   printliteral(cls[i])
  end
 end
end

function printvars(vars)
 print("[")
 for i in 1:length(vars)
  i != 1 && print(",")
  print(vars[i])
 end
 print("]")
end

function printclause(cid, core)
 print("$cid:")
 printvars(varsof(cid, core))

# if isempty(lidsof(cid,core)); print(".");
# else println(".") end
print(".")

 printlids(lidsof(cid,core), core)
 #printbody(cls.body)
end

printclauseof(lid, core) = printclause(cidof(lid, core),core)

function printclauses(core)
  for cid in sort(collect(keys(core.cdb)))
    printclause(cid, core)
println()
  end
end

function printcdb(cdb)
 if isempty(cdb)
  println("empty")
 else
  for cid in keys(cdb)
    println("$(cid): $(cdb[cid])")
  end
 end
end

function printldb(ldb)
 if isempty(ldb)
  println("empty")
 else
   for lid in keys(ldb)
     println("$(lid): $(ldb[lid].body)")
   end
 end
end

function printamap(amap)
 if isempty(keys(amap))
  println("empty")
 else
  for key in keys(amap)
    println("$key=>$(amap[key])")
  end
 end
end

function printstep(step)
  print("$(step.rid):<$(step.leftp):$(step.rightp)>=")
  printvars(step.sigma)
  print(":")
  printvars(step.rename[1])
  print("<-")
  printvars(step.rename[2])
end

function printproof(proof)
 for rid in sort(collect(keys(proof)))
   step = proof[rid]
   printstep(step)
   println()
 end
end

function printaproof0(rid, core)
 !(rid in keys(core.proof)) && return
 step=core.proof[rid]
 printaproof0(cidof(step.leftp,core), core)
 printaproof0(cidof(step.rightp,core), core)
 println("$rid=<$(cidof(step.leftp, core)):$(cidof(step.rightp, core))>")
end

function printns(shift)
  map(x->print(" "), 1:2shift)
end

function printaproof1(rid, core, shift=0)
  if rid in sort(collect(keys(core.proof)))
    step = core.proof[rid]
    printaproof1(cidof(step.leftp,core), core, shift+1)
    printaproof1(cidof(step.rightp,core), core, shift+1)

    println()
#    printns(shift)
    print("  ")
    print("<$(step.leftp):")
    print("$(step.rightp)>=")
    print_list(ovarsof(step.leftp, step.rightp, core))
    print("←")
    print_list(step.sigma)
    println()
#    printns(shift)
    printclause(rid, core)
  else
    println()
#    printns(shift)
    printclause(rid, core)
  end
end

function printproofs0(core)
 for rid in sort(contradictionsof(core))
  printaproof0(rid, core)
  println("\n---")
 end
end

function printproofs1(core)
 for rid in sort(contradictionsof(core))
  printaproof1(rid, core)
  println("\n---")
 end
end

function printmgu0(ovars, sigma, orig)
  for i in 1:length(ovars)
    if orig
      if ovars[i] != sigma[i]
        print("$(ovars[i])/$(sigma[i]); ")
      end
    else
      osig = deepcopy(sigma[i])
      if origof(ovars[i]) != origtermof(osig)
        print("$(origof(ovars[i]))/$(origtermof(osig));")
      end
    end
  end
end

function printmgu(rid, core, orig)
  if rid in keys(core.proof)
    step = core.proof[rid]
    printmgu(cidof(step.leftp,core), core, orig)
    printmgu(cidof(step.rightp,core), core, orig)

    println()
    print("$rid:")
    print("<$(step.leftp):")
    print("$(step.rightp)>=")
    printmgu0(ovarsof(step.leftp, step.rightp, core), step.sigma, orig)
  end
end

function printmgus(core, orig=false)
 for rid in contradictionsof(core)
  printmgu(rid, core, orig)
  println("\n---")
 end
end

function printproc(proc)
 for p in proc
  println("!$p")
 end
end

function printcore(core, showinfo=false)
println("name = $(core.name)")
println("max cid = $(core.maxcid)")
println("max rid = $(core.maxrid)")
showinfo && print_coreinfo(core)
#println("Psyms = $(core.allpsym)")
println("Proc")
printproc(core.proc)
println()
println("CDB")
 printcdb(core.cdb)
println()
println("LDB")
 printldb(core.ldb)
println()
println("LCMAP")
 printamap(core.lcmap)
println()
println("CLMAP")
 printamap(core.clmap)
println()
println("TEMPLATE(level0)")
 printtemplates0(core.level0,core)
println()
println("PROOFS")
 printproof(core.proof)
 println("\n-- end of core --")
println()
println("MEASURE")
 println("#trial   = $(core.trycnt[1])")
 println("#success = $(core.succnt[1])")
end

#### resolvent
function printliteral(lid, core)
 println("  $lid:$(literalof(lid,core).body)")
end

#### template

function printlid0(lid, core)
  print("$(lsymof(lid, core))")
end

function printlids0(lids, core)
  for lid in lids
   printlid0(lid, core)
  end
end

function printtemplate0(key, eq, core)
  print("$key= ")
  body=eq
  for i in 1:length(body)
    i!=1 && print("*")
    printlid0(body[i][3], core)
    print("[")
    printlids0(body[i][4],core)
    print("]($(body[i][1]))")
  end
end

function printtemplates0(eqs, core)
  for key in keys(eqs)
    eq = eqs[key]
    printtemplate0(key, eq, core)
    println()
  end
end

function printtemplate1(key, eq, core)
  print("$key= ")
  body=eq
  for i in 1:length(body)
    i!=1 && println("*")
    ll=body[i]
    printlid(ll[3], core)
    print("[")
    printlids(ll[4],core)
    print("]")
    print(" ($(ll[1]))")
  end
end

function printtemplates1(eqs,core)
  for key in keys(eqs)
    eq = eqs[key]
    printtemplate1(key, eq, core)
    println()
  end
end

##  print goal
function printgoal(goal, core)
  for lid in goal
    printliteral(lid, core)
  end
end

function printgoals(goals, core)
  for g in goals
    printgoal(g, core)
    println()
  end
end

## Analyzer
function print_list(set)
 print("[")
 for ei in 1:length(set)
  e = set[ei]
  if ei != 1; print(",") end
  print("$e")
 end
 print("]")
end

function print_set(set)
 print("[")
 for e in set
  print("$e ")
 end
 print("]")
end

function print_coreinfo(core)
 println("core info...")

 ci = analyze_sym(core)
 print("vars     = #$(length(ci.vsyms)):"); print_set(ci.vsyms);println()
 print("consts   = #$(length(ci.csyms)):"); print_set(ci.csyms);println()
 print("funcs    = #$(length(ci.fsyms)):"); print_set(ci.fsyms);println()
 print("preds    = #$(length(ci.psyms)):"); print_set(ci.psyms);println()
 print("clauses  = #$(length(keys(core.cdb))):"); print_list(collect(keys(core.cdb)));println()
 print("literals = #$(length(keys(core.ldb))):"); print_list(collect(keys(core.ldb)));println()

end

#### DYNAMIC Analysis

function printtracestep(vars, step)
  print("$(step.rid):<$(step.leftp):$(step.rightp)>=")
  printvars(step.sigma)
  print(" : ")
  printvars(step.rename[1])

  print("<-")
  printvars(step.rename[2])
end

function printprooftrace(rid, vars, proof)

 for rid in sort(collect(keys(proof)))
   step = proof[rid]
   printtracestep(vars, step)
   println()
 end
end

function printtrace(rid, vars, core)
 println("trace...$rid with $vars")
 step = core.proof[rid]
 printprooftrace(vars, core.proof)
end

function printatrace1(rid, vars, core)
  if rid in keys(core.proof)
    step = core.proof[rid]
    vars = printatrace1(cidof(step.leftp,core), vars, core)
    vars = printatrace1(cidof(step.rightp,core), vars, core)
    println();println("----")
    print("<$(step.leftp):")
    println("$(step.rightp)>")
    print("   unify:")
    print_termlist(ovarsof(step.leftp, step.rightp, core))
    print("←")
    print_termlist(step.sigma)
    println()
    vars = apply(ovarsof(step.leftp, step.rightp, core), vars, step.sigma)
    print("    VARS: "); print_termlist(vars); println()
    print("   rename:")
    print_termlist(step.rename[1])
    print("⇦")
    print_termlist(step.rename[2])
    println()
    vars = apply(step.rename[1], vars, step.rename[2])
    print("    VARS: "); print_termlist(vars);println()
    printclause(rid, core)
    println();println("----")
  else
    println()
    printclause(rid, core)
  end
  println()
  print("    VARS: ");print_termlist(vars)
  return vars
end

function traceofaproof(rid, vars, core)
  if rid in keys(core.proof)
    step = core.proof[rid]
    vars = traceofaproof(cidof(step.leftp,core), vars, core)
    vars = traceofaproof(cidof(step.rightp,core), vars, core)
    # print("<$(step.leftp):")
    # println("$(step.rightp)>")
    vars = apply(ovarsof(step.leftp, step.rightp, core), vars, step.sigma)
    vars = apply(step.rename[1], vars, step.rename[2])
  end
  return vars
end

function getBigvars(ci)
 bvars = []
 for x in ci.vsyms
  isuppercase(string(x)[1]) && push!(bvars, x)
 end
 bvars
end

function getVARS(core)
 ci = analyze_sym(core)
 getBigvars(ci)
end

function traceof(core)
  #vars = varsof(:C1, core)
  vars = getVARS(core)
  infos = []
  for rid in contradictionsof(core)
    inf = traceofaproof(rid, vars, core)
#    println("$vars = $inf")
    push!(infos, [rid, inf])
  end
  [vars, infos]
end

ext: jl

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# --- PRIVATE METHODS --------------------------------------------------------------------------------------- #
function _compute_profit_loss_at_expiration(asset::PSAbstractAsset, assetPriceValueArray::Array{Float64,1})::Array{Float64,1}

    # initialize -
    pl_value_array = Array{Float64,1}()

    # compute the intrinsic value -
    for (index, asset_price_value) in enumerate(assetPriceValueArray)
        result = intrinsic_value(asset,asset_price_value)
        if (isa(result.value,Exception) == true)
            return result
        end
        pl_value = result.value.pl_value

        # grab -
        push!(pl_value_array,pl_value)
    end

    # call -
    return pl_value_array
end

function _compute_profit_loss_at_expiration(asset::PSAbstractAsset, assetPriceStart::Float64, assetPriceStop::Float64;
    number_of_price_steps::Int64=1000)::Array{Float64,1}

    # initialize -
    profit_loss_array = Array{Float64,1}()

    # setup price range -
    asset_price_range = collect(range(assetPriceStart, assetPriceStop,length=number_of_price_steps))

    # return -
    return _compute_profit_loss_at_expiration(asset, asset_price_range)
end
# ----------------------------------------------------------------------------------------------------------- #

# --- PUBLIC METHODS ---------------------------------------------------------------------------------------- #
"""

    compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, assetPriceArray::Array{Float64,1})::(Union{PooksoftBase.PSResult{T}, Nothing} where T<:Any)

Compute the overall profit and loss (P/L) for a set option contracts with the same expiration date, and underlying asset.

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call contract models in this trade.
- `assetPriceArray::Array{Float64,1}`: An 1d array containing underlying asset prices to be used in the P/L calculation for this trade
"""
function compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, assetPriceArray::Array{Float64,1})::(Union{PSResult{T}, Nothing} where T<:Any)

    # TODO error checks -

    # get -
    number_of_price_steps = length(assetPriceArray)

    # initialize -
    number_of_assets = length(assetSet)
    assetProfitLossArray = zeros(number_of_price_steps, number_of_assets)
    tradeProfitLossArray = zeros(number_of_price_steps,2) # firstCol: price, secondCol: PL

    # process the set of assets -
    for (assetIndex, assetObject) in enumerate(assetSet)

        # process each asset type -
        profit_loss_array = _compute_profit_loss_at_expiration(assetObject, assetPriceArray)

        # add the PL for this asset to the overall array -
        for index = 1:number_of_price_steps
            assetProfitLossArray[index,assetIndex] = profit_loss_array[index]
        end
    end

    # ok, so the assetProfitLossArray holds the PL for each asset type in its cols -
    # to get the PL for the entire trade at a particular stock price, then sum across the cols -
    for price_index = 1:number_of_price_steps

        # grab col -
        tmp_row = assetProfitLossArray[price_index,:]

        # add -
        tradeProfitLossArray[price_index,1] = assetPriceArray[price_index]
        tradeProfitLossArray[price_index,2] = sum(tmp_row)
    end

    # return -
    return PSResult{Array{Float64,2}}(tradeProfitLossArray)
end

"""
    compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, assetPriceStart::Float64, assetPriceStop::Float64;
        number_of_price_steps::Int64=1000)::(Union{PooksoftBase.PSResult{T}, Nothing} where T<:Any)

Compute the overall profit and loss (P/L) for a set option contracts with the same expiration date, and underlying asset.

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call contract models in this trade.
- `assetPriceStart::Float64`: The start price for the underlying asset used to calculate the P/L values in this trade
- `assetPriceStop::Float64`: The start price for the underlying asset used to calculate the P/L values in this trade
- `number_of_price_steps::Int64=1000`: keyword arg describing the number of steps to take between the start and stop price. Default: 1000
"""
function compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, assetPriceStart::Float64, assetPriceStop::Float64;
    number_of_price_steps::Int64=1000)::(Union{PooksoftBase.PSResult{T}, Nothing} where T<:Any)

    # TODO: error checks -

    # setup price range -
    asset_price_array = collect(range(assetPriceStart, assetPriceStop,length=number_of_price_steps))

    # return -
    return compute_option_profit_and_loss_at_expiration(assetSet,asset_price_array)
end

"""
    compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, underlyingPriceRange::Tuple{Float64,Float64,Int64})::(Union{PooksoftBase.PSResult{T}, Nothing} where T<:Any)

Compute the overall profit and loss (P/L) for a set option contracts with the same expiration date, and underlying asset.

# Arguments
- `assetSet::Set{PSAbstractAsset}`: A set containing the put and call contract models in this trade.
- `underlyingPriceRange::Tuple{Float64,Float64,Int64}`: A tuple containing the price start, price stop and number of steps between the start and stop price to be used in the P/L calculation for this trade
"""
function compute_option_profit_and_loss_at_expiration(assetSet::Set{PSAbstractAsset}, underlyingPriceRange::Tuple{Float64,Float64,Int64})::(Union{PooksoftBase.PSResult{T}, Nothing} where T<:Any)

    # TODO: error checks -

    # setup price range -
    asset_price_array = collect(range(underlyingPriceRange[1], underlyingPriceRange[2],length=underlyingPriceRange[3]))

    # return -
    return compute_option_profit_and_loss_at_expiration(assetSet,asset_price_array)
end
# ----------------------------------------------------------------------------------------------------------- #


text:
function isone(a::AlgAssAbsOrdFracIdl)
  a = simplify!(a)
  return isone(denominator(a, false)) && isone(numerator(a, false))
end

###############################################################################
#
#  String I/O
#
###############################################################################

function show(io::IO, a::AlgAssAbsOrdFracIdl)
  print(io, "Fractional Ideal of ")
  show(IOContext(io, :compact => true), order(a))
  println(io, " with basis matrix ")
  print(io, basis_mat(a, false))
end

###############################################################################
#
#  Deepcopy
#
###############################################################################

function Base.deepcopy_internal(a::AlgAssAbsOrdFracIdl, dict::IdDict)
  b = typeof(a)(order(a), deepcopy(a.num), deepcopy(a.den))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :order && i != :num && i != :den
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

###############################################################################
#
#  Field access
#
###############################################################################

order(a::AlgAssAbsOrdFracIdl) = a.order

function numerator(a::AlgAssAbsOrdFracIdl, copy::Bool = true)
  if !isdefined(a, :num)
    if !isdefined(a, :basis_mat)
      error("No basis matrix defined")
    end
    a.num = ideal(order(a), numerator(basis_mat(a)))
  end
  if copy
    return deepcopy(a.num)
  else
    return a.num
  end
end

function denominator(a::AlgAssAbsOrdFracIdl, copy::Bool = true)
  if !isdefined(a, :den)
    if !isdefined(a, :basis_mat)
      error("No basis matrix defined")
    end
    a.den = denominator(basis_mat(a))
  end
  if copy
    return deepcopy(a.den)
  else
    return a.den
  end
end

function basis_mat(a::AlgAssAbsOrdFracIdl, copy::Bool = true)
  if !isdefined(a, :basis_mat)
    if !isdefined(a, :num) || !isdefined(a, :den)
      error("Numerator and/or denominator not defined")
    end
    a.basis_mat = FakeFmpqMat(basis_mat(a.num, copy = false), a.den)
  end
  if copy
    return deepcopy(a.basis_mat)
  else
    return a.basis_mat
  end
end

function basis_mat_inv(a::AlgAssAbsOrdFracIdl, copy::Bool = true)
  if !isdefined(a, :basis_mat_inv)
    a.basis_mat_inv = inv(basis_mat(a, false))
  end
  if copy
    return deepcopy(a.basis_mat_inv)
  else
    return a.basis_mat_inv
  end
end

################################################################################
#
#  Construction
#
################################################################################

function frac_ideal(O::AlgAssAbsOrd{S, T}, M::FakeFmpqMat, M_in_hnf::Bool = false) where {S, T}
  !M_in_hnf ? M = hnf(M, :lowerleft) : nothing

  return AlgAssAbsOrdFracIdl{S, T}(O, M)
end

function frac_ideal(O::AlgAssAbsOrd{S, T}, a::AlgAssAbsOrdIdl{S, T}, d::fmpz) where {S, T}
  return AlgAssAbsOrdFracIdl{S, T}(O, a, d)
end

function frac_ideal(O::AlgAssAbsOrd{S, T}, a::T) where {S, T}
  d = denominator(a, O)
  return frac_ideal(O, ideal(O, O(d*a)), d)
end

function frac_ideal_from_z_gens(O::AlgAssAbsOrd{S, T}, b::Vector{T}) where {S, T}
  d = degree(O)
  den = lcm([ denominator(bb) for bb in b ])
  num = ideal_from_z_gens(O, [ O(den*bb) for bb in b ])
  return frac_ideal(O, num, den)
end

################################################################################
#
#  Arithmetic
#
################################################################################

function +(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  d = lcm(denominator(a, false), denominator(b, false))
  ma = div(d, denominator(a, false))
  mb = div(d, denominator(b, false))
  return frac_ideal(order(a), numerator(a, false)*ma + numerator(b, false)*mb, d)
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  d = denominator(a, false)*denominator(b, false)
  return frac_ideal(order(a), numerator(a, false)*numerator(b, false), d)
end

^(A::AlgAssAbsOrdFracIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AlgAssAbsOrdFracIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

################################################################################
#
#  Ad hoc binary operations
#
################################################################################

function +(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  c = a*denominator(b, false) + numerator(b, false)
  return frac_ideal(order(a), c, denominator(b))
end

+(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T} = b + a

function *(a::AlgAssAbsOrdIdl{S, T}, b::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  return typeof(b)(order(a), a*numerator(b, false), denominator(b))
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, b::AlgAssAbsOrdIdl{S, T}) where {S, T}
  return typeof(a)(order(a), numerator(a, false)*b, denominator(a))
end

function *(x::T, a::AlgAssAbsOrdFracIdl{S, T}) where {S, T}
  O = order(a)
  return frac_ideal(O, x)*a
end

function *(a::AlgAssAbsOrdFracIdl{S, T}, x::T) where {S, T}
  O = order(a)
  return a*frac_ideal(O, x)
end

function *(x::T, a::AlgAssAbsOrdIdl{S, T}) where {S, T}
  O = order(a)
  return frac_ideal(O, x)*a
end

function *(a::AlgAssAbsOrdIdl{S, T}, x::T) where {S, T}
  O = order(a)
  return a*frac_ideal(O, x)
end

################################################################################
#
#  Equality
#
################################################################################

function ==(A::AlgAssAbsOrdFracIdl, B::AlgAssAbsOrdFracIdl)
  order(A) !== order(B) && return false
  return basis_mat(A, false) == basis_mat(B, false)
end

################################################################################
#
#  Simplification
#
################################################################################

function simplify!(a::AlgAssAbsOrdFracIdl)
  b = basis_mat(numerator(a, false), copy = false)
  g = gcd(denominator(a, false), content(b))

  if g != 1
    a.num = divexact(a.num, g)
    a.den = divexact(a.den, g)
  end

  return a
end

################################################################################
#
#  Random elements
#
################################################################################

function rand(a::AlgAssAbsOrdFracIdl, B::Int)
  z = rand(numerator(a, false), B)
  return fmpq(1, denominator(a, false))*elem_in_algebra(z, copy = false)
end


text:
using PowerSystems
const PSY = PowerSystems

############### Data Network ########################
include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusInverter.raw")
threebus_sys = System(threebus_file_dir, runchecks = false)
add_source_to_ref(threebus_sys)

############### Data devices ########################

function dyn_gen_second_order(generator)
    return PSY.DynamicGenerator(
        name = get_name(generator),
        ω_ref = 1.0, # ω_ref,
        machine = machine_oneDoneQ(), #machine
        shaft = shaft_no_damping(), #shaft
        avr = avr_type1(), #avr
        prime_mover = tg_none(), #tg
        pss = pss_none(), #pss
    )
end

function inv_case78(static_device)
    return DynamicInverter(
        name = get_name(static_device),
        ω_ref = 1.0, # ω_ref,
        converter = converter_high_power(), #converter
        outer_control = outer_control(), #outer control
        inner_control = inner_control(), #inner control voltage source
        dc_source = dc_source_lv(), #dc source
        freq_estimator = pll(), #pll
        filter = filt(), #filter
    )
end

for g in get_components(Generator, threebus_sys)
    if get_number(get_bus(g)) == 102
        case_gen = dyn_gen_second_order(g)
        add_component!(threebus_sys, case_gen, g)
    elseif get_number(get_bus(g)) == 103
        case_inv = inv_case78(g)
        add_component!(threebus_sys, case_inv, g)
    end
end

#Create Ybus_Fault
sys2 = System(threebus_file_dir, runchecks = false)
add_source_to_ref(sys2)
fault_branches = get_components(ACBranch, sys2)
for br in fault_branches
    if get_name(br) == "BUS 1-BUS 2-i_1"
        br.r = 3 * br.r
        br.x = 3 * br.x
        b_new = (from = br.b.from / 3, to = br.b.to / 3)
        br.b = b_new
    end
end
Ybus_fault = Ybus(sys2).data

for l in get_components(PSY.PowerLoad, threebus_sys)
    PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
end


text:
module QCMaterial

import PyCall: pyimport, PyNULL, PyVector
using Requires

# Refer to https://discourse.julialang.org/t/pyimport-works-from-repl-defined-module-but-not-in-my-package/43539.
const ofermion = PyNULL()
const ofpyscf = PyNULL()
const qulacs = PyNULL()
const pyutil = PyNULL()

function __init__()
    copy!(ofermion, pyimport("openfermion"))
    copy!(ofpyscf, pyimport("openfermionpyscf"))
    copy!(qulacs, pyimport("qulacs"))
    pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
    copy!(pyutil, pyimport("pyutil"))
    @require MPI="da04e1cc-30fd-572f-bb4f-1f8673147195" include("mpi.jl")
end

include("util.jl")
include("core.jl")
include("hamiltonian.jl")
include("uccsd.jl")
include("computation.jl")
include("hartree_fock.jl")
include("no_mpi.jl")
include("vqs.jl")

end


text:
@testset "cheb_diff" begin
    # Test the Chebyshev differentiation matrix implementation
    x2 = [0, 1//2, 1]
    D2 = [-3 4 -1; -1 0 1; 1 -4 3]
    x3 = [0, 1//4, 3//4, 1]
    D3 = [-38//6 8 -8//3 1; -2 2//3 2 -2//3; 2//3 -2 -2//3 2; -1 8//3 -8 38//6]
    # x4 = [1, 1/sqrt(2), 0, -1/sqrt(2), -1] # could work out symbolically if want to test BigFloat accuracy
    @test all(cheb_diff(Float64, 2) .≈ (x2, D2))
    @test all(cheb_diff(Float64, 3) .≈ (x3, D3))
    @test all(eltype.(cheb_diff(Float32, 3)) .=== (Float32, Float32))
end


text:
using Unicode
# urgh, doesn't seem good to use this global, but it works...
global all_paths = []

function read_map(input_file)
    # return a dictionary with nodes as keys, and list of connected nodes as values
    map_dict = Dict("start"=>[],
                    "end"=>[])
    for line in readlines(input_file)
        first,second = split(line,"-")
        if (!).(in(first, keys(map_dict)))
            map_dict[first] = []
        end
        push!(map_dict[first], second)
        # opposite way around
        if (!).(in(second, keys(map_dict)))
            map_dict[second] = []
        end
        push!(map_dict[second], first)
    end
    # by hand don't allow to go back from 'end' to anywhere else
    map_dict["end"] = []
    return map_dict
end

function is_cave_small(cave_name)
    return all(x->Unicode.islowercase(x),cave_name)
end

function can_i_visit_cave(cave, path)
    if (is_cave_small(cave) == false)
        return true
    elseif cave == "start"
        return false
    elseif (!).(in(cave, path))
        return true
    else
        # see if any other small caves have been visited twice
        for other_cave in path
            if all(x->Unicode.isuppercase(x),other_cave)
                continue
            end
            if count(x->x==other_cave, path) == 2
                return false
            end
        end
    end
    return true
end

function test_can_i_visit_cave()
    @assert can_i_visit_cave("A",["start","A","b"]) == true
    @assert can_i_visit_cave("b",["start","A","b"]) == true
    @assert can_i_visit_cave("start",["start","A","b"]) == false
    @assert can_i_visit_cave("c",["start","A","b"]) == true
    @assert can_i_visit_cave("c",["start","A","b","b"]) == true
    @assert can_i_visit_cave("c",["start","A","b","b","c"]) == false
end

test_can_i_visit_cave()

function make_path_branches(current_path, map_dict)
    current_cave = current_path[length(current_path)]
    if current_cave == "end"
        push!(all_paths,current_path)
        return
    end
    for possibility in map_dict[current_cave]
        path = deepcopy(current_path)
        if can_i_visit_cave(possibility, path)
            push!(path, possibility)
            make_path_branches(path, map_dict)
        end
    end
    return
end

test_map = read_map("input.txt")

make_path_branches(["start"], test_map)
println("length of all paths is ",length(all_paths))


text:
# Thermodynamic data of gas

"return the enthalpy of substance, '[KJ/mol]"
function GasEnthalpy(subName::String, Temp::Real)
    try
        Trange, coeff = gasThermoData[subName]
        if Temp > Trange[end] || Temp < Trange[1]
            printstyled(color=:red,
                @sprintf("%4.4f K is out of range, %4.2f K - %4.2f K \n",
                Temp, Trange[1], Trange[end]))
            return NaN
        else
            T = Temp*1E-3
            loc = findlast(x -> x<=Temp, Trange)
            return sum( [T, T^2/2, T^3/3, T^4/4, -1/T, 1, 0, -1] .* coeff[loc,:]) + coeff[loc, end]
        end
    catch
        printstyled(color=:red,@sprintf("%s not found in database.\n",subName))
        return NaN
    end
end

function GasEnthalpy()
    println(keys(gasThermoData))
end

" return the entropy of substance, [kJ/mol*K]"
function GasEntropy(subName::String, Temp::Real)
    try
        Trange, coeff = gasThermoData[subName]
        if Temp > Trange[end] || Temp < Trange[1]
            printstyled(color=:red,
                @sprintf("%4.4f K is out of range, %4.2f K - %4.2f K \n",
                Temp, Trange[1], Trange[end]))
            return NaN
        else
            T = Temp*1E-3
            loc = findlast(x -> x<=Temp, Trange)
            return sum( [log(T), T, T^2/2, T^3/3, -1/(2*T^2), 0, 1, 0] .* coeff[loc,:])*1E-3
        end
    catch
        printstyled(color=:red,@sprintf("%s not found in database.\n",subName))
        return NaN
    end
end

function GasEntropy()
    println(keys(gasThermoData))
end

"return the heat capacity of substance, kJ/mol*K"

function GasCp(subName::String, Temp::Real)
    try
        Trange, coeff = gasThermoData[subName]
        if Temp > Trange[end] || Temp < Trange[1]
            printstyled(color=:red,
                @sprintf("%4.4f K is out of range, %4.2f K - %4.2f K \n",
                Temp, Trange[1], Trange[end]))
            return NaN
        else
            T = Temp*1E-3
            loc = findlast(x -> x<=Temp, Trange)
            return sum( [1, T, T^2, T^3, 1/T^2, 0, 0, 0] .* coeff[loc,:])*1E-3
        end
    catch
        printstyled(color=:red,@sprintf("%s not found in database.\n",subName))
        return NaN
    end
end

function GasCp()
    println(keys(gasThermoData))
end


text:
using Optim
export solveDMD_withOptimBFGS

# BFGS solver
function solveDMD_withOptimBFGS(params;opts=Optim.Options())

    function fg!(f,gr,ar)
        if f != nothing && gr != nothing
            copyto!(params.ar,ar)
            return aBFuncGrad(gr,params)
        end
        if f == nothing && gr != nothing
            copyto!(params.ar,ar)
            aBGrad(gr,params)
            return
        end
        if f != nothing && gr == nothing
            copyto!(params.ar,ar)
            return aBFunc(params)
        end
        return
    end

    res = optimize(Optim.only_fg!(fg!),copy(params.ar),Optim.BFGS())
    return res

end


text:
using Kaleido
using Setfield
using Test

macro test_error(ex)
    quote
        let err = nothing
            try
                $(esc(ex))
            catch err
            end
            @test err isa Exception
            err
        end
    end
end


text:
using CallStack
using Test

const FILE = @__FILE__

foo() = nothing
bar() = foo()
baz() = bar()

buf = IOBuffer()
@callstack buf foo()
@test String(take!(buf)) ==
"""
   1  foo() in Main at $(FILE):6
"""

buf = IOBuffer()
@callstack buf bar()
@test String(take!(buf)) ==
"""
   1  bar() in Main at $(FILE):7
   2    foo() in Main at $(FILE):6
"""

buf = IOBuffer()
@callstack buf baz()
@test String(take!(buf)) ==
"""
   1  baz() in Main at $(FILE):8
   2    bar() in Main at $(FILE):7
   3      foo() in Main at $(FILE):6
"""


text:
include("common.jl")

example = "3,4,3,1,2"

struct FishyBusiness
    counts::Vector{Int64}
    FishyBusiness() = new(zeros(Int64, 9))
end

function Base.parse(::Type{FishyBusiness}, input::AbstractString)
    result = FishyBusiness()
    for value in split(input, ",")
        value = parse(Int64, value)
        i = value + 1
        result.counts[i] += 1
    end
    return result
end

function evolve!(fb::FishyBusiness)
    @inbounds begin
        num_new_lanternfish = first(fb.counts)

        # Bump up the zeros to 7.
        fb.counts[7 + 1] += fb.counts[0 + 1]
        fb.counts[0 + 1] = 0

        # Shift all other numbers down.
        for day in 1:8
            i = day + 1
            fb.counts[i - 1] = fb.counts[i]
        end

        # Now set the count at 8 days to be the number of new lanternfish.
        fb.counts[8 + 1] = num_new_lanternfish
    end
    return fb
end

num_fish(fb::FishyBusiness) = sum(fb.counts)

function num_fish(input::AbstractString, num_days::Integer)
    fb = parse(FishyBusiness, input)
    for _ in 1:num_days
        evolve!(fb)
    end
    return num_fish(fb)
end

get_ans(input::AbstractString) = num_fish(input, 80)

@assert get_ans(example) == 5934
@time get_ans(load_input(6))

# Part 2

get_ans2(input::AbstractString) = num_fish(input, 256)

@assert get_ans2(example) == 26984457539
get_ans2(load_input(6))


text:
# Is a macro useful for the helper functions in AG(), MAG(), etc?

module C
  macro makefunc(ex)
    return esc(quote
      wat() = $ex
    end)
  end
end

using .C

C.@makefunc(345)

wat() |> display


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiAppsV1RollingUpdateDeployment <: SwaggerModel
    maxSurge::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxSurge
    maxUnavailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxUnavailable

    function IoK8sApiAppsV1RollingUpdateDeployment(;maxSurge=nothing, maxUnavailable=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1RollingUpdateDeployment, Symbol("maxSurge"), maxSurge)
        setfield!(o, Symbol("maxSurge"), maxSurge)
        validate_property(IoK8sApiAppsV1RollingUpdateDeployment, Symbol("maxUnavailable"), maxUnavailable)
        setfield!(o, Symbol("maxUnavailable"), maxUnavailable)
        o
    end
end # type IoK8sApiAppsV1RollingUpdateDeployment

const _property_map_IoK8sApiAppsV1RollingUpdateDeployment = Dict{Symbol,Symbol}(Symbol("maxSurge")=>Symbol("maxSurge"), Symbol("maxUnavailable")=>Symbol("maxUnavailable"))
const _property_types_IoK8sApiAppsV1RollingUpdateDeployment = Dict{Symbol,String}(Symbol("maxSurge")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("maxUnavailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiAppsV1RollingUpdateDeployment }) = collect(keys(_property_map_IoK8sApiAppsV1RollingUpdateDeployment))
Swagger.property_type(::Type{ IoK8sApiAppsV1RollingUpdateDeployment }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiAppsV1RollingUpdateDeployment[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1RollingUpdateDeployment }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1RollingUpdateDeployment[property_name]

function check_required(o::IoK8sApiAppsV1RollingUpdateDeployment)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1RollingUpdateDeployment }, name::Symbol, val)
end


text:
const MatchPairs = Set{Union{Bool,Pair}}

# Catch all
match(x, y) = (x == y) === true ? MatchPairs() : MatchPairs(false)
match(x::Sym, y) = match(promote(x, y)...)
match(x, y::Sym) = match(promote(x, y)...)

match(x::Sym, y::Sym{<:Wild}) = MatchPairs(issubtag(tagof(x), tagof(y)) ? [y => x] : false)

function match(x::Sym, y::Sym)
    issubtag(tagof(x), tagof(y)) || return MatchPairs(false)
    gethead(x) === gethead(y) || return MatchPairs(false)
    xargs, yargs = getargs(x), getargs(y)
    length(xargs) == length(yargs) || return MatchPairs(false)
    m = mapreduce(x -> match(x[1], x[2]), union, zip(xargs, yargs))
    any(isfalse, m) && return MatchPairs(false)
    hashead(x, :symbol) && return MatchPairs([y => x])
    return m
end

function ismatch(pairs::MatchPairs)
    matches = Dict()
    for pair in pairs
        pair === true && continue
        pair === false && return false
        val = get(matches, pair[1], pair[2])
        isequal(val, pair[2]) || return false
        matches[pair[1]] = pair[2]
    end
    return true
end


text:
@testset "767.reorganize-string.jl" begin
    @test reorganize_string("vvvlo") == "vlvov" || reorganize_string("vvvlo") == "vovlv"
    @test reorganize_string("aab") == "aba"
    @test reorganize_string("aaba") == ""
end


text:
# Automatically generated test file
function exponential()

g = [
1.000000000000000,
1.284025416687741,
1.648721270700128,
2.117000016612675,
2.718281828459045,
3.490342957461841,
4.481689070338065,
5.754602676005730,
7.389056098930650,
9.487735836358526,
12.182493960703473,
15.642631884188171,
20.085536923187668,
25.790339917193062,
33.115451958692312,
42.521082000062783,
54.598150033144236]

ζ_vec = [50.9108, 15.0308, 4.43768, 1.31017, 0.386814, 0.114202, 0.0337169, 0.00995454, 0.00293897, 0.000867696]

I_sols = [
[1, 17],
[1, 12, 17],
[1, 9, 14, 17],
[1, 8, 12, 15, 17],
[1, 7, 11, 14, 16, 17],
[1, 6, 10, 13, 15, 16, 17],
[1, 6, 9, 11, 13, 15, 16, 17],
[1, 6, 9, 11, 13, 14, 15, 16, 17],
[1, 5, 8, 10, 12, 13, 14, 15, 16, 17],
[1, 5, 8, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 4, 7, 9, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 4, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 3, 5, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]]

f_sols = [
852.714249209732770,
50.910756467254942,
11.746481088717701,
4.139741013029379,
2.306336402589295,
1.321064892548748,
0.738536172471868,
0.369987216734444,
0.188919802507371,
0.103873014590135,
0.052932685417545,
0.022788524811403,
0.010314950639440,
0.003535039406415,
0.001084619914764,
0.000000000005457]

return g, ζ_vec, I_sols, f_sols

end


text:
## data layout

export DataLayout, dispose,
       byteorder, pointersize, intptr,
       sizeof, storage_size, abi_size,
       abi_alignment, frame_alignment, preferred_alignment,
       element_at, offsetof

# forward definition of DataLayout in src/module.jl
reftype(::Type{DataLayout}) = API.LLVMTargetDataRef

DataLayout(rep::String) = DataLayout(API.LLVMCreateTargetData(rep))

DataLayout(tm::TargetMachine) = DataLayout(API.LLVMCreateTargetDataLayout(ref(tm)))

function DataLayout(f::Core.Function, args...)
    data = DataLayout(args...)
    try
        f(data)
    finally
        dispose(data)
    end
end

dispose(data::DataLayout) = API.LLVMDisposeTargetData(ref(data))

Base.convert(::Type{String}, data::DataLayout) =
    unsafe_message(API.LLVMCopyStringRepOfTargetData(ref(data)))

function Base.show(io::IO, data::DataLayout)
    @printf(io, "DataLayout(%s)", convert(String, data))
end

byteorder(data::DataLayout) = API.LLVMByteOrder(ref(data))

pointersize(data::DataLayout) = API.LLVMPointerSize(ref(data))
pointersize(data::DataLayout, addrspace::Integer) =
    API.LLVMPointerSizeForAS(ref(data), addrspace)

intptr(data::DataLayout) = IntegerType(API.LLVMIntPtrType(ref(data)))
intptr(data::DataLayout, addrspace::Integer) =
    IntegerType(API.LLVMIntPtrTypeForAS(ref(data), addrspace))
intptr(data::DataLayout, ctx::Context) =
    IntegerType(API.LLVMIntPtrTypeInContext(ref(ctx), ref(data)))
intptr(data::DataLayout, addrspace::Integer, ctx::Context) =
    IntegerType(API.LLVMIntPtrTypeForASInContext(ref(ctx), ref(data), addrspace))

Base.sizeof(data::DataLayout, typ::LLVMType) =
    Int(API.LLVMSizeOfTypeInBits(ref(data), ref(typ)) / 8)
storage_size(data::DataLayout, typ::LLVMType) = API.LLVMStoreSizeOfType(ref(data), ref(typ))
abi_size(data::DataLayout, typ::LLVMType) = API.LLVMABISizeOfType(ref(data), ref(typ))

abi_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMABIAlignmentOfType(ref(data), ref(typ))
frame_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMCallFrameAlignmentOfType(ref(data), ref(typ))
preferred_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMPreferredAlignmentOfType(ref(data), ref(typ))
preferred_alignment(data::DataLayout, var::GlobalVariable) =
    API.LLVMPreferredAlignmentOfGlobal(ref(data), ref(var))

element_at(data::DataLayout, typ::StructType, offset::Integer) =
    API.LLVMElementAtOffset(ref(data), ref(typ), Culonglong(offset))

offsetof(data::DataLayout, typ::StructType, element::Integer) =
    API.LLVMOffsetOfElement(ref(data), ref(typ), element)


text:
# acc reachability script
# ENV["JULIA_DEBUG"] = Main # for debugging

using OVERTVerify
using LazySets
using Dates
using JLD2

controller = "acc_controller"
controller_filepath = "nnet_files/jmlr/"*controller*".nnet"
println("Controller is: ", controller)
query = OvertQuery(
    ACC,  # problem
    controller_filepath,    # network file
    Id(),      	# last layer activation layer Id()=linear, or ReLU()=relu
    "MIP",     	# query solver, "MIP" or "ReluPlex"
    55,        	# ntime
    0.1,       	# dt
    -1,        	# N_overt
    )

# x1,x2,x3 are lead vehicle variables
# x4,x5,x6 are ego vehicle variables
x_lead = [90,110]
v_lead = [32, 32.2]
gamma_lead = [0,0]
x_ego = [10,11]
v_ego = [30, 30.2]
gamma_ego = [0, 0]
var_list = [x_lead, v_lead, gamma_lead, x_ego, v_ego, gamma_ego]
input_set = Hyperrectangle(
    low=[v_range[1] for v_range in var_list],
    high=[v_range[2] for v_range in var_list]
    )

concretization_intervals = [20, 20, 15]
t1 = Dates.time()
concrete_state_sets, symbolic_state_sets, concrete_meas_sets, symbolic_meas_sets = symbolic_reachability_with_concretization(query, input_set, concretization_intervals)
t2 = Dates.time()
dt = (t2-t1)
print("elapsed time= $(dt) seconds")

# Intersect all sets with output constraint and see if
# reachable set is fully within safe set OR check to see if it ever intersects unsafe set
# they are equivalent
# We want the measurement to be greater than 10, always. So the unsafe set if <= 10
avoid_sets = [HalfSpace([1.], 10.)] # 1*y <= 10

reachable_meas_sets = clean_up_meas_sets(concrete_meas_sets, symbolic_meas_sets, concretization_intervals)

t1 = time()
safe, violations = check_avoid_set_intersection(reachable_meas_sets, input_set, avoid_sets)
dt_check = time() - t1

JLD2.@save "src/examples/jmlr/data/acc_reachability_data_55.jld2" query input_set concretization_intervals concrete_state_sets concrete_meas_sets symbolic_state_sets symbolic_meas_sets dt controller avoid_sets reachable_meas_sets safe violations dt_check


text:
"Calculate the Hamming Distance between two DNA strands"
function distance(a, b)
    length(a)==length(b) || throw(ArgumentError("diffrend strands length"))
    sum(collect(a).!=collect(b))
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# test/probmod.jl
# Testing problem modification
# Must be run as part of runtests.jl, as it needs a list of solvers.
#############################################################################
using JuMP, FactCheck

# If solvers not loaded, load them (i.e running just these tests)
!isdefined(:lp_solvers) && include("solvers.jl")

const TOL = 1e-4

facts("[probmod] Testing problem modification basics") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do

    # max 1.1x + 1.0y
    # st     x +    y <= 3
    #     0 <= x <= 3
    #     1 <= y <= 3
    # x* = 2, y* = 1
    m = Model(solver=solver)
    @variable(m, 0 <= x <= 3)
    @variable(m, 1 <= y <= 3)
    @objective(m, :Max, 1.1x + 1.0y)
    maincon = @constraint(m, x + y <= 3)
    @fact solve(m) --> :Optimal
    @fact m.internalModelLoaded --> true
    @fact getvalue(x) --> roughly(2.0, TOL)
    @fact getvalue(y) --> roughly(1.0, TOL)

    # Test adding a variable
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 3
    #     0 <= x <= 3
    #     1 <= y <= 3
    #     0 <= z <= 5
    # x* = 0, y* = 1, z* = 2
    @variable(m, 0 <= z <= 5, objective=100.0, inconstraints=[maincon], coefficients=[1.0])
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(0.0, TOL)
    @fact getvalue(y) --> roughly(1.0, TOL)
    @fact getvalue(z) --> roughly(2.0, TOL)
    @fact getobjectivevalue(m) --> roughly(201.0, TOL)

    # Test changing bounds
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 3
    #     0 <= x <= 3
    #     0 <= y <= 3
    #     0 <= z <= 2
    # x* = 1, y* = 0, z* = 2
    setlowerbound(y, 0.0)
    setupperbound(z, 2.0)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(1.0, TOL)
    @fact getvalue(y) --> roughly(0.0, TOL)
    @fact getvalue(z) --> roughly(2.0, TOL)
    m.internalModelLoaded = false
end
end
end

facts("[probmod] Testing problem modification part two") do
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 3
    #     0 <= x <= 3
    #     1 <= y <= 3
    #     0 <= z <= 0
    m = Model(solver=solver)
    @variable(m, 0 <= x <= 3)
    @variable(m, 0 <= y <= 3)
    @variable(m, 0 <= z <= 0)
    @objective(m, :Max, 1.1x + 1.0y + 100.0z)
    maincon = @constraint(m, x + y + z <= 3)
    @fact solve(m) --> :Optimal

    # Test changing problem type
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 3
    #     0 <= x <= 3
    #     0 <= y <= 3
    #     0 <= z <= 1.5, Integer
    # x* = 2, y* = 0, z* = 1
    setupperbound(z, 1.5)
    m.colCat[3] = :Int
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(2.0, TOL)
    @fact getvalue(y) --> roughly(0.0, TOL)
    @fact getvalue(z) --> roughly(1.0, TOL)

    # Test changing constraint bound (<= constraint)
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 2
    #     0 <= x <= 3
    #     0 <= y <= 3
    #     0 <= z <= 1.5, Integer
    # x* = 1, y* = 0, z* = 1
    JuMP.setRHS(maincon, 2.0)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(1.0, TOL)
    @fact getvalue(y) --> roughly(0.0, TOL)
    @fact getvalue(z) --> roughly(1.0, TOL)

    # Test adding a constraint
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 2
    #        x        +      z <= 0
    #     0 <= x <= 3
    #     0 <= y <= 3
    #     0 <= z <= 1.5, Integer
    # x* = 0, y* = 2, z* = 0
    xz0ref = @constraint(m, x + z <=0)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(0.0, TOL)
    @fact getvalue(y) --> roughly(2.0, TOL)
    @fact getvalue(z) --> roughly(0.0, TOL)

    # Test changing constraint bound (>= constraint)
    # max 1.1x + 1.0y + 100.0z
    # st     x +    y +      z <= 2
    #        x        +      z <= 1
    #        x +    y          >= 1
    #     0 <= x <= 3
    #     0 <= y <= 3
    #     0 <= z <= 1.5, Integer
    # x* = 1, y* = 0, z* = 1
    JuMP.setRHS(xz0ref, 2.0)
    xyg0ref = @constraint(m, x + y >= 0)
    @fact solve(m) --> :Optimal
    JuMP.setRHS(xyg0ref, 1)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(1.0, TOL)
    @fact getvalue(y) --> roughly(0.0, TOL)
    @fact getvalue(z) --> roughly(1.0, TOL)
end
end
end

facts("[probmod] Test adding a range constraint and modifying it") do
    m = Model()
    @variable(m, x)
    rangeref = @constraint(m, -10 <= x <= 10)
    @fact_throws JuMP.setRHS(rangeref, 11)
end

facts("[probmod] Test adding a 'decoupled' variable (#205)") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x >= 0)
    @objective(m, Min, x)
    solve(m)
    @variable(m, y >= 0)
    @constraint(m, x + y == 1)
    @objective(m, Min, 2x+y)
    solve(m)
    @fact getvalue(x) --> roughly(0.0, TOL)
    @fact getvalue(y) --> roughly(1.0, TOL)
end
end
end

facts("[probmod] Test JuMP.build") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x >= 0)
    @variable(m, y >= 0)
    @constraint(m, x + y == 1)
    @objective(m, Max, y)
    JuMP.build(m)
    @fact internalmodel(m) --> not(nothing)
    @fact m.internalModelLoaded --> true
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getvalue(x) --> roughly( 0.0, TOL)
    @fact getvalue(y) --> roughly( 1.0, TOL)
    @fact getobjectivevalue(m) --> roughly(1.0, TOL)
    @fact getdual(x)  --> roughly(-1.0, TOL)
    @fact getdual(y)  --> roughly( 0.0, TOL)
end
end
end

facts("[probmod] Test JuMP.build with MIP") do
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, x >= 0, Int)
    @variable(m, y, Bin)
    @constraint(m, x + y == 1)
    @objective(m, Max, y)
    JuMP.build(m)
    @fact internalmodel(m) --> not(nothing)
    @fact m.internalModelLoaded --> true
    stat = solve(m)
    @fact stat --> :Optimal
    @fact getvalue(x) --> roughly( 0.0, TOL)
    @fact getvalue(y) --> roughly( 1.0, TOL)
    @fact getobjectivevalue(m) --> roughly(1.0, TOL)
end
end
end

facts("[probmod] Test adding empty constraints") do
for solver in lp_solvers
context("With solver $(typeof(solver))") do
    m = Model(solver=solver)
    @variable(m, 0 <= x <= 9)
    @objective(m, Max, x)
    @constraint(m, x <= 5)
    @constraint(m, 0 <= 1)
    solve(m)
    @constraint(m, 0 <= 1)
    @fact solve(m) --> :Optimal
    @fact getvalue(x) --> roughly(5.0, TOL)
end
end
end

facts("[probmod] Test bound modification on binaries") do
for solver in ip_solvers
context("With solver $(typeof(solver))") do
    # Test semantics for modifying bounds on binary variables:
    # Variables should be restricted to the intersection of
    # {0,1} and their bounds.
    mod = Model(solver=solver)
    @variable(mod, x, Bin)
    @objective(mod, Max, x)
    solve(mod)
    @fact getvalue(x) --> roughly(1.0)
    setupperbound(x, 2.0)
    solve(mod)
    @fact getvalue(x) --> roughly(1.0)
    setupperbound(x, 0.0)
    solve(mod)
    @fact getvalue(x) --> roughly(0.0)
    # same thing, other direction
    mod = Model(solver=solver)
    @variable(mod, x, Bin)
    @objective(mod, Min, x)
    solve(mod)
    @fact getvalue(x) --> roughly(0.0)
    setlowerbound(x, -1.0)
    solve(mod)
    @fact getvalue(x) --> roughly(0.0)
    setlowerbound(x, 1.0)
    solve(mod)
    @fact getvalue(x) --> roughly(1.0)
end
end
end

facts("[probmod] Applicable regressions") do

function methods_test(solvername, solverobj, supp)
    mod = Model(solver=solverobj)
    @variable(mod, x >= 0)
    @constraint(mod, 2x == 2)
    solve(mod, suppress_warnings=true)
    internal_mod = internalmodel(mod)
    context(solvername) do
        for (it,(meth, args)) in enumerate(mpb_methods)
            if supp[it]
                @fact applicable(meth, internal_mod, args...) --> true
                @fact method_exists(meth, map(typeof, tuple(internal_mod, args...))) --> true
            end
        end
    end
end

const mpb_methods =
    [(MathProgBase.addquadconstr!, (Cint[1],Float64[1.0],Cint[1],Cint[1],Float64[1],'>',1.0)),
     (MathProgBase.setquadobjterms!, (Cint[1], Cint[1], Float64[1.0])),
     (MathProgBase.addconstr!,   ([1],[1.0],1.0,1.0)),
     (MathProgBase.addsos1!,     ([1],[1.0])),
     (MathProgBase.addsos2!,     ([1],[1.0])),
     (MathProgBase.addvar!,      ([1],[1.0],1.0,1.0,1.0)),
     (MathProgBase.setvarLB!,    ([1.0],)),
     (MathProgBase.setvarUB!,    ([1.0],)),
     (MathProgBase.setconstrLB!, ([1.0],)),
     (MathProgBase.setconstrUB!, ([1.0],)),
     (MathProgBase.setobj!,      ([1.0],)),
     (MathProgBase.setsense!,    (:Min,)),
     (MathProgBase.setvartype!,  ([:Cont],)),
     (MathProgBase.getinfeasibilityray, ()),
     (MathProgBase.getunboundedray, ()),
     (MathProgBase.getreducedcosts, ()),
     (MathProgBase.getconstrduals, ()),
     (MathProgBase.setwarmstart!, ([1.0]))]

if grb
    supp = (true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true)
    methods_test("Gurobi", Gurobi.GurobiSolver(OutputFlag=0), supp)
end
if cpx
    supp = (true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true)
    methods_test("CPLEX", CPLEX.CplexSolver(), supp)
end
if cbc
    supp = (false,false,true,false,false,true,true,true,true,true,true,true,true,true,true,true,true,false)
    methods_test("Clp", Clp.ClpSolver(), supp)
end
if cbc
    supp = (false,false,false,false,false,false,false,false,false,false,false,true,true,false,false,false,false,false)
    methods_test("Cbc", Cbc.CbcSolver(), supp)
end
if glp
    supp = (false,false,true,false,false,true,true,true,true,true,true,true,false,true,true,true,true,false)
    methods_test("GLPK LP", GLPKMathProgInterface.GLPKSolverLP(), supp)
    supp = (false,false,true,false,false, true,true,true,true,true,true,true,true,false,false,false,false,false)
    methods_test("GLPK MIP", GLPKMathProgInterface.GLPKSolverMIP(), supp)
end
if mos
    supp = (true,true,true,false,false,true,true,true,true,true,true,true,true,true,true,true,true,false)
    methods_test("Mosek", Mosek.MosekSolver(), supp)
end

end


text:
"""
Configurations for "watch" mode.
The configurations will only be active when used with [`report_and_watch_file`](@ref).

---
- `revise_all::Bool = true` \\
  Redirected to [`Revise.entr`](https://timholy.github.io/Revise.jl/stable/user_reference/#Revise.entr)'s `all` keyword argument.
  When set to `true`, JET will retrigger analysis as soon as code updates are detected in
    any module tracked by Revise.
  Currently when encountering `import/using` statements, JET won't perform analysis, but
    rather will just load the modules as usual execution (this also means Revise will track
    those modules).
  So if you're editing both files analyzed by JET and modules that are used within the files,
    this configuration should be enabled.
---
- `revise_modules = nothing` \\
  Redirected to [`Revise.entr`](https://timholy.github.io/Revise.jl/stable/user_reference/#Revise.entr)'s `modules` positional argument.
  If a iterator of `Module` is given, JET will retrigger analysis whenever code in `modules` updates.

  !!! tip
      This configuration is useful when your're also editing files that are not tracked by Revise,
      e.g. editing functions defined in `Base`:
      ```julia
      # re-performe analysis when you make a change to `Base`
      julia> report_and_watch_file(yourfile; revise_modules = [Base])
      ```
---
"""
struct WatchConfig
    # Revise configurations
    revise_all::Bool
    revise_modules
    @jetconfigurable WatchConfig(; revise_all::Bool = true,
                                   revise_modules   = nothing,
                                   ) =
        return new(revise_all,
                   revise_modules,
                   )
end

"""
    report_and_watch_file(filename::AbstractString;
                          jetconfigs...)

Watches `filename` and keeps re-triggering analysis with [`report_file`](@ref) on code update.
JET will try to analyze all the `include`d files reachable from `filename`, and it will
  re-trigger analysis if there is code update detected in any of the `include`d files.

This function internally uses [Revise.jl](https://timholy.github.io/Revise.jl/stable/) to
  track code updates. Revise also offers possibilities to track changes in files that are
  not directly analyzed by JET, or even changes in `Base` files. See [Watch Configurations](@ref)
  for more details.

See also: [`report_file`](@ref)
"""
function report_and_watch_file(args...; kwargs...)
    if @isdefined(Revise)
        _report_and_watch_file(args...; kwargs...)
    else
        init_revise!()
        @invokelatest _report_and_watch_file(args...; kwargs...)
    end
end

# HACK to avoid Revise loading overhead when just using `@report_call`, etc.
function init_revise!()
    @eval (@__MODULE__) using Revise
end

function _report_and_watch_file(filename::AbstractString; jetconfigs...)
    config = WatchConfig(; jetconfigs...)

    res = report_file(filename; jetconfigs...)
    display(res)
    included_files = res.res.included_files

    interrupted = false
    while !interrupted
        try
            Revise.entr(collect(included_files), config.revise_modules;
                        postpone = true, all = config.revise_all) do
                res = report_file(filename; jetconfigs...)
                display(res)
                next_included_files = res.res.included_files
                if any(∉(included_files), next_included_files)
                    # refresh watch files
                    throw(InsufficientWatches(next_included_files))
                end
                return nothing
            end
            interrupted = true # `InterruptException` was gracefully handled within `entr`, shutdown watch mode
        catch err
            # handle "expected" errors, keep running

            if isa(err, InsufficientWatches)
                included_files = err.included_files
                continue
            elseif isa(err, LoadError) ||
                   (isa(err, ErrorException) && startswith(err.msg, "lowering returned an error")) ||
                   isa(err, Revise.ReviseEvalException)
                continue

            # async errors
            elseif isa(err, CompositeException)
                errs = err.exceptions
                i = findfirst(e->isa(e, TaskFailedException), errs)
                if !isnothing(i)
                    tfe = errs[i]::TaskFailedException
                    res = tfe.task.result
                    if isa(res, InsufficientWatches)
                        included_files = res.included_files
                        continue
                    elseif isa(res, LoadError) ||
                           (isa(res, ErrorException) && startswith(res.msg, "lowering returned an error")) ||
                           isa(res, Revise.ReviseEvalException)
                        continue
                    end
                end
            end

            # fatal uncaught error happened in Revise.jl
            rethrow(err)
        end
    end
end

struct InsufficientWatches <: Exception
    included_files::Set{String}
end


text:
using WidgetsBase
using Test

@testset "Widgets.jl" begin
    s1 = Slider(1:100)
    @test s1[] == 1
    b = Button("hi")
    @test b[] == false
    t = TextField("Write!")
    @test t[] == "Write!"
    number_input = NumberInput(66.0)
    @test number_input[] == 66
    rslider = RangeSlider(1:100; value=[10, 80])
    @test rslider[] == [10, 80]

    checbkox = Checkbox(true)
    @test checbkox[] == true
end


text:
@inline function initialize!(integrator,cache::DP5ThreadedCache,f=integrator.f)
  integrator.kshortsize = 4
  integrator.k = [cache.update,cache.bspl,cache.dense_tmp3,cache.dense_tmp4]
  integrator.fsalfirst = cache.k1; integrator.fsallast = cache.k7
  f(integrator.t,integrator.uprev,integrator.fsalfirst) # Pre-start fsal
end

@inline function perform_step!(integrator,cache::DP5ThreadedCache,f=integrator.f)
  @unpack t,dt,uprev,u,k = integrator
  uidx = eachindex(integrator.uprev)
  @unpack a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,b1,b3,b4,b5,b6,b7,c1,c2,c3,c4,c5,c6 = cache.tab
  @unpack k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp = cache
  @unpack d1,d3,d4,d5,d6,d7 = cache.tab
  dp5threaded_loop1(dt,tmp,uprev,a21,k1,uidx)
  f(@muladd(t+c1*dt),tmp,k2)
  dp5threaded_loop2(dt,tmp,uprev,a31,k1,a32,k2,uidx)
  f(@muladd(t+c2*dt),tmp,k3)
  dp5threaded_loop3(dt,tmp,uprev,a41,k1,a42,k2,a43,k3,uidx)
  f(@muladd(t+c3*dt),tmp,k4)
  dp5threaded_loop4(dt,tmp,uprev,a51,k1,a52,k2,a53,k3,a54,k4,uidx)
  f(@muladd(t+c4*dt),tmp,k5)
  dp5threaded_loop5(dt,tmp,uprev,a61,k1,a62,k2,a63,k3,a64,k4,a65,k5,uidx)
  f(t+dt,tmp,k6)
  dp5threaded_loop6(dt,u,uprev,a71,k1,a73,k3,a74,k4,a75,k5,a76,k6,update,uidx)
  f(t+dt,u,integrator.fsallast)
  if integrator.opts.adaptive
    dp5threaded_adaptiveloop(dt,utilde,uprev,b1,k1,b3,k3,b4,k4,b5,k5,b6,k6,b7,k7,atmp,u,integrator.opts.abstol,integrator.opts.reltol,uidx)
    integrator.EEst = integrator.opts.internalnorm(atmp)
  end
  dp5threaded_denseloop(bspl,update,k1,k3,k4,k5,k6,k7,integrator.k,d1,d3,d4,d5,d6,d7,uidx)
  @pack integrator = t,dt,u,k
end

@noinline function dp5threaded_denseloop(bspl,update,k1,k3,k4,k5,k6,k7,k,d1,d3,d4,d5,d6,d7,uidx)
  Threads.@threads for i in uidx
    bspl[i] = k1[i] - update[i]
    k[3][i] = update[i] - k7[i] - bspl[i]
    k[4][i] = (d1*k1[i]+d3*k3[i]+d4*k4[i]+d5*k5[i]+d6*k6[i]+d7*k7[i])
  end
end

@noinline function dp5threaded_loop1(dt,tmp,uprev,a21,k1,uidx)
  a = dt*a21
  Threads.@threads for i in uidx
    tmp[i] = @muladd uprev[i]+a*k1[i]
  end
end

@noinline function dp5threaded_loop2(dt,tmp,uprev,a31,k1,a32,k2,uidx)
  Threads.@threads for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a31*k1[i]+a32*k2[i])
  end
end

@noinline function dp5threaded_loop3(dt,tmp,uprev,a41,k1,a42,k2,a43,k3,uidx)
  Threads.@threads for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a41*k1[i]+a42*k2[i]+a43*k3[i])
  end
end

@noinline function dp5threaded_loop4(dt,tmp,uprev,a51,k1,a52,k2,a53,k3,a54,k4,uidx)
  Threads.@threads for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i])
  end
end

@noinline function dp5threaded_loop5(dt,tmp,uprev,a61,k1,a62,k2,a63,k3,a64,k4,a65,k5,uidx)
  Threads.@threads for i in uidx
    tmp[i] = @muladd uprev[i]+dt*(a61*k1[i]+a62*k2[i]+a63*k3[i]+a64*k4[i]+a65*k5[i])
  end
end

@noinline function dp5threaded_loop6(dt,u,uprev,a71,k1,a73,k3,a74,k4,a75,k5,a76,k6,update,uidx)
  Threads.@threads for i in uidx
    update[i] = @muladd a71*k1[i]+a73*k3[i]+a74*k4[i]+a75*k5[i]+a76*k6[i]
    u[i] = @muladd uprev[i]+dt*update[i]
  end
end

@noinline function dp5threaded_adaptiveloop(dt,utilde,uprev,b1,k1,b3,k3,b4,k4,b5,k5,b6,k6,b7,k7,atmp,u,abstol,reltol,uidx)
  Threads.@threads for i in uidx
    utilde[i] = @muladd uprev[i] + dt*(b1*k1[i] + b3*k3[i] + b4*k4[i] + b5*k5[i] + b6*k6[i] + b7*k7[i])
    atmp[i] = ((utilde[i]-u[i])./@muladd(abstol+max(abs(uprev[i]),abs(u[i]))*reltol))
  end
end


text:
###############################################################################
# Stochiometry for a given reaction is a vector of pairs mapping species id to
# stochiometric coefficient.
###############################################################################

@inline @fastmath function evalrxrate(speciesvec::AbstractVector{T}, rxidx::S,
                              majump::MassActionJump{U,V,W})::R where {T,S,R,U <: AbstractVector{R},V,W}
    val = one(T)
    @inbounds for specstoch in majump.reactant_stoch[rxidx]
        specpop = speciesvec[specstoch[1]]
        val    *= specpop
        @inbounds for k = 2:specstoch[2]
            specpop -= one(specpop)
            val     *= specpop
        end
    end

    @inbounds return val * majump.scaled_rates[rxidx]
end

@inline @fastmath function executerx!(speciesvec::AbstractVector{T}, rxidx::S,
                                      majump::MassActionJump{U,V,W}) where {T,S,U,V,W}
    @inbounds net_stoch = majump.net_stoch[rxidx]
    @inbounds for specstoch in net_stoch
        speciesvec[specstoch[1]] += specstoch[2]
    end
    nothing
end

function scalerates!(unscaled_rates::AbstractVector{U}, stochmat::AbstractVector{V}) where {U,S,T,W <: Pair{S,T}, V <: AbstractVector{W}}
    @inbounds for i in eachindex(unscaled_rates)
        coef = one(T)
        @inbounds for specstoch in stochmat[i]
            coef *= factorial(specstoch[2])
        end
        unscaled_rates[i] /= coef
    end
    nothing
end

function scalerate(unscaled_rate::U, stochmat::AbstractVector{Pair{S,T}}) where {U <: Number, S, T}
    coef = one(T)
    @inbounds for specstoch in stochmat
        coef *= factorial(specstoch[2])
    end
    unscaled_rate /= coef
end

###############################################################################
# dependency graph when MassActionJump uses pairs to represent (species,stoich)
###############################################################################

# map from species to Set of reactions depending on that species
function spec_to_dep_rxs_map(numspec, ma_jumps::MassActionJump)

    numrxs = get_num_majumps(ma_jumps)

    # map from a species to reactions that depend on it
    spec_to_dep_rxs = [Set{Int}() for n = 1:numspec]
    for rx in 1:numrxs
        for (spec,stoch) in ma_jumps.reactant_stoch[rx]
            push!(spec_to_dep_rxs[spec], rx)
        end
    end

    spec_to_dep_rxs
end

# dependency graph is a map from a reaction to a vector of reactions
# that should depend on species it changes
function make_dependency_graph(numspec, ma_jumps::MassActionJump)

    numrxs          = get_num_majumps(ma_jumps)
    spec_to_dep_rxs = spec_to_dep_rxs_map(numspec, ma_jumps)

    # create map from rx to reactions depending on it
    dep_sets = [SortedSet{Int}() for n = 1:numrxs]
    for rx in 1:numrxs

        # rx changes spec, hence rxs depending on spec depend on rx
        for (spec,stoch) in ma_jumps.net_stoch[rx]
            for dependent_rx in spec_to_dep_rxs[spec]
                push!(dep_sets[rx], dependent_rx)
            end
        end
    end

    # convert to Vectors of Vectors
    dep_graph = Vector{Vector{Int}}(undef, numrxs)
    for rx = 1:numrxs
        dep_graph[rx] = [dep for dep in dep_sets[rx]]
    end

    dep_graph
end

# update dependency graph to make sure jumps depend on themselves
function add_self_dependencies!(dg)
    for (i,jump_deps) in enumerate(dg)
        if !any(y->isequal(y,i), jump_deps)
            push!(jump_deps, i)
            sort!(jump_deps)
        end
    end
end


text:
# Abstracts for algorithm parameters
mutable struct SPEA2 <: AbstractNSGA
    N::Int
    η_cr::Float64
    p_cr::Float64
    η_m::Float64
    p_m::Float64
    fitness::Vector{Float64}
end

"""
    SPEA2(;
        N = 100,
        η_cr = 20,
        p_cr = 0.9,
        η_m = 20,
        p_m = 1.0 / D,
        information = Information(),
        options = Options(),
    )

Parameters for the metaheuristic NSGA-II.

Parameters:

- `N` Population size.
- `η_cr`  η for the crossover.
- `p_cr` Crossover probability.
- `η_m`  η for the mutation operator.
- `p_m` Mutation probability (1/D for D-dimensional problem by default).

To use SPEA2, the output from the objective function should be a 3-touple
`(f::Vector, g::Vector, h::Vector)`, where `f` contains the objective functions,
`g` and `h` are inequality, equality constraints respectively.

A feasible solution is such that `g_i(x) ≤ 0 and h_j(x) = 0`.

```julia
using Metaheuristics

# Dimension
D = 2

# Objective function
f(x) = ( x, [sum(x.^2) - 1], [0.0] )

# bounds
bounds = [-1 -1;
           1  1.0
        ]

# define the parameters (use `SPEA2()` for using default parameters)
nsga2 = SPEA2(N = 100, p_cr = 0.85)

# optimize
status = optimize(f, bounds, nsga2)

# show results
display(status)
```

"""
function SPEA2(;
    N = 100,
    η_cr = 20,
    p_cr = 0.9,
    η_m = 20,
    p_m = -1,
    information = Information(),
    options = Options(),
)

    parameters = SPEA2(N, promote( Float64(η_cr), p_cr, η_m, p_m )...,[])
    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function update_state!(
    status::State,
    parameters::SPEA2,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )

    if isempty(parameters.fitness)
        parameters.fitness = compute_fitness(status.population)
    end

    Q = reproduction(status, parameters, problem)

    append!(status.population, create_solutions(Q, problem))

    # non-dominated sort, crowding distance, elitist removing
    environmental_selection!(status.population, parameters)
end

function environmental_selection!(population, parameters::SPEA2)
    Distance = compute_distances(population)
    fitness = compute_fitness(population, Distance)
    N = parameters.N

    next = fitness .< 1
    K = count(next)
    if K < N
        rank = sortperm(fitness)
        next[rank[1:N]] .= true
    elseif K > N
        del  = truncation(population[next], K-N, Distance[next,next])
        temp = findall(next)
        next[temp[del]] .= false
    end

    deleteat!(population, .!next)
    deleteat!(fitness, .!next)
    parameters.fitness = fitness

    return
end

function compute_distances(population)
    N = length(population)

    distances = zeros(N,N)
    for i in 1:N
        distances[i,i] = Inf
        for j in i+1:N
            distances[i,j] = norm(fval(population[i]) - fval(population[j]))
            distances[j,i] = distances[i,j]
        end
    end

    return distances
end

function truncation(population, K, distance = compute_distances(population))
    del = zeros(Bool, length(population))
    while count(del) < K
        remain   = findall(.!del)
        #temp     = sortperm(distance[remain,remain], dims = 2)
        # rank     = sortslicesperm(temp, dims=1)
        nn = argmin(distance[remain,remain])
        del[remain[nn.I[1]]] = true
    end

    return del
end

function compute_fitness(population, Distance = compute_distances(population))
    N = length(population)
    dominate = zeros(Bool,N, N)
    for i in 1:N
        for j in i+1:N
            k = compare(population[i], population[j])
            dominate[i,j] = k == 1
            dominate[j,i] = k == 2
        end
    end

    S = sum(dominate,dims=2)[:,1]
    R = [ sum(S[dominate[:,i]]) for i in 1:N]

    distance = sort(Distance,dims=2)

    D = 1 ./ (distance[:,floor(Int,sqrt(N))] .+ 2)
    return R + D
end