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

module TestTriangular

debug = false
using Test, LinearAlgebra, SparseArrays, Random
using LinearAlgebra: BlasFloat, errorbounds, full!, naivesub!, transpose!,
    UnitUpperTriangular, UnitLowerTriangular,
    mul!, rdiv!, rmul!, lmul!

debug && println("Triangular matrices")

n = 9
srand(123)

debug && println("Test basic type functionality")
@test_throws DimensionMismatch LowerTriangular(randn(5, 4))
@test LowerTriangular(randn(3, 3)) |> t -> [size(t, i) for i = 1:3] == [size(Matrix(t), i) for i = 1:3]

# The following test block tries to call all methods in base/linalg/triangular.jl in order for a combination of input element types. Keep the ordering when adding code.
for elty1 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat}, Int)
    # Begin loop for first Triangular matrix
    for (t1, uplo1) in ((UpperTriangular, :U),
                        (UnitUpperTriangular, :U),
                        (LowerTriangular, :L),
                        (UnitLowerTriangular, :L))

        # Construct test matrix
        A1 = t1(elty1 == Int ? rand(1:7, n, n) : convert(Matrix{elty1}, (elty1 <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> cholesky(t't).U |> t -> uplo1 == :U ? t : copy(t')))

        debug && println("elty1: $elty1, A1: $t1")

        # Convert
        @test convert(AbstractMatrix{elty1}, A1) == A1
        @test convert(Matrix, A1) == A1

        # full!
        @test full!(copy(A1)) == A1

        # similar
        @test isa(similar(A1), t1)
        @test eltype(similar(A1)) == elty1
        @test isa(similar(A1, Int), t1)
        @test eltype(similar(A1, Int)) == Int
        @test isa(similar(A1, (3,2)), Matrix{elty1})
        @test isa(similar(A1, Int, (3,2)), Matrix{Int})

        #copyto!
        simA1 = similar(A1)
        copyto!(simA1, A1)
        @test simA1 == A1

        # getindex
        let mA1 = Matrix(A1)
            # linear indexing
            for i in 1:length(A1)
                @test A1[i] == mA1[i]
            end
            # cartesian indexing
            for i in 1:size(A1, 1), j in 1:size(A1, 2)
                @test A1[i,j] == mA1[i,j]
            end
        end
        @test isa(A1[2:4,1], Vector)

        # setindex! (and copy)
        A1c = copy(A1)
        for i = 1:size(A1, 1)
            for j = 1:size(A1, 2)
                if uplo1 == :U
                    if i > j
                        A1c[i,j] = 0
                        @test_throws ArgumentError A1c[i,j] = 1
                    elseif i == j && t1 == UnitUpperTriangular
                        A1c[i,j] = 1
                        @test_throws ArgumentError A1c[i,j] = 0
                    else
                        A1c[i,j] = 0
                        @test A1c[i,j] == 0
                    end
                else
                    if i < j
                        A1c[i,j] = 0
                        @test_throws ArgumentError A1c[i,j] = 1
                    elseif i == j && t1 == UnitLowerTriangular
                        A1c[i,j] = 1
                        @test_throws ArgumentError A1c[i,j] = 0
                    else
                        A1c[i,j] = 0
                        @test A1c[i,j] == 0
                    end
                end
            end
        end

        # istril/istriu
        if uplo1 == :L
            @test istril(A1)
            @test !istriu(A1)
        else
            @test istriu(A1)
            @test !istril(A1)
        end

        #tril/triu
        if uplo1 == :L
            @test tril(A1,0)  == A1
            @test tril(A1,-1) == LowerTriangular(tril(Matrix(A1), -1))
            @test tril(A1,1)  == t1(tril(tril(Matrix(A1), 1)))
            @test_throws ArgumentError tril!(A1, -n - 2)
            @test_throws ArgumentError tril!(A1, n)
            @test triu(A1,0)  == t1(diagm(0 => diag(A1)))
            @test triu(A1,-1) == t1(tril(triu(A1.data,-1)))
            @test triu(A1,1)  == zeros(size(A1)) # or just @test iszero(triu(A1,1))?
            @test_throws ArgumentError triu!(A1, -n)
            @test_throws ArgumentError triu!(A1, n + 2)
        else
            @test triu(A1,0)  == A1
            @test triu(A1,1)  == UpperTriangular(triu(Matrix(A1), 1))
            @test triu(A1,-1) == t1(triu(triu(Matrix(A1), -1)))
            @test_throws ArgumentError triu!(A1, -n)
            @test_throws ArgumentError triu!(A1, n + 2)
            @test tril(A1,0)  == t1(diagm(0 => diag(A1)))
            @test tril(A1,1)  == t1(triu(tril(A1.data,1)))
            @test tril(A1,-1) == zeros(size(A1)) # or just @test iszero(tril(A1,-1))?
            @test_throws ArgumentError tril!(A1, -n - 2)
            @test_throws ArgumentError tril!(A1, n)
        end

        # factorize
        @test factorize(A1) == A1

        # [c]transpose[!] (test views as well, see issue #14317)
        let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
            # transpose
            @test copy(transpose(A1)) == transpose(Matrix(A1))
            @test copy(transpose(viewA1)) == transpose(Matrix(viewA1))
            # adjoint
            @test copy(A1') == Matrix(A1)'
            @test copy(viewA1') == Matrix(viewA1)'
            # transpose!
            @test transpose!(copy(A1)) == transpose(A1)
            @test transpose!(t1(view(copy(A1).data, vrange, vrange))) == transpose(viewA1)
            # adjoint!
            @test adjoint!(copy(A1)) == adjoint(A1)
            @test adjoint!(t1(view(copy(A1).data, vrange, vrange))) == adjoint(viewA1)
        end

        # diag
        @test diag(A1) == diag(Matrix(A1))

        # real
        @test real(A1) == real(Matrix(A1))
        @test imag(A1) == imag(Matrix(A1))
        @test abs.(A1) == abs.(Matrix(A1))

        # Unary operations
        @test -A1 == -Matrix(A1)

        # copy and copyto! (test views as well, see issue #14317)
        let vrange = 1:n-1, viewA1 = t1(view(A1.data, vrange, vrange))
            # copy
            @test copy(A1) == copy(Matrix(A1))
            @test copy(viewA1) == copy(Matrix(viewA1))
            # copyto!
            B = similar(A1)
            copyto!(B, A1)
            @test B == A1
            B = similar(copy(transpose(A1)))
            copyto!(B, copy(transpose(A1)))
            @test B == copy(transpose(A1))
            B = similar(viewA1)
            copyto!(B, viewA1)
            @test B == viewA1
            B = similar(copy(transpose(viewA1)))
            copyto!(B, copy(transpose(viewA1)))
            @test B == transpose(viewA1)
        end

        #exp/log
        if (elty1 == Float64 || elty1 == ComplexF64) && (t1 == UpperTriangular || t1 == LowerTriangular)
            @test exp(Matrix(log(A1))) ≈ A1
        end

        # scale
        if (t1 == UpperTriangular || t1 == LowerTriangular)
            unitt = istriu(A1) ? UnitUpperTriangular : UnitLowerTriangular
            if elty1 == Int
                cr = 2
            else
                cr = 0.5
            end
            ci = cr * im
            if elty1 <: Real
                A1tmp = copy(A1)
                rmul!(A1tmp, cr)
                @test A1tmp == cr*A1
                A1tmp = copy(A1)
                lmul!(cr, A1tmp)
                @test A1tmp == cr*A1
                A1tmp = copy(A1)
                A2tmp = unitt(A1)
                mul!(A1tmp, A2tmp, cr)
                @test A1tmp == cr * A2tmp
            else
                A1tmp = copy(A1)
                rmul!(A1tmp, ci)
                @test A1tmp == ci*A1
                A1tmp = copy(A1)
                lmul!(ci, A1tmp)
                @test A1tmp == ci*A1
                A1tmp = copy(A1)
                A2tmp = unitt(A1)
                mul!(A1tmp, A2tmp, ci)
                @test A1tmp == ci * A2tmp
            end
        end

        # Binary operations
        @test A1*0.5 == Matrix(A1)*0.5
        @test 0.5*A1 == 0.5*Matrix(A1)
        @test A1/0.5 == Matrix(A1)/0.5
        @test 0.5\A1 == 0.5\Matrix(A1)

        # inversion
        @test inv(A1) ≈ inv(lu(Matrix(A1)))
        inv(Matrix(A1)) # issue #11298
        @test isa(inv(A1), t1)
        # make sure the call to LAPACK works right
        if elty1 <: BlasFloat
            @test LinearAlgebra.inv!(copy(A1)) ≈ inv(lu(Matrix(A1)))
        end

        # Determinant
        @test det(A1) ≈ det(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
        @test logdet(A1) ≈ logdet(lu(Matrix(A1))) atol=sqrt(eps(real(float(one(elty1)))))*n*n
        lada, ladb = logabsdet(A1)
        flada, fladb = logabsdet(lu(Matrix(A1)))
        @test lada ≈ flada atol=sqrt(eps(real(float(one(elty1)))))*n*n
        @test ladb ≈ fladb atol=sqrt(eps(real(float(one(elty1)))))*n*n

        # Matrix square root
        @test sqrt(A1) |> t -> t*t ≈ A1

        # naivesub errors
        @test_throws DimensionMismatch naivesub!(A1,Vector{elty1}(undef,n+1))

        # eigenproblems
        if !(elty1 in (BigFloat, Complex{BigFloat})) # Not handled yet
            vals, vecs = eigen(A1)
            if (t1 == UpperTriangular || t1 == LowerTriangular) && elty1 != Int # Cannot really handle degenerate eigen space and Int matrices will probably have repeated eigenvalues.
                @test vecs*diagm(0 => vals)/vecs ≈ A1 atol=sqrt(eps(float(real(one(vals[1])))))*(opnorm(A1,Inf)*n)^2
            end
        end

        # Condition number tests - can be VERY approximate
        if elty1 <:BlasFloat
            for p in (1.0, Inf)
                @test cond(A1,p) ≈ cond(A1,p) atol=(cond(A1,p)+cond(A1,p))
            end
            @test cond(A1,2) == cond(Matrix(A1),2)
        end

        if !(elty1 in (BigFloat, Complex{BigFloat})) # Not implemented yet
            svd(A1)
            elty1 <: BlasFloat && svd!(copy(A1))
            svdvals(A1)
        end

        # Begin loop for second Triangular matrix
        for elty2 in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat}, Int)
            for (t2, uplo2) in ((UpperTriangular, :U),
                                (UnitUpperTriangular, :U),
                                (LowerTriangular, :L),
                                (UnitLowerTriangular, :L))

                debug && println("elty1: $elty1, A1: $t1, elty2: $elty2")

                A2 = t2(elty2 == Int ? rand(1:7, n, n) : convert(Matrix{elty2}, (elty2 <: Complex ? complex.(randn(n, n), randn(n, n)) : randn(n, n)) |> t -> cholesky(t't).U |> t -> uplo2 == :U ? t : copy(t')))

                # Convert
                if elty1 <: Real && !(elty2 <: Integer)
                    @test convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
                elseif elty2 <: Real && !(elty1 <: Integer)
                    @test_throws InexactError convert(AbstractMatrix{elty2}, A1) == t1(convert(Matrix{elty2}, A1.data))
                end

                # Binary operations
                @test A1 + A2 == Matrix(A1) + Matrix(A2)
                @test A1 - A2 == Matrix(A1) - Matrix(A2)

                # Triangular-Triangualar multiplication and division
                @test A1*A2 ≈ Matrix(A1)*Matrix(A2)
                @test transpose(A1)*A2 ≈ transpose(Matrix(A1))*Matrix(A2)
                @test A1'A2 ≈ Matrix(A1)'Matrix(A2)
                @test A1*transpose(A2) ≈ Matrix(A1)*transpose(Matrix(A2))
                @test A1*A2' ≈ Matrix(A1)*Matrix(A2)'
                @test transpose(A1)*transpose(A2) ≈ transpose(Matrix(A1))*transpose(Matrix(A2))
                @test A1'A2' ≈ Matrix(A1)'Matrix(A2)'
                @test A1/A2 ≈ Matrix(A1)/Matrix(A2)
                @test A1\A2 ≈ Matrix(A1)\Matrix(A2)
                offsizeA = Matrix{Float64}(I, n+1, n+1)
                @test_throws DimensionMismatch offsizeA / A2
                @test_throws DimensionMismatch offsizeA / transpose(A2)
                @test_throws DimensionMismatch offsizeA / A2'
                @test_throws DimensionMismatch offsizeA * A2
                @test_throws DimensionMismatch offsizeA * transpose(A2)
                @test_throws DimensionMismatch offsizeA * A2'
                @test_throws DimensionMismatch transpose(A2) * offsizeA
                @test_throws DimensionMismatch A2'  * offsizeA
                @test_throws DimensionMismatch A2   * offsizeA
            end
        end

        for eltyB in (Float32, Float64, BigFloat, ComplexF32, ComplexF64, Complex{BigFloat})
            B = convert(Matrix{eltyB}, (elty1 <: Complex ? real(A1) : A1)*fill(1., n, n))

            debug && println("elty1: $elty1, A1: $t1, B: $eltyB")

            if !(eltyB in (BigFloat, Complex{BigFloat})) # rand does not support BigFloat and Complex{BigFloat} as of Dec 2015
                Tri = Tridiagonal(rand(eltyB,n-1),rand(eltyB,n),rand(eltyB,n-1))
                @test lmul!(Tri,copy(A1)) ≈ Tri*Matrix(A1)
            end

            # Triangular-dense Matrix/vector multiplication
            @test A1*B[:,1] ≈ Matrix(A1)*B[:,1]
            @test A1*B ≈ Matrix(A1)*B
            @test transpose(A1)*B[:,1] ≈ transpose(Matrix(A1))*B[:,1]
            @test A1'B[:,1] ≈ Matrix(A1)'B[:,1]
            @test transpose(A1)*B ≈ transpose(Matrix(A1))*B
            @test A1'B ≈ Matrix(A1)'B
            @test A1*transpose(B) ≈ Matrix(A1)*transpose(B)
            @test A1*B' ≈ Matrix(A1)*B'
            @test B*A1 ≈ B*Matrix(A1)
            @test transpose(B[:,1])*A1 ≈ transpose(B[:,1])*Matrix(A1)
            @test B[:,1]'A1 ≈ B[:,1]'Matrix(A1)
            @test transpose(B)*A1 ≈ transpose(B)*Matrix(A1)
            @test B'A1 ≈ B'Matrix(A1)
            @test B*transpose(A1) ≈ B*transpose(Matrix(A1))
            @test B*A1' ≈ B*Matrix(A1)'
            @test transpose(B[:,1])*transpose(A1) ≈ transpose(B[:,1])*transpose(Matrix(A1))
            @test B[:,1]'A1' ≈ B[:,1]'Matrix(A1)'
            @test transpose(B)*transpose(A1) ≈ transpose(B)*transpose(Matrix(A1))
            @test B'A1' ≈ B'Matrix(A1)'

            if eltyB == elty1
                @test mul!(similar(B),A1,B)  ≈ A1*B
                @test mul!(similar(B), A1, adjoint(B)) ≈ A1*B'
                @test mul!(similar(B), A1, transpose(B)) ≈ A1*transpose(B)
                @test mul!(similar(B), adjoint(A1), B) ≈ A1'*B
                @test mul!(similar(B), transpose(A1), B) ≈ transpose(A1)*B
                # test also vector methods
                B1 = vec(B[1,:])
                @test mul!(similar(B1),A1,B1)  ≈ A1*B1
                @test mul!(similar(B1), adjoint(A1), B1) ≈ A1'*B1
                @test mul!(similar(B1), transpose(A1), B1) ≈ transpose(A1)*B1
            end
            #error handling
            Ann, Bmm, bm = A1, Matrix{eltyB}(undef, n+1, n+1), Vector{eltyB}(undef, n+1)
            @test_throws DimensionMismatch lmul!(Ann, bm)
            @test_throws DimensionMismatch rmul!(Bmm, Ann)
            @test_throws DimensionMismatch lmul!(transpose(Ann), bm)
            @test_throws DimensionMismatch lmul!(adjoint(Ann), bm)
            @test_throws DimensionMismatch rmul!(Bmm, adjoint(Ann))
            @test_throws DimensionMismatch rmul!(Bmm, transpose(Ann))

            # ... and division
            @test A1\B[:,1] ≈ Matrix(A1)\B[:,1]
            @test A1\B ≈ Matrix(A1)\B
            @test transpose(A1)\B[:,1] ≈ transpose(Matrix(A1))\B[:,1]
            @test A1'\B[:,1] ≈ Matrix(A1)'\B[:,1]
            @test transpose(A1)\B ≈ transpose(Matrix(A1))\B
            @test A1'\B ≈ Matrix(A1)'\B
            @test A1\transpose(B) ≈ Matrix(A1)\transpose(B)
            @test A1\B' ≈ Matrix(A1)\B'
            @test transpose(A1)\transpose(B) ≈ transpose(Matrix(A1))\transpose(B)
            @test A1'\B' ≈ Matrix(A1)'\B'
            Ann, bm = A1, Vector{elty1}(undef,n+1)
            @test_throws DimensionMismatch Ann\bm
            @test_throws DimensionMismatch Ann'\bm
            @test_throws DimensionMismatch transpose(Ann)\bm
            if t1 == UpperTriangular || t1 == LowerTriangular
                @test_throws LinearAlgebra.SingularException naivesub!(t1(zeros(elty1,n,n)),fill(eltyB(1),n))
            end
            @test B/A1 ≈ B/Matrix(A1)
            @test B/transpose(A1) ≈ B/transpose(Matrix(A1))
            @test B/A1' ≈ B/Matrix(A1)'
            @test transpose(B)/A1 ≈ transpose(B)/Matrix(A1)
            @test B'/A1 ≈ B'/Matrix(A1)
            @test transpose(B)/transpose(A1) ≈ transpose(B)/transpose(Matrix(A1))
            @test B'/A1' ≈ B'/Matrix(A1)'

            # Error bounds
            !(elty1 in (BigFloat, Complex{BigFloat})) && !(eltyB in (BigFloat, Complex{BigFloat})) && errorbounds(A1, A1\B, B)

        end
    end
end

# Matrix square root
Atn = UpperTriangular([-1 1 2; 0 -2 2; 0 0 -3])
Atp = UpperTriangular([1 1 2; 0 2 2; 0 0 3])
@test sqrt(Atn) |> t->t*t ≈ Atn
@test typeof(sqrt(Atn)[1,1]) <: Complex
@test sqrt(Atp) |> t->t*t ≈ Atp
@test typeof(sqrt(Atp)[1,1]) <: Real

Areal   = randn(n, n)/2
Aimg    = randn(n, n)/2
A2real  = randn(n, n)/2
A2img   = randn(n, n)/2

for eltya in (Float32, Float64, ComplexF32, ComplexF64, BigFloat, Int)
    A = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(Areal, Aimg) : Areal)
    # a2 = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(a2real, a2img) : a2real)
    εa = eps(abs(float(one(eltya))))

    for eltyb in (Float32, Float64, ComplexF32, ComplexF64)
        εb = eps(abs(float(one(eltyb))))
        ε = max(εa,εb)

        debug && println("\ntype of A: ", eltya, " type of b: ", eltyb, "\n")

        debug && println("Solve upper triangular system")
        Atri = UpperTriangular(lu(A).U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
        b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
        x = Matrix(Atri) \ b

        debug && println("Test error estimates")
        if eltya != BigFloat && eltyb != BigFloat
            for i = 1:2
                @test  norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
            end
        end
        debug && println("Test forward error [JIN 5705] if this is not a BigFloat")

        x = Atri \ b
        γ = n*ε/(1 - n*ε)
        if eltya != BigFloat
            bigA = big.(Atri)
            x̂ = fill(1., n, 2)
            for i = 1:size(b, 2)
                @test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
            end
        end

        debug && println("Test backward error [JIN 5705]")
        for i = 1:size(b, 2)
            @test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
        end

        debug && println("Solve lower triangular system")
        Atri = UpperTriangular(lu(A).U) |> t -> eltya <: Complex && eltyb <: Real ? real(t) : t # Here the triangular matrix can't be too badly conditioned
        b = convert(Matrix{eltyb}, Matrix(Atri)*fill(1., n, 2))
        x = Matrix(Atri)\b

        debug && println("Test error estimates")
        if eltya != BigFloat && eltyb != BigFloat
            for i = 1:2
                @test  norm(x[:,1] .- 1) <= errorbounds(UpperTriangular(A), x, b)[1][i]
            end
        end

        debug && println("Test forward error [JIN 5705] if this is not a BigFloat")
        b = (b0 = Atri*fill(1, n, 2); convert(Matrix{eltyb}, eltyb == Int ? trunc.(b0) : b0))
        x = Atri \ b
        γ = n*ε/(1 - n*ε)
        if eltya != BigFloat
            bigA = big.(Atri)
            x̂ = fill(1., n, 2)
            for i = 1:size(b, 2)
                @test norm(x̂[:,i] - x[:,i], Inf)/norm(x̂[:,i], Inf) <= condskeel(bigA, x̂[:,i])*γ/(1 - condskeel(bigA)*γ)
            end
        end

        debug && println("Test backward error [JIN 5705]")
        for i = 1:size(b, 2)
            @test norm(abs.(b[:,i] - Atri*x[:,i]), Inf) <= γ * norm(Atri, Inf) * norm(x[:,i], Inf)
        end
    end
end

# Issue 10742 and similar
@test istril(UpperTriangular(diagm(0 => [1,2,3,4])))
@test istriu(LowerTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(UpperTriangular(diagm(0 => [1,2,3,4])))
@test isdiag(LowerTriangular(diagm(0 => [1,2,3,4])))
@test !isdiag(UpperTriangular(rand(4, 4)))
@test !isdiag(LowerTriangular(rand(4, 4)))

# Test throwing in fallbacks for non BlasFloat/BlasComplex in A_rdiv_Bx!
let n = 5
    A = rand(Float16, n, n)
    B = rand(Float16, n-1, n-1)
    @test_throws DimensionMismatch rdiv!(A, LowerTriangular(B))
    @test_throws DimensionMismatch rdiv!(A, UpperTriangular(B))
    @test_throws DimensionMismatch rdiv!(A, UnitLowerTriangular(B))
    @test_throws DimensionMismatch rdiv!(A, UnitUpperTriangular(B))

    @test_throws DimensionMismatch rdiv!(A, adjoint(LowerTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, adjoint(UpperTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, adjoint(UnitLowerTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, adjoint(UnitUpperTriangular(B)))

    @test_throws DimensionMismatch rdiv!(A, transpose(LowerTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, transpose(UpperTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, transpose(UnitLowerTriangular(B)))
    @test_throws DimensionMismatch rdiv!(A, transpose(UnitUpperTriangular(B)))
end

# Test that UpperTriangular(LowerTriangular) throws. See #16201
@test_throws ArgumentError LowerTriangular(UpperTriangular(randn(3,3)))
@test_throws ArgumentError UpperTriangular(LowerTriangular(randn(3,3)))

# Issue 16196
@test UpperTriangular(Matrix(1.0I, 3, 3)) \ view(fill(1., 3), [1,2,3]) == fill(1., 3)

# dimensional correctness:
const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :TestHelpers) || @eval Main include(joinpath($(BASE_TEST_PATH), "TestHelpers.jl"))
using .Main.TestHelpers: Furlong
let A = UpperTriangular([Furlong(1) Furlong(4); Furlong(0) Furlong(1)])
    @test sqrt(A) == Furlong{1//2}.(UpperTriangular([1 2; 0 1]))
end

@testset "similar should preserve underlying storage type" begin
    local m, n = 4, 3
    sparsemat = sprand(m, m, 0.5)
    for TriType in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular)
        trisparsemat = TriType(sparsemat)
        @test isa(similar(trisparsemat), typeof(trisparsemat))
        @test isa(similar(trisparsemat, Float32), TriType{Float32,<:SparseMatrixCSC{Float32}})
        @test isa(similar(trisparsemat, (n, n)), typeof(sparsemat))
        @test isa(similar(trisparsemat, Float32, (n, n)), SparseMatrixCSC{Float32})
    end
end

@testset "special printing of Lower/UpperTriangular" begin
    @test occursin(r"3×3 (LinearAlgebra\.)?LowerTriangular{Int64,Array{Int64,2}}:\n 2  ⋅  ⋅\n 2  2  ⋅\n 2  2  2",
                   sprint(show, MIME"text/plain"(), LowerTriangular(2ones(Int64,3,3))))
    @test occursin(r"3×3 (LinearAlgebra\.)?UnitLowerTriangular{Int64,Array{Int64,2}}:\n 1  ⋅  ⋅\n 2  1  ⋅\n 2  2  1",
                   sprint(show, MIME"text/plain"(), UnitLowerTriangular(2ones(Int64,3,3))))
    @test occursin(r"3×3 (LinearAlgebra\.)?UpperTriangular{Int64,Array{Int64,2}}:\n 2  2  2\n ⋅  2  2\n ⋅  ⋅  2",
                   sprint(show, MIME"text/plain"(), UpperTriangular(2ones(Int64,3,3))))
    @test occursin(r"3×3 (LinearAlgebra\.)?UnitUpperTriangular{Int64,Array{Int64,2}}:\n 1  2  2\n ⋅  1  2\n ⋅  ⋅  1",
                   sprint(show, MIME"text/plain"(), UnitUpperTriangular(2ones(Int64,3,3))))
end

@testset "adjoint/transpose triangular/vector multiplication" begin
    for elty in (Float64, ComplexF64), trity in (UpperTriangular, LowerTriangular)
        A1 = trity(rand(elty, 1, 1))
        b1 = rand(elty, 1)
        A4 = trity(rand(elty, 4, 4))
        b4 = rand(elty, 4)
        @test A1 * b1' ≈ Matrix(A1) * b1'
        @test_throws DimensionMismatch A4 * b4'
        @test A1 * transpose(b1) ≈ Matrix(A1) * transpose(b1)
        @test_throws DimensionMismatch A4 * transpose(b4)
        @test A1' * b1' ≈ Matrix(A1') * b1'
        @test_throws DimensionMismatch A4' * b4'
        @test A1' * transpose(b1) ≈  Matrix(A1') * transpose(b1)
        @test_throws DimensionMismatch A4' * transpose(b4)
        @test transpose(A1) * transpose(b1) ≈  Matrix(transpose(A1)) * transpose(b1)
        @test_throws DimensionMismatch transpose(A4) * transpose(b4)
        @test transpose(A1) * b1' ≈ Matrix(transpose(A1)) * b1'
        @test_throws DimensionMismatch transpose(A4) * b4'
        @test b1' * transpose(A1) ≈ b1' * Matrix(transpose(A1))
        @test b4' * transpose(A4) ≈ b4' * Matrix(transpose(A4))
        @test transpose(b1) * A1' ≈ transpose(b1) * Matrix(A1')
        @test transpose(b4) * A4' ≈ transpose(b4) * Matrix(A4')
    end
end

end # module TestTriangular

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
@testset "closure" begin
    magic(x) = false
    sentinel = gensym("sentinel")
    @test magic(sentinel) == false

    # Getting closers to work means having a function created in the current scope
    patch = @patch magic(x) = x == sentinel
    apply(patch) do
        @test (@mock magic(sentinel)) == true
    end
end

# https://github.com/invenia/Mocking.jl/issues/56
@testset "kwarg default uses arg" begin
    function f end

    p = @patch function f(x; y=x)
        (x, x, y, y)
    end

    apply(p) do
        @test (@mock f(1)) == (1, 1, 1, 1)
        @test (@mock f(1, y=2)) == (1, 1, 2, 2)
    end
end

# https://github.com/invenia/Mocking.jl/issues/57
@testset "kwarg default uses symbol" begin
    function f end

    p = @patch function f(x; y=:foo)
        (x, x, y, y)
    end

    apply(p) do
        @test (@mock f(1)) == (1, 1, :foo, :foo)
        @test (@mock f(1, y=2)) == (1, 1, 2, 2)
    end
end

ext: jl


text:
#############################################################################
# ECOS.jl
# Wrapper around the ECOS solver https://github.com/ifa-ethz/ecos
# See http://github.com/jump-dev/ECOS.jl
#############################################################################
# ECOS.jl
# Contains the wrapper itself
#############################################################################

module ECOS

import ECOS_jll

using CEnum

include("gen/libecos_common.jl")
include("gen/libecos_api.jl")

"""
    unsafe_add_settings(p::Ptr{pwork}, options::Dict{Symbol})

Add settings to the ECOS model.
"""
function unsafe_add_settings(p::Ptr{pwork}, options::Dict{Symbol})
    problem = unsafe_load(p)::pwork
    stgs = unsafe_load(problem.stgs)::settings
    new_stgs = settings(
        [get(options, k, getfield(stgs, k)) for k in fieldnames(settings)]...,
    )
    unsafe_store!(problem.stgs, new_stgs)
    return
end

include("MOI_wrapper/MOI_wrapper.jl")

function __init__()
    global ecos = ECOS_jll.libecos
    libecos_version = VersionNumber(unsafe_string(ECOS_ver()))
    if libecos_version != v"2.0.8"
        error(
            "Current ECOS version installed is $(libecos_version), but we " *
            "require version 2.0.8",
        )
    end
    return
end

end # module


text:
module CodeGeneration

using CodeManager

function module_listeners(mod::Module)
  expr = Expr(:toplevel,
    :(using CodeManager),
    :(_module_listeners = Dict{Symbol, Module}()),
    :(addListener(mod::Module) = if !haskey(_module_listeners, mod) _module_listeners[Symbol(mod)] = mod end),
    :(callOnListeners(name::Symbol) = for (_, listener) in _module_listeners CodeManager.safe_call(listener, name) end),
    )
  Core.eval(mod,expr)
end
export module_listeners

end #CodeGeneration

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module Repr
using InteractiveUtils

import Term:
    truncate,
    escape_brackets,
    highlight,
    do_by_line,
    unescape_brackets,
    split_lines,
    term_theme

import ..Layout: vLine, rvstack, lvstack, Spacer, vstack, cvstack, hLine, pad
import ..Renderables: RenderableText, info, AbstractRenderable
import ..Consoles: console_width
import ..Panels: Panel, TextBox
import ..Style: apply_style
import ..Tprint: tprint
import ..Tables: Table
import ..TermMarkdown: parse_md

export @with_repr, termshow, install_term_repr

include("_repr.jl")
include("_inspect.jl")

function termshow(io::IO, obj)
    return print(
        io,
        Panel(
            repr_get_obj_fields_display(obj);
            fit = true,
            subtitle = escape_brackets(string(typeof(obj))),
            subtitle_justify = :right,
            width = 40,
            justify = :center,
            style = term_theme[].repr_panel_style,
            subtitle_style = term_theme[].repr_name_style,
        ),
    )
end

termshow(obj) = termshow(stdout, obj)

# ---------------------------------------------------------------------------- #
#                                     EXPR                                     #
# ---------------------------------------------------------------------------- #
function termshow(io::IO, e::Expr)
    content = repr_get_obj_fields_display(e)
    content =
        cvstack("{green}$(highlight(string(e))){/green}", hLine(content.measure.w), content)
    return print(
        io,
        Panel(
            content;
            fit = true,
            subtitle = escape_brackets(string(typeof(e))),
            subtitle_justify = :right,
            width = 40,
            justify = :center,
            style = term_theme[].repr_panel_style,
            subtitle_style = term_theme[].repr_name_style,
        ),
    )
end

# ---------------------------------------------------------------------------- #
#                                  DICTIONARY                                  #
# ---------------------------------------------------------------------------- #
function termshow(io::IO, obj::AbstractDict)
    short_string(x) = truncate(string(x), 30)
    # prepare text renderables
    k =
        RenderableText.(
            short_string.(keys(obj));
            style = term_theme[].repr_accent_style * " bold",
        )
    ktypes =
        RenderableText.(
            map(k -> "{{" * short_string(typeof(k)) * "}}", collect(keys(obj)));
            style = term_theme[].repr_type_style * " dim",
        )
    vals =
        RenderableText.(
            short_string.(values(obj));
            style = term_theme[].repr_values_style * " bold",
        )
    vtypes =
        RenderableText.(
            map(k -> "{{" * short_string(typeof(k)) * "}}", collect(values(obj)));
            style = term_theme[].repr_type_style * " dim",
        )

    # trim if too many
    arrows = if length(k) > 10
        k, ktypes, vals, vtypes = k[1:10], ktypes[1:10], vals[1:10], vtypes[1:10]

        push!(k, RenderableText("⋮"; style = term_theme[].repr_accent_style))
        push!(ktypes, RenderableText("⋮"; style = term_theme[].repr_type_style * " dim"))
        push!(vals, RenderableText("⋮"; style = term_theme[].repr_values_style))
        push!(vtypes, RenderableText("⋮"; style = term_theme[].repr_type_style * " dim"))

        vstack(RenderableText.(repeat(["=>"], length(k) - 1); style = "red bold")...)
    else
        vstack(RenderableText.(repeat(["=>"], length(k)); style = "red bold")...)
    end

    # prepare other renderables
    space = Spacer(1, length(k))
    line = vLine(length(k); style = "dim #7e9dd9")

    _keys_renderables = cvstack(ktypes...) * line * space * cvstack(k...)
    _values_renderables = cvstack(vals...) * space * line * cvstack(vtypes...)

    return print(
        io,
        Panel(
            _keys_renderables * space * arrows * space * _values_renderables;
            fit = true,
            title = escape_brackets(string(typeof(obj))),
            title_justify = :left,
            width = 40,
            justify = :center,
            style = term_theme[].repr_panel_style,
            title_style = term_theme[].repr_name_style,
            padding = (2, 2, 1, 1),
            subtitle = "{bold white}$(length(keys(obj))){/bold white}{default} items{/default}",
            subtitle_justify = :right,
        ),
    )
end

# ---------------------------------------------------------------------------- #
#                                ABSTRACT ARRAYS                               #
# ---------------------------------------------------------------------------- #
termshow(io::IO, mtx::AbstractMatrix) = print(
    io,
    repr_panel(
        mtx,
        matrix2content(mtx),
        "{bold white}$(size(mtx, 1)) × $(size(mtx, 2)){/bold white}{default} {/default}",
    ),
)

termshow(io::IO, vec::Union{Tuple,AbstractVector}) = print(
    io,
    repr_panel(
        vec,
        vec2content(vec),
        "{bold white}$(length(vec)){/bold white}{default} items{/default}";
        justify = :left,
    ),
)

function termshow(io::IO, arr::AbstractArray)
    I0 = CartesianIndices(size(arr)[3:end])
    I = I0[1:min(10, length(I0))]

    panels::Vector{Union{Panel,Spacer}} = []
    for (n, i) in enumerate(I)
        i_string = join([string(i) for i in Tuple(i)], ", ")
        push!(
            panels,
            Panel(
                matrix2content(arr[:, :, i]; max_w = 60, max_items = 25, max_D = 5);
                subtitle = "[:, :, $i_string]",
                subtitle_justify = :right,
                width = 22,
                style = "dim yellow",
                subtitle_style = "default",
                title = "($n)",
                title_style = "dim bright_blue",
            ),
        )
        push!(panels, Spacer(1, 2))
    end

    if length(I0) > length(I)
        push!(
            panels,
            Panel(
                "{dim bright_blue bold underline}$(length(I0) - length(I)){/dim bright_blue bold underline}{dim bright_blue} frames omitted{/dim bright_blue}";
                width = panels[end - 1].measure.w,
                style = "dim yellow",
            ),
        )
    end

    return print(
        io,
        repr_panel(
            arr,
            vstack(panels...),
            "{white}" * join(string.(size(arr)), " × ") * "{/white}",
        ),
    )
end

function termshow(io::IO, obj::DataType)
    ts = term_theme[].repr_type_style
    field_names = apply_style.(string.(fieldnames(obj)), term_theme[].repr_accent_style)
    field_types = apply_style.(map(f -> "::" * string(f), obj.types), ts)

    line = vLine(length(field_names); style = term_theme[].repr_name_style)
    space = Spacer(1, length(field_names))
    fields = rvstack(field_names...) * space * lvstack(string.(field_types)...)

    type_name = apply_style(string(obj), term_theme[].repr_name_style * " bold")
    sup = supertypes(obj)[2]
    type_name *= " {bright_blue dim}<: $sup{/bright_blue dim}"
    content =
        "    " * repr_panel(
            nothing,
            string(type_name / ("  " * line * fields)),
            nothing;
            fit = false,
            width = min(console_width() - 5, 80),
            justify = :center,
        )

    # get docstring
    doc, _ = get_docstring(obj)
    doc = parse_md(doc; width = min(100, console_width()))
    doc = split_lines(doc)
    if length(doc) > 100
        doc = [
            doc[1:min(100, length(doc))]...,
            "{dim bright_blue}$(length(doc)-100) lines omitted...{/dim bright_blue}",
        ]
    end
    doc = join(doc, "\n")

    print(io, content)
    print(io, hLine(console_width(), "Docstring"; style = "green dim", box = :HEAVY))
    tprint(io, doc)
end

function termshow(io::IO, fun::Function)
    # get methods
    _methods = split_lines(string(methods(fun)))
    N = length(_methods)

    _methods = length(_methods) > 1 ? _methods[2:min(11, N)] : []
    _methods = map(m -> join(split(join(split(m, "]")[2:end]), " in ")[1]), _methods)
    _methods = map(
        m -> replace(
            m,
            string(fun) => "{bold #a5c6d9}$(string(fun)){/bold #a5c6d9}";
            count = 1,
        ),
        _methods,
    )
    counts = RenderableText.("(" .* string.(1:length(_methods)) .* ") "; style = "bold dim")
    length(_methods) < N - 1 && push!(
        _methods,
        "\n{bold dim bright_blue}$(N - length(_methods)-1){/bold dim bright_blue}{dim bright_blue} methods omitted...{/dim bright_blue}",
    )
    methods_contents = if N > 1
        rvstack(counts...) * lvstack(RenderableText.(highlight.(_methods))...)
    else
        split_lines(string(methods(fun)))[1]
    end

    panel =
        "       " * repr_panel(
            nothing,
            methods_contents,
            "{white bold}$(N-1){/white bold} methods",
            title = "Function: {bold bright_blue}$(string(fun)){/bold bright_blue}",
            width = console_width() - 12,
        )

    # get docstring
    doc, _ = get_docstring(fun)
    doc = parse_md(doc; width = console_width() - 8)
    doc = split_lines(doc)
    if length(doc) > 100
        doc = [
            doc[1:min(100, length(doc))]...,
            "{dim bright_blue}$(length(doc)-100) lines omitted...{/dim bright_blue}",
        ]
    end
    print(io, panel)
    print(io, hLine(console_width(), "Docstring"; style = "green dim", box = :HEAVY))
    print(io, "   " * RenderableText(join(doc, "\n"); width = console_width() - 8))
end

# ---------------------------------------------------------------------------- #
#                                 INSTALL REPR                                 #
# ---------------------------------------------------------------------------- #
function install_term_repr()
    @eval begin
        import Term: termshow

        Base.show(io::IO, ::MIME"text/plain", num::Number) =
            tprint(io, string(num); highlight = true)

        Base.show(io::IO, num::Number) = tprint(io, string(num); highlight = true)

        Base.show(io::IO, ::MIME"text/plain", obj::AbstractDict) = termshow(io, obj)

        Base.show(io::IO, ::MIME"text/plain", obj::Union{AbstractArray,AbstractMatrix}) =
            termshow(io, obj)

        Base.show(io::IO, ::MIME"text/plain", fun::Function) = termshow(io, fun)

        Base.show(io::IO, ::MIME"text/plain", obj::DataType) = termshow(io, obj)

        Base.show(io::IO, ::MIME"text/plain", expr::Expr) = termshow(io, expr)
    end
end

# ---------------------------------------------------------------------------- #
#                                   WITH REPR                                  #
# ---------------------------------------------------------------------------- #

"""
    with_repr(typedef::Expr)

Function for the macro @with_repr which creates a `Base.show` method for a type.

The `show` method shows the field names/types for the
type and the values of the fields.

# Example
```
@with_repr struct TestStruct2
    x::Int
    name::String
    y
end
```
"""
function with_repr(typedef::Expr)
    tn = typename(typedef) # the name of the type
    showfn = :(Base.show(io::IO, ::MIME"text/plain", obj::$tn) = termshow(io, obj))

    quote
        $typedef
        $showfn
    end
end

"""
with_repr(typedef::Expr)

Function for the macro @with_repr which creates a `Base.show` method for a type.

The `show` method shows the field names/types for the
type and the values of the fields.

# Example
```
@with_repr struct TestStruct2
x::Int
name::String
y
end
```
"""
macro with_repr(typedef)
    return esc(with_repr(typedef))
end

end


text:
# example script for running ChhabraJensen (direct) method

using DataFrames, CSV, GLM, Plots, Plots.PlotMeasures, Polynomials, Polynomials.PolyCompat

include("mfdfa1.jl")

dataFile = DataFrame(CSV.File("exampleTS.csv", header = false))

dataSeries = dataFile[:,1] |> vec

# check for negative points
minimum(dataSeries)

# parameters
q = LinRange(-5,5,101)
exponents = LinRange(4,10,19)
scale = round.(2 .^exponents)
m = 1

mf_results = MFDFA(dataSeries, scale, q, m)

# Plots

# mf-spectrum plot
mf_spectrum_plot = plot(mf_results["hq"], mf_results["Dq"],seriestype=:scatter, legend=false, xlabel = "hq", ylabel = "Dq", title=string("mfw: ", round(mf_results["mfw"];digits=3)))

# q-order plot
Hq_q_plot = plot(mf_results["q"], mf_results["Hq"], seriestype=:scatter, legend=false, xlabel = "q", ylabel = "Hq")

# combined plots
plot(Hq_q_plot, mf_spectrum_plot, size = (1000,500), margin = 25px)


text:
"""
    CoulombTransformation{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}

Coulomb transformation tensor with even spacing for t-points.

If `k` is zero this tensor solves Poisson equation. If nonzero it works as
Helmholz equation Greens function.

# Fields
- `elementgrid::SMatrix{NG,NE}` : Gauss points in elements
- `t::SVector{NT}{Float64}`     : t-coordinates
- `wt::Vector{T}`               : Weights for t-integration
- `w::SVector{NG}{Float64}`     : Weights for position coordinate integration
- `tmin::Float64`               : t-integration lower limit  (default=0.)
- `tmax::Float64`               : t-integration upper limit  (default=25.)
- `k::T`                        : Helmholz equation constant (default=0.)

# Contruction
    CoulombTransformation(ceg::CubicElementGrid, nt::Int; tmin=0., tmax=25., k=0.)

# Example
```jldoctest
julia> ceg = CubicElementGrid(5, 4, 16)
Cubic elements grid with 4^3 elements with 16^3 Gauss points

julia> CoulombTransformation(ceg, 48)
Coulomb transformation tensor for 48 t-points, tmin=0.0 tmax=25.0 and k=0.0

julia> ct = CoulombTransformation(ceg, 48; tmin=4., tmax=30., k=0.1)
Coulomb transformation tensor for 48 t-points, tmin=4.0 tmax=30.0 and k=0.1

julia> typeof(ct) <: AbstractArray{Float64,5}
true

julia> ct[3,3,1,1,5] ≈ 0.0158700408
true
```
"""
struct CoulombTransformation{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}
    elementgrid::SMatrix{NG,NE}
    t::SVector{NT}{Float64}
    wt::Vector{T}
    w::SVector{NG}{Float64}
    tmin::Float64
    tmax::Float64
    k::T
    function CoulombTransformation(ceg::CubicElementGrid, nt::Int; tmin=0., tmax=25., k=0.)
        t, wt = gausspoints(nt; elementsize=(tmin, tmax))
        grid = grid1d(ceg)
        s = size(grid)
        wt = wt .* exp.(-(k^2)./(4t.^2))
        new{eltype(wt), nt, s[2], s[1]}(grid, t, wt, ceg.w, tmin, tmax, k)
    end
end

"""
    CoulombTransformationLog{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}

Coulomb transformation tensor with logarithmic spacing for t-points.

If `k` is zero this tensor solves Poisson equation. If nonzero it works as
Helmholz equation Greens function.

# Fields
- `elementgrid::SMatrix{NG,NE}` : Gauss points in elements
- `t::SVector{NT}{Float64}`     : t-coordinates
- `wt::Vector{T}`               : Weights for t-integration
- `w::SVector{NG}{Float64}`     : Weights for position coordinate integration
- `tmin::Float64`               : t-integration lower limit
- `tmax::Float64`               : t-integration upper limit
- `k::T`                        : Helmholz equation constant

# Construction
    CoulombTransformationLog(ceg::CubicElementGrid, nt::Int; kwargs...)

## Keywords
- `tmin=0.`
- `tmax=25.`
- `k=0.`

# Example
```jldoctest
julia> ceg = CubicElementGrid(5, 4, 16)
Cubic elements grid with 4^3 elements with 16^3 Gauss points

julia> CoulombTransformationLog(ceg, 48)
Logarithmic Coulomb transformation tensor for 48 t-points, tmin=0.0, tmax=25.0 and k=0.0

julia> ct = CoulombTransformationLog(ceg, 48; tmin=4., tmax=30., k=0.1)
Logarithmic Coulomb transformation tensor for 48 t-points, tmin=4.0, tmax=30.0 and k=0.1

julia> typeof(ct) <: AbstractArray{Float64,5}
true

julia> ct[3,3,1,1,5] ≈ 0.0594740698
true
```
"""
struct CoulombTransformationLog{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}
    elementgrid::SMatrix{NG,NE}
    t::SVector{NT}{Float64}
    wt::Vector{T}
    w::SVector{NG}{Float64}
    tmin::Float64
    tmax::Float64
    k::T
    function CoulombTransformationLog(ceg::CubicElementGrid, nt::Int;
                                      tmin=0., tmax=25., k=0.)
        s, ws = gausspoints(nt; elementsize=(log(tmin+1e-12), log(tmax)))
        t = exp.(s)
        wt = ws .* t
        grid = grid1d(ceg)
        ss = size(grid)
        wt = wt .* exp.(-(k^2)./(4t.^2))
        new{eltype(wt), nt, ss[2], ss[1]}(grid, t, wt, ceg.w, tmin, tmax, k)
    end
end

"""
    CoulombTransformationLocal{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}

Coulomb transformation tensor with even spacing for t-points and t-integration
local average correction.

If `k` is zero this tensor solves Poisson equation. If nonzero it works as
Helmholz equation Greens function.

# Fields
- `elementgrid::SMatrix{NG,NE}` : Gauss points in elements
- `t::SVector{NT}{Float64}`     : t-coordinates
- `wt::Vector{T}`               : Weights for t-integration
- `w::SVector{NG}{Float64}`     : Weights for position coordinate integration
- `tmin::Float64`               : t-integration lower limit
- `tmax::Float64`               : t-integration upper limit
- `k::T`                        : Helmholz equation constant
- `δ::Float64`                  : Correction constant for t-integration ∈ ]0,1]
- `δp::SMatrix{NG,NE}`          : Upper average points for t-integration correction
- `δm::SMatrix{NG,NE}`          : Lower average points for t-integration correction

# Construction
    CoulombTransformationLocal(ceg::CubicElementGrid, nt::Int; kwargs...)

## Keywords
- `tmin=0.`
- `tmax=25.`
- `δ=0.25`
- `k=0.`

# Example
```jldoctest
julia> ceg = CubicElementGrid(5, 4, 16)
Cubic elements grid with 4^3 elements with 16^3 Gauss points

julia> CoulombTransformationLocal(ceg, 48)
Coulomb transformation tensor with local correction for 48 t-points, tmin=0.0, tmax=25.0 and k=0.0

julia> ct = CoulombTransformationLocal(ceg, 48; tmin=4., tmax=30., k=0.1)
Coulomb transformation tensor with local correction for 48 t-points, tmin=4.0, tmax=30.0 and k=0.1

julia> typeof(ct) <: AbstractArray{Float64,5}
true

julia> ct[3,3,1,1,5] ≈ 0.0591078365
true
```
"""
struct CoulombTransformationLocal{T, NT, NE, NG} <: AbstractCoulombTransformationLocal{NT, NE, NG}
    elementgrid::SMatrix{NG,NE}
    t::SVector{NT}{Float64}
    wt::Vector{T}
    w::SVector{NG}{Float64}
    tmin::Float64
    tmax::Float64
    k::T
    δ::Float64
    δp::SMatrix{NG,NE}
    δm::SMatrix{NG,NE}
    function CoulombTransformationLocal(ceg::CubicElementGrid, nt::Int;
                                   tmin=0., tmax=25., δ=0.25, k=0.)
        @assert 0 < δ <= 1
        @assert 0 <= tmin < tmax
        t, wt = gausspoints(nt; elementsize=(tmin, tmax))
        grid = grid1d(ceg)
        s = size(grid)
        δp = zeros(s...)
        δm = zeros(s...)
        δm[2:end,:] = grid[1:end-1,:] .- grid[2:end,:]
        δp[1:end-1,:] = grid[2:end,:] .- grid[1:end-1,:]
        for j in axes(grid, 2)
            δm[1,j] =  austrip(get_1d_element(ceg.elements, j).low) - grid[1,j]
            δp[end,j] = austrip(get_1d_element(ceg.elements, j).high) - grid[end,j]
        end
        wt = wt .* exp.(-(k^2)./(4t.^2))
        new{eltype(wt), nt, s[2], s[1]}(grid, t, wt, ceg.w, tmin, tmax, k, δ, δ.*δp, δ.*δm)
    end
end

"""
    CoulombTransformationLogLocal{T, NT, NE, NG} <: AbstractCoulombTransformationSingle{NT, NE, NG}

Coulomb transformation tensor with logarithmic spacing for t-points and
t-integration local average correction.

If `k` is zero this tensor solves Poisson equation. If nonzero it works as
Helmholz equation Greens function.

# Fields
- `elementgrid::SMatrix{NG,NE}` : Gauss points in elements
- `t::SVector{NT}{Float64}`     : t-coordinates
- `wt::Vector{T}`               : Weights for t-integration
- `w::SVector{NG}{Float64}`     : Weights for position coordinate integration
- `tmin::Float64`               : t-integration lower limit
- `tmax::Float64`               : t-integration upper limit
- `k::T`                        : Helmholz equation constant
- `δ::Float64`                  : Correction constant for t-integration ∈ ]0,1]
- `δp::SMatrix{NG,NE}`          : Upper average points for t-integration correction
- `δm::SMatrix{NG,NE}`          : Lower average points for t-integration correction

# Construction
    CoulombTransformationLogLocal(ceg::CubicElementGrid, nt::Int; kwargs...)

## Keywords
- `tmin=0.`
- `tmax=25.`
- `δ=0.25`
- `k=0.`

# Example
```jldoctest
julia> ceg = CubicElementGrid(5, 4, 16)
Cubic elements grid with 4^3 elements with 16^3 Gauss points

julia> CoulombTransformationLogLocal(ceg, 48)
Logarithmic Coulomb transformation tensor with local correction for 48 t-points, tmin=0.0, tmax=25.0 and k=0.0

julia> ct = CoulombTransformationLogLocal(ceg, 48; tmin=4., tmax=30., k=0.1)
Logarithmic Coulomb transformation tensor with local correction for 48 t-points, tmin=4.0, tmax=30.0 and k=0.1

julia> typeof(ct) <: AbstractArray{Float64,5}
true

julia> ct[3,3,1,1,5] ≈ 0.0591078365
true
```
"""
struct CoulombTransformationLogLocal{T, NT, NE, NG} <: AbstractCoulombTransformationLocal{NT, NE, NG}
    elementgrid::SMatrix{NG,NE}
    t::SVector{NT}{Float64}
    wt::Vector{T}
    w::SVector{NG}{Float64}
    tmin::Float64
    tmax::Float64
    k::T
    δ::Float64
    δp::SMatrix{NG,NE}
    δm::SMatrix{NG,NE}
    function CoulombTransformationLogLocal(ceg::CubicElementGrid, nt::Int;
                                   tmin=0., tmax=25., δ=0.25, k=0.)
        @assert 0 < δ <= 1
        @assert 0 <= tmin < tmax
        s, ws = gausspoints(nt; elementsize=(log(tmin+1e-12), log(tmax)))
        t = exp.(s)
        wt = ws .* t
        t, wt = gausspoints(nt; elementsize=(tmin, tmax))
        grid = grid1d(ceg)
        ss = size(grid)
        δp = zeros(ss...)
        δm = zeros(ss...)
        δm[2:end,:] = grid[1:end-1,:] .- grid[2:end,:]
        δp[1:end-1,:] = grid[2:end,:] .- grid[1:end-1,:]
        for j in 1:ss[2]
            δm[1,j] =  austrip(get_1d_element(ceg.elements, j).low) - grid[1,j]
            δp[end,j] = austrip(get_1d_element(ceg.elements, j).high) - grid[end,j]
        end
        wt = wt .* exp.(-(k^2)./(4t.^2))
        new{eltype(wt), nt, ss[2], ss[1]}(grid, t, wt, ceg.w, tmin, tmax, k, δ, δ.*δp, δ.*δm)
    end
end

"""
    CoulombTransformationCombination{T, NE, NG} <: AbstractCoulombTransformationCombination

Coulomb transformation tensor that is composed of different types of tensors
with different t-integration ranges.

This tensor is intended to to be used to control how t-integration is performed.

# Fields
- `tensors::Vector{AbstractCoulombTransformation}` : Coulomb transformation subtensors
- `ref::Vector{Int}`  : index of tensor wich is referenced at given t-point index
- `t::Vector{Float64}`  : t-integration points
- `wt::Vector{T}`  : weights for t-integration
- `tmin::Float64`               : t-integration lower limit
- `tmax::Float64`               : t-integration upper limit
- `tindex::Vector{Int}`         : Index of t in subtensor
- `k::T`                        : Helmholz equation constant

# Construction
    CoulombTransformationCombination(t::AbstractCoulombTransformation...)

# Example
```jldoctest
julia> ceg = CubicElementGrid(5, 4, 16)
Cubic elements grid with 4^3 elements with 16^3 Gauss points

julia> ct1 = CoulombTransformationLog(ceg, 48)
Logarithmic Coulomb transformation tensor for 48 t-points, tmin=0.0, tmax=25.0 and k=0.0

julia> ct2 = CoulombTransformationLogLocal(ceg, 48; tmin=25, tmax=50)
Logarithmic Coulomb transformation tensor with local correction for 48 t-points, tmin=4.0, tmax=30.0 and k=0.1

julia> cc = CoulombTransformationCombination(ct1)
Coulomb transformation combination tensor for 48 t-points, tmin=0.0, tmax=25.0 and k=0.0

julia> push!(cc, ct2)
Coulomb transformation combination tensor for 96 t-points, tmin=0.0, tmax=50.0 and k=0.0

julia> CoulombTransformationCombination(ct1, ct2)
Coulomb transformation combination tensor for 96 t-points, tmin=0.0, tmax=50.0 and k=0.0
```
"""
mutable struct CoulombTransformationCombination{T, NE, NG} <: AbstractCoulombTransformationCombination
    tensors::Vector{AbstractCoulombTransformation}
    ref::Vector{Int}
    t::Vector{Float64}
    wt::Vector{T}
    tmin::Float64
    tmax::Float64
    tindex::Vector{Int}
    k::T
    function CoulombTransformationCombination(ct::AbstractCoulombTransformationSingle)
        s = size(ct)
        new{eltype(ct.wt), s[3], s[2]}([ct], ones(s[end]), ct.t, ct.wt, ct.tmin, ct.tmax, 1:length(ct.t), ct.k)
    end
end

function CoulombTransformationCombination(t::AbstractCoulombTransformationSingle...)
    ct = CoulombTransformationCombination(t[1])
    for x in t[2:end]
        push!(ct,x)
    end
    return ct
end

function loglocalct(ceg::CubicElementGrid, nt::Int; δ=0.25, tmax=300, tboundary=20, k=0.)
    @warn "loglocalct is depricated use optimal_coulomb_tranformation"
    @assert 0 < tboundary < tmax
    s, ws = gausspoints(nt; elementsize=(log(1e-12), log(tmax)))
    t = exp.(s)
    l = length(t[t .< tboundary])
    ct1 = CoulombTransformationLog(ceg, l; tmax=tboundary, k=k)
    ct2 = CoulombTransformationLogLocal(ceg, nt-l; tmin=tboundary, tmax=tmax, δ=δ, k=k)
    return CoulombTransformationCombination(ct1,ct2)
end

"""
    optimal_coulomb_tranformation(args...; kwargs...)

Returns Coulomb transformation tensor with optimal pre tested parameters.

# Args
- `ceg::CubicElementGrid`  :  grid where transformation is done
- `nt::Int`                :  number of t-points

# Keywords
- `δ=0.25`          : parameter for local average correction
- `tmax=700`        : maximum t value
- `tboundary=20`    : point after which local average correction is applied
- `k=0.`            : constant for Helmholz equation
"""
function optimal_coulomb_tranformation(ceg::CubicElementGrid, nt::Int; δ=0.25, tmax=700, tboundary=20, k=0.)
    @assert 0 < tboundary < tmax
    s, ws = gausspoints(nt; elementsize=(log(1e-12), log(tmax)))
    t = exp.(s)
    l = length(t[t .< tboundary])
    ct1 = CoulombTransformationLog(ceg, l; tmax=tboundary, k=k)
    ct2 = CoulombTransformationLogLocal(ceg, nt-l; tmin=tboundary, tmax=tmax, δ=δ, k=k)
    return CoulombTransformationCombination(ct1,ct2)
end

(ct::AbstractCoulombTransformation)(r,t) = exp(-(t*r)^2)
(ct::AbstractCoulombTransformationSingle)(r,p::Int) = exp(-ct.t[p]^2*r^2)

function (ct::AbstractCoulombTransformationLocal)(r,t)
    rmax = (r+ct.δ)*t
    rmin = (r-ct.δ)*t
    return 0.5*√π*erf(rmin, rmax)/(rmax-rmin)
end

function (ct::AbstractCoulombTransformationLocal)(r,p::Int)
    rmax = (r+ct.δ)*ct.t[p]
    rmin = (r-ct.δ)*ct.t[p]
    return 0.5*√π*erf(rmin, rmax)/(rmax-rmin)
end

function (ct::CoulombTransformationCombination)(r,p::Int)
    i = ct.ref[p]
    return ct.tensors[i](r, ct.t[p])
end

function Base.size(ct::AbstractCoulombTransformationSingle)
    ng, ne = size(ct.elementgrid)
    nt = length(ct.t)
    return (ng,ng,ne,ne,nt)
end

function Base.size(ct::CoulombTransformationCombination{T,NE,NG}) where {T,NE,NG}
    nt = length(ct.ref)
    return (NG,NG,NE,NE,nt)
end

function Base.getindex(ct::AbstractCoulombTransformationSingle, α::Int, β::Int, I::Int, J::Int, p::Int)
    r = ct.elementgrid[α,I] - ct.elementgrid[β,J]
    return ct.w[β]*ct(r, p)
end

function Base.getindex(ct::AbstractCoulombTransformationLocal, α::Int, β::Int, I::Int, J::Int, p::Int)
    r = abs(ct.elementgrid[α,I] - ct.elementgrid[β,J])
    δ = min(abs(ct.δp[α,I]-ct.δm[β,J]), abs(ct.δm[α,I]-ct.δp[β,J]))
    rmax = (r+δ)*ct.t[p]
    rmin = (r-δ)*ct.t[p]
    return ct.w[β]*0.5*√π*erf(rmin, rmax)/(rmax-rmin)
end

function Base.getindex(ct::CoulombTransformationCombination, α::Int, β::Int, I::Int, J::Int, p::Int)
    i = ct.ref[p]
    return ct.tensors[i][α, β, I, J, ct.tindex[p]]
end

function  Base.push!(ctc::CoulombTransformationCombination{T,NE,NG},
            tt::AbstractCoulombTransformationSingle{NT,NE,NG}) where {T,NT,NE,NG}
    @assert ctc.tensors[end].tmax <= tt.tmin
    @assert ctc.k == tt.k
    push!(ctc.tensors,tt)
    ctc.ref = vcat(ctc.ref, ones(NT).*length(ctc.tensors))
    ctc.t = vcat( ctc.t, tt.t)
    ctc.wt = vcat( ctc.wt, tt.wt)
    ctc.tmax = tt.tmax
    ctc.tindex = vcat( ctc.tindex, 1:length(tt.t))
    return ctc
end

function Base.show(io::IO, ::MIME"text/plain", cct::CoulombTransformation)
    s = size(cct)
    print(io, "Coulomb transformation tensor for $(s[end]) t-points,"*
              " tmin=$(cct.tmin), tmax=$(cct.tmax) and k=$(cct.k)"
    )
end

function Base.show(io::IO, ::MIME"text/plain", cct::CoulombTransformationLog)
    s = size(cct)
    print(io, "Logarithmic Coulomb transformation tensor for $(s[end]) t-points,"*
              " tmin=$(cct.tmin), tmax=$(cct.tmax) and k=$(cct.k)"
    )
end

function Base.show(io::IO, ::MIME"text/plain", cct::CoulombTransformationLocal)
    s = size(cct)
    print(io, "Coulomb transformation tensor with local correction for $(s[end])"*
              " t-points, tmin=$(cct.tmin), tmax=$(cct.tmax) and k=$(cct.k)"
    )
end

function Base.show(io::IO, ::MIME"text/plain", cct::CoulombTransformationCombination)
    s = size(cct)
    print(io, "Coulomb transformation combination tensor for $(s[end])"*
              " t-points, tmin=$(cct.tmin), tmax=$(cct.tmax) and k=$(cct.k)"
    )
end

function Base.show(io::IO, ::MIME"text/plain", cct::CoulombTransformationLogLocal)
    s = size(cct)
    print(io, "Logarithmic Coulomb transformation tensor with local correction for $(s[end])"*
              " t-points, tmin=$(cct.tmin), tmax=$(cct.tmax) and k=$(cct.k)"
    )
end

## Integration weight tensor

"""
    ω_tensor(ceg::CubicElementGrid) -> Matrix{Float64}

Return a tensor that contains integration constants to perform integral over
element grid.
"""
function ω_tensor(ceg::CubicElementGrid)
    n = size(ceg)[end]
    l = length(ceg.w)
    ω = reshape( repeat(ceg.w, n), l, n )
end

## Help functions to create charge density

function density_tensor(grid::AbstractArray, r::AbstractVector, a)
    return [ exp(-a*sum( (x-r).^2 )) for x in grid ]
end

density_tensor(grid; r=SVector(0.,0.,0.), a=1.) = density_tensor(grid, r, a)

## Derivative tensor

abstract type AbstractDerivativeTensor <: AbtractTransformationTensor{2} end

"""
    DerivativeTensor{NG} <: AbstractDerivativeTensor

Tensor that performs derivate over [`CubicElementGrid`](@ref).

# Fields
- `gauss_points::SVector{NG}{Float64}`  : point for integration
- `values::Matrix{Float64}`             : tensor values

# Construction
    DerivativeTensor(ceg::CubicElementGrid)
"""
struct DerivativeTensor{NG} <: AbstractDerivativeTensor
    gauss_points::SVector{NG}{Float64}
    values::Matrix{Float64}
    function DerivativeTensor(ceg::CubicElementGrid)
        function _lpoly(x)
            # Legendre polynomials (without prefactor)
            l = length(x)
            out = []
            for i in 1:l
                z = vcat(x[1:i-1],x[i+1:end])
                push!(out, fromroots(z))
            end
            return out
        end
        function _pre_a(x)
            # Legendre plynomial prefactors
            l = length(x)
            a = ones(l)
            for i in 1:l
                for j in 1:l
                    if j != i
                        a[i] *= x[i] - x[j]
                    end
                end
            end
            return a.^-1
        end
        # Derivative is done by analytically deriving Gauss-Legendre basis
        # functions and then calculating values over the derivatives.
        # f(x) = ∑aᵢpᵢ(x) → f'(x) = ∑aᵢpᵢ'(x)
        x = BigFloat.(ceg.gpoints)
        xmin=-0.5*elementsize(ceg.elements)
        xmax=-xmin
        grid = grid1d(ceg)
        ng, ne = size(grid)
        a=_pre_a(x)
        p=_lpoly(x)
        pd = derivative.(p)
        ϕ = zeros(length(x), length(pd))
        for i in eachindex(pd)
            ϕ[:,i] = a[i]*pd[i].(x)
        end
        new{ng}(x, ϕ)
    end
end

function Base.getindex(dt::DerivativeTensor, i::Int, j::Int)
    return dt.values[i,j]
end

Base.size(dt::DerivativeTensor) = size(dt.values)

function Base.show(io::IO, ::MIME"text/plain", dt::DerivativeTensor)
    print(io, "Derivative tensor size=$(size(dt))")
end

function kinetic_energy(ceg::CubicElementGrid, dt, ψ)
    @warn "kinetic_energy is deprecated use operators instead"
    tmp = similar(ψ)

    @tensor tmp[i,j,k,I,J,K] = dt[i,l] * ψ[l,j,k,I,J,K]
    ex = 0.5*integrate(tmp, ceg, tmp)

    @tensor tmp[i,j,k,I,J,K] = dt[j,l] * ψ[i,l,k,I,J,K]
    ey = 0.5*integrate(tmp, ceg, tmp)

    @tensor tmp[i,j,k,I,J,K] = dt[k,l] * ψ[i,j,l,I,J,K]
    ez = 0.5*integrate(tmp, ceg, tmp)

    return ex+ey+ez
end

## Momentum tensor

struct MomentumTensor{NG} <: AbstractArray{ComplexF64, 3}
    dt::DerivativeTensor{NG}
    function MomentumTensor(dt::DerivativeTensor{T}) where T
        @warn "MomentumTensor is deprecated use operators instead"
        new{T}(dt)
    end
end

function MomentumTensor(ceg::CubicElementGrid)
    return MomentumTensor(DerivativeTensor(ceg))
end

Base.size(mt::MomentumTensor) = (size(mt.dt)...,3)

Base.getindex(mt::MomentumTensor,i::Int,j::Int,xyz::Int) = -ComplexF64(0,mt.dt[i,j])

function Base.show(io::IO, ::MIME"text/plain", mt::MomentumTensor)
    print(io, "Momentum tensor size=$(size(mt))")
end

## Nuclear potential tensors

abstract type AbstractNuclearPotential <: AbtractTransformationTensor{3} end
abstract type AbstractNuclearPotentialSingle{T} <: AbstractNuclearPotential end
abstract type AbstractNuclearPotentialCombination{T} <: AbstractNuclearPotential end

"""
    NuclearPotentialTensor{T}

Tensor for nuclear potential in 1D. To get 3D you need 3 of these tensors and
integrate t-coordinate.

Easy way is to use [`PotentialTensor`](@ref) to get 3D tensor.

# Fields
- `elementgrid::Matrix{T}`   : Grid definition, collums are elements, rows Gauss points
- `t::Vector{Float64}`       : t-coordinate
- `wt::Vector{Float64}`      : Integration weights for t
- `tmin::Float64`            : Minimum t-value
- `tmax::Float64`            : Maximum t-value
- `r::T`                     : Nuclear coordinate

# Creation
    NuclearPotentialTensor(r, ceg::CubicElementGrid, nt::Int; tmin=0, tmax=30)

- `r`                       : Nuclear coordinate (1D)
- `ceg::CubicElementGrid`   : Grid for the potential
- `nt::Int`                 : Number of t-points
"""
struct NuclearPotentialTensor{T} <: AbstractNuclearPotentialSingle{T}
    elementgrid::Matrix{T}
    t::Vector{Float64}
    wt::Vector{Float64}
    tmin::Float64
    tmax::Float64
    r::T
    function NuclearPotentialTensor(r, ceg::CubicElementGrid, nt::Int; tmin=0, tmax=30)
        @assert dimension(r) == NoDims || dimension(r) == dimension(u"m")
        rt = austrip(r)
        t, wt = gausspoints(nt; elementsize=(tmin, tmax))
        grid = Array(grid1d(ceg)) .- rt
        new{eltype(grid)}(grid, t, wt, tmin, tmax, rt)
    end
end

"""
    NuclearPotentialTensorLog{T}

Tensor for nuclear potential in 1D. To get 3D you need 3 of these tensors and
integrate t-coordinate. This version has logarithmic spacing for t-points.

Easy way is to use [`PotentialTensor`](@ref) to get 3D tensor.

# Fields
- `elementgrid::Matrix{T}`   : Grid definition, collums are elements, rows Gauss points
- `t::Vector{Float64}`       : t-coordinate
- `wt::Vector{Float64}`      : Integration weights for t
- `tmin::Float64`            : Minimum t-value
- `tmax::Float64`            : Maximum t-value
- `r::T`                     : Nuclear coordinate

# Creation
    NuclearPotentialTensorLog(r, ceg::CubicElementGrid, nt::Int; tmin=0, tmax=30)

- `r`                       : Nuclear coordinate (1D)
- `ceg::CubicElementGrid`   : Grid for the potential
- `nt::Int`                 : Number of t-points
"""
struct NuclearPotentialTensorLog{T} <: AbstractNuclearPotentialSingle{T}
    elementgrid::Matrix{T}
    t::Vector{Float64}
    wt::Vector{Float64}
    tmin::Float64
    tmax::Float64
    r::T
    function NuclearPotentialTensorLog(r, ceg::CubicElementGrid, nt::Int; tmin=0, tmax=30)
        @assert dimension(r) == NoDims || dimension(r) == dimension(u"m")
        rt = austrip(r)
        s, ws = gausspoints(nt; elementsize=(log(tmin+1e-12), log(tmax)))
        t = exp.(s)
        wt = ws .* t
        grid = Array(grid1d(ceg)) .- rt
        new{eltype(grid)}(grid, t, wt, tmin, tmax, rt)
    end
end

"""
    NuclearPotentialTensorLogLocal{T}

Tensor for nuclear potential in 1D. To get 3D you need 3 of these tensors and
integrate t-coordinate. This version has logarithmic spacing for t-points and
local average correction.

Easy way is to use [`PotentialTensor`](@ref) to get 3D tensor.

# Fields
- `elementgrid::Matrix{T}`   : Grid definition, collums are elements, rows Gauss points
- `t::Vector{Float64}`       : t-coordinate
- `wt::Vector{Float64}`      : Integration weights for t
- `tmin::Float64`            : Minimum t-value
- `tmax::Float64`            : Maximum t-value
- `r::T`                     : Nuclear coordinate
- `δ::Float64`               : Correction parameter ∈[0,1]
- `δp::Matrix{T}`            : Upper limits for correction
- `δm::Matrix{T}`            : Lower limits for correction

# Creation
    NuclearPotentialTensorLogLocal(r, ceg::CubicElementGrid, nt::Int; tmin=0, tmax=30, δ=0.25)

- `r`                       : Nuclear coordinate (1D)
- `ceg::CubicElementGrid`   : Grid for the potential
- `nt::Int`                 : Number of t-points
"""
struct NuclearPotentialTensorLogLocal{T} <: AbstractNuclearPotentialSingle{T}
    elementgrid::Matrix{T}
    t::Vector{Float64}
    wt::Vector{Float64}
    tmin::Float64
    tmax::Float64
    r::T
    δ::Float64
    δp::Matrix{T}
    δm::Matrix{T}
    function NuclearPotentialTensorLogLocal(r, ceg::CubicElementGrid, nt::Int;
                                           tmin=0, tmax=30, δ=0.25)
        @assert 0 < δ <= 1
        @assert 0 <= tmin < tmax
        @assert dimension(r) == NoDims || dimension(r) == dimension(u"m")
        rt = austrip(r)
        s, ws = gausspoints(nt; elementsize=(log(tmin+1e-12), log(tmax)))
        t = exp.(s)
        wt = ws .* t
        t, wt = gausspoints(nt; elementsize=(tmin, tmax))
        grid = Array(grid1d(ceg)) .- rt
        ss = size(grid)

        # get range around points [x+δm, x+δp]
        # and calculate average on that range (later on)
        δp = zeros(ss...)
        δm = zeros(ss...)
        δm[2:end,:] = grid[1:end-1,:] .- grid[2:end,:]
        δp[1:end-1,:] = grid[2:end,:] .- grid[1:end-1,:]
        for j in 1:ss[2]
            δm[1,j] = ceg.elements[j].low - grid[1,j]
            δp[end,j] = ceg.elements[j].high - grid[end,j]
        end
        new{eltype(δm)}(grid, t, wt, tmin, tmax, rt, δ, δ.*δp, δ.*δm)
    end
end

"""
    NuclearPotentialTensorGaussian{T}

Tensor for nuclear potential in 1D. To get 3D you need 3 of these tensors and
integrate t-coordinate. This version has logarithmic spacing for t-points
with gaussian nuclear density.

Easy way is to use [`PotentialTensor`](@ref) to get 3D tensor.

# Fields
- `elementgrid::Matrix{T}`   : Grid definition, collums are elements, rows Gauss points
- `t::Vector{Float64}`       : t-coordinate
- `wt::Vector{Float64}`      : Integration weights for t
- `tmin::Float64`            : Minimum t-value
- `tmax::Float64`            : Maximum t-value
- `r::T`                     : Nuclear coordinate
- `β::Float64`               : Width of Gaussian exp(-βr²)

# Creation
    NuclearPotentialTensorGaussian(r, ceg::CubicElementGrid, nt; β=100, tmin=0, tmax=20)

- `r`                       : Nuclear coordinate (1D)
- `ceg::CubicElementGrid`   : Grid for the potential
- `nt::Int`                 : Number of t-points
"""
struct NuclearPotentialTensorGaussian{T} <: AbstractNuclearPotentialSingle{T}
    elementgrid::Matrix{T}
    t::Vector{Float64}
    wt::Vector{Float64}
    tmin::Float64
    tmax::Float64
    r::T
    β::Float64
    function NuclearPotentialTensorGaussian(r, ceg::CubicElementGrid, nt, β; tmin=0, tmax=20)
        @assert dimension(r) == NoDims || dimension(r) == dimension(u"m")
        rt = austrip(r)
        s, ws = gausspoints(nt; elementsize=(log(tmin+1e-12), log(tmax)))
        t = exp.(s)
        wt = ws .* t
        grid = Array(grid1d(ceg)) .- rt
        new{eltype(grid)}(grid, t, wt, tmin, tmax, rt, β)
    end
end

function NuclearPotentialTensorGaussian(r, ceg::CubicElementGrid, nt; σ=0.1, tmin=0, tmax=20)
    min_d = elementsize(ceg.elements) / size(ceg)[1] |> austrip
    β = 0.5*(σ*min_d)^-2
    return NuclearPotentialTensorGaussian(r, ceg, nt, β; tmin=tmin, tmax=tmax)
end

"""
    NuclearPotentialTensorCombination{T}

Combine different tensors to bigger one. Allows combining normal, logarithmic
and local correction ones to one tensor.

# Fields
- `subtensors::Vector{AbstractNuclearPotentialSingle{T}}`  :  combination of these
- `ref::Vector{Int}`    :  transforms t-index to subtensor index
- `index::Vector{Int}`  :  transforms t-index to subtensors t-index
- `wt::Vector{Float64}` :  t-integration weight
- `tmin::Float64`       :  minimum t-value
- `tmax::Float64`       :  maximum t-value
- `r::T`                :  nuclear coordinate

# Creation
    NuclearPotentialTensorCombination(npt::AbstractNuclearPotentialSingle{T}...)
"""
mutable struct NuclearPotentialTensorCombination{T} <: AbstractNuclearPotentialCombination{T}
    subtensors::Vector{AbstractNuclearPotentialSingle{T}}
    ref::Vector{Int}
    index::Vector{Int}
    t::Vector{Float64}
    wt::Vector{Float64}
    tmin::Float64
    tmax::Float64
    r::T
    function NuclearPotentialTensorCombination(npt::AbstractNuclearPotentialSingle)
        nt = size(npt)[end]
        new{eltype(npt.elementgrid)}([npt],
                   ones(nt),
                   1:nt,
                   deepcopy(npt.t),
                   deepcopy(npt.wt),
                   npt.tmin,
                   npt.tmax,
                   npt.r
                   )
    end
end

function NuclearPotentialTensorCombination(npt::AbstractNuclearPotentialSingle{T}...) where {T}
    tmp = NuclearPotentialTensorCombination(npt[1])
    for x in npt[2:end]
        push!(tmp, x)
    end
    return tmp
end

function Base.show(io::IO, ::MIME"text/plain", npt::AbstractNuclearPotential)
    print(io, "Nuclear potential tensor for $(length(npt.t)) t-points,"*
              " tmin=$(npt.tmin), tmax=$(npt.tmax)"
    )
end

Base.size(npt::AbstractNuclearPotential) = (size(npt.elementgrid)..., length(npt.t))

function Base.size(npct::NuclearPotentialTensorCombination)
    return (size(npct.subtensors[1])[1:2]..., length(npct.t))
end

function Base.getindex(npt::AbstractNuclearPotential, i::Int, j::Int, p::Int)
    r = npt.elementgrid[i,j]
    return exp(-npt.t[p]^2 * r^2)
end

function Base.getindex(npt::NuclearPotentialTensorLogLocal, i::Int, j::Int, p::Int)
    r = npt.elementgrid[i,j]
    δ = abs(npt.δp[i,j] - npt.δm[i,j])
    rmax = (r+δ) * npt.t[p]
    rmin = (r-δ) * npt.t[p]
    return 0.5 * √π * erf(rmin, rmax) / (rmax-rmin)
end

function Base.getindex(npt::NuclearPotentialTensorGaussian, i::Int, j::Int, p::Int)
    r = npt.elementgrid[i,j]
    t = npt.t[p]
    a = npt.β
    q = a * t^2 / (t^2 + a)
    return sqrt(a/(t^2+a)) * exp(-q * r^2)
end

function Base.getindex(npct::NuclearPotentialTensorCombination, i::Int, j::Int, p::Int)
    s = npct.ref[p]
    t = npct.index[p]
    return npct.subtensors[s][i,j,t]
end

function Base.push!(nptc::NuclearPotentialTensorCombination{T},
                    npt::AbstractNuclearPotentialSingle{T}) where {T}
    @assert nptc.r == npt.r
    @assert nptc.tmax <= npt.tmin
    push!(nptc.subtensors, npt)
    n = size(npt)[end]
    append!(nptc.ref, (length(nptc.subtensors))*ones(n))
    append!(nptc.t, npt.t)
    append!(nptc.wt, npt.wt)
    append!(nptc.index, 1:n)
    nptc.tmax = npt.tmax
    return nptc
end

"""
    PotentialTensor{Tx,Ty,Tz}

Used to create nuclear potential.

# Creation
    PotentialTensor(npx, npy, npz)

`npx`, `npy` and `npz` are nuclear potential tensors.
To get the best result usually they should be [`NuclearPotentialTensorCombination`](@ref).

# Usage
Create tensor and then call `Array` on it to get nuclear potential.

"""
struct PotentialTensor{Tx<:AbstractNuclearPotential,
                       Ty<:AbstractNuclearPotential,
                       Tz<:AbstractNuclearPotential}
    x::Tx
    y::Ty
    z::Tz
    function PotentialTensor(npx, npy, npz)
        @assert size(npx)[end] == size(npy)[end] == size(npz)[end]
        @assert npx.tmin == npy.tmin == npx.tmin
        @assert npx.tmax == npy.tmax == npx.tmax
        new{typeof(npx), typeof(npy), typeof(npz)}(npx, npy, npz)
    end
end

Base.show(io::IO, ::MIME"text/plain", pt::PotentialTensor) = print(io, "PotentialTensor")

function Base.size(pt::PotentialTensor)
    s1 = size(pt.x)
    s2 = size(pt.y)
    s3 = size(pt.z)
    return s1[1], s2[1], s3[1], s1[2], s2[2], s3[2]
end

function Base.getindex(pt::PotentialTensor, i::Int, j::Int, k::Int, I::Int, J::Int, K::Int)
    out = pt.x[i,I,:]
    out .*= pt.y[j,J,:]
    out .* pt.z[k,K,:]
    out .* pt.x.wt
    return out
end

function Base.Array(pt::PotentialTensor)
    x = Array(pt.x)
    y = Array(pt.y)
    z = Array(pt.z)
    @tullio out[i,j,k,I,J,K] := x[i,I,t] * y[j,J,t] * z[k,K,t] * pt.x.wt[t]
    return 2/sqrt(π) .* out
end

"""
    nuclear_potential(ceg::CubicElementGrid, q, r)
    nuclear_potential(ceg::CubicElementGrid, aname::String, r)

Gives nuclear potential for given nuclear charge or atom name.
"""
function nuclear_potential(ceg::CubicElementGrid, q, r)

    ncx = optimal_nuclear_tensor(ceg, r[1])
    ncy = optimal_nuclear_tensor(ceg, r[2])
    ncz = optimal_nuclear_tensor(ceg, r[3])

    return nuclear_potential(ceg, q, ncx, ncy, ncz)
end

function optimal_nuclear_tensor(ceg, x)
    np = NuclearPotentialTensor(x, ceg, 64; tmin=0, tmax=70)
    npl = NuclearPotentialTensorLog(x, ceg, 16; tmin=70, tmax=120)
    npll = NuclearPotentialTensorGaussian(x, ceg, 16;  σ=0.1, tmin=120, tmax=300)
    nplll = NuclearPotentialTensorGaussian(x, ceg, 32; σ=0.1, tmin=300, tmax=20000)
    return NuclearPotentialTensorCombination(np, npl, npll, nplll)
end

function nuclear_potential(ceg,
                           q,
                           npx::AbstractNuclearPotential,
                           npy::AbstractNuclearPotential,
                           npz::AbstractNuclearPotential)
    @assert dimension(q) == dimension(u"C") || dimension(q) == NoDims
    pt = PotentialTensor(npx, npy, npz)
    qau = austrip(q)
    return (qau*u"hartree/e_au") * ScalarOperator(ceg, Array(pt))
end

function nuclear_potential(ceg::CubicElementGrid, aname::String, r)
    if length(aname) <= 2    # is atomic symbol
        e = elements[Symbol(aname)]
    else
        e = elements[aname]
    end
    @debug "element number $(e.number)"
    return nuclear_potential(ceg, e.number, r)
end


text:
using DelimitedFiles

input = joinpath(@__DIR__, "input")
entries = readdlm(input, Int)

function twosum(entries)
    for i in 1:(length(entries))
        for j in (i + 1):(length(entries))
            if entries[i] + entries[j] == 2020
                return entries[i] * entries[j]
            end
        end
    end

    return 0
end

function threesum(entries)
    for i in 1:(length(entries))
        for j in (i + 1):(length(entries))
            for k in (j + 1):(length(entries))
                if entries[i] + entries[j] + entries[k] == 2020
                    return entries[i] * entries[j] * entries[k]
                end
            end
        end
    end

    return 0
end

p1 = twosum(entries)
p2 = threesum(entries)

@assert(p1 == 876459)
@assert(p2 == 116168640)

println("-----------------------------------------------------------------------")
println("report repair -- part one :: $p1")
println("report repair -- part two :: $p2")
println("-----------------------------------------------------------------------")


text:
function is_valid_unstructured_mesh(dataset::_T,
                                    repeat_cells::Bool
                                   ) where {_T<:AbstractVTKUnstructuredData}
    T = repr(_T)
    _point_coords, _cell_types, _cell_connectivity, _point_data, _cell_data = dataset.point_coords,
        dataset.cell_types, dataset.cell_connectivity, dataset.point_data, dataset.cell_data

    in(dim(dataset), [2,3]) || return false, "$T: Only 2 and 3 dimensional mesh data are allowed."

    length(_cell_types) == length(_cell_connectivity) || return false, "$T: The cell types vector and the cell connectivity vector have unequal lengths."
    length(_cell_connectivity) == 0 || (_min = get_lowest_index(_cell_connectivity)) >= 1 || return false, "$T: The lowest point index used in the cell connectivity vector is $_min < 1. Please consider using the add_point_id_offset!(cell_connectivity, $(1-_min)) function to make the lowest index 1."
    length(_cell_connectivity) == 0 || get_highest_index(_cell_connectivity) <= num_of_points(dataset) || return false, "$T: The highest point index used in the cell connectivity vector is more than the number of points passed."

    _keys = collect(keys(_point_data))
    length(_keys) == 0 ||
    num_of_points(dataset) == (var_dim(dataset, _keys[1], "Point") == 1 ? length(_point_data[_keys[1]]) : size(_point_data[_keys[1]], 2)) &&
    begin
        _out = true
        _var_dim1 = var_dim(dataset, _keys[1], "Point")
        for m in _keys[2:end]
            _var_dimi = var_dim(dataset, m, "Point")
            (_out = (_var_dim1 == 1 ? length(_point_data[_keys[1]]) : size(_point_data[_keys[1]], 2)) ==
                (_var_dimi == 1 ? length(_point_data[m]) : size(_point_data[m], 2))) || break
        end
        _out
    end || return false, "$T: The number of points defined is not consistent with the number of points in at least one of the point data fields."

    _keys = collect(keys(_cell_data))
    length(_keys) == 0 ||
    num_of_cells(dataset) == (var_dim(dataset, _keys[1], "Cell") == 1 ? length(_cell_data[_keys[1]]) : size(_cell_data[_keys[1]], 2)) &&
    begin
        _out = true
        _var_dim1 = var_dim(dataset, _keys[1], "Cell")
        for m in _keys[2:end]
            _var_dimi = var_dim(dataset, m, "Cell")
            (_out = (_var_dim1 == 1 ? length(_cell_data[_keys[1]]) : size(_cell_data[_keys[1]], 2)) ==
                (_var_dimi == 1 ? length(_cell_data[m]) : size(_cell_data[m], 2))) || break
        end
        _out
    end || return false, "$T: The number of cells defined is not consistent with the number of cells in at least one of the cell data fields."

    begin
        _out = true
        for i in 1:length(_cell_connectivity)
            (_out = VTK_CELL_TYPE[cell_type(dataset, i)].nodes == -1 ||
                length(_cell_connectivity[i]) == VTK_CELL_TYPE[cell_type(dataset, i)].nodes) || break
        end
        _out
    end || return false, "$T: The number of points in at least one cell does not match its corresponding cell type."

    if !repeat_cells
        begin
            _out = true
            for i in 1:length(_cell_connectivity)
                for j in i+1:length(_cell_connectivity)
                    (_out = !similar_cells(_cell_connectivity[i], cell_type(dataset, i),
                        _cell_connectivity[j], cell_type(dataset, j))) || break
                end
            end
            _out
        end || return false, "$T: Repeat cells are not allowed."
    end

    return true, ""
end

function is_valid(dataset::VTKUnstructuredData; repeat_cells=false)
    return is_valid_unstructured_mesh(dataset, repeat_cells)
end

function is_valid(dataset::VTKPolyData; repeat_cells=false)
    T = "VTKPolyData"
    _cell_types, _cell_connectivity = dataset.cell_types, dataset.cell_connectivity

    length(_cell_types) == 0 || begin
        _out = true
        poly_cells = [POINT_CELLS; LINE_CELLS; FACE_CELLS]
        for i in 1:length(_cell_connectivity)
            (_out = in(cell_type(dataset, i), poly_cells)) || break
        end
        _out
    end || return false, "$T: At least one of the cells is a volume cell. $T can only have point, line or face cells, consider VTKUnstructuredData."

    return is_valid_unstructured_mesh(dataset, repeat_cells)
end

function is_valid(dataset::VTKStructuredData)
    T = "VTKStructuredData"
    _dim = dim(dataset)
    in(_dim, [2,3]) && _dim == length(size(dataset.point_coords)[2:end]) || return false, "$T: Only 2 and 3 dimensional mesh data are allowed. First index iterates over the coordinate components, and the remaining indices iterate over the extents of the mesh."

    _extents = extents(dataset)
    all([_extents[i] > 0 for i in 1:length(_extents)]) || return false, "$T: Extent must be at least 1 in each dimensions."

    _keys = collect(keys(dataset.point_data))
    length(_keys) == 0 || begin
        _out = true
        for i in _keys
            var_extents = size(dataset.point_data[i])
            (_out = var_extents == _extents) || (_out = var_extents[2:end] == _extents) || break
        end
        _out
    end || return false, "$T: The point extents of the mesh and one of the point data fields are not consistent."

    cell_extents = (_extents .- 1)
    _keys = collect(keys(dataset.cell_data))
    length(_keys) == 0 || begin
        _out = true
        for i in _keys
            var_extents = size(dataset.cell_data[i])
            (_out = var_extents == cell_extents) || (_out = var_extents[2:end] == cell_extents) || break
        end
        _out
    end || return false, "$T: The cell extents of the mesh and one of the cell data fields are not consistent."

    return true, ""
end

function is_valid(dataset::VTKRectilinearData{S}) where {S}
    T = "VTKRectilinearData"
    in(length(dataset.point_coords), [2,3]) || return false, "$T: Only 2 and 3 dimensional mesh data are allowed."

    !(S <: Real) || begin
        _out = true
        for i in length(dataset.point_coords), j in 1:length(dataset.point_coords[i])-1
            (_out = dataset.point_coords[i][j] < dataset.point_coords[i][j+1]) || break
        end
        _out
    end || return false, "$T: Coordinate vectors should be in ascending order."

    _extents = extents(dataset)
    all([_extents[i] > 0 for i in 1:length(_extents)]) || return false, "$T: Extent must be at least 1 in each dimensions."
    _keys = collect(keys(dataset.point_data))
    length(_keys) == 0 || begin
        _out = true
        for i in _keys
            var_extents = size(dataset.point_data[i])
            (_out = var_extents == _extents) || (_out = var_extents[2:end] == _extents) || break
        end
        _out
    end || return false, "$T: The point extents of the mesh and one of the point data fields are not consistent."

    _cell_extents = cell_extents(dataset)
    _keys = collect(keys(dataset.cell_data))
    length(_keys) == 0 || begin
        _out = true
        for i in _keys
            var_extents = size(dataset.cell_data[i])
            (_out = var_extents == _cell_extents) || (_out = var_extents[2:end] == _cell_extents) || break
        end
        _out
    end || return false, "$T: The cell extents of the mesh and one of the cell data fields are not consistent."

    return true, ""
end

function is_valid(dataset::VTKUniformRectilinearData{S}) where {S}
    T = "VTKUniformRectilinearData"
    !(S <: Real) || begin
        _out = true
        for i in 1:length(dataset.spacing)
            (_out = dataset.spacing[i] > 0) || return false, "Spacing must be positive for all dimensions."
        end
        _out
    end
    _extents = extents(dataset)
    all([_extents[i] > 0 for i in 1:length(_extents)]) || return false, "$T: Extent must be at least 1 in each dimensions."

    in(length(_extents), [2,3]) || return false, "$T: Only 2 and 3 dimensional mesh data are allowed."
    length(_extents) == length(dataset.spacing) == length(dataset.origin) || return false, "$T: The length of \"origin\", the length of \"spacing\", and the length of \"extents\" are not consistent."

    _keys = collect(keys(dataset.point_data))
    length(_keys) == 0 || begin
        _out = true
        for m in _keys
            var_extents = size(dataset.point_data[m])
            (_out = var_extents == _extents) || (_out = var_extents[2:end] == _extents) || break
        end
        _out
    end || return false, "$T: The point extents of the mesh and one of the point data fields are not consistent."

    cell_extents = _extents .- 1
    _keys = collect(keys(dataset.cell_data))
    length(_keys) == 0 || begin
        _out = true
        for m in _keys
            var_extents = size(dataset.cell_data[m])
            (_out = var_extents == cell_extents) || (_out = var_extents[2:end] == cell_extents) || break
        end
        _out
    end || return false, "$T: The cell extents of the mesh and one of the cell data fields are not consistent."

    return true, ""
end

function is_valid(dataset::VTKMultiblockData)
    for block in dataset
        valid, _error = is_valid(block)
        valid ? continue : return false, _error
    end
    return true, ""
end

function is_valid(dataset::VTKTimeSeriesData)
    T = "VTKTimeSeriesData"
    timemarkers = dataset.timemarkers
    data = dataset.data
    begin
        _out = true
        for i in 1:length(timemarkers)-1
            (_out = timemarkers[i] < timemarkers[i+1]) || break
        end
        _out
    end || return false, "$T: Time markers should be in ascending order."

    length(timemarkers) == length(data) || return false, "$T: The number of time markers and the length of data must be equal."

    for block in dataset
        valid, _error = is_valid(block)
        valid ? continue : return false, _error
    end
    return true, ""
end