text:
using LazyArrays, FillArrays, LinearAlgebra, StaticArrays, ArrayLayouts, Test, Base64
import LazyArrays: MemoryLayout, DenseColumnMajor, materialize!, call, paddeddata,
                    MulAdd, Applied, ApplyLayout, DefaultApplyStyle, sub_materialize, resizedata!,
                    CachedVector, ApplyLayout, arguments

@testset "concat" begin
    @testset "Vcat" begin
        @testset "Vector" begin
            A = @inferred(Vcat(Vector(1:10), Vector(1:20)))
            @test eltype(A) == Int
            @test @inferred(axes(A)) == (Base.OneTo(30),)
            @test @inferred(A[5]) == A[15] == 5
            @test_throws BoundsError A[31]
            @test reverse(A) == Vcat(Vector(reverse(1:20)), Vector(reverse(1:10)))
            b = Array{Int}(undef, 31)
            @test_throws DimensionMismatch copyto!(b, A)
            b = Array{Int}(undef, 30)
            @test @allocated(copyto!(b, A)) == 0
            @test b == vcat(A.args...)
            @test copy(A) isa Vcat
            @test copy(A) == A
            @test copy(A) !== A
            @test vec(A) === A
            @test A' == transpose(A) == Vector(A)'
            @test permutedims(A) == permutedims(Vector(A))

            A = @inferred(Vcat(1:10, 1:20))
            @test @inferred(length(A)) == 30
            @test @inferred(A[5]) == A[15] == 5
            @test_throws BoundsError A[31]
            @test reverse(A) == Vcat(reverse(1:20), reverse(1:10))
            b = Array{Int}(undef, 31)
            @test_throws DimensionMismatch copyto!(b, A)
            b = Array{Int}(undef, 30)
            copyto!(b, A)
            @test_broken @allocated(copyto!(b, A)) == 0
            @test @allocated(copyto!(b, A)) ≤ 200
            @test b == vcat(A.args...)
            @test copy(A) === A
            @test vec(A) === A
            @test A' == transpose(A) == Vector(A)'
            @test A' == Hcat((1:10)', (1:20)')
            @test transpose(A) == Hcat(transpose(1:10), transpose(1:20))
            @test permutedims(A) == permutedims(Vector(A))

            @test map(copy,A) isa Vcat
            @test Applied(A)[3] == 3
        end

        @testset "Matrix" begin
            A = Vcat(randn(2,10), randn(4,10))
            @test @inferred(length(A)) == 60
            @test @inferred(size(A)) == (6,10)
            @test_throws BoundsError A[61]
            @test_throws BoundsError A[7,1]
            b = Array{Float64}(undef, 7,10)
            @test_throws DimensionMismatch copyto!(b, A)
            b = Array{Float64}(undef, 6,10)
            @test_broken @allocated(copyto!(b, A)) == 0
            @test @allocated(copyto!(b, A)) ≤ 200
            @test b == vcat(A.args...)
            @test copy(A) isa Vcat
            @test copy(A) == A
            @test copy(A) !== A
            @test vec(A) == vec(Matrix(A))
            @test A' == transpose(A) == Matrix(A)'
            @test permutedims(A) == permutedims(Matrix(A))
            @test_throws BoundsError A[7,2] = 6
            @test Applied(A)[1,3] == A[1,3]

            A = Vcat(randn(2,10).+im.*randn(2,10), randn(4,10).+im.*randn(4,10))
            @test eltype(A) == ComplexF64
            @test @inferred(length(A)) == 60
            @test @inferred(size(A)) == (6,10)
            @test_throws BoundsError A[61]
            @test_throws BoundsError A[7,1]
            b = Array{ComplexF64}(undef, 7,10)
            @test_throws DimensionMismatch copyto!(b, A)
            b = Array{ComplexF64}(undef, 6,10)
            @test_broken @allocated(copyto!(b, A)) == 0
            @test @allocated(copyto!(b, A)) ≤ 200
            @test b == vcat(A.args...)
            @test copy(A) isa Vcat
            @test copy(A) == A
            @test copy(A) !== A
            @test vec(A) == vec(Matrix(A))
            @test A' == Matrix(A)'
            @test transpose(A) == transpose(Matrix(A))
            @test permutedims(A) == permutedims(Matrix(A))

            @testset "indexing" begin
                A = Vcat(randn(2,10), randn(4,10))
                @test A[2,1:5] == Matrix(A)[2,1:5]
                @test A[2,:] == Matrix(A)[2,:]
                @test A[1:5,2] == Matrix(A)[1:5,2]
                @test A[:,2] == Matrix(A)[:,2]
                @test A[:,:] == A[1:6,:] == A[:,1:10] == A[1:6,1:10] == A
            end
        end

        @testset "etc" begin
            @test Vcat() isa Vcat{Any,1,Tuple{}}

            A = Vcat(1,zeros(3,1))
            @test A isa AbstractMatrix
            @test A[1,1] == 1.0
            @test A[2,1] == 0.0
            @test axes(A) == (Base.OneTo(4),Base.OneTo(1))
            @test permutedims(A) == permutedims(Matrix(A))
        end

        @testset "Vcat adjoints of vectors" begin
            # This special case was added to catch fast paths but
            # could be removed
            v = Vcat((1:5)', (2:6)')
            @test copyto!(Matrix{Float64}(undef,2,5), v) == Matrix(v) == [(1:5)'; (2:6)']
        end

        @testset "adjoint sub" begin
            @test arguments(view(Vcat(1,1:10)',1,:)) == (Fill(1,1),1:10)
        end

        @testset "Empty Vcat" begin
            @test @inferred(Vcat{Int}([1])) == [1]
            @test @inferred(Vcat{Int}()) == Int[]
        end

        @testset "in" begin
            @test 1 in Vcat(1, 1:10_000_000_000)
            @test 100_000_000 in Vcat(1, 1:10_000_000_000)
        end

        @testset "convert" begin
            for T in (Float32, Float64, ComplexF32, ComplexF64)
                Z = Vcat(zero(T),Zeros{T}(10))
                @test convert(AbstractArray,Z) ≡ Z
                @test convert(AbstractArray{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                @test convert(AbstractVector{T},Z) ≡ AbstractVector{T}(Z) ≡ Z
            end
        end

         @testset "Any/All" begin
            @test all(Vcat(true, Fill(true,100_000_000)))
            @test any(Vcat(false, Fill(true,100_000_000)))
            @test all(iseven, Vcat(2, Fill(4,100_000_000)))
            @test any(iseven, Vcat(2, Fill(1,100_000_000)))
            @test_throws TypeError all(Vcat(1))
            @test_throws TypeError any(Vcat(1))
        end

        @testset "isbitsunion #45" begin
            @test copyto!(Vector{Vector{Int}}(undef,6), Vcat([[1], [2], [3]], [[1], [2], [3]])) ==
                [[1], [2], [3], [1], [2], [3]]

            a = Vcat{Union{Float64,UInt8}}([1.0], [UInt8(1)])
            @test Base.isbitsunion(eltype(a))
            r = Vector{Union{Float64,UInt8}}(undef,2)
            @test copyto!(r, a) == a
            @test r == a
            @test copyto!(Vector{Float64}(undef,2), a) == [1.0,1.0]
        end

         @testset "maximum/minimum Vcat" begin
            x = Vcat(1:2, [1,1,1,1,1], 3)
            @test maximum(x) == 3
            @test minimum(x) == 1
        end

        @testset "searchsorted" begin
            a = Vcat(1:1_000_000, [10_000_000_000,12_000_000_000], 14_000_000_000)
            b = Vcat(1, 3:1_000_000, [2_000_000, 3_000_000])
            @test @inferred(searchsortedfirst(a, 6_000_000_001)) == 1_000_001
            @test @inferred(searchsortedlast(a, 2)) == 2
            @test @inferred(searchsortedfirst(b, 5)) == 4
            @test @inferred(searchsortedlast(b, 1)) == 1
        end

        @testset "args with hcat and view" begin
            A = Vcat(fill(2.0,1,10),ApplyArray(hcat, Zeros(1), fill(3.0,1,9)))
            @test arguments(view(A,:,10)) == ([2.0], [3.0])
        end

        @testset "union" begin
            a = Vcat([1,3,4],5:7)
            b = Vcat([1,3,4],5:7)
            union(a,b)
        end

        @testset "==" begin
            A = Vcat([1,2],[0])
            B = Vcat([1,2],[0])
            C = Vcat([1],[2,0])
            @test A == B == C == [1,2,0]
            @test A ≠ [1,2,4]
        end

        @testset "resizedata!" begin
            # allow emulating a cached Vector
            a = Vcat([1,2], Zeros(8))
            @test resizedata!(a, 2) ≡ a
            @test_throws BoundsError resizedata!(a,3)
        end

        @testset "Axpy" begin
            a = Vcat([1.,2],Zeros(1_000_000))
            b = Vcat([1.,2],Zeros(1_000_000))
            BLAS.axpy!(2.0, a, b)
            @test b[1:10] == [3; 6; zeros(8)]
            BLAS.axpy!(2.0, view(a,:), b)
            @test b[1:10] == [5; 10; zeros(8)]
        end

        @testset "l/rmul!" begin
            a = Vcat([1.,2],Zeros(1_000_000))
            @test ArrayLayouts.lmul!(2,a) ≡ a
            @test a[1:10] == [2; 4; zeros(8)]
            @test ArrayLayouts.rmul!(a,2) ≡ a
            @test a[1:10] == [4; 8; zeros(8)]
        end

        @testset "Dot" begin
            a = Vcat([1,2],Zeros(1_000_000))
            b = Vcat([1,2,3],Zeros(1_000_000))
            @test @inferred(dot(a,b)) ≡ 5.0
            @test @inferred(dot(a,1:1_000_002)) ≡ @inferred(dot(1:1_000_002,a)) ≡ 5.0
        end

        @testset "search" begin
            a = Vcat([1,2], 5:100)
            v = Vector(a)
            @test searchsortedfirst(a, 0) ≡ searchsortedfirst(v, 0) ≡ 1
            @test searchsortedfirst(a, 2) ≡ searchsortedfirst(v, 2) ≡ 2
            @test searchsortedfirst(a, 4) ≡ searchsortedfirst(v, 4) ≡ 3
            @test searchsortedfirst(a, 50) ≡ searchsortedfirst(v, 50) ≡ 48
            @test searchsortedfirst(a, 101) ≡ searchsortedfirst(v, 101) ≡ 99
            @test searchsortedlast(a, 0) ≡ searchsortedlast(v, 0) ≡ 0
            @test searchsortedlast(a, 2) ≡ searchsortedlast(v, 2) ≡ 2
            @test searchsortedlast(a, 4) ≡ searchsortedlast(v, 4) ≡ 2
            @test searchsortedlast(a, 50) ≡ searchsortedlast(v, 50) ≡ 48
            @test searchsortedlast(a, 101) ≡ searchsortedlast(v, 101) ≡ 98
            @test searchsorted(a, 0) ≡ searchsorted(v, 0) ≡ 1:0
            @test searchsorted(a, 2) ≡ searchsorted(v, 2) ≡ 2:2
            @test searchsorted(a, 4) ≡ searchsorted(v, 4) ≡ 3:2
            @test searchsorted(a, 50) ≡ searchsorted(v, 50) ≡ 48:48
            @test searchsorted(a, 101) ≡ searchsorted(v, 101) ≡ 99:98
        end

        @testset "print" begin
            @test Base.replace_in_print_matrix(Vcat(1:3,Zeros(10)), 4, 1, "0.0") == " ⋅ "
        end

        @testset "projection" begin
            a = Vcat(Ones(5), Zeros(5))
            b = randn(10)
            @test colsupport(a .* b) ≡ Base.OneTo(5)
            @test a .* b isa Vcat
            @test Diagonal(a) * b isa Vcat
            @test b .* a isa Vcat
            @test a .* b == Diagonal(a) * b == b .* a
        end
    end

    @testset "Hcat" begin
        A = @inferred(Hcat(1:10, 2:11))
        @test_throws BoundsError A[1,3]
        @test_throws BoundsError A[11,1]

        @test @inferred(call(A)) == hcat
        @test @inferred(size(A)) == (10,2)
        @test @inferred(A[5]) == @inferred(A[5,1]) == 5
        @test @inferred(A[11]) == @inferred(A[1,2]) == 2
        b = Array{Int}(undef, 11, 2)
        @test_throws DimensionMismatch copyto!(b, A)
        b = Array{Int}(undef, 10, 2)
        @test_broken @allocated(copyto!(b, A)) == 0
        @test @allocated(copyto!(b, A)) ≤ 200
        @test b == hcat(A.args...)
        @test copy(A) === A
        @test vec(A) == vec(Matrix(A))
        @test vec(A) === Vcat(1:10,2:11)
        @test A' == Matrix(A)'
        @test A' === Vcat((1:10)', (2:11)')

        A = Hcat(Vector(1:10), Vector(2:11))
        b = Array{Int}(undef, 10, 2)
        copyto!(b, A)
        @test b == hcat(A.args...)
        @test @allocated(copyto!(b, A)) == 0
        @test @allocated(copyto!(b, A)) ≤ 100
        @test copy(A) isa Hcat
        @test copy(A) == A
        @test copy(A) !== A
        @test vec(A) == vec(Matrix(A))
        @test vec(A) === Vcat(A.args...)
        @test A' == Matrix(A)'
        @test_throws BoundsError A[11,1] = 5
        @test_throws BoundsError A[5,3] = 5

        A = @inferred(Hcat(1, zeros(1,5)))
        @test A == hcat(1, zeros(1,5))
        @test vec(A) == vec(Matrix(A))
        @test A' == Matrix(A)'

        A = @inferred(Hcat(Vector(1:10), randn(10, 2)))
        b = Array{Float64}(undef, 10, 3)
        copyto!(b, A)
        @test b == hcat(A.args...)
        @test @allocated(copyto!(b, A)) == 0
        @test vec(A) == vec(Matrix(A))

        A = Hcat(randn(5).+im.*randn(5), randn(5,2).+im.*randn(5,2))
        b = Array{ComplexF64}(undef, 5, 3)
        copyto!(b, A)
        @test b == hcat(A.args...)
        @test @allocated(copyto!(b, A)) == 0
        @test vec(A) == vec(Matrix(A))
        @test A' == Matrix(A)'
        @test transpose(A) == transpose(Matrix(A))

        @testset "getindex bug" begin
            A = randn(3,3)
            H = Hcat(A,A)
            @test H[1,1] == applied(hcat,A,A)[1,1] == A[1,1]
        end

        @testset "adjoint vec / permutediims" begin
            @test vec(Hcat([1,2]', 3)) == 1:3
            @test permutedims(Hcat([1,2]', 3)) == reshape(1:3,3,1)
        end

        @testset "indexing" begin
            A = Hcat(randn(2,2), randn(2,3))
            @test A[2,1:5] == A[2,:] == Matrix(A)[2,:]
            @test A[1:2,2] == A[:,2] == Matrix(A)[:,2]
            @test A[:,2] == Matrix(A)[:,2]
            @test A[:,:] == A[1:2,:] == A[:,1:5] == A[1:2,1:5] == A
        end

        @testset "Hcat getindex" begin
            A = Hcat(1, (1:10)')
            @test A[1,:] isa Vcat{<:Any,1}
            @test A[1,:][1:10] == A[1,1:10]
        end
    end

    @testset "Hcat/Vcat adjoints" begin
        v = Vcat(1, 1:5)
        h = Hcat(1:5, 2:6)
        @test v' isa Adjoint
        @test transpose(v) isa Transpose
        @test MemoryLayout(v') isa DualLayout{ApplyLayout{typeof(hcat)}}
        @test MemoryLayout(transpose(v)) isa DualLayout{ApplyLayout{typeof(hcat)}}
        @test MemoryLayout(Adjoint(h)) isa ApplyLayout{typeof(vcat)}
        @test MemoryLayout(Transpose(h)) isa ApplyLayout{typeof(vcat)}
        @test copy(v') ≡ v'
        @test copy(Adjoint(h)) ≡ h'
        @test copy(transpose(v)) ≡ transpose(v)
        @test copy(Transpose(h)) ≡ transpose(h)

        @test arguments(v') ≡ (1, (1:5)')
        @test arguments(Adjoint(h)) ≡ ((1:5)', (2:6)')
        @test arguments(transpose(v)) ≡ (1, transpose(1:5))
        @test arguments(Transpose(h)) ≡ (transpose(1:5), transpose(2:6))

        @test v'*(1:6) ≡ 71

        A = Vcat(rand(2,3), randn(3,3))
        @test A' isa Hcat
        @test transpose(A) isa Hcat
        @test MemoryLayout(Adjoint(A)) isa ApplyLayout{typeof(hcat)}
        @test MemoryLayout(Transpose(A)) isa ApplyLayout{typeof(hcat)}
        @test Hcat(Adjoint(A)) == Hcat(Transpose(A)) == transpose(A) == A'
    end

    @testset "Hvcat" begin
        A = ApplyArray(hvcat, 2, randn(5,5), randn(5,6), randn(6,5), randn(6,6))
        @test eltype(A) == Float64
        @test A == hvcat(A.args...)
        @test_throws BoundsError A[2,12]
        @test_throws BoundsError A[12,2]
        @test A[[1,6],3] == Matrix(A)[[1,6],3]
        @test copyto!(similar(A), A) == A

        B = ApplyArray(hvcat, (3,2,1), randn(5,2), ones(5,3), randn(5,6), randn(6,5), randn(6,6), randn(2,11))
        @test eltype(B) == Float64
        @test B == hvcat(B.args...)
        @test_throws BoundsError A[2,14]
        @test_throws BoundsError A[14,2]
        @test copyto!(similar(B), B) == B

        A = ApplyArray(hvcat, 2, 1, 2, 3, 4)
        @test A == copyto!(similar(A), A) == [1 2; 3 4]
        @test_throws DimensionMismatch copyto!(zeros(Int,1,1), A)
        B = ApplyArray(hvcat, 3, 1, 2, 3, 4)
        @test_throws ArgumentError copyto!(similar(B), B)
        B = ApplyArray(hvcat, (1,2), 1, 2)
        @test_throws ArgumentError copyto!(zeros(Int,1,2), B)

        V = ApplyArray(hvcat, 2, [1,2], [3,4], [5,6], [7,8])
        @test V == [1 3; 2 4; 5 7; 6 8]

        W = ApplyArray(hvcat, 2, [1,2], [3], [5,6], [7,8])
        @test_throws ArgumentError copyto!(similar(W), W)
    end

    @testset "DefaultApplyStyle" begin
        v = Applied{DefaultApplyStyle}(vcat, (1, zeros(3)))
        @test v[1] == 1
        v = Applied{DefaultApplyStyle}(vcat, (1, zeros(3,1)))
        @test v[1,1] == 1
        H = Applied{DefaultApplyStyle}(hcat, (1, zeros(1,3)))
        @test H[1,1] == 1
    end

    @testset "setindex!" begin
        x = randn(5)
        y = randn(6)
        A = Vcat(x, y, 3)
        A[1] = 1
        @test A[1] == x[1] == 1
        A[6] = 2
        @test A[6] == y[1] == 2
        @test_throws MethodError A[12] = 3
        @test_throws BoundsError A[13] = 3

        x = randn(2,2); y = randn(3,2)
        A = Vcat(x,y)
        A[1,1] = 1
        @test A[1,1] == x[1,1] == 1
        A[3,1] = 2
        @test A[3,1] == y[1,1] == 2
        A[6] = 3
        @test A[1,2] == x[1,2] == 3

        x = randn(2,2); y = randn(2,3)
        B = Hcat(x,y)
        B[1,1] = 1
        @test B[1,1] == x[1,1] == 1
        B[1,3] = 2
        @test B[1,3] == y[1,1] == 2
    end

    @testset "fill!" begin
        A = Vcat([1,2,3],[4,5,6])
        fill!(A,2)
        @test A == fill(2,6)

        A = Vcat(2,[4,5,6])
        @test fill!(A,2) == fill(2,4)
        @test_throws ArgumentError fill!(A,3)

        A = Hcat([1,2,3],[4,5,6])
        fill!(A,2)
        @test A == fill(2,3,2)
    end

    @testset "Mul" begin
        A = Hcat([1.0 2.0],[3.0 4.0])
        B = Vcat([1.0,2.0],[3.0,4.0])

        @test MemoryLayout(typeof(A)) isa ApplyLayout{typeof(hcat)}
        @test MemoryLayout(typeof(B)) isa ApplyLayout{typeof(vcat)}
        @test A*B == Matrix(A)*Vector(B) == mul!(Vector{Float64}(undef,1),A,B) == (Vector{Float64}(undef,1) .= @~ A*B)
        @test materialize!(MulAdd(1.1,A,B,2.2,[5.0])) == 1.1*Matrix(A)*Vector(B)+2.2*[5.0]

        @test B * A ≈ Array(B) * Array(A) ≈ Vcat([1 2]', [3 4]') * A

        @test B * BroadcastArray(exp, [1 2]) ≈ B * exp.([1 2])

        A = Hcat([1.0 2.0; 3 4],[3.0 4.0; 5 6])
        B = Vcat([1.0,2.0],[3.0,4.0])
        @test MemoryLayout(typeof(A)) isa ApplyLayout{typeof(hcat)}
        @test MemoryLayout(typeof(B)) isa ApplyLayout{typeof(vcat)}
        @test A*B == Matrix(A)*Vector(B) == mul!(Vector{Float64}(undef,2),A,B) == (Vector{Float64}(undef,2) .= @~ A*B)
        @test materialize!(MulAdd(1.1,A,B,2.2,[5.0,6])) ≈ 1.1*Matrix(A)*Vector(B)+2.2*[5.0,6]

        A = Hcat([1.0 2.0; 3 4],[3.0 4.0; 5 6])
        B = Vcat([1.0 2.0; 3 4],[3.0 4.0; 5 6])
        @test MemoryLayout(typeof(A)) isa ApplyLayout{typeof(hcat)}
        @test MemoryLayout(typeof(B)) isa ApplyLayout{typeof(vcat)}
        @test A*B == Matrix(A)*Matrix(B) == mul!(Matrix{Float64}(undef,2,2),A,B) == (Matrix{Float64}(undef,2,2) .= @~ A*B)
        @test materialize!(MulAdd(1.1,A,B,2.2,[5.0 6; 7 8])) ≈ 1.1*Matrix(A)*Matrix(B)+2.2*[5.0 6; 7 8]
    end

    @testset "broadcast" begin
        x = Vcat(1:2, [1,1,1,1,1], 3)
        y = 1:8
        f = (x,y) -> cos(x*y)
        @test f.(x,y) isa Vcat
        @test @inferred(broadcast(f,x,y)) == f.(Vector(x), Vector(y))

        @test (x .+ y) isa Vcat
        @test (x .+ y).args[1] isa AbstractRange
        @test (x .+ y).args[end] isa Int

        z = Vcat(1:2, [1,1,1,1,1], 3)
        @test (x .+ z) isa BroadcastArray
        @test (x + z) isa BroadcastArray
        @test Vector( x .+ z) == Vector( x + z) == Vector(x) + Vector(z) == z .+ x

        @testset "Lazy mixed with Static treats as Lazy" begin
            s = SVector(1,2,3,4,5,6,7,8)
            @test f.(x , s) isa Vcat
            @test f.(x , s) == f.(Vector(x), Vector(s))
        end

        @testset "special cased" begin
            @test Vcat(1, Ones(5))  + Vcat(2, Fill(2.0,5)) ≡ Vcat(3, Fill(3.0,5))
            @test Vcat(SVector(1,2,3), Ones(5))  + Vcat(SVector(4,5,6), Fill(2.0,5)) ≡ Vcat(SVector(5,7,9), Fill(3.0,5))
            @test Vcat([1,2,3],Fill(1,7)) .* Zeros(10) ≡ Zeros(10) .* Vcat([1,2,3],Fill(1,7)) ≡ Zeros(10)
        end

        @testset "2-arg" begin
            a = Vcat([1,2], 1:3)
            b = Vcat([1], 1:4)
            @test a+b == b+a
        end

        H = Hcat(1, zeros(1,10))
        @test H/2 isa Hcat
        @test 2\H isa Hcat
        @test H./Ref(2) isa Hcat
        @test Ref(2).\H isa Hcat
        @test H/2  == H./Ref(2) == 2\H == Ref(2) .\ H == [1/2 zeros(1,10)]
    end

    @testset "norm" begin
        for a in (Vcat(1,2,Fill(5,3)), Hcat([1,2],randn(2,2)), Vcat(1,Float64[])),
            p in (-Inf, 0, 0.1, 1, 2, 3, Inf)
            @test norm(a,p) ≈ norm(Array(a),p)
        end
    end

    @testset "SubV/Hcat" begin
        A = Vcat(1,[2,3], Fill(5,10))
        V = view(A,3:5)
        @test MemoryLayout(V) isa ApplyLayout{typeof(vcat)}
        @inferred(arguments(V))
        @test arguments(V)[1] ≡ Fill(1,0)
        @test A[parentindices(V)...] == copy(V) == Array(A)[parentindices(V)...]

        A = Vcat((1:100)', Zeros(1,100),Fill(1,2,100))
        V = view(A,:,3:5)
        @test MemoryLayout(V) isa ApplyLayout{typeof(vcat)}
        @test A[parentindices(V)...] == copy(V) == Array(A)[parentindices(V)...]
        V = view(A,2:3,3:5)
        @test MemoryLayout(V) isa ApplyLayout{typeof(vcat)}
        @test A[parentindices(V)...] == copy(V) == Array(A)[parentindices(V)...]
    end

    @testset "col/rowsupport" begin
        H = Hcat(Diagonal([1,2,3]), Zeros(3,3), Diagonal([1,2,3]))
        V = Vcat(Diagonal([1,2,3]), Zeros(3,3), Diagonal([1,2,3]))
        @test colsupport(H,2) == rowsupport(V,2) == 2
        @test colsupport(H,4) == rowsupport(V,4) == 1:0
        @test colsupport(H,8) == rowsupport(V,8) == 2
        @test colsupport(H,10) == rowsupport(V,10)== 1:0
        @test rowsupport(H,1) == colsupport(V,1) == 1:7
        @test rowsupport(H,2) == colsupport(V,2) == 2:8
        @test colsupport(H,3:4) == rowsupport(V,3:4) == Base.OneTo(3)
        @test rowsupport(H,2:3) == colsupport(V,2:3) == 2:9
    end

    @testset "print" begin
        H = Hcat(Diagonal([1,2,3]), Zeros(3,3))
        V = Vcat(Diagonal([1,2,3]), Zeros(3,3))
        @test stringmime("text/plain", H) == "hcat(3×3 Diagonal{$Int, Vector{$Int}}, 3×3 Zeros{Float64}):\n 1.0   ⋅    ⋅    ⋅    ⋅    ⋅ \n  ⋅   2.0   ⋅    ⋅    ⋅    ⋅ \n  ⋅    ⋅   3.0   ⋅    ⋅    ⋅ "
        @test stringmime("text/plain", V) == "vcat(3×3 Diagonal{$Int, Vector{$Int}}, 3×3 Zeros{Float64}):\n 1.0   ⋅    ⋅ \n  ⋅   2.0   ⋅ \n  ⋅    ⋅   3.0\n  ⋅    ⋅    ⋅ \n  ⋅    ⋅    ⋅ \n  ⋅    ⋅    ⋅ "
        v = Vcat(1, Zeros(3))
        @test colsupport(v,1) == 1:1
        @test stringmime("text/plain", v) == "vcat($Int, 3-element Zeros{Float64}):\n 1.0\n  ⋅ \n  ⋅ \n  ⋅ "
        A = Vcat(Ones{Int}(1,3), Diagonal(1:3))
        @test stringmime("text/plain", A) == "vcat(1×3 Ones{$Int}, 3×3 Diagonal{$Int, UnitRange{$Int}}):\n 1  1  1\n 1  ⋅  ⋅\n ⋅  2  ⋅\n ⋅  ⋅  3"
    end
end


text:
import Base.@__doc__
import Base.Meta.isexpr
import Base.isidentifier

"""
    number_of_arguments(method::Function)

    number_of_arguments(method::Method)

Gives the number of arguments taken by `method`. If a generic function
is supplied, it must be associated with exactly one method.
"""
function number_of_arguments(method::Function)
    fmethods = methods(method).ms
    @assert length(fmethods) == 1

    number_of_arguments(first(fmethods))
end

number_of_arguments(method::Method) = length(method.sig.parameters) - 1

"""
    argument_names(method::Function)

Gives the original names of all arguments taken by `method`. Requires
that `method` is the only method belonging to its generic function.
"""
function argument_names(method::Function)
    n_args = number_of_arguments(method)

    argtypes = repeat([Any]; outer=n_args)
    lowered_code = first(code_lowered(method, argtypes))

    arguments = filter(isidentifier, lowered_code.slotnames[2:end])
    @assert length(arguments) == number_of_arguments(method)

    arguments
end

"""
    @partially_applicable f(x, args...; kwargs...) = ()

    @partially_applicable function g(x, args...; kwargs...) end

Pretty hacky. When applied to a function definition `f` of a function with
more than one argument, defines two methods on the function: one equivalent
to that defined by `f`, and one that takes one fewer argument than `f`,
amounting to all but the first argument, and partially applies them,
returning the resulting single-parameter

# Examples

```jldoctest
julia> @partially_applicable f(x, a, b) = a*x + b
f (generic function with 2 methods)

julia> f(28, 38, 273) == 28 |> f(38, 273)
true

julia> @partially_applicable g(x; kwargs...) = (println(kwargs); x*5)
g (generic function with 2 methods)

julia> g(2, sortaneat=:yes) == 2 |> g(sortaneat=:yes)
Any[(:sortaneat,:yes)]
Any[(:sortaneat,:yes)]
true
```
"""
macro partially_applicable(func)
    func_original = copy(func)

    is_parameter(x) = isa(x, Symbol) || isexpr(x, :(::))

    (var_index, var) =
        [(i, x) for (i, x) in enumerate(func.args[1].args) if is_parameter(x)][2]

    deleteat!(func.args[1].args, var_index)

    if isa(var, Expr)
        var = var.args[1]
    end

    func.args[2] = Expr(:(->), var, func.args[2])

    quote
        @__doc__ $(esc(func_original))
                 $(esc(func))
    end
end

ext: jl


text:
#=struct VmapFunc
    fun
    in_axes
    out_axes
    obj
end

_ndims(x::AbstractArray) = ndims(x)
_ndims(x::AbstractArray{<:AbstractArray}) = map(ndims, x)

_max_ndims(x::AbstractArray{<:Number}) = ndims(x)
_max_ndims(x::AbstractArray) = maximum(map(ndims, x))
_max_ndims(x::Union{Tuple,NamedTuple,Dict})

mapover(f, iselement, x) = iselement(x) ? f(x) : map(e -> mapover(f, iselement, e), x)
mapover(f, iselement, x::Union{Dict, NamedTuple}) = iselement(x) ? f(x) : map(e -> mapover(f, iselement, e), values(x))

mapreduceover(f, op, iselement, x) = iselement(x) ? f(x) : reduce(op, map(e -> mapreduceover(f, op, iselement, e), x))
mapreduceover(f, op, iselement, x::Union{Dict, NamedTuple}) = iselement(x) ? f(x) : reduce(op, map(e -> mapreduceover(f, op, iselement, e), values(x)))

_isleaf(x) = x isa Number || x isa AbstractArray{<:Number}

(f::Vmapfunc)(args...; kwargs...)=begin
    isnothing(f.obj) || f.obj(args...; kwargs...)
    in_axes = f.in_axes
    out_axes = f.out_axes

    if in_axes isa Colon
        py_in_axes = 0
    elseif in_axes isa Integer
        maxdim = mapreduceover(ndims, max, _isleaf, args)
        py_in_axes = maxdim - in_axes
    elseif in_axes isa Tuple

    end
    if f._in
end=#

"""
    vmap(f, in_axes, out_axes)

Vectorizing map. Creates a function which maps fun over argument axes.

Parameters:
    fun (Callable) – Function to be mapped over additional axes.
    in_axes – A nonnegative integer, None, or (nested) standard Python container
        (tuple/list/dict) thereof specifying which input array axes to map over.
        If each positional argument to fun is an array, then in_axes can be a nonnegative
        integer, a None, or a tuple of integers and Nones with length equal to the number
        of positional arguments to fun. An integer or None indicates which array
        axis to map over for all arguments (with None indicating not to map any axis),
        and a tuple indicates which axis to map for each corresponding positional argument.
        If the positional arguments to fun are container types, the corresponding element
        of in_axes can itself be a matching container, so that distinct array axes can be
        mapped for different container elements. in_axes must be a container tree prefix of
        the positional argument tuple passed to fun.
        At least one positional argument must have in_axes not None. The sizes of the
        mapped input axes for all mapped positional arguments must all be equal.

    out_axes – A nonnegative integer, None, or (nested) standard Python container
        (tuple/list/dict) thereof indicating where the mapped axis should appear in
        the output. All outputs with a mapped axis must have a non-None out_axes specification.
Return type:
        Callable

Returns:
    Batched/vectorized version of fun with arguments that correspond to those of fun,
    but with extra array axes at positions indicated by in_axes, and a return value
    that corresponds to that of fun, but with extra array axes at positions indicated by out_axes.
"""
#vmap(f; in_axes, out_axes=:) = vmap(f, in_axes, out_axes=out_axes)
vmap(f; in_axes_fromend, out_axes_fromend=0) = vmap_endorder(f, in_axes_fromend, out_axes=out_axes)

function vmap(f, in_axes; out_axes=:)
    throw(ErrorException("Please use the argument in_axes_fromend, and write the batched dimensions
    starting from the last (0), second last (1)... etc.
    This is required because arrays are stored transposed."))
end

function vmap_endorder(f, in_axes=0, out_axes=0)
    # Change from Julia's 1 based to Python's 0 based indexing
    return jax.vmap(f, in_axes=in_axes, out_axes=out_axes)
end


text:
#=
PTE:
- Julia version:
- Author: Dan
- Date: 2021-06-03
=#
import Base.show
export Element, getElement
"""
	Element
Represents an Element. [`Atom`](@ref) holds objects of `Element`.
"""
struct Element
    name_                   ::String
    symbol_                 ::String
    group_                  ::Int64
    period_                 ::Int64
    atomic_number_          ::Int64
    atomic_weight_          ::Float64
    atomic_radius_          ::Float64
    covalent_radius_        ::Float64
    van_der_waals_radius_   ::Float64
    is_metal_               ::Bool
    electronegativity_      ::Float64
	properties_				::Vector{Tuple{String,UInt8}}
    Element() = new()

	Element(name_                   ::String,
			symbol_                 ::String,
			group_                  ::Int64,
			period_                 ::Int64,
			atomic_number_          ::Int64,
			atomic_weight_          ::Float64,
			atomic_radius_          ::Float64,
			covalent_radius_        ::Float64,
			van_der_waals_radius_   ::Float64,
			is_metal_               ::Bool,
			electronegativity_      ::Float64) = begin

				 new(	name_,
						symbol_,
						group_,
						period_,
						atomic_number_,
						atomic_weight_,
						atomic_radius_,
						covalent_radius_,
						van_der_waals_radius_,
						is_metal_,
						electronegativity_,
						Vector{Tuple{String,UInt8}}()
				)
			end
end

Base.show(io::IO, el::Element) = print(io, el.symbol_)

getElement(at::AbstractAtom) = at.element_

"Periodic Table of Elements"
const elements_ = Element[
#			 name			symbol	group	period	number	weight		atomicradius	cov.radius	vdw.radius	metal	en
    Element("Actinium",		"Ac",	3,	 	7,	 	89,		227.0278,	1.88,			0.0,		0.0,		true,	1.3	),
	Element("Aluminum", 	"Al",   13,    	3,      13,     26.981539,  1.43,  			1.25, 		2.05, 		true,   1.61),
	Element("Americium",	"Am",   0,    	7,      95,     243.0614,   1.73,  			0.0,  		0.0,  		true,   1.3	),
	Element("Antimony", 	"Sb",   15,    	5,      51,     121.76,     1.82,  			1.41, 		2.2,  		false,  2.05),
	Element("Argon",    	"Ar",   18,    	3,      18,     39.948,     1.74,  			0.0,  		1.91, 		false,  0.0	),
	Element("Arsenic",  	"As",   15,   	4,      33,     74.92159,   1.25,  			1.21, 		2.0,  		false,  2.18),
	Element("Astatine", 	"At",   17,   	6,      85,     209.9871,   0.0,   			0.0,  		0.0,  		false,  1.96),
	Element("Barium",   	"Ba",   2,    	6,      56,     137.327,    2.17,  			1.98, 		0.0,  		true,   0.89),
	Element("Berkelium",	"Bk",   0,   	7,      97,     247.0703,   1.70,  			0.0,  		0.0,  		true,   1.3	),
	Element("Beryllium",	"Be",   2,   	2,       4,     9.012182,  	1.13,  			0.89, 		0.0,  		true,   1.57),
	Element("Bismuth",  	"Bi",   15,     6,      83,     208.98037,  1.55,  			1.52, 		2.4,  		true,   2.0	),
	Element("Bohrium",  	"Bh",   7,    	7,     107,     262.12,     0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Boron",    	"B",    13,    	2,       5,     10.811,     0.83,  			0.88, 		2.08, 		false,  2.04),
	Element("Bromine",  	"Br",   17,    	4,      35,     79.904,     0.0,   			1.14, 		1.95, 		false,  2.96),
	Element("Cadmium",  	"Cd",   12,    	5,      48,     112.411,    1.49,  			1.41, 		0.0,  		true,   1.69),
	Element("Caesium",  	"Cs",   1,    	6,      55,     132.90543,  2.654, 			2.35, 		2.62, 		true,   0.79),
	Element("Calcium",  	"Ca",   2,    	4,      20,     40.078,     1.97,  			1.74, 		0.0,  		true,   1.0	),
	Element("Californium",	"Cf",   0,    	7,      98,     251.0796,   1.69,  			0.0,  		0.0,  		true,   1.3	),
	Element("Carbon",  		"C",    14,    	2,       6,     12.011,     0.77,  			0.77, 		1.85, 		false,  2.55),
	Element("Cerium",  		"Ce",   0,    	6,      58,     140.115,    1.825, 			1.65, 		0.0,  		true,   1.12),
	Element("Chlorine",		"Cl",   17,    	3,      17,     35.4527,    0.0,   			0.99, 		1.81, 		false,  3.16),
	Element("Chromium",		"Cr",   6,    	4,      24,     51.9961,    1.25,  			0.0,  		0.0,  		true,   1.66),
	Element("Cobalt",  		"Co",   9,    	4,      27,     58.9332,    1.25,  			1.16, 		0.0,  		true,   1.88),
	Element("Copper",  		"Cu",   11,    	4,      29,     63.546,     1.28,  			1.17, 		0.0, 		true,    1.9),
	Element("Curium",  		"Cm",   0,    	7,      96,     247.0703,   1.74,  			0.0,  		0.0,  		true,   1.3	),
	Element("Dubnium", 		"Db",   4,    	7,     104,     261.11,     0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Dysprosium",	"Dy",   0,    	6,      66,     162.5,      1.77,  			1.59, 		0.0,  		true,   1.23),
	Element("Einsteinium",	"Es",   0,    	7,      99,     252.083,    2.03,  			0.0,  		0.0,  		true,   1.3	),
	Element("Erbium",   	"Er",   0,    	6,      68,     167.26,     1.76,  			1.57, 		0.0,  		true,   1.25),
	Element("Europium", 	"Eu",   0,    	6,      63,     151.965,    2.04,  			1.85, 		0.0,  		true,   1.2	),
	Element("Fermium",  	"Fm",   0,    	7,     100,     257.0951,   0.0,   			0.0,  		0.0,  		false,  1.3	),
	Element("Fluorine", 	"F",    17,    	2,       9,     18.9984032, 0.709, 			0.58, 		1.35, 		false,  3.98),
	Element("Francium", 	"Fr",   1,    	7,      87,     223.0197,   2.7,   			0.0,  		0.0,  		true,   0.7	),
	Element("Gadolinium",	"Gd",   0,    	6,      64,     157.25,     1.8,   			1.61, 		0.0,  		true,   0.94),
	Element("Gallium",  	"Ga",   13,    	4,      31,     69.723,     1.22,  			1.25, 		0.0,  		true,   1.81),
	Element("Germanium",	"Ge",   14,    	4,      32,     72.61,      1.23,  			1.22, 		0.0,  		false,  2.01),
	Element("Gold",     	"Au",   11,    	6,      79,     196.96654,  1.44,  			1.34, 		0.0,  		true,   2.0	),
	Element("Hafnium",  	"Hf",   4,    	6,      72,     178.49,     1.56,  			1.44, 		0.0,  		true,   1.5	),
	Element("Hahnium",  	"Hn",   8,    	7,     108,     0.0,       	0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Helium",   	"He",   18,    	1,       2,     4.002602,  	1.28,  			0.0,  		1.22, 		false,  0.0	),
	Element("Holmium",  	"Ho",   0,    	6,      67,     164.93032,  1.77,  			1.58, 		0.0,  		true,   1.24),
	Element("Hydrogen", 	"H",    1,    	1,       1,     1.00797,   	0.78,  			0.3,  		1.2,  		false,  2.2	),	#aka Hydrogenium
	Element("Indium",   	"In",   13,    	5,      49,     114.818,    1.63,  			1.5,  		0.0,  		true,   1.78),
	Element("Iodine",   	"I",    17,    	5,      53,     126.90447,  0.0,   			1.33, 		2.15, 		false,  2.66),	#aka Jod
	Element("Iridium",  	"Ir",   9,    	6,      77,     192.217,    1.36,  			1.26, 		0.0,  		true,   2.28),
	Element("Iron",     	"Fe",   8,    	4,      26,     55.845,     1.24,  			1.16, 		0.0,  		true,   1.83),	#aka Ferrum
	Element("Joliotium",	"Jl",  	5,    	7,     105,     262.114,    0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Krypton",  	"Kr",   18,    	4,      36,     83.80,      0.0,   			1.89, 		1.98, 		false,  0.0	),
	Element("Lanthanum",	"La",   3,    	6,      57,     138.9055,   1.88,  			1.69, 		0.0,  		true,   1.1	),
	Element("Lawrencium",	"Lr",  	3,    	7,     103,     262.11,     0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Lead",     	"Pb",   14,    	6,      82,     207.2,      1.75,  			1.54, 		0.0,  		true,   2.02),	#aka Plumbum
	Element("Lithium",  	"Li",   1,    	2,       3,     6.941,     	1.52,  			1.23, 		0.0,  		true,   0.98),
	Element("Lutetium", 	"Lu",   3,    	6,      71,     174.967,    1.72,  			1.56, 		0.0,  		true,   1.3	),
	Element("Magnesium",	"Mg",   2,    	3,      12,     24.30506,   1.6,   			1.36, 		0.0,  		true,   1.31),
	Element("Manganese",	"Mn",   7,    	4,      25,     54.93805,   1.24,  			1.77, 		0.0,  		true,   1.55),	#aka Mangan
	Element("Meitnerium",	"Mt",  	9,    	7,     109,     0.0,       	0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Mendelevium",	"Md",  	0,    	7,     101,     258.1,      0.0,   			0.0,  		0.0,  		true,   1.3	),
	Element("Mercury",  	"Hg",   12,    	6,      80,     200.59,     1.60,  			1.44, 		0.0,  		true,   1.8	),	#aka Hydrargyrum
	Element("Molybdenum",   "Mo",   6,    	5,      42,     95.94,      1.36,  			1.29, 		0.0,  		true,   2.16),
	Element("Neodymium",    "Nd",   0,    	6,      60,     144.24,     1.82,  			1.64, 		0.0,  		true,   1.14),
	Element("Neon",         "Ne",   18,    	2,      10,     20.1797,    0.0,   			0.0,  		1.6,  		false,  0.0	),
	Element("Neptunium",    "Np",  	0,    	7,      93,     237.0482,   1.5,   			0.0,  		0.0,  		true,   1.28),
	Element("Nickel",       "Ni",   10,    	4,      28,     58.6934,    1.25,  			1.15, 		0.0,  		true,   1.91),
	Element("Niobium",      "Nb",   5,    	5,      41,     92.90638,   1.43,  			1.34, 		0.0,  		true,   1.6	),
	Element("Nitrogen", 	"N",    15,    	2,       7,     14.00674,   0.71,  			0.7,  		1.54, 		false,  3.04),	#aka Nitrogenium
	Element("Nobelium",     "No",  	0,    	7,     102,     259.1009,   0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Osmium",       "Os",   8,    	6,      76,     190.23,     1.35,  			1.26, 		0.0,  		true,   2.2	),
	Element("Oxygen",     	"O",    16,    	2,       8,     15.9994,    0.6,   			0.66, 		1.4,  		false,  3.44),	#aka Oxygenium
	Element("Palladium",    "Pd",   10,    	5,      46,     106.42,     1.38,  			1.28, 		0.0,  		true,   2.2	),
	Element("Phosphorus",   "P",    15,    	3,      15,     30.973762,  1.15,  			1.10, 		1.9,  		false,  2.19),
	Element("Platinum",     "Pt",   10,    	6,      78,     195.08,     1.38,  			1.29, 		0.0,  		true,   2.54),
	Element("Plutonium",    "Pu",  	7,    	0,      94,     244.0642,   0.0,   			0.0,  		0.0,  		true,   1.3	),
	Element("Polonium",     "Po", 	16,    	6,      84,     208.9824,   1.67,  			1.53, 		0.0,  		false,  2.2	),
	Element("Potassium",   	"K",    1,    	4,      19,     39.0983,    2.27,  			2.03, 		2.31, 		true,   0.82),	#aka Kalium
	Element("Praseodymium", "Pr",   0,    	6,      59,     140.90765,  1.83,  			1.65, 		0.0,  		true,   1.13),
	Element("Promethium",   "Pm",  	0,    	6,      61,     144.9127,   1.81,  			0.0,  		0.0,  		true,   0.94),
	Element("Protactinium", "Pa",   0,    	7,      91,     231.03588,  1.61,  			0.0,  		0.0,  		true,   1.38),
	Element("Radium",       "Ra",	2,    	7,      88,     226.0254,   2.23,  			0.0,  		0.0,  		true,   0.89),
	Element("Radon",        "Rn",	18,    	6,      86,     222.0176,   0.0,   			0.0,  		0.0,  		false,  0.7	),
	Element("Rhenium",      "Re",   7,    	6,      75,     186.207,    1.37,  			1.28, 		0.0,  		true,   2.2	),
	Element("Rhodium",      "Rh",   9,    	5,      45,     102.9055,   1.34,  			1.25, 		0.0,  		true,   2.28),
	Element("Rubidium",     "Rb",   1,    	5,      37,     85.4678,    1.475, 			0.0,  		2.44, 		true,   0.82),
	Element("Ruthenium",    "Ru",   8,    	5,      44,     101.07,     1.34,  			1.24, 		0.0,  		true,   2.2	),
	Element("Rutherfordium","Rf",   6,    	7,     106,     263.118,    0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Samarium",     "Sm",   0,    	6,      62,     150.36,     1.8,   			1.66, 		0.0,  		true,   1.17),
	Element("Scandium",     "Sc",   3,    	4,      21,     44.95591,   1.61,  			1.44, 		0.0,  		true,   1.36),
	Element("Selenium",     "Se",   16,    	4,      34,     78.96,      2.15,  			1.17, 		2.0,  		true,   2.55),
	Element("Silicon",     	"Si",   14,    	3,      14,     28.0855,    1.17,  			1.17, 		2.0,  		false,  1.9	),	#aka Silicium
	Element("Silver",      	"Ag",   11,    	5,      47,     107.8682,   1.44,  			1.34, 		0.0,  		true,   1.93),	#aka Argentum
	Element("Sodium",      	"Na",   1,    	3,      11,     22.989768,  1.54,  			0.0,  		2.31, 		true,   0.93),	#aka Natrium
	Element("Strontium",    "Sr",   2,    	5,      38,     87.62,      2.15,  			1.92, 		0.0,  		true,   0.95),	#aka
	Element("Sulphur",      "S",    16,    	3,      16,     32.066,     1.04,  			1.04, 		1.85, 		false,  2.58),	#aka Sulfur
	Element("Tantalum",     "Ta",   5,    	6,      73,     180.9479,   1.43,  			1.34, 		0.0,  		true,   2.36),
	Element("Technetium",   "Tc",   7,    	5,      43,     98.9072,    1.36,  			0.0,  		0.0,  		true,   1.9	),
	Element("Tellurium",    "Te",   16,    	5,      52,     127.6,      1.43,  			1.37, 		2.2,  		false,  2.1	),
	Element("Terbium",      "Tb",   0,    	6,      65,     158.92534,  1.78,  			1.59, 		0.0,  		true,   1.22),
	Element("Thallium",     "Tl",  	13,    	6,      81,     204.3833,   1.7,   			1.55, 		0.0,  		true,   2.33),
	Element("Thorium",      "Th",   0,    	7,      90,     232.0381,   1.80,  			0.0,  		0.0,  		true,   0.0	),
	Element("Thulium",      "Tm",   0,    	6,      69,     168.93421,  1.75,  			1.56, 		0.0,  		true,   0.96),
	Element("Tin",          "Sn",  	14,    	5,      50,     118.71,     1.41,  			1.4,  		2.0,  		true,   1.96),	#aka Stannum
	Element("Titanium",     "Ti",   4,    	4,      22,     47.867,     1.45,  			1.32, 		0.0,  		true,   1.54),
	Element("Tungsten",     "W",    6,    	6,      74,     183.84,     1.37,  			1.3,  		0.0,  		true,   1.9	),	#aka Wolfram
	Element("Ununbium",     "Uub",  12,    	7,     112,     0.0,       	0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Ununnilium",   "Uun",  10,    	7,     110,     0.0,       	0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Ununninium",   "Uuu",  11,    	7,     111,     0.0,       	0.0,   			0.0,  		0.0,  		true,   0.0	),
	Element("Uranium",      "U",    0,    	7,      92,     238.0289,   1.54,  			0.0,  		0.0,  		true,   1.26),
	Element("Vanadium",     "V",    5,    	4,      23,     50.9415,    1.32,  			0.0,  		0.0,  		true,   1.63),
	Element("Xenon",        "Xe",   18,    	5,      54,     131.29,     2.18,  			2.09, 		2.16, 		false,  2.6	),
	Element("Ytterbium",    "Yb",   0,    	6,      70,     173.04,     1.94,  			1.7,  		0.0,  		true,   1.27),
	Element("Yttrium",      "Y",    3,    	5,      39,     88.90585,   1.81,  			1.62, 		0.0,  		true,   1.22),
	Element("Zinc",         "Zn",   12,    	4,      30,     65.39,      1.33,  			1.25, 		0.0,  		true,   1.65),	#aka Zincum
	Element("Zirconium",    "Zr",   4,    	5,      40,     91.224,     1.6,   			1.45, 		0.0,  		true,   1.3	)
]

const symbol_to_element_ = Dict{String,Element}(
	el.symbol_ => el for el in elements_
)
symbol_to_element_["Cb"] = symbol_to_element_["Nb"]	#Niobium has 2 names and 2 symbols

Element(symbol::String) = symbol_to_element_[capitalize(symbol)]


text:
module ModuleConstraints
using MultipleViewGeometry.ModuleTypes, MultipleViewGeometry.ModuleMathAliases, MultipleViewGeometry.ModuleOperators
using StaticArrays, LinearAlgebra
export satisfy, EpipolarConstraint, Constraint

include("constraints.jl")
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^3*(a+b*tan(c+d*x^2)), x) == :(1//4*a*x^4+1//4*im*b*x^4+-1//2*b*d^-1*x^2*log(1+exp(2*im*(c+d*x^2)))+1//4*im*b*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^2))))
@test integrate(x*(a+b*tan(c+d*x^2)), x) == :((1/2)*a*x^2+-1//2*b*d^-1*log(cos(c+d*x^2)))
@test integrate(x^3*(a+b*tan(c+d*x^2))^2, x) == :(-1//4*b^2*x^4+1//4*a^2*x^4+(1/2)*b^2*d^-2*log(cos(c+d*x^2))+(1/2)*im*a*b*x^4+(1/2)*b^2*d^-1*x^2*tan(c+d*x^2)+(1/2)*im*a*b*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^2)))+-1*a*b*d^-1*x^2*log(1+exp(2*im*(c+d*x^2))))
@test integrate(x*(a+b*tan(c+d*x^2))^2, x) == :((1/2)*x^2*(a^2+-1*b^2)+(1/2)*b^2*d^-1*tan(c+d*x^2)+-1*a*b*d^-1*log(cos(c+d*x^2)))
@test integrate(x^3*(a+b*tan(c+d*x^2))^-1, x) == :(x^4*(4a+4*im*b)^-1+b*d^-1*x^2*(2*a^2+2*b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^2)))+-1*im*b*d^-2*(4*a^2+4*b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^2))))
@test integrate(x*(a+b*tan(c+d*x^2))^-1, x) == :(a*x^2*(2*a^2+2*b^2)^-1+b*d^-1*(2*a^2+2*b^2)^-1*log(a*cos(c+d*x^2)+b*sin(c+d*x^2)))
@test integrate(x^3*(a+b*tan(c+d*x^2))^-2, x) == :(-1*x^4*(4*a^2+4*b^2)^-1+(1/2)*b*d^-2*(a^2+b^2)^-2*(b+2*a*d*x^2)*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^2)))+-1*b*d^-1*x^2*(a+b*tan(c+d*x^2))^-1*(2*a^2+2*b^2)^-1+1//8*a^-1*d^-2*(a+im*b)^-1*(b+2*a*d*x^2)^2*(a^2+b^2)^-1+-1//2*im*a*b*d^-2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^2))))
@test integrate(x*(a+b*tan(c+d*x^2))^-2, x) == :((1/2)*x^2*(a^2+b^2)^-2*(a^2+-1*b^2)+-1*b*d^-1*(a+b*tan(c+d*x^2))^-1*(2*a^2+2*b^2)^-1+a*b*d^-1*(a^2+b^2)^-2*log(a*cos(c+d*x^2)+b*sin(c+d*x^2)))
@test integrate(x^3*(a+b*tan(c+d*x^(1/2))), x) == :(1//4*a*x^4+1//4*im*b*x^4+-21*b*d^-3*x^5//2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-2*b*d^-1*x^7//2*log(1+exp(2*im*(c+d*x^(1/2))))+105*b*d^-5*x^3//2*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+-315//2*b*d^-7*x^(1/2)*Polylogarithms.polylog(7,-1*exp(2*im*(c+d*x^(1/2))))+-315//4*im*b*d^-8*Polylogarithms.polylog(8,-1*exp(2*im*(c+d*x^(1/2))))+7*im*b*d^-2*x^3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+-105//2*im*b*d^-4*x^2*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+315//2*im*b*x*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*tan(c+d*x^(1/2))), x) == :(1//3*a*x^3+1//3*im*b*x^3+-10*b*d^-3*x^3//2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-2*b*d^-1*x^5//2*log(1+exp(2*im*(c+d*x^(1/2))))+15*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+15//2*im*b*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^(1/2))))+-15*im*b*x*d^-4*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+5*im*b*d^-2*x^2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*tan(c+d*x^(1/2))), x) == :((1/2)*a*x^2+(1/2)*im*b*x^2+-3*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-2*b*d^-1*x^3//2*log(1+exp(2*im*(c+d*x^(1/2))))+-3//2*im*b*d^-4*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+3*im*b*x*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate(a+b*tan(c+d*x^(1/2)), x) == :(a*x+im*b*x+im*b*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+-2*b*d^-1*x^(1/2)*log(1+exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*tan(c+d*x^(1/2)))^2, x) == :(-1//3*b^2*x^3+1//3*a^2*x^3+-15*b^2*d^-6*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+-2*im*b^2*d^-1*x^5//2+2*b^2*d^-1*x^5//2*tan(c+d*x^(1/2))+10*b^2*d^-2*x^2*log(1+exp(2*im*(c+d*x^(1/2))))+30*x*b^2*d^-4*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+2//3*im*a*b*x^3+-20*im*b^2*d^-3*x^3//2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+-20*a*b*d^-3*x^3//2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-4*a*b*d^-1*x^5//2*log(1+exp(2*im*(c+d*x^(1/2))))+15*im*a*b*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^(1/2))))+30*im*b^2*d^-5*x^(1/2)*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+30*a*b*d^-5*x^(1/2)*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^(1/2))))+-30*im*a*b*x*d^-4*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+10*im*a*b*d^-2*x^2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*tan(c+d*x^(1/2)))^2, x) == :((1/2)*a^2*x^2+-1//2*b^2*x^2+3*b^2*d^-4*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+im*a*b*x^2+-2*im*b^2*d^-1*x^3//2+2*b^2*d^-1*x^3//2*tan(c+d*x^(1/2))+6*x*b^2*d^-2*log(1+exp(2*im*(c+d*x^(1/2))))+-6*im*b^2*d^-3*x^(1/2)*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2))))+-6*a*b*d^-3*x^(1/2)*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^(1/2))))+-4*a*b*d^-1*x^3//2*log(1+exp(2*im*(c+d*x^(1/2))))+-3*im*a*b*d^-4*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^(1/2))))+6*im*a*b*x*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate((a+b*tan(c+d*x^(1/2)))^2, x) == :(x*a^2+-1*x*b^2+2*b^2*d^-2*log(cos(c+d*x^(1/2)))+2*im*a*b*x+2*b^2*d^-1*x^(1/2)*tan(c+d*x^(1/2))+-4*a*b*d^-1*x^(1/2)*log(1+exp(2*im*(c+d*x^(1/2))))+2*im*a*b*d^-2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^3*(a+b*tan(c+d*x^(1/2)))^-1, x) == :(x^4*(4a+4*im*b)^-1+-105*b*d^-5*x^3//2*(a^2+b^2)^-1*Polylogarithms.polylog(5,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+2*b*d^-1*x^7//2*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+21*b*d^-3*x^5//2*(a^2+b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+315*im*b*d^-8*(4*a^2+4*b^2)^-1*Polylogarithms.polylog(8,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+315*b*d^-7*x^(1/2)*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(7,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-315*im*b*x*d^-6*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-7*im*b*d^-2*x^3*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+105*im*b*d^-4*x^2*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*tan(c+d*x^(1/2)))^-1, x) == :(x^3*(3a+3*im*b)^-1+-15*im*b*d^-6*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-15*b*d^-5*x^(1/2)*(a^2+b^2)^-1*Polylogarithms.polylog(5,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+2*b*d^-1*x^5//2*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+10*b*d^-3*x^3//2*(a^2+b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-5*im*b*d^-2*x^2*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+15*im*b*x*d^-4*(a^2+b^2)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*tan(c+d*x^(1/2)))^-1, x) == :(x^2*(2a+2*im*b)^-1+2*b*d^-1*x^3//2*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+3*im*b*d^-4*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+3*b*d^-3*x^(1/2)*(a^2+b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-3*im*b*x*d^-2*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2)))))
@test integrate((a+b*tan(c+d*x^(1/2)))^-1, x) == :(x*(a+im*b)^-1+-1*im*b*d^-2*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+2*b*d^-1*x^(1/2)*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*tan(c+d*x^(1/2)))^-2, x) == :(1//3*x^3*(a+-1*im*b)^-2+-4//3*b^2*x^3*(a^2+b^2)^-2+-15*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-15*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+4*b*x^3*(a+-1*im*b)^-2*(-3b+3*im*a)^-1+-10*b^2*d^-2*x^2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-4*im*b^2*d^-1*x^5//2*(a^2+b^2)^-2+10*b^2*d^-2*x^2*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+15*b*d^-6*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+30*x*b^2*d^-4*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+30*x*b^2*d^-4*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-30*b*x*d^-4*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-30*b*d^-5*x^(1/2)*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-20*im*b^2*d^-3*x^3//2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-20*im*b^2*d^-3*x^3//2*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-4*im*b^2*d^-1*x^5//2*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+4*b*d^-1*x^5//2*(a+im*b)^-1*(a+-1*im*b)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+4*b^2*d^-1*x^5//2*(a+im*b)^-1*(b+im*a)^-2*(-1b+im*a+(b+im*a)*exp(2*im*(c+d*x^(1/2))))^-1+10*b*d^-2*x^2*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+20*b*d^-3*x^3//2*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+30*im*b^2*d^-5*x^(1/2)*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+30*im*b^2*d^-5*x^(1/2)*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x*(a+b*tan(c+d*x^(1/2)))^-2, x) == :((1/2)*x^2*(a+-1*im*b)^-2+-2*b^2*x^2*(a^2+b^2)^-2+2*b*x^2*(a+-1*im*b)^-2*(-1b+im*a)^-1+3*b^2*d^-4*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+3*b^2*d^-4*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-6*x*b^2*d^-2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-4*im*b^2*d^-1*x^3//2*(a^2+b^2)^-2+-3*b*d^-4*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+6*x*b^2*d^-2*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-6*im*b^2*d^-3*x^(1/2)*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-6*im*b^2*d^-3*x^(1/2)*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+-4*im*b^2*d^-1*x^3//2*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+4*b*d^-1*x^3//2*(a+im*b)^-1*(a+-1*im*b)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+4*b^2*d^-1*x^3//2*(a+im*b)^-1*(b+im*a)^-2*(-1b+im*a+(b+im*a)*exp(2*im*(c+d*x^(1/2))))^-1+6*b*x*d^-2*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2))))+6*b*d^-3*x^(1/2)*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^(1/2)))))
@test integrate((a+b*tan(c+d*x^(1/2)))^-2, x) == :(-1*x*(a^2+b^2)^-1+(1/2)*a^-1*d^-2*(a+im*b)^-1*(b+2*a*d*x^(1/2))^2*(a^2+b^2)^-1+-2*b*d^-1*x^(1/2)*(a+b*tan(c+d*x^(1/2)))^-1*(a^2+b^2)^-1+2*b*d^-2*(a^2+b^2)^-2*(b+2*a*d*x^(1/2))*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2))))+-2*im*a*b*d^-2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^(1/2)))))
@test integrate(x^2*(a+b*tan(c+d*x^1//3)), x) == :(1//3*a*x^3+1//3*im*b*x^3+945//2*b*d^-9*Polylogarithms.polylog(9,-1*exp(2*im*(c+d*x^1//3)))+-945*b*d^-7*x^2//3*Polylogarithms.polylog(7,-1*exp(2*im*(c+d*x^1//3)))+-42*b*d^-3*x^2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-3*b*d^-1*x^8//3*log(1+exp(2*im*(c+d*x^1//3)))+315*b*d^-5*x^4//3*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+-945*im*b*d^-8*x^1//3*Polylogarithms.polylog(8,-1*exp(2*im*(c+d*x^1//3)))+-126*im*b*d^-4*x^5//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+12*im*b*d^-2*x^7//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3)))+630*im*b*x*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^1//3))))
@test integrate(x*(a+b*tan(c+d*x^1//3)), x) == :((1/2)*a*x^2+(1/2)*im*b*x^2+-15*b*x*d^-3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-3*b*d^-1*x^5//3*log(1+exp(2*im*(c+d*x^1//3)))+45//2*b*d^-5*x^1//3*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+45//4*im*b*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^1//3)))+-45//2*im*b*d^-4*x^2//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+15//2*im*b*d^-2*x^4//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3))))
@test integrate(a+b*tan(c+d*x^1//3), x) == :(a*x+im*b*x+-3//2*b*d^-3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-3*b*d^-1*x^2//3*log(1+exp(2*im*(c+d*x^1//3)))+3*im*b*d^-2*x^1//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3))))
@test integrate(x^2*(a+b*tan(c+d*x^1//3))^2, x) == :(-1//3*b^2*x^3+1//3*a^2*x^3+-1260*x*b^2*d^-6*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+-3*im*b^2*d^-1*x^8//3+3*b^2*d^-1*x^8//3*tan(c+d*x^1//3)+24*b^2*d^-2*x^7//3*log(1+exp(2*im*(c+d*x^1//3)))+252*b^2*d^-4*x^5//3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+945*im*b^2*d^-9*Polylogarithms.polylog(8,-1*exp(2*im*(c+d*x^1//3)))+945*a*b*d^-9*Polylogarithms.polylog(9,-1*exp(2*im*(c+d*x^1//3)))+1890*b^2*d^-8*x^1//3*Polylogarithms.polylog(7,-1*exp(2*im*(c+d*x^1//3)))+2//3*im*a*b*x^3+-1890*im*b^2*d^-7*x^2//3*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^1//3)))+-1890*a*b*d^-7*x^2//3*Polylogarithms.polylog(7,-1*exp(2*im*(c+d*x^1//3)))+-84*im*b^2*d^-3*x^2*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3)))+-84*a*b*d^-3*x^2*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-6*a*b*d^-1*x^8//3*log(1+exp(2*im*(c+d*x^1//3)))+630*im*b^2*d^-5*x^4//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+630*a*b*d^-5*x^4//3*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+-1890*im*a*b*d^-8*x^1//3*Polylogarithms.polylog(8,-1*exp(2*im*(c+d*x^1//3)))+-252*im*a*b*d^-4*x^5//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+24*im*a*b*d^-2*x^7//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3)))+1260*im*a*b*x*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^1//3))))
@test integrate(x*(a+b*tan(c+d*x^1//3))^2, x) == :((1/2)*a^2*x^2+-1//2*b^2*x^2+-45//2*b^2*d^-6*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+im*a*b*x^2+-3*im*b^2*d^-1*x^5//3+3*b^2*d^-1*x^5//3*tan(c+d*x^1//3)+15*b^2*d^-2*x^4//3*log(1+exp(2*im*(c+d*x^1//3)))+45*b^2*d^-4*x^2//3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-30*im*x*b^2*d^-3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3)))+-30*a*b*x*d^-3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+-6*a*b*d^-1*x^5//3*log(1+exp(2*im*(c+d*x^1//3)))+45*im*b^2*d^-5*x^1//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+45*a*b*d^-5*x^1//3*Polylogarithms.polylog(5,-1*exp(2*im*(c+d*x^1//3)))+45//2*im*a*b*d^-6*Polylogarithms.polylog(6,-1*exp(2*im*(c+d*x^1//3)))+-45*im*a*b*d^-4*x^2//3*Polylogarithms.polylog(4,-1*exp(2*im*(c+d*x^1//3)))+15*im*a*b*d^-2*x^4//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3))))
@test integrate((a+b*tan(c+d*x^1//3))^2, x) == :(x*a^2+-1*x*b^2+-3*im*b^2*d^-1*x^2//3+-3*im*b^2*d^-3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3)))+-3*a*b*d^-3*Polylogarithms.polylog(3,-1*exp(2*im*(c+d*x^1//3)))+2*im*a*b*x+3*b^2*d^-1*x^2//3*tan(c+d*x^1//3)+6*b^2*d^-2*x^1//3*log(1+exp(2*im*(c+d*x^1//3)))+-6*a*b*d^-1*x^2//3*log(1+exp(2*im*(c+d*x^1//3)))+6*im*a*b*d^-2*x^1//3*Polylogarithms.polylog(2,-1*exp(2*im*(c+d*x^1//3))))
@test integrate(x^2*(a+b*tan(c+d*x^1//3))^-1, x) == :(x^3*(3a+3*im*b)^-1+-945*b*d^-9*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(9,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-315*b*d^-5*x^4//3*(a^2+b^2)^-1*Polylogarithms.polylog(5,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+3*b*d^-1*x^8//3*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+42*b*d^-3*x^2*(a^2+b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+945*b*d^-7*x^2//3*(a^2+b^2)^-1*Polylogarithms.polylog(7,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-630*im*b*x*d^-6*(a^2+b^2)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-12*im*b*d^-2*x^7//3*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+126*im*b*d^-4*x^5//3*(a^2+b^2)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+945*im*b*d^-8*x^1//3*(a^2+b^2)^-1*Polylogarithms.polylog(8,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3))))
@test integrate(x*(a+b*tan(c+d*x^1//3))^-1, x) == :(x^2*(2a+2*im*b)^-1+-45*im*b*d^-6*(4*a^2+4*b^2)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-45*b*d^-5*x^1//3*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(5,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+3*b*d^-1*x^5//3*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+15*b*x*d^-3*(a^2+b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-15*im*b*d^-2*x^4//3*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+45*im*b*d^-4*x^2//3*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3))))
@test integrate((a+b*tan(c+d*x^1//3))^-1, x) == :(x*(a+im*b)^-1+3*b*d^-3*(2*a^2+2*b^2)^-1*Polylogarithms.polylog(3,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+3*b*d^-1*x^2//3*(a^2+b^2)^-1*log(1+(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3)))+-3*im*b*d^-2*x^1//3*(a^2+b^2)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-2*(a^2+b^2)*exp(2*im*(c+d*x^1//3))))
@test integrate(x^2*(a+b*tan(c+d*x^1//3))^-2, x) == :(1//3*x^3*(a+-1*im*b)^-2+-4//3*b^2*x^3*(a^2+b^2)^-2+4*b*x^3*(a+-1*im*b)^-2*(-3b+3*im*a)^-1+-1260*x*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-1260*x*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-945*b*d^-9*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(9,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-24*b^2*d^-2*x^7//3*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-6*im*b^2*d^-1*x^8//3*(a^2+b^2)^-2+24*b^2*d^-2*x^7//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+252*b^2*d^-4*x^5//3*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+252*b^2*d^-4*x^5//3*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+945*im*b^2*d^-9*(a^2+b^2)^-2*Polylogarithms.polylog(8,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+945*im*b^2*d^-9*(a^2+b^2)^-2*Polylogarithms.polylog(9,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+1890*b^2*d^-8*x^1//3*(a^2+b^2)^-2*Polylogarithms.polylog(7,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+1890*b^2*d^-8*x^1//3*(a^2+b^2)^-2*Polylogarithms.polylog(8,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-1890*im*b^2*d^-7*x^2//3*(a^2+b^2)^-2*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-1890*im*b^2*d^-7*x^2//3*(a^2+b^2)^-2*Polylogarithms.polylog(7,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-1890*b*d^-8*x^1//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(8,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-630*b*d^-5*x^4//3*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-252*b*d^-4*x^5//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-84*im*b^2*d^-3*x^2*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-84*im*b^2*d^-3*x^2*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-6*im*b^2*d^-1*x^8//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b*d^-1*x^8//3*(a+im*b)^-1*(a+-1*im*b)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b^2*d^-1*x^8//3*(a+im*b)^-1*(b+im*a)^-2*(-1b+im*a+(b+im*a)*exp(2*im*(c+d*x^1//3)))^-1+24*b*d^-2*x^7//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+84*b*d^-3*x^2*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+630*im*b^2*d^-5*x^4//3*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+630*im*b^2*d^-5*x^4//3*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+1260*b*x*d^-6*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+1890*b*d^-7*x^2//3*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(7,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3))))
@test integrate(x*(a+b*tan(c+d*x^1//3))^-2, x) == :((1/2)*x^2*(a+-1*im*b)^-2+-2*b^2*x^2*(a^2+b^2)^-2+2*b*x^2*(a+-1*im*b)^-2*(-1b+im*a)^-1+-45//2*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-45//2*b^2*d^-6*(a^2+b^2)^-2*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-15*b^2*d^-2*x^4//3*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-6*im*b^2*d^-1*x^5//3*(a^2+b^2)^-2+15*b^2*d^-2*x^4//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+45*b*d^-6*(a+-1*im*b)^-2*(-2b+2*im*a)^-1*Polylogarithms.polylog(6,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+45*b^2*d^-4*x^2//3*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+45*b^2*d^-4*x^2//3*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-45*b*d^-5*x^1//3*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-45*b*d^-4*x^2//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-30*im*x*b^2*d^-3*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-30*im*x*b^2*d^-3*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-6*im*b^2*d^-1*x^5//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b*d^-1*x^5//3*(a+im*b)^-1*(a+-1*im*b)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b^2*d^-1*x^5//3*(a+im*b)^-1*(b+im*a)^-2*(-1b+im*a+(b+im*a)*exp(2*im*(c+d*x^1//3)))^-1+15*b*d^-2*x^4//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+30*b*x*d^-3*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+45*im*b^2*d^-5*x^1//3*(a^2+b^2)^-2*Polylogarithms.polylog(4,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+45*im*b^2*d^-5*x^1//3*(a^2+b^2)^-2*Polylogarithms.polylog(5,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3))))
@test integrate((a+b*tan(c+d*x^1//3))^-2, x) == :(x*(a+-1*im*b)^-2+-4*x*b^2*(a^2+b^2)^-2+4*b*x*(a+-1*im*b)^-2*(-1b+im*a)^-1+-6*im*b^2*d^-1*x^2//3*(a^2+b^2)^-2+-6*b^2*d^-2*x^1//3*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-3*im*b^2*d^-3*(a^2+b^2)^-2*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-3*im*b^2*d^-3*(a^2+b^2)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+3*b*d^-3*(a+im*b)^-1*(a+-1*im*b)^-2*Polylogarithms.polylog(3,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b^2*d^-2*x^1//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+-6*im*b^2*d^-1*x^2//3*(a^2+b^2)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b*d^-1*x^2//3*(a+im*b)^-1*(a+-1*im*b)^-2*log(1+(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b*d^-2*x^1//3*(a+-1*im*b)^-2*(-1b+im*a)^-1*Polylogarithms.polylog(2,-1*(a+im*b)^-1*(a+-1*im*b)*exp(2*im*(c+d*x^1//3)))+6*b^2*d^-1*x^2//3*(a+im*b)^-1*(b+im*a)^-2*(-1b+im*a+(b+im*a)*exp(2*im*(c+d*x^1//3)))^-1)

max_stars_repo_path: test/4 Trig functions/4.3 Tangent/4.3.11 (e x)^m (a+b tan(c+d x^n))^p.jl


text:
# Computes the function cel(kc,p,a,b) from Bulirsch (1969).
function cel_bulirsch(k2::T,p::T,a::T,b::T) where {T <: Real}
@assert (k2 <= 1.0)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 != 1.0
  kc = sqrt(1.0-k2)
else
  kc = eps(k2)
end
# Initialize values:
ee = kc; m=1.0
if p > 0.0
  p = sqrt(p); b /= p
else
  q=k2; g=1.0-p; f = g-k2
  q *= (b-a*p); p=sqrt(f/g); a=(a-b)/g
  b = -q/(g*g*p)+a*p
end
# Compute recursion:
f=a; a += b/p; g=ee/p; b += f*g; b +=b; p +=g; g=m; m += kc
while abs(g-kc) > g*ca
  kc=sqrt(ee); kc += kc; ee = kc*m
  f=a; a += b/p; g=ee/p; b += f*g; b +=b; p +=g; g=m; m += kc
end
return pi/2*(a*m+b)/(m*(m+p))
end

# Version called with kc (this is to improve precision of computation):
function cel_bulirsch(k2::T,kc::T,p::T,a::T,b::T) where {T <: Real}
@assert (k2 <= 1.0)
ca = sqrt(eps(k2))
# Avoid undefined k2=1 case:
if k2 == 1.0
  kc = eps(k2)
end
# Initialize values:
ee = kc; m=1.0
if p > 0.0
  p = sqrt(p); b /= p
else
  q=k2; g=1.0-p; f = g-k2
  q *= (b-a*p); p=sqrt(f/g); a=(a-b)/g
  b = -q/(g*g*p)+a*p
end
# Compute recursion:
f=a; a += b/p; g=ee/p; b += f*g; b +=b; p +=g; g=m; m += kc
while abs(g-kc) > g*ca
  kc=sqrt(ee); kc += kc; ee = kc*m
  f=a; a += b/p; g=ee/p; b += f*g; b +=b; p +=g; g=m; m += kc
end
return pi/2*(a*m+b)/(m*(m+p))
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
module ItemGraphs

import LightGraphs: AbstractGraph, DiGraph, add_edge!, add_vertex!, nv,
    dijkstra_shortest_paths, enumerate_paths, has_vertex, has_path

export ItemGraph, add_vertex!, add_edge!, items, has_vertex, has_path,
    ItemGraphException, edgeitem, edgeitems

struct ItemGraphException <: Exception
    msg::String
end
Base.show(io::IO, ex::ItemGraphException) = print(io, ex.msg)

"""
    ItemGraph{T,S}([graph]; lazy=false)

Create a new ItemGraph with items of type `T` based on a `graph` from `LightGraphs`.
Items of type `S` can be assigned to the edges of the graph, this type can be omitted
and will be `Float64` by default.
`graph` can be omitted as well and a `DiGraph` will be used by default.
If `lazy` is set to true, the paths between items will be computed on-demand when
[`items`](@ref) is called.
Otherwise all paths are recomputed whenever a new edge is inserted.
"""
struct ItemGraph{T,S,G}
    graph::G
    ids::Dict{Int,T}
    items::Dict{T,Int}
    edges::Dict{T,Dict{T,S}}
    paths::Dict{T,Dict{T,Vector{T}}}
    lazy::Bool
    ItemGraph{T,S}(graph::G; lazy = false) where {T, S, G <: AbstractGraph} = new{T,S,G}(graph,
        Dict{Int,T}(),
        Dict{T,Int}(),
        Dict{T,Dict{T,S}}(),
        Dict{T,Dict{T,Vector{T}}}(),
        lazy,)
end

ItemGraph{T}(; kwargs...) where {T} = ItemGraph{T,Float64}(DiGraph(); kwargs...)
ItemGraph{T, S}(; kwargs...) where {T, S} = ItemGraph{T,S}(DiGraph(); kwargs...)

getid(graph::ItemGraph{T}, item::T) where {T} = graph.items[item]
getitem(graph::ItemGraph, id) = graph.ids[id]
"""
    has_vertex(graph, item)

Return true if `item` is contained in `graph`.
"""
has_vertex(graph::ItemGraph, item) = haskey(graph.items, item)
"""
    has_path(graph, from, to)

Return true if there is a path between `from` and `to` in `graph`.
"""
has_path(graph::ItemGraph, from, to) = has_path(graph.graph, getid(graph, from), getid(graph, to))

"""
    add_vertex!(graph, item)

Add `item` to `graph`.
"""
function add_vertex!(graph::ItemGraph{T}, item::T) where T
    has_vertex(graph, item) && return

    add_vertex!(graph.graph)
    id = nv(graph.graph)
    push!(graph.ids, id => item)
    push!(graph.items, item => id)
    nothing
end

"""
    add_edge!(graph, from, to, [item])

Add an edge between `from` and `to` to `graph`. Optionally assign an item
to the edge.
"""
function add_edge!(graph::ItemGraph{T,S}, from::T, to::T,
        item::S=zero(S)) where {T, S}
    add_vertex!(graph, from)
    add_vertex!(graph, to)
    add_edge!(graph.graph, getid(graph, from), getid(graph, to))

    !graph.lazy && calculate_paths!(graph)

    edges = get!(graph.edges, from, Dict{Int,S}())
    push!(edges, to => item)
    nothing
end

"""
    edgeitem(graph, from, to)

Get the item assigned to the edge between `from` and `to`.
"""
function edgeitem(graph::ItemGraph{T}, from::T, to::T) where T
    graph.edges[from][to]
end

"""
    edgeitems(graph, from, to)

Get all items assigned to the edges between `from` and `to`.
"""
function edgeitems(graph::ItemGraph{T,S}, from::T, to::T) where {T, S}
    path = calculate_path(graph, from, to)
    edges = Vector{S}(undef, length(path) - 1)
    for i in eachindex(edges)
        edges[i] = graph.edges[path[i]][path[i+1]]
    end
    edges
end

"""
    items(graph, from, to)

Get the items on the path between and including `from` and `to`.
Will throw an `ItemGraphsException` if either `from` or `to` are not a part of `graph` or if there
is no path between them.
"""
function items(graph, from, to)
    for vertex in (from, to)
        !has_vertex(graph, vertex) && throw(ItemGraphException("There is no vertex '$vertex'."))
    end

    !has_path(graph, from, to) && throw(ItemGraphException("There is no path between '$from' and '$to'."))

    calculate_path(graph, from, to)
end

function calculate_paths!(graph::ItemGraph{T}) where T
    for (origin, oid) in graph.items
        d = dijkstra_shortest_paths(graph.graph, oid)
        for (target, tid) in graph.items
            origin == target && continue
            path = getitem.(Ref(graph), enumerate_paths(d, tid))
            paths = get!(graph.paths, origin, Dict{T,Vector{T}}())
            push!(paths, target => path)
        end
    end
end

function calculate_path(graph::ItemGraph{T}, from, to) where T
    !graph.lazy && return graph.paths[from][to]

    origin = getid(graph, from)
    target = getid(graph, to)

    getitem.(Ref(graph), enumerate_paths(dijkstra_shortest_paths(graph.graph, origin), target))::Vector{T}
end

end # module


text:
module MosimoBase

using Reexport: @reexport
@reexport using LinearAlgebra, Statistics, StaticArrays, Random, DataFrames
using MacroTools
using Random: AbstractRNG, GLOBAL_RNG

export Nullable
include("./bases/bases.jl")

export Vector2, Vector3, MosiVector,
    Vector2s, Vector3s,
    AbstractVector2s, AbstractVector3s,
    flatten, unflatten,
    to_matrix, ndims
include("./bases/vectors.jl")

export read_vector_all,
    read_vectors, read_vectors_all
include("./utils/io.jl")

export MosiRNG,
    make_seed, new_rng, restore_rng, rng_state
include("./utils/rng.jl")

export period_check, update_periods!,
    default_distance, pbc_distance, pbc_box,
    original_vector, original_vectors
include("./utils/pbc.jl")

export polar_to_cartesian, cartesian_to_polar,
    random_2d_direction,
    random_point_on_circle, random_point_on_disk
include("./utils/2d_math.jl")

export spherical_to_cartesian, cartesian_to_spherical,
    ρ_hat, θ_hat, φ_hat,
    ∂Ω_∂ρ, ∂Ω_∂θ, ∂Ω_∂φ
include("./utils/spherical_coordinates.jl")

export random_3d_direction, random_direction_on_plane,
    rotate_around, decompose_vector
include("./utils/3d_math.jl")

export orthogonalize, projection_matrix
include("./utils/projection_matrix.jl")

export MosiSystem, Molecule, MolecularSystem, ConfigurationSystem,
    natoms, positions, velocities, periods, has_pbc, original_positions,
    distance_function, update_periods!
include("./types/system.jl")

export MosiModel, UnknownModel,
    name, is_3d, vectype, constraints, constraint_gradients,
    potential_energy_function, force_function, potential_energy_gradients
include("./types/model.jl")

export potential_energy
include("./types/energy.jl")

export SimulationState, SimulationTape, SimulationSetup, SimulationResult, SimulationError,
    SimpleTape,
    system, times, mosi_model, get_configuration_func, init_state, tape, observables
include("./types/simulation.jl")

export TapeFiles, MultiFileMemoryMapTape,
    has_vs, update_tape, configuration
include("./types/tape_files.jl")

export Ensemble, CanonicalEnsemble, MicrocanonicalEnsemble, PotentialEnergyLandscapeEnsemble
include("./types/ensembles.jl")

end


text:
using SpikingNetworks
parseParameters"""
τ=10
N=1000
v0=[20*i/N for i=1:N]"""
parseEquations"""
dv/dt=(v0-v)/τ"""
net=init"""
N=N
v=0*ones(N)"""
spike"""
v.>10
5"""
reset"v[i]=0"

net.dt=0.1
@time run!(net,1000/net.dt)
t,neuron=rasterPlot(net)
using MatlabPlot
mplot(SpikingNetwork.v0,map(x->length(x),net.tSpike))


text:
module StateSpaceCPG

using Optim, Distributions, LinearAlgebra, StaticArrays, Dates, Printf, StatsBase

import Base: size, show

export statespace, kfas, diagnostics

include("prints.jl")
include("structures.jl")
include("utils.jl")
include("model_unknowns.jl")
include("models.jl")
include("estimation.jl")
include("opt_methods.jl")
include("kalman_utils.jl")
include("univariate_kalman.jl")
include("kalman.jl")
include("sqrt_kalman.jl")
include("forecast.jl")
include("diagnostics.jl")

"""
    statespace(model; filter_type = KalmanFilter, opt_method = LBFGS(), verbose = 1)

Estimate the pre-specified state-space model. The function will only estimate the entries that are declared `NaN`. If there are no NaNs in the `model` it will
only perform the filter and smoother computations.
"""
function statespace(model::StateSpaceModel{T}; filter_type::DataType = KalmanFilter{T},
                    opt_method::AbstractOptimizationMethod = LBFGS(model, 3), verbose::Int = 1) where T

    if !(verbose in [0, 1, 2, 3])
        @warn("Incorrect verbose value input (should be 0, 1, 2 or 3): switching to default value 1")
        verbose = 1
    end

    print_header(verbose)

    # Maximum likelihood estimation
    estimate_statespace!(model, filter_type, opt_method; verbose = verbose)

    # Kalman filter and smoothing
    filter_output, smoothed_state = kfas(model, filter_type)

    print_bottom(verbose)

    return StateSpace{T}(model, filter_output, smoothed_state, filter_type, opt_method)
end

"""
    kfas(model::StateSpaceModel{T}, filter_type::DataType) where T

Perform Kalman filter and smoother according to the chosen `filter_type`.
"""
function kfas(model::StateSpaceModel{T}, filter_type::DataType) where T
    error(filter_type , " not implemented") # Returns an error if it cannot
                                            # find a specialized kfas
end

end # end module

max_stars_repo_name: PaulMainwood/StateSpaceCPG

max_issues_repo_name: PaulMainwood/StateSpaceCPG

max_forks_repo_name: PaulMainwood/StateSpaceCPG


text:
function hello()
    return "Hello, World!";
end

max_stars_repo_name: AbstractGeek/exercism-expeditions

max_issues_repo_name: AbstractGeek/exercism-expeditions

max_forks_repo_name: AbstractGeek/exercism-expeditions


text:
# include("pascaltriangle.jl") # define pascaltriangle(m,n)

function in2bEne(jj::Int64, Msize::Int64, Np::Int64, matp::Matrix{Int64})

    # this function return the index of Fock state at the index of many body state
    indfk = copy(jj)

    # matp = zeros(Int,Msize+1,Np+1);
    # pascaltriangle!(Msize,Np,matp) # the size is Msize+1 times Np+1
    # note the indices are m+1 and n+1 for N^m_n

    # unlike Ponomarev's way, I count the index of Fock state in the backwards
    # so that "1" means |Np,0,0,...>.
    indfk = matp[Msize+1,Np+1] + 1 - indfk # the indices are m+1 and n+1 for N^m_n

    if indfk > matp[Msize+1,Np+1] # the indices are m+1 and n+1 for N^m_n
       error("indfk larger than the maximum of N^M_N.")
    end

    indfk = indfk - 1
    # vecmb = spzeros(Int64,Msize+1) #spzeros(Int64,Msize+1)

    indM = Msize-1
    indN = Np

    Eneminumhalf = 0

    while true

          if indN == 0
             break
          end

          if indfk >= matp[indM+1,indN+1]

             indfk = indfk - matp[indM+1,indN+1] # the indices are m+1 and n+1 for N^m_n
             # vecmb[Msize-indM] = vecmb[Msize-indM] + 1
             indN = indN - 1

             Eneminumhalf = Eneminumhalf + ceil(Int64,(Msize-indM)/2)-1

          else
             indM = indM - 1
          end

    end

    # the last element for the coefficient
    # vecmb[Msize+1] = 1

    # return vecmb
    return Eneminumhalf #+1/2

end

max_stars_repo_name: leios/SOCED

max_issues_repo_name: leios/SOCED

max_forks_repo_name: leios/SOCED


text:
"""
Aluminum plate with force at the center. One quarter is modeled.

The plate is excited with the transducer at the initial time (i.e. distributed
loading mimicking a concentrated force is applied).

Reference:  Ostachowicz W, Kudela P, Krawczuk M, Zak A. ch. 5: i-xii; John Wiley & Sons, Ltd . 2012
"""
module plate_expl_damp_explore_examples

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.AssemblyModule
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using SymRCM
using VisualStructures: default_layout_3d, plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using PlotlyJS
using Gnuplot;
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection
using ThreadedSparseCSR
using UnicodePlots
using PGFPlotsX
using InteractiveUtils
using BenchmarkTools
using FinEtools.MeshExportModule.VTKWrite: vtkwritecollection, vtkwrite

# Aluminium Alloy
const E = 68*phun("GPa");
const nu = 0.33;
const rho = 2660.0*phun("kg/m^3")
const d1 = 1000*phun("mm")
const dD = 106*phun("mm")
const dF = dD + 197*phun("mm")
const thickness = 1.0*phun("mm") # the book says 1.0
const D = E * thickness^3 / 12 / (1 - nu^2)
const cp = (omega) -> sqrt(omega) * (D / rho / thickness)^(1/4)
const cg = (omega) -> 2*cp(omega)
const omegad = 2*pi*1000*phun("Hz")
const ksi = 100.0/2/omegad
const carrier_frequency = 35*phun("kilo*Hz")
const modulation_frequency = 7*phun("kilo*Hz")
const totalforce = 1*phun("N")
const forcepatchradius = 3*phun("mm")
const forcedensity = totalforce/(pi*forcepatchradius^2)
const tend = 1.0*phun("milli*s")
const visualize = true
const visualizeclear = !true
const visualizevtk = !true
const distributedforce = !true

const color = "black"
const algorithm = :cd
const color = "red"
const algorithm = :predcorr

function cd_loop!(M, K, ksi, U0, V0, tend, dt, force!, peek)
    U = similar(U0)
    V = similar(U0)
    A = similar(U0)
    F = similar(U0)
    E = similar(U0)
    C = similar(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invMC = similar(U0)
    invMC .= 1.0 ./ (vec(diag(M)) .+ (dt/2) .* C)
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end
    dt2_2 = ((dt^2)/2)
    dt_2 = (dt/2)

    t = 0.0
    @. U = U0; @. V = V0; # Initial Conditions
    A .= invMC .* force!(F, t); # Compute initial acceleration
    peek(0, U, V, t)
    for step in 1:nsteps
        t = t + dt
        @. U += dt*V + dt2_2*A;   # Displacement update
        force!(F, t);  # External forces are computed
        ThreadedSparseCSR.bmul!(E, K, U)  # Evaluate elastic restoring forces
        @. F -= E + C * (V + dt_2 * A)  # Calculate total force
        @. V += dt_2 * A; # Update the velocity, part one: add old acceleration
        @. A = invMC * F # Compute the acceleration
        @. V += dt_2 * A; # Update the velocity, part two: add new acceleration
        peek(step, U, V, t)
    end
end

function predcorr_loop!(M, K, ksi, U0, V0, tend, dt, force!, peek)
    U = similar(U0)
    V = similar(U0)
    A = similar(U0)
    F = similar(U0)
    E = similar(U0)
    C = similar(U0)
    C .= (ksi*2*omegad) .* vec(diag(M))
    invM = similar(U0)
    invM .= 1.0 ./ vec(diag(M))
    nsteps = Int64(round(tend/dt))
    if nsteps*dt < tend
        dt = tend / (nsteps+1)
    end
    dt2_2 = ((dt^2)/2)
    dt_2 = (dt/2)

    t = 0.0
    # Initial Conditions
    @. U = U0; @. V = V0
    A .= invM .* force!(F, t);
    peek(0, U, V, t)
    for step in 1:nsteps
        t = t + dt
        @. U += dt*V + dt2_2*A; # Displacement update
        @. V += dt_2 * A # Predict the velocity
        force!(F, t); # External loading
        F .-= ThreadedSparseCSR.bmul!(E, K, U) # Add elastic restoring forces
        @. F -= C * V # Add damping force
        @. A = invM * F # Compute the acceleration
        @. V += dt_2 * A # Correct the velocity
        peek(step, U, V, t)
    end
end

function _execute(nref = 2, nthr = 0)
    tolerance = d1/nref/100

    xs = sort(unique(vcat(
        linearspace(0.0, dD, 2),
        linearspace(dD, dF, 3),
        linearspace(dF, d1/2, 3),
        )))
    ys = collect(xs)

    fens, fes = T3blockx(xs, ys); # Mesh
    for r in 1:nref
        fens, fes = T3refine(fens, fes)
    end

    fens.xyz = xyz3(fens)

    @info "Mesh $(count(fens)) nodes, $(count(fes)) elements"

    # vtkwrite("plate_expl.vtu", fens, fes)
    # vtkwrite("plate_expl-boundary.vtu", fens, bfes)

    # Renumber the nodes
    femm = FEMMBase(IntegDomain(fes, TriRule(1)))
    C = connectionmatrix(femm, count(fens))
    perm = symrcm(C)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), mater)
    # Set up
    femm.drilling_stiffness_scale = 1.0

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Symmetry EBC's
    lx = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    ly = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2, 4, 6]
        setebc!(dchi, ly, true, i)
    end
    for i in [1, 5, 6]
        setebc!(dchi, lx, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi, perm);
    # numberdofs!(dchi);

    # Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    SM = FinEtoolsFlexStructures.AssemblyModule
    K = FEMMShellT3FFModule.stiffness(femm, SM.SysmatAssemblerSparseCSRSymm(0.0), geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, SysmatAssemblerSparseDiag(), geom0, dchi);

    # Solve
    function pwr(K, M)
        invM = fill(0.0, size(M, 1))
        invM .= 1.0 ./ (vec(diag(M)))
        v = rand(size(M, 1))
        w = fill(0.0, size(M, 1))
        for i in 1:30
            ThreadedSparseCSR.bmul!(w, K, v)
            wn = norm(w)
            w .*= (1.0/wn)
            v .= invM .* w
            vn = norm(v)
            v .*= (1.0/vn)
        end
        sqrt((v' * (K * v)) / (v' * M * v))
    end
    @time omega_max = pwr(K, M)
    @show omega_max = max(omega_max, 20*2*pi*carrier_frequency)
    # @show dt = Float64(0.9* 2/omega_max) * (sqrt(1+ksi^2) - ksi)
    @show dt = Float64(0.99 * 2/omega_max)

    U0 = gathersysvec(dchi)
    V0 = deepcopy(U0)

    points = Dict(
    "A" => selectnode(fens; nearestto=[0.0 0.0 0.0])[1],
    "B" => selectnode(fens; nearestto=[dD 0.0 0.0])[1],
    "C" => selectnode(fens; nearestto=[dF 0.0 0.0])[1],
    "D" => selectnode(fens; nearestto=[dD dD 0.0])[1],
    "E" => selectnode(fens; nearestto=[dD dF 0.0])[1],
    "F" => selectnode(fens; nearestto=[dF dF 0.0])[1],
    )

    pointdofs = Dict()
    for k in keys(points)
        pointdofs[k] = dchi.dofnums[points[k], 3]
    end

    function computetrac!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        dx = XYZ[1] - fens.xyz[points["A"], 1]
        dy = XYZ[2] - fens.xyz[points["A"], 2]
        dz = XYZ[3] - fens.xyz[points["A"], 3]
        d = sqrt(dx^2+dy^2+dz^2)
        forceout .= 0.0
        if d < forcepatchradius
            forceout[3] = -forcedensity
        end
        return forceout
    end

    if distributedforce
        lfemm = FEMMBase(IntegDomain(fes, TriRule(6)))
        fi = ForceIntensity(FFlt, 6, computetrac!);
        Fmag = distribloads(lfemm, geom0, dchi, fi, 2);
    else
        Fmag = fill(0.0, dchi.nfreedofs)
        Fmag[pointdofs["A"]] = -totalforce/4 # only a quarter of the plate is modeled
    end

    @show sum(Fmag), -totalforce/4

    function force!(F, t)
        mul = 0.0
        if t <= 1/modulation_frequency
            mul = 0.5 * (1 - cos(2*pi*modulation_frequency*t)) * sin(2*pi*carrier_frequency*t)
        end
        F .= mul .* Fmag
        return F
    end

    nsteps = Int(round(tend/dt))
    pointvelocities = Dict()
    for k in keys(points)
        pointvelocities[k] = fill(0.0, nsteps+1)
    end
    displacements = []
    nbtw = Int(round(nsteps/100))

    peek(step, U, V, t) = begin
        for k in keys(points)
            b = pointvelocities[k]
            b[step+1] = U[pointdofs[k]]
        end
        if rem(step+1, nbtw) == 0
            push!(displacements, deepcopy(U))
        end
        nothing
    end

    @info "$nsteps steps"

    nth = (nthr == 0 ? Base.Threads.nthreads() : nthr)
    @info "$nth threads used"

    if algorithm == :cd
        cd_loop!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek)
    else
        predcorr_loop!(M, K, ksi, U0, V0, nsteps*dt, dt, force!, peek)
    end

    savecsv("plate_expl_$nref", t = collect(0.0:dt:(nsteps*dt))./phun("milli*s"),
        v = pointvelocities["C"]./phun("mm/s"))

    if visualize
        # @gp  "set terminal windows 0 "  :-
        if visualizeclear
            @gp "clear"
        end

        # for k in sort(string.(keys(points)))
        #     @gp  :- collect(0.0:dt:(nsteps*dt)) pointvelocities[k] " lw 2 lc rgb '$color' with lines title 'Deflection at $(string(k))' "  :-
        # end

        @gp  :- "set xrange [0:1.0] " :- # in milliseconds
        @gp  :- "set yrange [-1.1e-5:1.1e-5] " :-
        for (i, k) in enumerate(["C"])
            @gp  :-  collect(0.0:dt:(nsteps*dt))./phun("milli*s") pointvelocities[k]./phun("mm/s") " lw 2 lc rgb '$color' with lines title 'Velocity at $(string(k))' "  :-
            # @gp  :- collect(0.0:dt:(nsteps*dt))./phun("milli*s") pointvelocities[k] " lw 2 lc rgb '$color' pt $i with lp title 'Velocity at $(string(k))' "  :-
        end

        @gp  :- "set tics font \"Times-Roman\"" :-
        @gp  :- "set xlabel 'Time [ms]' font \"Times-Roman\"" :-
        @gp  :- "set ylabel 'Velocity [mm/s]' font \"Times-Roman\"" :-
        @gp  :- "set title 'Free-floating center-loaded plate' font \"Times-Roman\""

        objects = []

        for (i, k) in enumerate(["C"])
            xdata = collect(0.0:dt:(nsteps*dt))./phun("milli*s")
            ydata = pointvelocities[k]./phun("mm/s")
            @pgf p = PGFPlotsX.Plot(
            {
            color = "black",
            line_width  = 1.0
            },
            Coordinates([v for v in  zip(xdata, ydata)])
            )
            push!(objects, p)
            push!(objects, LegendEntry("Point $k"))
        end

        @pgf ax = Axis(
        {
        xlabel = "Time [ms]",
        ylabel = "Velocity [mm/s]",
        ymin = -1.1e-5,
        ymax = 1.1e-5,
        xmin = 0.0,
        xmax = 1.0,
        xmode = "linear",
        ymode = "linear",
        yminorgrids = "true",
        grid = "both",
        legend_style = {
        at = Coordinate(0.5, 1.05),
        anchor = "south",
        legend_columns = -1
        },
        },
        objects...
        )

        display(ax)
        pgfsave("plate_expl_$nref.pdf", ax)

        # Visualization
        if visualizevtk
            @info "Dumping visualization"
            times = Float64[]
            vectors = []
            for i in 1:length(displacements)
                scattersysvec!(dchi, displacements[i])
                push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
                push!(times, i*dt*nbtw)
            end
            vtkwritecollection("plate_expl_$nref", fens, fes, times; vectors = vectors)
        end
    end

    return true
end

function test(nrefs = [5], nthr = 0)
    @info "Cracked Plate, nrefs = $nrefs: parallel CSR"
    results = []
    for nref in nrefs
        push!(results, _execute(nref, nthr))
    end
    return results
end

function allrun(nrefs = [5,], nthr = 0)
    println("#####################################################")
    println("# test ")
    test(nrefs, nthr)

    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

allrun()

end # module
nothing


text:
function bracket_sign_change(f′, a, b; k=2)
    if a > b; a,b = b,a; end # ensure a < b
    center, half_width = (b+a)/2, (b-a)/2
    while f′(a)*f′(b) > 0
        half_width *= k
        a = center - half_width
        b = center + half_width
    end
    return (a,b)
end

df(x) = 10*cos(10*x)-3*sin(3*x)

as,bs = bracket_sign_change(df,3.4,3.8)
as,df(as),bs,df(bs)