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

using Test

import Base.LinAlg.BlasInt

@test_throws ArgumentError Base.LinAlg.LAPACK.chkuplo('Z')
@test_throws ArgumentError Base.LinAlg.LAPACK.chkside('Z')
@test_throws ArgumentError Base.LinAlg.LAPACK.chkdiag('Z')
@test_throws ArgumentError Base.LinAlg.LAPACK.chktrans('Z')

@testset "syevr" begin
    guardsrand(123) do
        Ainit = randn(5,5)
        @testset for elty in (Float32, Float64, Complex64, Complex128)
            if elty == Complex64 || elty == Complex128
                A = complex.(Ainit, Ainit)
            else
                A = Ainit
            end
            A = convert(Array{elty, 2}, A)
            Asym = A'A
            vals, Z = LAPACK.syevr!('V', copy(Asym))
            @test Z*(Diagonal(vals)*Z') ≈ Asym
            @test all(vals .> 0.0)
            @test LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[vals .< 1.0]
            @test LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[4:5]
            @test vals ≈ LAPACK.syev!('N','U',copy(Asym))

            @test_throws DimensionMismatch LAPACK.sygvd!(1,'V','U',copy(Asym),ones(elty,6,6))
        end
    end
end

@testset "gglse" begin
    let
        @testset for elty in (Float32, Float64, Complex64, Complex128)
            A = convert(Array{elty, 2}, [1 1 1 1; 1 3 1 1; 1 -1 3 1; 1 1 1 3; 1 1 1 -1])
            c = convert(Array{elty, 1}, [2, 1, 6, 3, 1])
            B = convert(Array{elty, 2}, [1 1 1 -1; 1 -1 1 1; 1 1 -1 1])
            d = convert(Array{elty, 1}, [1, 3, -1])
            @test LAPACK.gglse!(A, c, B, d)[1] ≈ convert(Array{elty}, [0.5, -0.5, 1.5, 0.5])
        end
    end
end

@testset "gebrd, bdsqr, throw for bdsdc" begin
    let
        n = 10
        @testset for elty in (Float32, Float64)
            d, e = convert(Vector{elty}, randn(n)), convert(Vector{elty}, randn(n - 1))
            U, Vt, C = eye(elty, n), eye(elty, n), eye(elty, n)
            s, _ = LAPACK.bdsqr!('U', copy(d), copy(e), Vt, U, C)
            @test Array(Bidiagonal(d, e, :U)) ≈ U*Diagonal(s)*Vt

            @test_throws ArgumentError LAPACK.bdsqr!('A', d, e, Vt, U, C)
            @test_throws DimensionMismatch LAPACK.bdsqr!('U', d, [e; 1], Vt, U, C)
            @test_throws DimensionMismatch LAPACK.bdsqr!('U', d, e, Vt[1:end - 1, :], U, C)
            @test_throws DimensionMismatch LAPACK.bdsqr!('U', d, e, Vt, U[:,1:end - 1], C)
            @test_throws DimensionMismatch LAPACK.bdsqr!('U', d, e, Vt, U, C[1:end - 1, :])

            @test_throws ArgumentError LAPACK.bdsdc!('U','Z',d,e)

            A = rand(elty,n,n)
            B = copy(A)
            B, d, e, tauq, taup = LAPACK.gebrd!(B)
            U, Vt, C = eye(elty, n), eye(elty, n), eye(elty, n)
            s, _ = LAPACK.bdsqr!('U',d,e[1:n-1],Vt, U, C)
            @test s ≈ svdvals(A)
        end
    end
end

@testset "Issue #7886" begin
    let
        x, r = LAPACK.gelsy!([0 1; 0 2; 0 3.], [2, 4, 6.])
        @test x ≈ [0,2]
        @test r == 1
    end
end

@testset "geqrt(3)" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = copy(A)
        C,T = LAPACK.geqrt!(A,zeros(elty,10,10))
        D,S = LAPACK.geqrt3!(A,zeros(elty,10,10))
        @test C ≈ D
    end
end

@testset "gbtrf and gbtrs" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        d = rand(elty,6)
        dl = rand(elty,5)
        du = rand(elty,5)
        dl2 = rand(elty,4)
        AB = zeros(elty,6,6)
        AB[6,1:4] = dl2
        AB[5,1:5] = dl
        AB[4,:] = d
        AB[3,2:6] = du
        AB,ipiv = LAPACK.gbtrf!(2,1,6,AB)
        C = rand(elty,6,6)
        D = copy(C)
        D = LAPACK.gbtrs!('N',2,1,6,AB,ipiv,D)
        A = diagm(-2 => dl2, -1 => dl, 0 => d, 1 => du)
        @test A\C ≈ D
        @test_throws DimensionMismatch LAPACK.gbtrs!('N',2,1,6,AB,ipiv,ones(elty,7,6))
        @test_throws Base.LinAlg.LAPACKException LAPACK.gbtrf!(2,1,6,zeros(AB))
    end
end

@testset "geqp3, geqrt error handling" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        @test_throws DimensionMismatch LAPACK.geqlf!(A,zeros(elty,11))
        @test_throws DimensionMismatch LAPACK.gelqf!(A,zeros(elty,11))
        @test_throws DimensionMismatch LAPACK.geqp3!(A,ones(Base.LinAlg.BlasInt,11),zeros(elty,10))
        @test_throws DimensionMismatch LAPACK.geqp3!(A,ones(Base.LinAlg.BlasInt,10),zeros(elty,11))
        @test_throws ArgumentError LAPACK.geqrt!(A,zeros(elty,11,10))
        @test_throws DimensionMismatch LAPACK.geqrt3!(A,zeros(elty,11,10))
        @test_throws DimensionMismatch LAPACK.geqrt3!(ones(elty,10,11),zeros(elty,11,11))
        @test_throws DimensionMismatch LAPACK.geqrf!(A,zeros(elty,11))
        @test_throws DimensionMismatch LAPACK.gerqf!(A,zeros(elty,11))
    end
end

@testset "gels, gesv, getrs, getri error handling" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = rand(elty,11,11)
        @test_throws DimensionMismatch LAPACK.gels!('N',A,B)
        @test_throws DimensionMismatch LAPACK.gels!('T',A,B)
        @test_throws DimensionMismatch LAPACK.gesv!(A,B)
        @test_throws DimensionMismatch LAPACK.getrs!('N',A,ones(Base.LinAlg.BlasInt,10),B)
        @test_throws DimensionMismatch LAPACK.getrs!('T',A,ones(Base.LinAlg.BlasInt,10),B)
        @test_throws DimensionMismatch LAPACK.getri!(A,ones(Base.LinAlg.BlasInt,11))
    end
end

@testset "gelsy, gelsd" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = rand(elty,10,10)
        C, j = LAPACK.gelsd!(copy(A),copy(B))
        D, k = LAPACK.gelsy!(copy(A),copy(B))
        @test C ≈ D
        @test_throws DimensionMismatch LAPACK.gelsd!(A,rand(elty,12,10))
        @test_throws DimensionMismatch LAPACK.gelsy!(A,rand(elty,12,10))
    end
end

@testset "gglse errors" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        @test_throws DimensionMismatch LAPACK.gglse!(A,zeros(elty,10),rand(elty,12,11),zeros(elty,12))
        @test_throws DimensionMismatch LAPACK.gglse!(A,zeros(elty,11),rand(elty,10,10),zeros(elty,10))
        @test_throws DimensionMismatch LAPACK.gglse!(A,zeros(elty,10),rand(elty,10,10),zeros(elty,11))
    end
end

@testset "gesvd, ggsvd" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,5)
        U,S,V = svd(A)
        lU,lS,lVt = LAPACK.gesvd!('S','S',A)
        @test U ≈ lU
        @test S ≈ lS
        @test V' ≈ lVt
        B = rand(elty,10,10)
        # xggsvd3 replaced xggsvd in LAPACK 3.6.0
        if LAPACK.version() < v"3.6.0"
            @test_throws DimensionMismatch LAPACK.ggsvd!('S','S','S',A,B)
        else
            @test_throws DimensionMismatch LAPACK.ggsvd3!('S','S','S',A,B)
        end
    end
end

@testset "geevx, ggev errors" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = rand(elty,10,10)
        @test_throws ArgumentError LAPACK.geevx!('M','N','N','N',A)
        @test_throws ArgumentError LAPACK.geevx!('N','Z','N','N',A)
        @test_throws ArgumentError LAPACK.geevx!('N','N','Z','N',A)
        @test_throws ArgumentError LAPACK.geevx!('N','N','N','Z',A)
        @test_throws ArgumentError LAPACK.ggev!('N','B',A,B)
        @test_throws ArgumentError LAPACK.ggev!('B','N',A,B)
        @test_throws DimensionMismatch LAPACK.ggev!('N','N',A,zeros(elty,12,12))
    end
end

@testset "gebal/gebak" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10) * Diagonal(exp10.(linspace(-10,10,10)))
        B = copy(A)
        ilo, ihi, scale = LAPACK.gebal!('S',B)
        Bvs = eigvecs(B)
        Avs = eigvecs(A)
        Bvs = LAPACK.gebak!('S','R',ilo,ihi,scale,Bvs)
        @test norm(diff(Avs ./ Bvs)) < 100 * eps(abs(float(one(elty))))
    end
end

@testset "gels" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        X = rand(elty,10)
        B,Y,z = LAPACK.gels!('N',copy(A),copy(X))
        @test A\X ≈ Y
    end
end

@testset "getrf/getri" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        iA = inv(A)
        A, ipiv = LAPACK.getrf!(A)
        A = LAPACK.getri!(A, ipiv)
        @test A ≈ iA
    end
end

@testset "geev" begin
    # complex is easier for now
    @testset for elty in (Complex64, Complex128)
        A = rand(elty,10,10)
        Aw, Avl, Avr = LAPACK.geev!('N','V',copy(A))
        fA = eigfact(A)
        @test fA[:values] ≈ Aw
        @test fA[:vectors] ≈ Avr
    end
end

@testset "gtsv" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        du = rand(elty,9)
        d  = rand(elty,10)
        dl = rand(elty,9)
        b  = rand(elty,10)
        c = Tridiagonal(dl,d,du) \ b
        b = LAPACK.gtsv!(dl,d,du,b)
        @test b ≈ c
        @test_throws DimensionMismatch LAPACK.gtsv!(zeros(elty,11),d,du,b)
        @test_throws DimensionMismatch LAPACK.gtsv!(dl,d,zeros(elty,11),b)
        @test_throws DimensionMismatch LAPACK.gtsv!(dl,d,du,zeros(elty,11))
        @test LAPACK.gtsv!(elty[],elty[],elty[],elty[]) == elty[]
    end
end

@testset "gttrs,gttrf errors" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        du = rand(elty,9)
        d  = rand(elty,10)
        dl = rand(elty,9)
        b  = rand(elty,10)
        @test_throws DimensionMismatch LAPACK.gttrf!(zeros(elty,11),d,du)
        @test_throws DimensionMismatch LAPACK.gttrf!(dl,d,zeros(elty,11))
        @test_throws DimensionMismatch LAPACK.gttrs!('N',zeros(elty,11),d,du,ones(elty,9),zeros(Base.LinAlg.BlasInt,10),b)
        @test_throws DimensionMismatch LAPACK.gttrs!('N',dl,d,zeros(elty,11),ones(elty,9),zeros(Base.LinAlg.BlasInt,10),b)
        @test_throws DimensionMismatch LAPACK.gttrs!('N',dl,d,du,ones(elty,9),zeros(Base.LinAlg.BlasInt,10),zeros(elty,11))
        A = lufact(Tridiagonal(dl,d,du))
        b  = rand(elty,10,5)
        c = copy(b)
        dl,d,du,du2,ipiv = LAPACK.gttrf!(dl,d,du)
        c = LAPACK.gttrs!('N',dl,d,du,du2,ipiv,c)
        @test A\b ≈ c
    end
end

@testset "orglq and friends errors" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        A,tau = LAPACK.gelqf!(A)
        @test_throws DimensionMismatch LAPACK.orglq!(A,tau,11)
        @test_throws DimensionMismatch LAPACK.ormlq!('R','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormlq!('L','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormlq!('R','N',A,zeros(elty,11),rand(elty,10,10))
        @test_throws DimensionMismatch LAPACK.ormlq!('L','N',A,zeros(elty,11),rand(elty,10,10))

        B = copy(A)
        C = LAPACK.orglq!(B,tau)
        @test LAPACK.ormlq!('R','N',A,tau,eye(elty,10)) ≈ C

        A = rand(elty,10,10)
        A,tau = LAPACK.geqrf!(A)
        @test_throws DimensionMismatch LAPACK.orgqr!(A,tau,11)
        B = copy(A)
        @test LAPACK.orgqr!(B,tau) ≈ LAPACK.ormqr!('R','N',A,tau,eye(elty,10))
        @test_throws DimensionMismatch LAPACK.ormqr!('R','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormqr!('L','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormqr!('R','N',A,zeros(elty,11),rand(elty,10,10))
        @test_throws DimensionMismatch LAPACK.ormqr!('L','N',A,zeros(elty,11),rand(elty,10,10))

        A = rand(elty,10,10)
        A,tau = LAPACK.geqlf!(A)
        @test_throws DimensionMismatch LAPACK.orgql!(A,tau,11)
        B = copy(A)
        @test LAPACK.orgql!(B,tau) ≈ LAPACK.ormql!('R','N',A,tau,eye(elty,10))
        @test_throws DimensionMismatch LAPACK.ormql!('R','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormql!('L','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormql!('R','N',A,zeros(elty,11),rand(elty,10,10))
        @test_throws DimensionMismatch LAPACK.ormql!('L','N',A,zeros(elty,11),rand(elty,10,10))

        A = rand(elty,10,10)
        A,tau = LAPACK.gerqf!(A)
        @test_throws DimensionMismatch LAPACK.orgrq!(A,tau,11)
        B = copy(A)
        @test LAPACK.orgrq!(B,tau) ≈ LAPACK.ormrq!('R','N',A,tau,eye(elty,10))
        @test_throws DimensionMismatch LAPACK.ormrq!('R','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormrq!('L','N',A,tau,rand(elty,11,11))
        @test_throws DimensionMismatch LAPACK.ormrq!('R','N',A,zeros(elty,11),rand(elty,10,10))
        @test_throws DimensionMismatch LAPACK.ormrq!('L','N',A,zeros(elty,11),rand(elty,10,10))

        A = rand(elty,10,11)
        Q = copy(A)
        Q,tau = LAPACK.gerqf!(Q)
        R = triu(Q[:,2:11])
        LAPACK.orgrq!(Q,tau)
        @test Q*Q' ≈ eye(elty,10)
        @test R*Q ≈ A
        @test_throws DimensionMismatch LAPACK.orgrq!(zeros(elty,11,10),zeros(elty,10))

        C = rand(elty,10,10)
        V = rand(elty,10,10)
        T = zeros(elty,10,11)
        @test_throws DimensionMismatch LAPACK.gemqrt!('L','N',V,T,C)
        @test_throws DimensionMismatch LAPACK.gemqrt!('R','N',V,T,C)

        C = rand(elty,10,10)
        V = rand(elty,11,10)
        T = zeros(elty,10,10)
        @test_throws DimensionMismatch LAPACK.gemqrt!('R','N',V,T,C)
        @test_throws DimensionMismatch LAPACK.gemqrt!('L','N',V,T,C)

        # test size(T) = (nb,k) ensures 1 <= nb <= k
        T = zeros(elty,10,10)
        V = rand(elty,5,10)
        @test_throws DimensionMismatch LAPACK.gemqrt!('L','N',V,T,C)
        C = rand(elty,10,10)
        V = rand(elty,10,10)
        T = zeros(elty,11,10)
        @test_throws DimensionMismatch LAPACK.gemqrt!('R','N',V,T,C)

        @test_throws DimensionMismatch LAPACK.orghr!(1, 10, C, zeros(elty,11))
    end
end

@testset "sytri, sytrs, and sytrf" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        A = A + A.' #symmetric!
        B = copy(A)
        B,ipiv = LAPACK.sytrf!('U',B)
        @test triu(inv(A)) ≈ triu(LAPACK.sytri!('U',B,ipiv))
        @test_throws DimensionMismatch LAPACK.sytrs!('U',B,ipiv,rand(elty,11,5))
        @test LAPACK.sytrf!('U',zeros(elty,0,0)) == (zeros(elty,0,0),zeros(BlasInt,0))
    end

    # Rook-pivoting variants
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        A = A + A.' #symmetric!
        B = copy(A)
        B,ipiv = LAPACK.sytrf_rook!('U',B)
        @test triu(inv(A)) ≈ triu(LAPACK.sytri_rook!('U',B,ipiv))
        @test_throws DimensionMismatch LAPACK.sytrs_rook!('U',B,ipiv,rand(elty,11,5))
        @test LAPACK.sytrf_rook!('U',zeros(elty,0,0)) == (zeros(elty,0,0),zeros(BlasInt,0))
        A = rand(elty,10,10)
        A = A + A.' #symmetric!
        b = rand(elty,10)
        c = A \ b
        b,A = LAPACK.sysv_rook!('U',A,b)
        @test b ≈ c
        @test_throws DimensionMismatch LAPACK.sysv_rook!('U',A,rand(elty,11))

        # syconvf_rook error handling
        # way argument is wrong
        @test_throws ArgumentError LAPACK.syconvf_rook!('U', 'U', A, rand(BlasInt, 10))
        # ipiv has wrong length
        @test_throws ArgumentError LAPACK.syconvf_rook!('U', 'R', A, rand(BlasInt, 9))
        # e has wrong length
        @test_throws ArgumentError LAPACK.syconvf_rook!('U', 'R', A, rand(BlasInt, 10), rand(elty, 9))
    end
end

@testset "hetrf, hetrs" begin
    @testset for elty in (Complex64, Complex128)
        A = rand(elty,10,10)
        A = A + A' #hermitian!
        B = copy(A)
        B,ipiv = LAPACK.hetrf!('U',B)
        @test_throws DimensionMismatch LAPACK.hetrs!('U',B,ipiv,rand(elty,11,5))
        @test_throws DimensionMismatch LAPACK.hetrs_rook!('U',B,ipiv,rand(elty,11,5))
    end
end

@testset "stev, stebz, stein, stegr" begin
    @testset for elty in (Float32, Float64)
        d = rand(elty,10)
        e = rand(elty,9)
        @test_throws DimensionMismatch LAPACK.stev!('U',d,rand(elty,10))
        @test_throws DimensionMismatch LAPACK.stebz!('A','B',zero(elty),zero(elty),0,0,-1.,d,rand(elty,10))
        @test_throws DimensionMismatch LAPACK.stegr!('N','A',d,rand(elty,10),zero(elty),zero(elty),0,0)
        @test_throws DimensionMismatch LAPACK.stein!(d,zeros(elty,10),zeros(d),zeros(BlasInt,10),zeros(BlasInt,10))
        @test_throws DimensionMismatch LAPACK.stein!(d,e,zeros(elty,11),zeros(BlasInt,10),zeros(BlasInt,10))
    end
end

@testset "trtri & trtrs" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        A = triu(A)
        B = copy(A)
        @test inv(A) ≈ LAPACK.trtri!('U','N',B)
        @test_throws DimensionMismatch LAPACK.trtrs!('U','N','N',B,zeros(elty,11,10))
    end
end

@testset "tgsen, tzrzf, & trsyl" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        Z = zeros(elty,10,10)
        @test_throws DimensionMismatch LAPACK.tgsen!(zeros(BlasInt,10),Z,zeros(elty,11,11),Z,Z)
        @test_throws DimensionMismatch LAPACK.tgsen!(zeros(BlasInt,10),Z,Z,zeros(elty,11,11),Z)
        @test_throws DimensionMismatch LAPACK.tgsen!(zeros(BlasInt,10),Z,Z,Z,zeros(elty,11,11))
        @test_throws DimensionMismatch LAPACK.trsyl!('N','N',Z,Z,zeros(elty,11,11))
        @test_throws DimensionMismatch LAPACK.tzrzf!(zeros(elty,10,5))
    end
end

@testset "sysv" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        guardsrand(123) do
            A = rand(elty,10,10)
            A = A + A.' #symmetric!
            b = rand(elty,10)
            c = A \ b
            b,A = LAPACK.sysv!('U',A,b)
            @test b ≈ c
            @test_throws DimensionMismatch LAPACK.sysv!('U',A,rand(elty,11))
        end
    end
end

@testset "hesv" begin
    @testset for elty in (Complex64, Complex128)
        A = rand(elty,10,10)
        A = A + A' #hermitian!
        b = rand(elty,10)
        c = A \ b
        b,A = LAPACK.hesv!('U',A,b)
        @test b ≈ c
        @test_throws DimensionMismatch LAPACK.hesv!('U',A,rand(elty,11))
        A = rand(elty,10,10)
        A = A + A' #hermitian!
        b = rand(elty,10)
        c = A \ b
        b,A = LAPACK.hesv_rook!('U',A,b)
        @test b ≈ c
        @test_throws DimensionMismatch LAPACK.hesv_rook!('U',A,rand(elty,11))
    end
end

@testset "ptsv" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        dv = ones(elty,10)
        ev = zeros(elty,9)
        rdv = real(dv)
        A = SymTridiagonal(dv,ev)
        if elty <: Complex
            A = Tridiagonal(conj(ev),dv,ev)
        end
        B = rand(elty,10,10)
        C = copy(B)
        @test A\B ≈ LAPACK.ptsv!(rdv,ev,C)
        @test_throws DimensionMismatch LAPACK.ptsv!(rdv,ones(elty,10),C)
        @test_throws DimensionMismatch LAPACK.ptsv!(rdv,ev,ones(elty,11,11))
    end
end

@testset "pttrf and pttrs" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        dv = ones(elty,10)
        ev = zeros(elty,9)
        rdv = real(dv)
        A = SymTridiagonal(dv,ev)
        if elty <: Complex
            A = Tridiagonal(conj(ev),dv,ev)
        end
        rdv,ev = LAPACK.pttrf!(rdv,ev)
        @test_throws DimensionMismatch LAPACK.pttrf!(rdv,dv)
        B = rand(elty,10,10)
        C = copy(B)
        if elty <: Complex
            @test A\B ≈ LAPACK.pttrs!('U',rdv,ev,C)
            @test_throws DimensionMismatch LAPACK.pttrs!('U',rdv,ones(elty,10),C)
            @test_throws DimensionMismatch LAPACK.pttrs!('U',rdv,ev,rand(elty,11,11))
        else
            @test A\B ≈ LAPACK.pttrs!(rdv,ev,C)
            @test_throws DimensionMismatch LAPACK.pttrs!(rdv,ones(elty,10),C)
            @test_throws DimensionMismatch LAPACK.pttrs!(rdv,ev,rand(elty,11,11))
        end
    end
end

@testset "posv and some errors for friends" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)/100
        A += real(diagm(0 => 10*real(rand(elty,10))))
        if elty <: Complex
            A = A + A'
        else
            A = A + A.'
        end
        B = rand(elty,10,10)
        D = copy(A)
        C = copy(B)
        D,C = LAPACK.posv!('U',D,C)
        @test A\B ≈ C
        @test_throws DimensionMismatch LAPACK.posv!('U',D,ones(elty,12,12))
        @test_throws DimensionMismatch LAPACK.potrs!('U',D,ones(elty,12,12))

        @test LAPACK.potrs!('U',zeros(elty,0,0),ones(elty,0)) == ones(elty,0)
    end
end

@testset "gesvx" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = rand(elty,10,5)
        C = copy(A)
        D = copy(B)
        X, rcond, f, b, r = LAPACK.gesvx!(C,D)
        @test X ≈ A\B
    end
end

@testset "gees, gges error throwing" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        A = rand(elty,10,10)
        B = rand(elty,11,11)
        @test_throws DimensionMismatch LAPACK.gges!('V','V',A,B)
    end
end

@testset "trrfs & trevc" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        T = triu(rand(elty,10,10))
        S = copy(T)
        select = zeros(Base.LinAlg.BlasInt,10)
        select[1] = 1
        select,Vr = LAPACK.trevc!('R','S',select,copy(T))
        @test Vr ≈ eigvecs(S)[:,1]
        select = zeros(Base.LinAlg.BlasInt,10)
        select[1] = 1
        select,Vl = LAPACK.trevc!('L','S',select,copy(T))
        select = zeros(Base.LinAlg.BlasInt,10)
        select[1] = 1
        select,Vln,Vrn = LAPACK.trevc!('B','S',select,copy(T))
        @test Vrn ≈ eigvecs(S)[:,1]
        @test Vln ≈ Vl
        @test_throws ArgumentError LAPACK.trevc!('V','S',select,copy(T))
        @test_throws DimensionMismatch LAPACK.trrfs!('U','N','N',T,rand(elty,10,10),rand(elty,10,11))
    end
end

@testset "laic1" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        @test_throws DimensionMismatch LAPACK.laic1!(1,rand(elty,10),real(rand(elty)),rand(elty,11),rand(elty))
    end
end

@testset "trsen" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        for job in ('N', 'E', 'V', 'B')
            for c in ('V', 'N')
                A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
                T,Q,d = schur(A)
                s, sep = Base.LinAlg.LAPACK.trsen!(job,c,Array{LinAlg.BlasInt}([0,1,0,0]),T,Q)[4:5]
                @test d[1] ≈ T[2,2]
                @test d[2] ≈ T[1,1]
                if c == 'V'
                    @test  Q*T*Q' ≈ A
                end
                if job == 'N' || job == 'V'
                    @test iszero(s)
                else
                    @test s ≈ 0.8080423 atol=1e-6
                end
                if job == 'N' || job == 'E'
                    @test iszero(sep)
                else
                    @test sep ≈ 2. atol=3e-1
                end
            end
        end
    end
end

@testset "trexc" begin
    @testset for elty in (Float32, Float64, Complex64, Complex128)
        for c in ('V', 'N')
            A = convert(Matrix{elty}, [7 2 2 1; 1 5 2 0; 0 3 9 4; 1 1 1 4])
            T,Q,d = schur(A)
            Base.LinAlg.LAPACK.trexc!(c,LinAlg.BlasInt(1),LinAlg.BlasInt(2),T,Q)
            @test d[1] ≈ T[2,2]
            @test d[2] ≈ T[1,1]
            if c == 'V'
                @test Q*T*Q' ≈ A
            end
        end
    end
end

@testset "Julia vs LAPACK" begin
    # Test our own linear algebra functionality against LAPACK
    @testset for elty in (Float32, Float64, Complex{Float32}, Complex{Float64})
        for nn in (5,10,15)
            if elty <: Real
                A = convert(Matrix{elty}, randn(10,nn))
            else
                A = convert(Matrix{elty}, complex.(randn(10,nn),randn(10,nn)))
            end    ## LU (only equal for real because LAPACK uses different absolute value when choosing permutations)
            if elty <: Real
                FJulia  = Base.LinAlg.generic_lufact!(copy(A))
                FLAPACK = Base.LinAlg.LAPACK.getrf!(copy(A))
                @test FJulia.factors ≈ FLAPACK[1]
                @test FJulia.ipiv ≈ FLAPACK[2]
                @test FJulia.info ≈ FLAPACK[3]
            end

            ## QR
            FJulia  = LinAlg.qrfactUnblocked!(copy(A))
            FLAPACK = Base.LinAlg.LAPACK.geqrf!(copy(A))
            @test FJulia.factors ≈ FLAPACK[1]
            @test FJulia.τ ≈ FLAPACK[2]
        end
    end
end

# Issue 13976
let A = [NaN 0.0 NaN; 0 0 0; NaN 0 NaN]
    @test_throws ArgumentError exp(A)
end

# Issue 14065 (and 14220)
let A = [NaN NaN; NaN NaN]
    @test_throws ArgumentError eigfact(A)
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
function parse_line(line::String)::Tuple{Int,Int,Char,SubString}
    lower_bound::Int = 1
    upper_bound::Int = 1
    char::Char = 'a'

    first_idx::Int = 1
    for (idx, c) in enumerate(line)
        if c == '-'
            lower_bound = parse(Int, SubString(line, first_idx, idx - 1))
            first_idx = idx + 1
        elseif c == ' '
            upper_bound = parse(Int, SubString(line, first_idx, idx - 1))
            char = line[idx + 1]
            first_idx = idx + 4
            break
        end
    end
    password = SubString(line, first_idx, lastindex(line))
    return lower_bound, upper_bound, char, password
end

function count_char(char::Char, string::SubString)::Int
    counter = 0
    for c in string
        if c == char
            counter += 1
        end
    end
    return counter
end

function run(s)
    # Your code here
    count = 0
    for line in readlines(IOBuffer(s))
        (lower_bound, upper_bound, char, password) = parse_line(line)
        if lower_bound <= count_char(char, password) <= upper_bound
            count += 1
        end
    end
    return count
end

#########################################

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()

ext: jl


text:
@testset "Rating Kernel Tests" begin
    preference = Persa.Preference([1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0])
    @test string(preference) == "[1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0]"
    @test size(preference) == 9
    @test eltype(preference) <: Real
    @test Persa.recommendation(preference) == 4.0

    preference = Persa.Preference([1, 2, 3, 4, 5])
    @test string(preference) == "[1, 2, 3, 4, 5]"
    @test size(preference) == 5
    @test eltype(preference) <: Int
    @test Persa.recommendation(preference) == 4.0

    @test Persa.value(Persa.MissingRating()) === missing
    @test Persa.value(Persa.MissingRating{Int}()) === missing

    @test Persa.value(Persa.Rating(4, preference)) == 4
    @test Persa.value(Persa.Rating(6, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.Rating(-1, preference)) == Persa.minimum(preference)

    @test Persa.value(Persa.PredictRating(4, preference)) == 4
    @test Persa.value(Persa.PredictRating(6, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.PredictRating(-1, preference)) == Persa.minimum(preference)
    @test Persa.value(Persa.PredictRating(5.5, preference)) == Persa.maximum(preference)
    @test Persa.value(Persa.PredictRating(-1.1, preference)) == Persa.minimum(preference)
    @test Persa.value(Persa.PredictRating(4.1, preference)) == 4.1

    predict = Persa.PredictRating(4.5, preference)
    label = Persa.Rating(4, preference)

    rating = Persa.Rating(4, preference)
    userPreference = Persa.UserPreference(1, 1, rating)
    @test userPreference[1] == 1
    @test userPreference[2] == 1
    @test userPreference[3] == 4
    @test isnan(userPreference[3]) == false
    @test ismissing(userPreference[3]) == false

    @test string(rating) == "4"

    io = IOBuffer()
    show(io, rating)
    str = String(take!(io))
    @test str == "Rating: 4"

    @test string(userPreference) == "(user: 1, item: 1, rating: 4)"

    rating = Persa.PredictRating(4.5, preference)
    userPreference = Persa.UserPreference(2, 3, rating)
    @test userPreference[1] == 2
    @test userPreference[2] == 3
    @test userPreference[3] == 4.5
    @test !isnan(userPreference[3]) == true
    @test !ismissing(userPreference[3]) == true

    io = IOBuffer()
    show(io, rating)
    str = String(take!(io))
    @test str == "Rating: 4.5 (4)"

    @test string(rating) == "4.5 (4)"
    @test string(userPreference) == "(user: 2, item: 3, rating: 4.5 (4))"

    rating = Persa.MissingRating()
    userPreference = Persa.UserPreference(2, 3, rating)
    @test userPreference[1] == 2
    @test userPreference[2] == 3
    @test isnan(userPreference[3]) == true
    @test ismissing(userPreference[3]) == true

    io = IOBuffer()
    show(io, rating)
    str = String(take!(io))
    @test str == "Rating: missing"

    @test string(rating) == "missing"
    @test string(userPreference) == "(user: 2, item: 3, rating: missing)"

    @testset "Operations with Rating" begin
        @testset "Operations with Missing Rating" begin
            rating = Persa.MissingRating{Int}()

            predictRating = Persa.PredictRating(4.2, preference)
            databaseRating = Persa.Rating(4, preference)

            @test_throws Exception  rating + 1
            @test_throws Exception  1 + rating
            @test_throws Exception  rating + 1.1
            @test_throws Exception  1.1 + rating
            @test_throws Exception  predictRating + rating
            @test_throws Exception  databaseRating + rating
            @test_throws Exception  rating + predictRating
            @test_throws Exception  rating + databaseRating

            @test_throws Exception  rating - 1
            @test_throws Exception  1 - rating
            @test_throws Exception  rating - 1.1
            @test_throws Exception  1.1 - rating
            @test_throws Exception  predictRating - rating
            @test_throws Exception  databaseRating - rating
            @test_throws Exception  rating - predictRating
            @test_throws Exception  rating - databaseRating

            @test_throws Exception rating * 1
            @test_throws Exception 1 * rating
            @test_throws Exception rating * 1.1
            @test_throws Exception 1.1 * rating
            @test_throws Exception  predictRating * rating
            @test_throws Exception  databaseRating * rating
            @test_throws Exception  rating * predictRating
            @test_throws Exception  rating * databaseRating

            @test_throws Exception rating / 1
            @test_throws Exception 1 / rating
            @test_throws Exception rating / 1.1
            @test_throws Exception 1.1 / rating
            @test_throws Exception  predictRating / rating
            @test_throws Exception  databaseRating / rating
            @test_throws Exception  rating / predictRating
            @test_throws Exception  rating / databaseRating

            @test_throws Exception  rating == 1
            @test_throws Exception  1 == rating
            @test_throws Exception  rating == 1.1
            @test_throws Exception  1.1 == rating
            @test_throws Exception  predictRating == rating
            @test_throws Exception  databaseRating == rating
            @test_throws Exception  rating == predictRating
            @test_throws Exception  rating == databaseRating

            @test_throws Exception  rating != 1
            @test_throws Exception  1 != rating
            @test_throws Exception  rating != 1.1
            @test_throws Exception  1.1 != rating
            @test_throws Exception  predictRating != rating
            @test_throws Exception  databaseRating != rating
            @test_throws Exception  rating != predictRating
            @test_throws Exception  rating != databaseRating

            @test_throws Exception  rating < 1
            @test_throws Exception  1 < rating
            @test_throws Exception  rating < 1.1
            @test_throws Exception  1.1 < rating
            @test_throws Exception  predictRating < rating
            @test_throws Exception  databaseRating < rating
            @test_throws Exception  rating < predictRating
            @test_throws Exception  rating < databaseRating

            @test_throws Exception  rating > 1
            @test_throws Exception  1 > rating
            @test_throws Exception  rating > 1.1
            @test_throws Exception  1.1 > rating
            @test_throws Exception  predictRating > rating
            @test_throws Exception  databaseRating > rating
            @test_throws Exception  rating > predictRating
            @test_throws Exception  rating > databaseRating

            @test_throws Exception  rating <= 1
            @test_throws Exception  1 <= rating
            @test_throws Exception  rating <= 1.1
            @test_throws Exception  1.1 <= rating
            @test_throws Exception  predictRating <= rating
            @test_throws Exception  databaseRating <= rating
            @test_throws Exception  rating <= predictRating
            @test_throws Exception  rating <= databaseRating

            @test_throws Exception  rating >= 1
            @test_throws Exception  1 >= rating
            @test_throws Exception  rating >= 1.1
            @test_throws Exception  1.1 >= rating
            @test_throws Exception  predictRating >= rating
            @test_throws Exception  databaseRating >= rating
            @test_throws Exception  rating >= predictRating
            @test_throws Exception  rating >= databaseRating

            @test_throws Exception  rating + rating
            @test_throws Exception  rating - rating
            @test_throws Exception  rating * rating
            @test_throws Exception  rating / rating

            @test_throws Exception  rating == rating
            @test_throws Exception  rating != rating
            @test_throws Exception  rating > rating
            @test_throws Exception  rating < rating
            @test_throws Exception  rating >= rating
            @test_throws Exception  rating <= rating
        end

        @testset "Operations with Rating" begin
            preference = Persa.Preference([1, 2, 3, 4, 5])
            value = 1
            increment = 0.2
            rating = Persa.Rating(value, preference)

            functions = [maximum, minimum]
            operations = [+, -, *, /]

            for o in operations
                for f in functions
                    @test o(rating, f(preference)) == o(value, f(preference))
                    @test o(f(preference), rating) == o(f(preference), value)
                end

                @test o(rating, rating) == o(value, value)
            end

            @test (rating > value) == false
            @test (rating >= value) == true
            @test (rating < value) == false
            @test (rating <= value) == true
            @test (rating == value) == true
            @test (rating != value) == false

            @test (value > rating) == false
            @test (value >= rating) == true
            @test (value < rating) == false
            @test (value <= rating) == true
            @test (value == rating) == true
            @test (value != rating) == false

            @test (rating > rating) == false
            @test (rating < rating) == false
            @test (rating >= rating) == true
            @test (rating <= rating) == true
            @test (rating == rating) == true
            @test (rating != rating) == false
        end

        @testset "Operations with Ratings Array" begin
            preference = Persa.Preference([1, 2, 3, 4, 5])
            values = [1, 2, 3, 4, 5]
            ratings = Persa.convert(values, preference)
            i = 1

            @test length(values) == length(ratings)
            for i = 1:length(values)
                @test values[i] == ratings[i]
            end
        end

        @testset "Rerating" begin
            preference = Persa.Preference([1, 2, 3, 4, 5])
            rating = Persa.Rating(1, preference)
            rerating = Persa.Rating(2, preference)

            @test rating[1] == 1
            @test rerating[1] == 2
            @test_throws BoundsError rating[2]
            @test_throws BoundsError rerating[2]

            Persa.rerating(rerating, rating)

            @test rerating[1] == 2
            @test rerating[2] == 1
            @test_throws BoundsError rerating[3]
        end
    end
end


text:
# This file was generated, do not modify it. # hide
Xs, ys = source.((X, y));

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
function chebyshev_polynomial_derivative(order::S,x::T) where {T<:AbstractFloat,S<:Integer}

  polynomial    = Array{T}(undef,1,order+1)
  poly_deriv    = Array{T}(undef,1,order+1)
  polynomial[1] = one(T)
  poly_deriv[1] = zero(T)

  for i = 2:order+1
    if i == 2
      polynomial[i] = x
      poly_deriv[i] = one(T)
    else
      polynomial[i] = 2*x*polynomial[i-1]-polynomial[i-2]
      poly_deriv[i] = 2*polynomial[i-1]+2*x*poly_deriv[i-1]-poly_deriv[i-2]
    end
  end

  return poly_deriv

end

function chebyshev_polynomial_derivative(order::S,x::Array{T,1}) where {T<:AbstractFloat,S<:Integer}

  poly_deriv = Array{T,2}(undef,length(x),order+1)

  for i = 1:length(x)
    poly_deriv[i,:] = chebyshev_polynomial_derivative(order,x[i])
  end

  return poly_deriv

end


text:
include("./rest.jl")

"""
    get_dataset(tableName::String)

Returns the entire dataset.
It is not recommended to retrieve datasets with more than 100k rows using this method.
For large datasets, please use the 'space_time' method and retrieve the data in smaller chunks.
Note that this method does not return the dataset metadata.
Use the 'get_dataset_metadata' method to get the dataset metadata.

Examples
≡≡≡≡≡≡≡≡≡≡
get_dataset("tblHOT_LAVA")

"""
function get_dataset(tableName::String)
    datasetID = get_dataset_ID(tableName)
    maxRow = 2000000
    api = API()
    status, response = query(api, "SELECT JSON_stats FROM tblDataset_Stats WHERE Dataset_ID=$datasetID")
    jDict = JSON.parse(response.JSON_stats[1])
    rows = convert(Int, jDict["lat"]["count"])
    rows > maxRow && error("""
        The requested dataset has $rows records.
        It is not recommended to retrieve datasets with more than $maxRow rows using this method.
        For large datasets, please use the 'space_time' method and retrieve the data in smaller chunks.
        """)
    status, response = query(api, "SELECT * FROM $tableName")
    return response
end

"""
    _subset(spName::String,
            table::String,
            variable::String,
            dt1::String,
            dt2::String,
            lat1::Float64,
            lat2::Float64,
            lon1::Float64,
            lon2::Float64,
            depth1::Float64,
            depth2::Float64
           )

Returns a subset of data according to space-time constraints.
This method is meant to be used internally.

Examples
≡≡≡≡≡≡≡≡≡≡
_subset("uspSpaceTime", "tblAltimetry_REP", "sla", "2016-01-01", "2016-01-01", 30., 31., -160., -159., 0., 0.);

"""
function _subset(spName::String,
                table::String,
                variable::String,
                dt1::String,
                dt2::String,
                lat1::Float64,
                lat2::Float64,
                lon1::Float64,
                lon2::Float64,
                depth1::Float64,
                depth2::Float64
                )
    statement = "EXEC $spName '$table', '$variable', '$dt1', '$dt2', $lat1, $lat2, $lon1, $lon2, $depth1, $depth2";
    api = API();
    status, response = query(api, statement);
    return response
end

"""
    _interval_to_uspName(interval::String)

Returns a timeseries-based stored procedure name according to the specified interval.
This method is meant to be used internally.

Examples
≡≡≡≡≡≡≡≡≡≡
_interval_to_uspName("week");

"""
function _interval_to_uspName(interval::String)
    usp = "";
    if interval == ""
        usp = "uspTimeSeries"
    elseif interval in ["w", "week", "weekly"]
        usp = "uspWeekly"
    elseif interval in ["m", "month", "monthly"]
        usp = "uspMonthly"
    elseif interval in ["q", "s", "season", "seasonal", "seasonality", "quarterly"]
        usp = "uspQuarterly"
    elseif interval in ["a", "y", "year", "yearly", "annual"]
        usp = "uspAnnual"
    else
        error("Invalid interval: $interval")
    end
    return usp
end

"""
    space_time(;
               table::String,
               variable::String,
               dt1::String,
               dt2::String,
               lat1::Number,
               lat2::Number,
               lon1::Number,
               lon2::Number,
               depth1::Number,
               depth2::Number
              )

Returns a subset of data according to the specified space-time constraints (dt1, dt2, lat1, lat2, lon1, lon2, depth1, depth2).
The results are ordered by time, lat, lon, and depth (if exists), respectively.

Examples
≡≡≡≡≡≡≡≡≡≡
space_time(
           table="tblArgoMerge_REP",
           variable="argo_merge_salinity_adj",
           dt1="2015-05-01",
           dt2="2015-05-30",
           lat1=28.1,
           lat2=35.4,
           lon1=-71.3,
           lon2=-50,
           depth1=0,
           depth2=100
           )

"""
function space_time(;
                    table::String,
                    variable::String,
                    dt1::String,
                    dt2::String,
                    lat1::Number,
                    lat2::Number,
                    lon1::Number,
                    lon2::Number,
                    depth1::Number,
                    depth2::Number
                    )
    return _subset("uspSpaceTime",
                   table,
                   variable,
                   dt1,
                   dt2,
                   float(lat1),
                   float(lat2),
                   float(lon1),
                   float(lon2),
                   float(depth1),
                   float(depth2)
                   )
end

"""
    time_series(;
               table::String,
               variable::String,
               dt1::String,
               dt2::String,
               lat1::Number,
               lat2::Number,
               lon1::Number,
               lon2::Number,
               depth1::Number,
               depth2::Number
              )

Returns a subset of data according to the specified space-time constraints (dt1, dt2, lat1, lat2, lon1, lon2, depth1, depth2).
The returned data subset is aggregated by time: at each time interval, the mean and standard deviation of the variable values within the space-time constraints are computed.
The sequence of these values construct the timeseries.
The timeseries data can be binned weekly, monthly, quarterly, or annually, if the interval parameter is set (this feature is not applicable to climatological datasets).
The resulted timeseries is returned in the form of a dataframe ordered by time.

Examples
≡≡≡≡≡≡≡≡≡≡
time_series(
           table="tblAltimetry_REP",
           variable="adt",
           dt1="2016-01-01",
           dt2="2018-01-01",
           lat1=33,
           lat2=35,
           lon1=-160,
           lon2=-159,
           depth1=0,
           depth2=0,
           interval="monthly"
           )

"""
function time_series(;
                    table::String,
                    variable::String,
                    dt1::String,
                    dt2::String,
                    lat1::Number,
                    lat2::Number,
                    lon1::Number,
                    lon2::Number,
                    depth1::Number,
                    depth2::Number,
                    interval::String=""
                    )
    usp = _interval_to_uspName(interval);
    if (usp != "uspTimeSeries") && (is_climatology(table))
        error(
              """
              Custom binning (monthly, weekly, ...) is not suppoerted for climatological data sets.
              Table $table represents a climatological data set.
              """
            );
    end
    return _subset(usp,
                   table,
                   variable,
                   dt1,
                   dt2,
                   float(lat1),
                   float(lat2),
                   float(lon1),
                   float(lon2),
                   float(depth1),
                   float(depth2)
                   )
end

"""
    depth_profile(;
                 table::String,
                 variable::String,
                 dt1::String,
                 dt2::String,
                 lat1::Number,
                 lat2::Number,
                 lon1::Number,
                 lon2::Number,
                 depth1::Number,
                 depth2::Number
                )

Returns a subset of data according to the specified space-time constraints (dt1, dt2, lat1, lat2, lon1, lon2, depth1, depth2).
The returned data subset is aggregated by depth: at each depth level the mean and standard deviation of the variable values within the space-time constraints are computed.
The sequence of these values construct the depth profile.
The resulting depth profile is returned in the form of a Pandas dataframe ordered by depth.

Examples
≡≡≡≡≡≡≡≡≡≡
depth_profile(
             table="tblPisces_NRT",
             variable="CHL",
             dt1="2016-04-30",
             dt2="2016-04-30",
             lat1=20,
             lat2=24,
             lon1=-170,
             lon2=-150,
             depth1=0,
             depth2=500
             )

"""
function depth_profile(;
                      table::String,
                      variable::String,
                      dt1::String,
                      dt2::String,
                      lat1::Number,
                      lat2::Number,
                      lon1::Number,
                      lon2::Number,
                      depth1::Number,
                      depth2::Number
                      )
    return _subset("uspDepthProfile",
                   table,
                   variable,
                   dt1,
                   dt2,
                   float(lat1),
                   float(lat2),
                   float(lon1),
                   float(lon2),
                   float(depth1),
                   float(depth2)
                   )
end

"""
    section(;
            table::String,
            variable::String,
            dt1::String,
            dt2::String,
            lat1::Number,
            lat2::Number,
            lon1::Number,
            lon2::Number,
            depth1::Number,
            depth2::Number
           )

Returns a subset of data according to the specified space-time constraints.
The results are ordered by time, lat, lon, and depth.

Examples
≡≡≡≡≡≡≡≡≡≡
section(
        table="tblPisces_NRT",
        variable="NO3",
        dt1="2016-04-30",
        dt2="2016-04-30",
        lat1=10,
        lat2=50,
        lon1=-159,
        lon2=-158,
        depth1=0,
        depth2=500
        )

"""
function section(;
                table::String,
                variable::String,
                dt1::String,
                dt2::String,
                lat1::Number,
                lat2::Number,
                lon1::Number,
                lon2::Number,
                depth1::Number,
                depth2::Number
                )
    return _subset("uspSectionMap",
                   table,
                   variable,
                   dt1,
                   dt2,
                   float(lat1),
                   float(lat2),
                   float(lon1),
                   float(lon2),
                   float(depth1),
                   float(depth2)
                   )
end

"""
    match(;
          sourceTable::String,
          sourceVariable::String,
          targetTables::Array,
          targetVariables::Array,
          dt1::String,
          dt2::String,
          lat1::Number,
          lat2::Number,
          lon1::Number,
          lon2::Number,
          depth1::Number,
          depth2::Number,
          timeTolerance::Array,
          latTolerance::Array,
          lonTolerance::Array,
          depthTolerance::Array
          )

Colocalizes the source variable (from source table) with the target variable (from target table).
The tolerance parameters set the matching boundaries between the source and target data sets.
Returns a dataframe containing the source variable joined with the target variable.

# Arguments
- `sourceTable::String`: table name of the source data set.
- `sourceVariable::String`: the source variable. The target variables are matched (colocalized) with this variable.
- `targetTables::Array`: table names of the target data sets to be matched with the source data.
- `targetVariables::Array`: variable names to be matched with the source variable.
- `dt1::String`: start date or datetime.
- `dt2::String`: end date or datetime.
- `lat1::Number`: start latitude [degree N]. This parameter sets the lower bound of the meridional cut. Note latitude ranges from -90 to 90.
- `lat2::Number`: end latitude [degree N]. This parameter sets the upper bound of the meridional cut. Note latitude ranges from -90 to 90.
- `lon1::Number`: start longitude [degree E]. This parameter sets the lower bound of the zonal cut. Note latitude ranges from -180 to 180.
- `lon2::Number`: end longitude [degree E]. This parameter sets the upper bound of the zonal cut. Note latitude ranges from -180 to 180.
- `depth1::Number`: start depth [m]. This parameter sets the lower bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `depth2::Number`: end depth [m]. This parameter sets the upper bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `timeTolerance::Number`: integer array of temporal tolerance values between pairs of source and target datasets. The size and order of values in this array should match those of targetTables. This parameter is in day units except when the target variable represents monthly climatology data in which case it is in month units. Note fractional values are not supported in the current version.
- `latTolerance::Number`: float array of spatial tolerance values in meridional direction [deg] between pairs of source and target data sets.
- `lonTolerance::Number`: float array of spatial tolerance values in zonal direction [deg] between pairs of source and target data sets.
- `depthTolerance::Number`: float array of spatial tolerance values in vertical direction [m] between pairs of source and target data sets.

Examples
≡≡≡≡≡≡≡≡≡≡
The source variable in this example is particulate pseudo cobalamin (Me_PseudoCobalamin_Particulate_pM) measured by
Ingalls lab during the KM1315 cruise. This variable is colocalized with one target variabele, picoprokaryote concentration,
from Darwin model.

match(
    sourceTable="tblKM1314_Cobalmins",
    sourceVariable="Me_PseudoCobalamin_Particulate_pM",
    targetTables=["tblDarwin_Phytoplankton"],
    targetVariables=["picoprokaryote"],
    dt1="2013-08-11",
    dt2="2013-09-05",
    lat1=22.25,
    lat2=450.25,
    lon1=-159.25,
    lon2=-127.75,
    depth1=-5,
    depth2=305,
    timeTolerance=[1],
    latTolerance=[0.25],
    lonTolerance=[0.25],
    depthTolerance=[5]
    );
"""
function match(;
               sourceTable::String,
               sourceVariable::String,
               targetTables::Array,
               targetVariables::Array,
               dt1::String,
               dt2::String,
               lat1::Number,
               lat2::Number,
               lon1::Number,
               lon2::Number,
               depth1::Number,
               depth2::Number,
               timeTolerance::Array,
               latTolerance::Array,
               lonTolerance::Array,
               depthTolerance::Array
               )
    return compile(
                   spname="uspMatch",
                   sourceTable=sourceTable,
                   sourceVariable=sourceVariable,
                   targetTables=targetTables,
                   targetVariables=targetVariables,
                   dt1=dt1,
                   dt2=dt2,
                   lat1=lat1,
                   lat2=lat2,
                   lon1=lon1,
                   lon2=lon2,
                   depth1=depth1,
                   depth2=depth2,
                   timeTolerance=timeTolerance,
                   latTolerance=latTolerance,
                   lonTolerance=lonTolerance,
                   depthTolerance=depthTolerance
                  )
end

"""
    along_track(;
                cruise::String,
                targetTables::Array,
                targetVariables::Array,
                depth1::Number,
                depth2::Number,
                timeTolerance::Array,
                latTolerance::Array,
                lonTolerance::Array,
                depthTolerance::Array
                )

Takes a cruise name and colocalizes the cruise track with the specified variable(s).

# Arguments
- `cruise::String`: cruise name.
- `targetTables::Array`: table names of the target data sets to be matched with the source data.
- `targetVariables::Array`: variable names to be matched with the source variable.
- `depth1::Number`: start depth [m]. This parameter sets the lower bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `depth2::Number`: end depth [m]. This parameter sets the upper bound of the vertical cut. Note depth is a positive number (it is 0 at surface and grows towards ocean floor).
- `timeTolerance::Number`: integer array of temporal tolerance values between pairs of source and target datasets. The size and order of values in this array should match those of targetTables. This parameter is in day units except when the target variable represents monthly climatology data in which case it is in month units. Note fractional values are not supported in the current version.
- `latTolerance::Number`: float array of spatial tolerance values in meridional direction [deg] between pairs of source and target data sets.
- `lonTolerance::Number`: float array of spatial tolerance values in zonal direction [deg] between pairs of source and target data sets.
- `depthTolerance::Number`: float array of spatial tolerance values in vertical direction [m] between pairs of source and target data sets.

Examples
≡≡≡≡≡≡≡≡≡≡
This example demonstrates how to colocalize the "gradients_1" cruise (official name: KOK1606) with 2 target variables:
"prochloro_abundance" from underway seaflow dataset "PO4" from Darwin climatology model.

along_track(
            cruise="gradients_1",
            targetTables=["tblSeaFlow", "tblDarwin_Nutrient_Climatology"],
            targetVariables=["prochloro_abundance", "PO4_darwin_clim"],
            depth1=0,
            depth2=5,
            timeTolerance=[0, 0],
            latTolerance=[0.01, 0.25],
            lonTolerance=[0.01, 0.25],
            depthTolerance=[0, 5]
            );
"""
function along_track(;
                    cruise::String,
                    targetTables::Array,
                    targetVariables::Array,
                    depth1::Number,
                    depth2::Number,
                    timeTolerance::Array,
                    latTolerance::Array,
                    lonTolerance::Array,
                    depthTolerance::Array
                    )
    df = cruise_bounds(cruise);
    return match(
                 sourceTable="tblCruise_Trajectory",
                 sourceVariable=string(df.ID[1]),
                 targetTables=targetTables,
                 targetVariables=targetVariables,
                 dt1=df.dt1[1],
                 dt2=df.dt2[1],
                 lat1=df.lat1[1],
                 lat2=df.lat2[1],
                 lon1=df.lon1[1],
                 lon2=df.lon2[1],
                 depth1=depth1,
                 depth2=depth2,
                 timeTolerance=timeTolerance,
                 latTolerance=latTolerance,
                 lonTolerance=lonTolerance,
                 depthTolerance=depthTolerance
                 )
end


text:
# # How to port NiLang to Zygote
using NiLang, NiLang.AD, Zygote

# Let's start from the Julia native implementation of `norm2` function.
function norm2(x::AbstractArray{T}) where T
    out = zero(T)
    for i=1:length(x)
        @inbounds out += x[i]^2
    end
    return out
end

# Zygote is able to generate correct gradients, but much slower than the original program.
using BenchmarkTools
x = randn(1000);
original_grad = norm2'(x)
@benchmark norm2'($x) seconds=1

# The orignal program is
@benchmark norm2($x) seconds=1

# Then we have the reversible implementation
@i function r_norm2(out::T, x::AbstractArray{T}) where T
    for i=1:length(x)
        @inbounds out += x[i]^2
    end
end

# The gradient generated by NiLang is much faster, which is comparable to the forward program
@benchmark (~r_norm2)(GVar($(norm2(x)), 1.0), $(GVar(x))) seconds=1

# to enjoy the speed of `NiLang` in `Zygote`, just bind the adjoint rule
Zygote.@adjoint function norm2(x::AbstractArray{T}) where T
    out = norm2(x)
    out, δy -> (grad((~r_norm2)(GVar(out, δy), GVar(x))[2]),)
end
@assert norm2'(x) ≈ original_grad

# See, much faster
@benchmark norm2'(x) seconds=1


text:
using QuadGK: quadgk
# -------------------------------------------------------------------------
# Quick version of functions (for testing)
function q_linearupperbound(xstar, gllstar, b, x, v)
    a = max(dot(-gllstar, v), 0.) + norm(x - xstar) * b
    tauprime = -a / b + sqrt((a / b)^2 + 2randexp() / b)
    (tauprime, a + b * tauprime)
end

# Benchmark dirty version (for testing)
function q_sampletaugaussian(mu,C,x,v)
    r = randexp()
    vciv = dot(v, C \ v)
    vcixmu = dot(v,C \ (x - mu))
    t = 0
    if vcixmu < 0 ## can be optimised because we now apriori
        t= -vcixmu / vciv + sqrt(2r / vciv)
    else
        t= -vcixmu / vciv + sqrt((vcixmu / vciv)^2 + 2r / vciv)
    end
    t
end
# ------------------------------------------------------------------------------

Random.seed!(123)

### test NextEvent capsule (trivial test for code cov)

tau = 0.1
dobounce = (g, v) -> -dot(g, v) < 1.0
flipindex = 4
ne = NextEvent(tau; dobounce=dobounce, flipindex=flipindex)
g  = randn(4)
v  = randn(4)
@test ne.tau == tau && ne.dobounce(g, v) == dobounce(g, v) &&
        ne.flipindex == flipindex
ne2 = NextEvent(tau)
@test ne2.tau == tau && ne2.dobounce(g, v) &&
        ne2.flipindex == -1
ne3 = NextEvent(tau; dobounce=dobounce)
@test ne3.tau == tau && ne3.dobounce(g, v) == dobounce(g, v) &&
        ne3.flipindex == -1

p = 5

gllstar = randn(p)
xstar = randn(p)
x = randn(p)
v = randn(p)
g = randn(p)
b = randn()

lb = LinearBound(gllstar, xstar, b)

## check that q_linearupperbound and nextevent_bps match.
Random.seed!(10); (a1, a2) = q_linearupperbound(gllstar, xstar, b, x, v)
acc = rand() < -dot(g, v) / a2
Random.seed!(10); bounce = nextevent_bps(lb, x, v)

@test abs(a1 - bounce.tau) <= 1e-12 && bounce.dobounce(g,v) == acc

## check that q_sampletaugaussian and nextevent_bps match.

Random.seed!(123)

p = 50
C1 = randn(p, p)
C1 *= C1'
mu = randn(p)

mvg = MvGaussianStandard(mu, C1)

x = randn(p)
v = randn(p)
g = randn(p)

Random.seed!(12); tauq = q_sampletaugaussian(mu, C1, x, v)
Random.seed!(12); bounce = nextevent_bps(mvg, x, v)

@test abs(tauq - bounce.tau) <= 1e-12 && bounce.dobounce(g, v)

## additional tests

nedef = NextEvent(0.5)
@test nedef.tau == 0.5 && nedef.dobounce(randn(), randn()) &&
        nedef.flipindex == -1

Random.seed!(555); a = nextevent_zz(mvg, x, v)
Random.seed!(555)
u1 = mvg.prec * x - mvg.precmu
u2 = mvg.prec * v
taus = zeros(p)
for i ∈ 1:p
    ai = u1[i] * v[i]
    bi = u2[i] * v[i]
    ci = ai ./ bi
    ei = max(0.0, ci) * ci
    taus[i] = -ci + sqrt(ei + 2randexp() / abs(bi))
end
tau, flipindex = findmin(taus)

@test a.tau == tau && a.flipindex == flipindex

## testing nextevent_bps for PMF

d = 3
r = 3
s = 0.1

pmfg = PMFGaussian(r, s, d)

n  = 25
xu = [randn(d) for i ∈ 1:n]
xv = [randn(d) for i ∈ 1:n]
wu = [randn(d) for i ∈ 1:n]
wv = [randn(d) for i ∈ 1:n]

Random.seed!(542)
tau = [nextevent_bps(pmfg, [xu[i]; xv[i]], [wu[i]; wv[i]]).tau for i ∈ 1:n]

Random.seed!(542)
ri = [randexp() for i ∈ 1:100]

# it must collide
p1(xu, xv, wu, wv) = PDSampler.Poly(
    [dot(xu,xv)-pmfg.r, (dot(xu, wv) + dot(xv, wu)), dot(wu, wv)])
p2(xu, xv, wu, wv) = PDSampler.Poly(
    [(dot(xu, wv) + dot(xv, wu)), 2dot(wu, wv)])
E(xu,xv,wu,wv) = p1(xu, xv, wu, wv) * p2(xu, xv, wu, wv)

χ(i) = t-> max(0.0, p1(xu[i], xv[i], wu[i], wv[i])(t) *
                    p2(xu[i], xv[i], wu[i], wv[i])(t))

### testing roots
t0(i) = -0.5(dot(xu[i], wv[i]) + dot(xv[i], wu[i])) / dot(wu[i], wv[i])
Δ(i) = t0(i)^2 - (dot(xu[i], xv[i]) - pmfg.r) / dot(wu[i], wv[i])

@test maximum(
        (Δ(i) > 0) ? maximum(map(χ(i), t0(i) .+ sqrt(Δ(i)) * [-1., 0., 1.])) :
                     χ(i)(t0(i))
        for i in 1:n) <= 1e-12

tm(i) = t0(i) - sqrt(abs(Δ(i)))
tp(i) = t0(i) + sqrt(abs(Δ(i)))

@test maximum(quadgk(χ(i), 0.0, tau[i])[1] - ri[i] for i ∈ 1:n) <= 1e-6


text:
using Test, DynamicPPL, Setfield
using DynamicPPL:
    leafcontext,
    setleafcontext,
    childcontext,
    setchildcontext,
    AbstractContext,
    NodeTrait,
    IsLeaf,
    IsParent,
    PointwiseLikelihoodContext,
    contextual_isassumption,
    ConditionContext,
    hasvalue,
    getvalue,
    hasvalue_nested,
    getvalue_nested

# Dummy context to test nested behaviors.
struct ParentContext{C<:AbstractContext} <: AbstractContext
    context::C
end
ParentContext() = ParentContext(DefaultContext())
DynamicPPL.NodeTrait(::ParentContext) = DynamicPPL.IsParent()
DynamicPPL.childcontext(context::ParentContext) = context.context
DynamicPPL.setchildcontext(::ParentContext, child) = ParentContext(child)
Base.show(io::IO, c::ParentContext) = print(io, "ParentContext(", childcontext(c), ")")

# TODO: Should we maybe put this in DPPL itself?
function Base.iterate(context::AbstractContext)
    if NodeTrait(context) isa IsLeaf
        return nothing
    end

    return context, context
end
function Base.iterate(_::AbstractContext, context::AbstractContext)
    return _iterate(NodeTrait(context), context)
end
_iterate(::IsLeaf, context) = nothing
function _iterate(::IsParent, context)
    child = childcontext(context)
    return child, child
end

Base.IteratorSize(::Type{<:AbstractContext}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{<:AbstractContext}) = Base.EltypeUnknown()

"""
    remove_prefix(vn::VarName)

Return `vn` but now with the prefix removed.
"""
function remove_prefix(vn::VarName)
    return VarName{Symbol(split(string(vn), string(DynamicPPL.PREFIX_SEPARATOR))[end])}(
        getlens(vn)
    )
end

"""
    varnames(vn::VarName, val)

Return iterator over all varnames that are represented by `vn` on `val`,
e.g. `varnames(@varname(x), rand(2))` results in an iterator over `[@varname(x[1]), @varname(x[2])]`.
"""
varnames(vn::VarName, val::Real) = [vn]
function varnames(vn::VarName, val::AbstractArray{<:Union{Real,Missing}})
    return (
        VarName(vn, getlens(vn) ∘ Setfield.IndexLens(Tuple(I))) for
        I in CartesianIndices(val)
    )
end
function varnames(vn::VarName, val::AbstractArray)
    return Iterators.flatten(
        varnames(VarName(vn, getlens(vn) ∘ Setfield.IndexLens(Tuple(I))), val[I]) for
        I in CartesianIndices(val)
    )
end

@testset "contexts.jl" begin
    child_contexts = [DefaultContext(), PriorContext(), LikelihoodContext()]

    parent_contexts = [
        ParentContext(DefaultContext()),
        SamplingContext(),
        MiniBatchContext(DefaultContext(), 0.0),
        PrefixContext{:x}(DefaultContext()),
        PointwiseLikelihoodContext(),
        ConditionContext((x=1.0,)),
        ConditionContext((x=1.0,), ParentContext(ConditionContext((y=2.0,)))),
        ConditionContext((x=1.0,), PrefixContext{:a}(ConditionContext((var"a.y"=2.0,)))),
        ConditionContext((x=[1.0, missing],)),
    ]

    contexts = vcat(child_contexts, parent_contexts)

    @testset "NodeTrait" begin
        @testset "$context" for context in contexts
            # Every `context` should have a `NodeTrait`.
            @test NodeTrait(context) isa NodeTrait
        end
    end

    @testset "leafcontext" begin
        @testset "$context" for context in child_contexts
            @test leafcontext(context) === context
        end

        @testset "$context" for context in parent_contexts
            @test NodeTrait(leafcontext(context)) isa IsLeaf
        end
    end

    @testset "setleafcontext" begin
        @testset "$context" for context in child_contexts
            # Setting to itself should return itself.
            @test setleafcontext(context, context) === context

            # Setting to a different context should return that context.
            new_leaf = context isa DefaultContext ? PriorContext() : DefaultContext()
            @test setleafcontext(context, new_leaf) === new_leaf

            # Also works for parent contexts.
            new_leaf = ParentContext(context)
            @test setleafcontext(context, new_leaf) === new_leaf
        end

        @testset "$context" for context in parent_contexts
            # Leaf contexts.
            new_leaf =
                leafcontext(context) isa DefaultContext ? PriorContext() : DefaultContext()
            @test leafcontext(setleafcontext(context, new_leaf)) === new_leaf

            # Setting parent contexts as "leaf" means that the new leaf should be
            # the leaf of the parent context we just set as the leaf.
            new_leaf = ParentContext((
                leafcontext(context) isa DefaultContext ? PriorContext() : DefaultContext()
            ))
            @test leafcontext(setleafcontext(context, new_leaf)) === leafcontext(new_leaf)
        end
    end

    # `IsParent` interface.
    @testset "childcontext" begin
        @testset "$context" for context in parent_contexts
            @test childcontext(context) isa AbstractContext
        end
    end

    @testset "setchildcontext" begin
        @testset "nested contexts" begin
            # Both of the following should result in the same context.
            context1 = ParentContext(ParentContext(ParentContext()))
            context2 = setchildcontext(
                ParentContext(), setchildcontext(ParentContext(), ParentContext())
            )
            @test context1 === context2
        end

        @testset "$context" for context in parent_contexts
            # Setting the child context to a leaf should now change the `leafcontext` accordingly.
            new_leaf =
                leafcontext(context) isa DefaultContext ? PriorContext() : DefaultContext()
            new_context = setchildcontext(context, new_leaf)
            @test childcontext(new_context) === leafcontext(new_context) === new_leaf
        end
    end

    @testset "contextual_isassumption" begin
        @testset "$context" for context in contexts
            # Any `context` should return `true` by default.
            @test contextual_isassumption(context, VarName{gensym(:x)}())

            if any(Base.Fix2(isa, ConditionContext), context)
                # We have a `ConditionContext` among us.
                # Let's first extract the conditioned variables.
                conditioned_values = DynamicPPL.conditioned(context)

                for (sym, val) in pairs(conditioned_values)
                    vn = VarName{sym}()

                    # We need to drop the prefix of `var` since in `contextual_isassumption`
                    # it will be threaded through the `PrefixContext` before it reaches
                    # `ConditionContext` with the conditioned variable.
                    vn_without_prefix = remove_prefix(vn)

                    # Let's check elementwise.
                    for vn_child in varnames(vn_without_prefix, val)
                        if get(val, getlens(vn_child)) === missing
                            @test contextual_isassumption(context, vn_child)
                        else
                            @test !contextual_isassumption(context, vn_child)
                        end
                    end
                end
            end
        end
    end

    @testset "getvalue_nested & hasvalue_nested" begin
        @testset "$context" for context in contexts
            fake_vn = VarName{gensym(:x)}()
            @test !hasvalue_nested(context, fake_vn)
            @test_throws ErrorException getvalue_nested(context, fake_vn)

            if any(Base.Fix2(isa, ConditionContext), context)
                # `ConditionContext` specific.

                # Let's first extract the conditioned variables.
                conditioned_values = DynamicPPL.conditioned(context)

                for (sym, val) in pairs(conditioned_values)
                    vn = VarName{sym}()

                    # We need to drop the prefix of `var` since in `contextual_isassumption`
                    # it will be threaded through the `PrefixContext` before it reaches
                    # `ConditionContext` with the conditioned variable.
                    vn_without_prefix = remove_prefix(vn)

                    for vn_child in varnames(vn_without_prefix, val)
                        # `vn_child` should be in `context`.
                        @test hasvalue_nested(context, vn_child)
                        # Value should be the same as extracted above.
                        @test getvalue_nested(context, vn_child) ===
                              get(val, getlens(vn_child))
                    end
                end
            end
        end
    end

    @testset "Evaluation" begin
        @testset "$context" for context in contexts
            # Just making sure that we can actually sample with each of the contexts.
            @test (gdemo_default(SamplingContext(context)); true)
        end
    end

    @testset "PrefixContext" begin
        ctx = @inferred PrefixContext{:f}(
            PrefixContext{:e}(
                PrefixContext{:d}(
                    PrefixContext{:c}(
                        PrefixContext{:b}(PrefixContext{:a}(DefaultContext()))
                    ),
                ),
            ),
        )
        vn = VarName{:x}()
        vn_prefixed = @inferred DynamicPPL.prefix(ctx, vn)
        @test DynamicPPL.getsym(vn_prefixed) == Symbol("a.b.c.d.e.f.x")
        @test getlens(vn_prefixed) === getlens(vn)

        vn = VarName{:x}(((1,),))
        vn_prefixed = @inferred DynamicPPL.prefix(ctx, vn)
        @test DynamicPPL.getsym(vn_prefixed) == Symbol("a.b.c.d.e.f.x")
        @test getlens(vn_prefixed) === getlens(vn)
    end
end


text:
function charybdis_deck(cur_T_bar::AbstractSymbol=symbols(:T_bar); cur_kwargs...)
  cur_dict = Dict(
    :deck => :charybdis,
    :T_bar => cur_T_bar,
    :is_pulsed => false,
    :H => 1.7,
    :Q => 25.0,
    :epsilon => 0.3,
    :kappa_95 => 1.8,
    :delta_95 => 0.35,
    :nu_n => 0.4,
    :nu_T => 1.1,
    :Z_eff => 1.75,
    :f_D => 0.9,
    :A => 2.5,
    :l_i => 0.557854,
    :N_G => 0.9,
    :eta_T => 0.4,
    :eta_RF => 0.5,
    :tau_FT => 1.6e9,
    :C_ejima => 0.3,
    :max_beta_N => 0.028,
    :max_q_95 => 2.5,
    :max_P_E => 3000.0,
    :max_P_W => 3.0,
    :max_q_DV => 10.0,
    :B_CS => 12.0,
    :sigma_CS => 600,
    :sigma_TF => 88888888,
    :J_CS => 50,
    :J_TF => 88888888,
    :k => 0.4,
    :eta_CD => 0.29282,
    :wave_theta => 135.0,
    :eta_LH => 0.6
  )

  merge!(cur_dict, Dict(cur_kwargs))

  BaseReactor(cur_dict)
end


text:
# test RESTful API first
#include("test_restapi_read.jl")
#include("test_restapi_write.jl")

using Base.Test
using BOSSArrays

ba = BOSSArray( collectionName  = "jingpengw_test",
                experimentName  = "test",
                channelName     = "em")

a = rand(UInt8, 200,200,1)

# the boss will overwrite zeros in the database, so make sure that all the value was overwritten by avoid zeros.
for i in eachindex(a)
    if a[i]== UInt8(0)
        a[i] = UInt8(1)
    end
end

@testset "test IO of BOSSArray" begin
    ba[10001:10200, 10001:10200, 101] = a

    b = ba[10001:10200, 10001:10200, 101]

    #@show a
    #@show b

    @test all(a.==b)
end
# @show arr

# using Images
# # img = Image(arr)
# using ImageView
# using FixedPointNumbers
# ImageView.imshow( reinterpret(N0f8, arr) )
#
# # ImageView.imshow(rand(UInt8, 200,200))
#
#
# sleep(60)


text:
import StatsBase

x, y = rand(500), rand(500)

####################
# Single binnings
####################
single_binning = RectangularBinning(0.5)

@test causality(x, y, TransferOperatorGridTest(binning = single_binning, ηs = 5)) isa Array{Float64,0}
@test causality(x, y, TransferOperatorGridTest(binning = single_binning, ηs = 2:10)) isa Vector

####################################################################################################
# Multiple binnings (summarised to a single value per η using the `binning_summary_statistic`).
####################################################################################################
multiple_binnings = [RectangularBinning(n) for n in 3:5]

causality(x, y, TransferOperatorGridTest(binning = multiple_binnings,
    binning_summary_statistic = StatsBase.mean, ηs = 1)) isa Array{Float64,0}

causality(x, y, TransferOperatorGridTest(binning = multiple_binnings,
    binning_summary_statistic = StatsBase.mean, ηs = 2:10)) isa Vector


text:
using PyTensorFlow
using PyCall
using LinearAlgebra
using PyPlot
using DelimitedFiles
using Random
Random.seed!(233)

np = pyimport("numpy")
include("poisson_op.jl")
include("laplacian_op.jl")
include("sat_op.jl")

const K_CONST =  9.869232667160130e-16 * 86400 * 1e3

mutable struct Ctx
  m; n; h; NT; Δt; Z; X; ρw; ρo;
  μw; μo; K; g; ϕ; qw; qo; sw0
end

function tfCtxGen(m,n,h,NT,Δt,Z,X,ρw,ρo,μw,μo,K,g,ϕ,qw,qo,sw0,ifTrue)
  tf_h = constant(h)
  # tf_NT = constant(NT)
  tf_Δt = constant(Δt)
  tf_Z = constant(Z)
  tf_X= constant(X)
  tf_ρw = constant(ρw)
  tf_ρo = constant(ρo)
  tf_μw = constant(μw)
  tf_μo = constant(μo)
  # tf_K = isa(K,Array) ? Variable(K) : K
  if ifTrue
    tf_K = constant(K)
  else
    tf_K = Variable(K)
  end
  tf_g = constant(g)
  # tf_ϕ = Variable(ϕ)
  tf_ϕ = constant(ϕ)
  tf_qw = constant(qw)
  tf_qo = constant(qo)
  tf_sw0 = constant(sw0)
  return Ctx(m,n,tf_h,NT,tf_Δt,tf_Z,tf_X,tf_ρw,tf_ρo,tf_μw,tf_μo,tf_K,tf_g,tf_ϕ,tf_qw,tf_qo,tf_sw0)
end

function Krw(Sw)
    return Sw ^ 1.5
end

function Kro(So)
    return So ^1.5
end

function ave_normal(quantity, m, n)
    aa = sum(quantity)
    return aa/(m*n)
end

# variables : sw, u, v, p
# (time dependent) parameters: qw, qo, ϕ
function onestep(sw, p, m, n, h, Δt, Z, ρw, ρo, μw, μo, K, g, ϕ, qw, qo)
    # step 1: update p
    # λw = Krw(sw)/μw
    # λo = Kro(1-sw)/μo
    λw = sw.*sw/μw
    λo = (1-sw).*(1-sw)/μo
    λ = λw + λo
    q = qw + qo + λw/(λo+1e-16).*qo
    # q = qw + qo
    potential_c = (ρw - ρo)*g .* Z

    # Step 1: implicit potential
    Θ = upwlap_op(K * K_CONST, λo, potential_c, h, constant(0.0))

    load_normal = (Θ+q/ALPHA) - ave_normal(Θ+q/ALPHA, m, n)

    # p = poisson_op(λ.*K* K_CONST, load_normal, h, constant(0.0), constant(1))
    p = upwps_op(K * K_CONST, λ, load_normal, p, h, constant(0.0), constant(0)) # potential p = pw - ρw*g*h

    # step 2: implicit transport
    sw = sat_op(sw, p, K * K_CONST, ϕ, qw, qo, μw, μo, sw, Δt, h)
    return sw, p
end

"""
impes(tf_ctx)
Solve the two phase flow equation.
`qw` and `qo` -- `NT x m x n` numerical array, `qw[i,:,:]` the corresponding value of qw at i*Δt
`sw0` and `p0` -- initial value for `sw` and `p`. `m x n` numerical array.
"""
function imseq(tf_ctx)
    ta_sw, ta_p = TensorArray(NT+1), TensorArray(NT+1)
    ta_sw = write(ta_sw, 1, tf_ctx.sw0)
    ta_p = write(ta_p, 1, constant(zeros(tf_ctx.m, tf_ctx.n)))
    i = constant(1, dtype=Int32)
    function condition(i, tas...)
        i <= tf_ctx.NT
    end
    function body(i, tas...)
        ta_sw, ta_p = tas
        sw, p = onestep(read(ta_sw, i), read(ta_p, i), tf_ctx.m, tf_ctx.n, tf_ctx.h, tf_ctx.Δt, tf_ctx.Z, tf_ctx.ρw, tf_ctx.ρo, tf_ctx.μw, tf_ctx.μo, tf_ctx.K, tf_ctx.g, tf_ctx.ϕ, tf_ctx.qw[i], tf_ctx.qo[i])
        ta_sw = write(ta_sw, i+1, sw)
        ta_p = write(ta_p, i+1, p)
        i+1, ta_sw, ta_p
    end

    _, ta_sw, ta_p = while_loop(condition, body, [i; ta_sw; ta_p;])
    out_sw, out_p = stack(ta_sw), stack(ta_p)
end

function vis(val, args...;kwargs...)
    close("all")
    ns = Int64.(round.(LinRange(1,size(val,1),9)))
    for i = 1:9
        subplot(330+i)
        imshow(val[ns[i],:,:], args...;kwargs...)
        colorbar()
    end
end


text:
using FMM3D
using Test
using Random
using LinearAlgebra

@testset "testing Helmholtz utilities" begin

    n = 1000
    nt = 1100
    sources = randn(3,n)
    targets = randn(3,nt)
    charges = randn(n) + im*randn(n)
    dipvecs = randn(3,n) + im*randn(3,n)
    zk = 2.1

    vals1c = h3ddir(zk,sources,targets,charges=charges,
                    pgt=2,thresh=1e-16)

    vals1cs = h3ddir(zk,sources,sources,charges=charges,
                    pgt=2,thresh=1e-16)

    vals1d = h3ddir(zk,sources,targets,dipvecs=dipvecs,
                    pgt=2,thresh=1e-16)

    vals1ds = h3ddir(zk,sources,sources,dipvecs=dipvecs,
                    pgt=2,thresh=1e-16)

    vals1dc = h3ddir(zk,sources,targets,charges=charges,
                     dipvecs=dipvecs,pgt=2,thresh=1e-16)

    vals1dcs = h3ddir(zk,sources,sources,charges=charges,
                     dipvecs=dipvecs,pgt=2,thresh=1e-16)

    eps = 1e-12
    vals2c = hfmm3d(eps,zk,sources,targets=targets,charges=charges,
                    pg=2,pgt=2)
    vals2d = hfmm3d(eps,zk,sources,targets=targets,dipvecs=dipvecs,
                    pg=2,pgt=2)
    vals2dc = hfmm3d(eps,zk,sources,targets=targets,charges=charges,
                     dipvecs=dipvecs,pg=2,pgt=2)

    @test norm(vals2c.pot-vals1cs.pottarg)/norm(vals1cs.pottarg) < eps
    @test norm(vals2c.grad-vals1cs.gradtarg)/norm(vals1cs.gradtarg) < eps
    @test norm(vals2c.pottarg-vals1c.pottarg)/norm(vals1c.pottarg) < eps
    @test norm(vals2c.gradtarg-vals1c.gradtarg)/norm(vals1c.gradtarg) < eps

    @test norm(vals2d.pot-vals1ds.pottarg)/norm(vals1ds.pottarg) < eps
    @test norm(vals2d.grad-vals1ds.gradtarg)/norm(vals1ds.gradtarg) < eps
    @test norm(vals2d.pottarg-vals1d.pottarg)/norm(vals1d.pottarg) < eps
    @test norm(vals2d.gradtarg-vals1d.gradtarg)/norm(vals1d.gradtarg) < eps

    @test norm(vals2dc.pot-vals1dcs.pottarg)/norm(vals1dcs.pottarg) < eps
    @test norm(vals2dc.grad-vals1dcs.gradtarg)/norm(vals1dcs.gradtarg) < eps
    @test norm(vals2dc.pottarg-vals1dc.pottarg)/norm(vals1dc.pottarg) < eps
    @test norm(vals2dc.gradtarg-vals1dc.gradtarg)/norm(vals1dc.gradtarg) < eps

end

@testset "testing Laplace utilities" begin

    n = 1000
    nt = 1100
    sources = randn(3,n)
    targets = randn(3,nt)
    charges = randn(n)
    dipvecs = randn(3,n)

    vals1c = l3ddir(sources,targets,charges=charges,
                    pgt=2,thresh=1e-16)

    vals1cs = l3ddir(sources,sources,charges=charges,
                    pgt=2,thresh=1e-16)

    vals1d = l3ddir(sources,targets,dipvecs=dipvecs,
                    pgt=2,thresh=1e-16)

    vals1ds = l3ddir(sources,sources,dipvecs=dipvecs,
                    pgt=2,thresh=1e-16)

    vals1dc = l3ddir(sources,targets,charges=charges,
                     dipvecs=dipvecs,pgt=2,thresh=1e-16)

    vals1dcs = l3ddir(sources,sources,charges=charges,
                     dipvecs=dipvecs,pgt=2,thresh=1e-16)

    eps = 1e-12
    vals2c = lfmm3d(eps,sources,targets=targets,charges=charges,
                    pg=2,pgt=2)
    vals2d = lfmm3d(eps,sources,targets=targets,dipvecs=dipvecs,
                    pg=2,pgt=2)
    vals2dc = lfmm3d(eps,sources,targets=targets,charges=charges,
                     dipvecs=dipvecs,pg=2,pgt=2)

    @test norm(vals2c.pot-vals1cs.pottarg)/norm(vals1cs.pottarg) < eps
    @test norm(vals2c.grad-vals1cs.gradtarg)/norm(vals1cs.gradtarg) < eps
    @test norm(vals2c.pottarg-vals1c.pottarg)/norm(vals1c.pottarg) < eps
    @test norm(vals2c.gradtarg-vals1c.gradtarg)/norm(vals1c.gradtarg) < eps

    @test norm(vals2d.pot-vals1ds.pottarg)/norm(vals1ds.pottarg) < eps
    @test norm(vals2d.grad-vals1ds.gradtarg)/norm(vals1ds.gradtarg) < eps
    @test norm(vals2d.pottarg-vals1d.pottarg)/norm(vals1d.pottarg) < eps
    @test norm(vals2d.gradtarg-vals1d.gradtarg)/norm(vals1d.gradtarg) < eps

    @test norm(vals2dc.pot-vals1dcs.pottarg)/norm(vals1dcs.pottarg) < eps
    @test norm(vals2dc.grad-vals1dcs.gradtarg)/norm(vals1dcs.gradtarg) < eps
    @test norm(vals2dc.pottarg-vals1dc.pottarg)/norm(vals1dc.pottarg) < eps
    @test norm(vals2dc.gradtarg-vals1dc.gradtarg)/norm(vals1dc.gradtarg) < eps

end

@testset "testing Stokes utilities" begin

    n = 1000
    nt = 1100
    sources = randn(3,n)
    targets = randn(3,nt)
    stoklet = randn(3,n)
    strslet = randn(3,n)
    strsvec = randn(3,n)

    vals1c = st3ddir(sources,targets,stoklet=stoklet,
                    ppregt=3,thresh=1e-16)

    vals1cs = st3ddir(sources,sources,stoklet=stoklet,
                    ppregt=3,thresh=1e-16)

    vals1d = st3ddir(sources,targets,strslet=strslet,strsvec=strsvec,
                    ppregt=3,thresh=1e-16)

    vals1ds = st3ddir(sources,sources,strslet=strslet,strsvec=strsvec,
                    ppregt=3,thresh=1e-16)

    vals1dc = st3ddir(sources,targets,stoklet=stoklet,
                     strslet=strslet,strsvec=strsvec,ppregt=3,thresh=1e-16)

    vals1dcs = st3ddir(sources,sources,stoklet=stoklet,
                     strslet=strslet,strsvec=strsvec,ppregt=3,thresh=1e-16)

    eps = 1e-12
    vals2c = stfmm3d(eps,sources,targets=targets,stoklet=stoklet,
                    ppreg=3,ppregt=3)
    vals2d = stfmm3d(eps,sources,targets=targets,strslet=strslet,strsvec=strsvec,
                    ppreg=3,ppregt=3)
    vals2dc = stfmm3d(eps,sources,targets=targets,stoklet=stoklet,
                     strslet=strslet,strsvec=strsvec,ppreg=3,ppregt=3)

    @test norm(vals2c.pot-vals1cs.pottarg)/norm(vals1cs.pottarg) < eps
    @test norm(vals2c.pre-vals1cs.pretarg)/norm(vals1cs.pretarg) < eps
    @test norm(vals2c.grad-vals1cs.gradtarg)/norm(vals1cs.gradtarg) < eps
    @test norm(vals2c.pottarg-vals1c.pottarg)/norm(vals1c.pottarg) < eps
    @test norm(vals2c.gradtarg-vals1c.gradtarg)/norm(vals1c.gradtarg) < eps
    @test norm(vals2c.pretarg-vals1c.pretarg)/norm(vals1c.pretarg) < eps

    @test norm(vals2d.pot-vals1ds.pottarg)/norm(vals1ds.pottarg) < eps
    @test norm(vals2d.pre-vals1ds.pretarg)/norm(vals1ds.pretarg) < eps
    @test norm(vals2d.grad-vals1ds.gradtarg)/norm(vals1ds.gradtarg) < eps
    @test norm(vals2d.pottarg-vals1d.pottarg)/norm(vals1d.pottarg) < eps
    @test norm(vals2d.pretarg-vals1d.pretarg)/norm(vals1d.pretarg) < eps
    @test norm(vals2d.gradtarg-vals1d.gradtarg)/norm(vals1d.gradtarg) < eps

    @test norm(vals2dc.pot-vals1dcs.pottarg)/norm(vals1dcs.pottarg) < eps
    @test norm(vals2dc.pre-vals1dcs.pretarg)/norm(vals1dcs.pretarg) < eps
    @test norm(vals2dc.grad-vals1dcs.gradtarg)/norm(vals1dcs.gradtarg) < eps
    @test norm(vals2dc.pottarg-vals1dc.pottarg)/norm(vals1dc.pottarg) < eps
    @test norm(vals2dc.pretarg-vals1dc.pretarg)/norm(vals1dc.pretarg) < eps
    @test norm(vals2dc.gradtarg-vals1dc.gradtarg)/norm(vals1dc.gradtarg) < eps

end

@testset "testing electromagnetics utilities" begin

    n = 10
    nt = 11
    sources = randn(3,n)
    targets = randn(3,nt)
    h_current = randn(3,n) + im*randn(3,n)
    e_current = randn(3,n) + im*randn(3,n)
    e_charge = randn(n) + im*randn(n)
    zk = 2.1

    vals1A = em3ddir(zk,sources,targets,h_current=h_current,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1As = em3ddir(zk,sources,sources,h_current=h_current,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1B = em3ddir(zk,sources,targets,e_current=e_current,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1Bs = em3ddir(zk,sources,sources,e_current=e_current,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1lambda = em3ddir(zk,sources,targets,e_charge=e_charge,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1lambdas = em3ddir(zk,sources,sources,e_charge=e_charge,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1all = em3ddir(zk,sources,targets,h_current=h_current,e_current=e_current,e_charge=e_charge,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    vals1alls = em3ddir(zk,sources,sources,h_current=h_current,e_current=e_current,e_charge=e_charge,
                     ifEtarg=true,ifdivEtarg=true,
                     ifcurlEtarg=true,thresh=1e-16)

    eps = 1e-12
    vals2A = emfmm3d(eps,zk,sources,targets=targets,h_current=h_current,
                    ifE=true,ifdivE=true,ifcurlE=true,
                    ifEtarg=true,ifdivEtarg=true,ifcurlEtarg=true)

    vals2B = emfmm3d(eps,zk,sources,targets=targets,e_current=e_current,
                    ifE=true,ifdivE=true,ifcurlE=true,
                    ifEtarg=true,ifdivEtarg=true,ifcurlEtarg=true)

    vals2lambda = emfmm3d(eps,zk,sources,targets=targets,e_charge=e_charge,
                    ifE=true,ifdivE=true,ifcurlE=true,
                    ifEtarg=true,ifdivEtarg=true,ifcurlEtarg=true)

    vals2all = emfmm3d(eps,zk,sources,targets=targets,h_current=h_current,e_current=e_current,e_charge=e_charge,
                    ifE=true,ifdivE=true,ifcurlE=true,
                    ifEtarg=true,ifdivEtarg=true,ifcurlEtarg=true)

    function absrelerr(vex,v)
        return norm(vex-v)/max(norm(vex),1)
    end

    @test absrelerr(vals1As.Etarg,vals2A.E) < eps
    @test absrelerr(vals1A.Etarg,vals2A.Etarg) < eps
    @test absrelerr(vals1As.divEtarg,vals2A.divE) < eps
    @test absrelerr(vals1A.divEtarg,vals2A.divEtarg) < eps
    @test absrelerr(vals1As.curlEtarg,vals2A.curlE) < eps
    @test absrelerr(vals1A.curlEtarg,vals2A.curlEtarg) < eps

    @test absrelerr(vals1Bs.Etarg,vals2B.E) < eps
    @test absrelerr(vals1B.Etarg,vals2B.Etarg) < eps
    @test absrelerr(vals1Bs.divEtarg,vals2B.divE) < eps
    @test absrelerr(vals1B.divEtarg,vals2B.divEtarg) < eps
    @test absrelerr(vals1Bs.curlEtarg,vals2B.curlE) < eps
    @test absrelerr(vals1B.curlEtarg,vals2B.curlEtarg) < eps

    @test absrelerr(vals1lambdas.Etarg,vals2lambda.E) < eps
    @test absrelerr(vals1lambda.Etarg,vals2lambda.Etarg) < eps
    @test absrelerr(vals1lambdas.divEtarg,vals2lambda.divE) < eps
    @test absrelerr(vals1lambda.divEtarg,vals2lambda.divEtarg) < eps
    @test absrelerr(vals1lambdas.curlEtarg,vals2lambda.curlE) < eps
    @test absrelerr(vals1lambda.curlEtarg,vals2lambda.curlEtarg) < eps

    @test absrelerr(vals1alls.Etarg,vals2all.E) < eps
    @test absrelerr(vals1all.Etarg,vals2all.Etarg) < eps
    @test absrelerr(vals1alls.divEtarg,vals2all.divE) < eps
    @test absrelerr(vals1all.divEtarg,vals2all.divEtarg) < eps
    @test absrelerr(vals1alls.curlEtarg,vals2all.curlE) < eps
    @test absrelerr(vals1all.curlEtarg,vals2all.curlEtarg) < eps

end

@testset "testing lower level routines" begin

    # besseljs3d

    nterms = 10
    z = 1.1 + im*1.2
    scale = 1.3
    ifder = 1

    fj10 = (-3.5183264829466616750769540*(1e-9) +
        im*8.88237983492960538163695769*(1e-9))

    fjs, fjder = besseljs3d(nterms,z,scale=scale,ifder=ifder)

    @test (abs(fj10-fjs[11]*(scale^10))/abs(fj10) < 1e-10)

end


text:
function innerjoin; end

const ⨝ = innerjoin

# TODO make signature compatible with theta-joins, etc.
function innerjoin(t1::Table{n1}, t2::Table{n2}) where {n1, n2}
    all_names = _union(n1, n2)
    same_names = _interesect(Val(n1), Val(n2))
    new_names = _tempnames(Val(same_names))

    @assert length(same_names) === 1 # TODO fix `Equals`...

    t2_renamed = Rename(same_names, new_names)(t2)
    t3 = filter(Equals{(same_names..., new_names...)}, t1 × t2_renamed)
    return Project{(all_names,)}()(t3)
end

# TODO how do "outer" joins work in this framework?


text:
# Opus error codes
const OPUS_OK               =  0
const OPUS_BAD_ARG          = -1
const OPUS_BUFFER_TOO_SMALL = -2
const OPUS_INTERNAL_ERROR   = -3
const OPUS_INVALID_PACKET   = -4
const OPUS_UNIMPLEMENTED    = -5
const OPUS_INVALID_STATE    = -6
const OPUS_ALLOC_FAIL       = -7

const OPUS_ERROR_MESSAGE_STRS = Dict(
    0  => "No Error",
    -1 => "Bad Argument",
    -2 => "Buffer Too Small",
    -3 => "Internal Error",
    -4 => "Invalid Packet",
    -5 => "Unimplemented",
    -6 => "Invalid State",
    -7 => "Allocation Failure"
)

const OPUS_APPLICATION_VOIP                = 2048
const OPUS_APPLICATION_AUDIO               = 2049
const OPUS_APPLICATION_RESTRICTED_LOWDELAY = 2051

immutable OpusHead
    # Magic bytes "OpusHead" (0x646165487375704f)
    opus_head::UInt64
    # Should always be equal to one
    version::UInt8
    # Number of channels, must be greater than zero
    channels::UInt8
    # This is the number of samples (at 48 kHz) to discard from the decoder
    # output when starting playback, and also the number to subtract from a
    # page's granule position to calculate its PCM sample position.
    # NOTE: This is currently completely ignored in Opus.jl
    preskip::UInt16
    # Samplerate of input stream (Let's face it, it's always 48 KHz)
    samplerate::UInt32
    # Output gain that should be applied
    output_gain::UInt16
    # Channel mapping family, we always assume this is zero
    channel_map_family::UInt8
end

OpusHead() = OpusHead(0x646165487375704f, 1, 1, 312, 48000, 0, 0)
OpusHead(samplerate, channels) = OpusHead(0x646165487375704f, 1, channels, 312, samplerate, 0, 0)

function OpusHead(io::IO)
    magic = read(io, UInt64)
    if magic != 0x646165487375704f
        error("Input packet is not an OpusHead!, magic is $(magic)")
    end
    version = read(io, UInt8)
    channels = read(io, UInt8)
    preskip = read(io, UInt16)
    samplerate = read(io, UInt32)
    output_gain = read(io, UInt16)
    channel_map_family = read(io, UInt8)
    return OpusHead( magic, version, channels, preskip, samplerate, output_gain, channel_map_family)
end
OpusHead(data::Vector{UInt8}) = OpusHead(IOBuffer(data))

function write(io::IO, x::OpusHead)
    write(io, x.opus_head)
    write(io, x.version)
    write(io, x.channels)
    write(io, x.preskip)
    write(io, x.samplerate)
    write(io, x.output_gain)
    write(io, x.channel_map_family)
end

function convert(::Type{Vector{UInt8}}, x::OpusHead)
    io = IOBuffer()
    write(io, x)
    seekstart(io)
    return readbytes(io)
end

function show(io::IO, x::OpusHead)
    write(io, "OpusHead packet ($(x.samplerate)Hz, $(x.channels) channels)")
end

type OpusTags
    # Magic signature "OpusTags" (0x736761547375704f)
    opus_tags::UInt64

    vendor_string::AbstractString
    tags::Vector{AbstractString}
end
OpusTags() = OpusTags(0x736761547375704f, "Opus.jl", AbstractString["encoder=Opus.jl"])

function read_opus_tag(io::IO)
    # First, read in a length
    len = read(io, UInt32)

    # Next, read the string
    return bytestring(readbytes(io, len))
end

function write_opus_tag(io::IO, tag::AbstractString)
    # First, write out the length
    write(io, UInt32(length(tag)))

    # Next, write out the tag itself
    write(io, tag)
end

function OpusTags(io::IO)
    magic = read(io, UInt64)
    if magic != 0x736761547375704f
        error("Input packet is not an OpusTags!, magic is $(magic)")
    end
    # First, read the vendor string
    vendor_string = read_opus_tag(io)

    # Next, read how many tags we've got
    num_tags = read(io, UInt32)

    # Read all the tags in, one after another
    tags = [read_opus_tag(io) for idx in 1:num_tags]

    return OpusTags(magic, vendor_string, tags)
end
OpusTags(data::Vector{UInt8}) = OpusTags(IOBuffer(data))

function write(io::IO, x::OpusTags)
    write(io, x.opus_tags)
    write_opus_tag(io, x.vendor_string)

    write(io, UInt32(length(x.tags)))
    for tagidx = 1:length(x.tags)
        write_opus_tag(io, x.tags[tagidx])
    end
end

function convert(::Type{Vector{UInt8}}, x::OpusTags)
    io = IOBuffer()
    write(io, x)
    seekstart(io)
    return readbytes(io)
end

function show(io::IO, x::OpusTags)
    write(io, "OpusTags packet\n")
    write(io, "  Vendor: $(x.vendor_string)\n")
    write(io, "  Tags:")
    for tag in x.tags
        write(io, "\n")
        write(io, "    $tag")
    end
end

function is_header_packet(packet::Vector{UInt8})
    # Check if it's an OpusHead or OpusTags packet
    if length(packet) > 8
        magic = bytestring(packet[1:8])
        if magic == "OpusHead" || magic == "OpusTags"
            return true
        end
    end

    return false
end


text:
# This file contains the definition for the Lines struct.

"""
calc_lineparams(lines_df::DataFrame)
	Calculates line parameters with regard to voltage level, number of circuits,
	and Thermal Reliability Margin (TRM).
	Returns reactance (p.u.), resistance (p.u.), and pmax (MW)
"""
function calc_lineparams(lines_df::DataFrame, Sbase::Number, TRM::Number)

	zbase(voltage::Number) = (voltage*1E3)^2 / (Sbase*1E6)

	# L = Symbol.(lines_df[!, 1])
	L = lines_df[!, 1]
	zbase_vlevel = map(zbase, lines_df[!, :voltage])	# Calculate reactance in p.u. given voltage level

	reactance = Dict(zip(L, lines_df[!, :reactance] ./ lines_df[!, :circuits] ./ zbase_vlevel))
	resistance = Dict(zip(L,  lines_df[!, :resistance] ./ lines_df[!, :circuits] ./ zbase_vlevel))

	pmax = Dict(zip(L, lines_df[!, :pmax] .* lines_df[!, :circuits] .* (1-TRM) ))
	return reactance, resistance, pmax
end

mutable struct Lines
    id
    from
    to
	voltage
	reactance
	resistance
	pmax

    function Lines(lines_df::DataFrame, Sbase::Number, TRM::Number)

		# L = Symbol.(lines_df[!, 1])
		L = lines_df[!, 1]
		reactance, resistance, pmax = calc_lineparams(lines_df, Sbase, TRM)                    # U = unit

		function lines_df_to_dict_with_id(df::DataFrame)
		    Dict(col[1] => Dict(zip(df[!, 1], col[2])) for col in eachcol(df, true))
		end

		lines_dict = lines_df_to_dict_with_id(lines_df)

        return new(L,
                   lines_dict[:from],
                   lines_dict[:to],
				   lines_dict[:voltage],
				   reactance,
				   resistance,
				   pmax
				   )
    end
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
function solve_op_heur(op::OrienteeringProblem; rng=nothing)
    openset = [1:length(op)]
    total_dist = 0.0
    @assert op.start != op.stop # code may not work for this case
    choice = -1
    path = [op.start]
    current = op.start
    while current != op.stop
        future = SimpleOP(op.r[openset],
                          [],
                          op.distance_limit-total_dist,
                          findfirst(openset,path[end]),
                          findfirst(openset,op.stop),
                          op.distances[openset,openset])
        choice = next_node_heur(future, rng=rng)
        choice = openset[choice]
        total_dist += op.distances[path[end],choice]
        @show push!(path,choice)
        splice!(openset,current)
        current = choice
    end
    return path
end

# S-algorithm from Tsiligirides 1984
function next_node_heur(op::OrienteeringProblem; rng=nothing)
    if rng==nothing
        rng = MersenneTwister(rand(Uint))
    end

    @assert op.start != op.stop

    is_feasible(i) = op.distances[op.start,i]+op.distances[i,op.stop]<=op.distance_limit
    feasible = nothing
    # try
    #     feasible = filter(is_feasible, [1:op.start-1, op.start+1:length(op)])
    # catch e
    #     @bp
    # end
    feasible = filter(is_feasible, [1:op.start-1, op.start+1:length(op)])
    @assert length(feasible) > 0

    # Tsiligiride's notation

    T_max = op.distance_limit
    T = 0.0
    S = op.r
    C = op.distances
    NPTS = length(op)
    LAST = op.start

    # r = 2.0
    r = 1.0
    a = 0.0

    # DEL(j) = sum([S[i]/C[j,i] for i in feasible])
    # E(j) = a*(T_max - T - C(LAST,j) - C(j,op.stop))*DEL(j)
    E(j) = 0.0 # because a = 0.0 was used

    A = zeros(length(op))
    for j in 1:length(A)
        if is_feasible(j) && j != op.start
            A[j] = ((S[j] + E(j))/C[op.start,j])^r
        end
    end

    sum_A = sum(A)
    P = [A[j]/sum_A for j in 1:length(op)]

    rn = rand(rng)
    accum = P[1]
    decision = 1
    while accum < rn
        decision += 1
        accum += P[decision]
    end
    return decision
end