text:
# Unit tests for Distances

struct FooDist <: PreMetric end # Julia 1.0 Compat: struct definition must be put in global scope

@testset "result_type" begin
    foodist(a, b) = a + b
    (::FooDist)(a, b) = foodist(a, b)
    for (Ta, Tb) in [
        (Int, Int),
        (Int, Float64),
        (Float32, Float32),
        (Float32, Float64),
    ]
        A, B = rand(Ta, 2, 3), rand(Tb, 2, 3)
        @test result_type(FooDist(), A, B) == result_type(FooDist(), Ta, Tb)
        @test result_type(foodist, A, B) == result_type(foodist, Ta, Tb) == typeof(foodist(oneunit(Ta), oneunit(Tb)))

        a, b = rand(Ta), rand(Tb)
        @test result_type(FooDist(), a, b) == result_type(FooDist(), Ta, Tb)
        @test result_type(foodist, a, b) == result_type(foodist, Ta, Tb) == typeof(foodist(oneunit(Ta), oneunit(Tb)))
    end
end

function test_metricity(dist, x, y, z)
    @testset "Test metricity of $(typeof(dist))" begin
        @test dist(x, y) == evaluate(dist, x, y)

        dxy = dist(x, y)
        dxz = dist(x, z)
        dyz = dist(y, z)
        if isa(dist, PreMetric)
            # Unfortunately small non-zero numbers (~10^-16) are appearing
            # in our tests due to accumulating floating point rounding errors.
            # We either need to allow small errors in our tests or change the
            # way we do accumulations...
            @test dist(x, x) + one(eltype(x)) ≈ one(eltype(x))
            @test dist(y, y) + one(eltype(y)) ≈ one(eltype(y))
            @test dist(z, z) + one(eltype(z)) ≈ one(eltype(z))
            @test dxy ≥ zero(eltype(x))
            @test dxz ≥ zero(eltype(x))
            @test dyz ≥ zero(eltype(x))
        end
        if isa(dist, SemiMetric)
            @test dxy ≈ dist(y, x)
            @test dxz ≈ dist(z, x)
            @test dyz ≈ dist(y, z)
        else # Not symmetric, so more PreMetric tests
            @test dist(y, x) ≥ zero(eltype(x))
            @test dist(z, x) ≥ zero(eltype(x))
            @test dist(z, y) ≥ zero(eltype(x))
        end
        if isa(dist, Metric)
            # Again we have small rounding errors in accumulations
            @test dxz ≤ dxy + dyz || dxz ≈ dxy + dyz
            dyx = dist(y, x)
            @test dyz ≤ dyx + dxz || dyz ≈ dyx + dxz
            dzy = dist(z, y)
            @test dxy ≤ dxz + dzy || dxy ≈ dxz + dzy
        end
    end
end

@testset "PreMetric, SemiMetric, Metric on $T" for T in (Float64, F64)
    Random.seed!(123)
    n = 100
    x = rand(T, n)
    y = rand(T, n)
    z = rand(T, n)

    test_metricity(SqEuclidean(), x, y, z)
    test_metricity(Euclidean(), x, y, z)
    test_metricity(Cityblock(), x, y, z)
    test_metricity(TotalVariation(), x, y, z)
    test_metricity(Chebyshev(), x, y, z)
    test_metricity(Minkowski(2.5), x, y, z)

    test_metricity(CosineDist(), x, y, z)
    test_metricity(CorrDist(), x, y, z)

    test_metricity(ChiSqDist(), x, y, z)

    test_metricity(Jaccard(), x, y, z)
    test_metricity(SpanNormDist(), x, y, z)

    test_metricity(BhattacharyyaDist(), x, y, z)
    test_metricity(HellingerDist(), x, y, z)
    test_metricity(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), x, y, z);

    x₁ = rand(T, 2)
    x₂ = rand(T, 2)
    x₃ = rand(T, 2)

    test_metricity(Haversine(6371.0), x₁, x₂, x₃)
    test_metricity(SphericalAngle(), x₁, x₂, x₃)

    k = rand(1:3, n)
    l = rand(1:3, n)
    m = rand(1:3, n)

    test_metricity(Hamming(), k, l, m)

    a = rand(Bool, n)
    b = rand(Bool, n)
    c = rand(Bool, n)

    test_metricity(RogersTanimoto(), a, b, c)
    test_metricity(BrayCurtis(), a, b, c)
    test_metricity(Jaccard(), a, b, c)

    for w in (rand(T, n), (rand(T, n)...,))
        test_metricity(PeriodicEuclidean(w), x, y, z)
        test_metricity(WeightedSqEuclidean(w), x, y, z)
        test_metricity(WeightedEuclidean(w), x, y, z)
        test_metricity(WeightedCityblock(w), x, y, z)
        test_metricity(WeightedMinkowski(w, 2.5), x, y, z)
        test_metricity(WeightedHamming(w), a, b, c)
    end

    Q = rand(T, n, n)
    Q = Q * Q'  # make sure Q is positive-definite

    test_metricity(SqMahalanobis(Q), x, y, z)
    test_metricity(Mahalanobis(Q), x, y, z)

    p = rand(T, n)
    q = rand(T, n)
    r = rand(T, n)
    p[p .< median(p)] .= 0
    p /= sum(p)
    q /= sum(q)
    r /= sum(r)

    test_metricity(KLDivergence(), p, q, r)
    test_metricity(RenyiDivergence(0.0), p, q, r)
    test_metricity(RenyiDivergence(1.0), p, q, r)
    test_metricity(RenyiDivergence(Inf), p, q, r)
    test_metricity(RenyiDivergence(0.5), p, q, r)
    test_metricity(RenyiDivergence(2), p, q, r)
    test_metricity(RenyiDivergence(10), p, q, r)
    test_metricity(JSDivergence(), p, q, r)
end

@testset "individual metrics" begin
    a = 1
    b = 2
    @test sqeuclidean(a, b) === 1
    @test euclidean(a, b) === 1.0
    @test jaccard(a, b) === 0.5
    @test cityblock(a, b) === 1
    @test totalvariation(a, b) === 0.5
    @test chebyshev(a, b) == 1.0
    @test braycurtis(a, b) === 1/3
    @test minkowski(a, b, 2) == 1.0
    @test hamming(a, b) === 1
    @test hamming("martha", "marhta") === 2
    @test hamming("es an ", " vs an") === 6
    @test hamming("", "") === 0
    @test peuclidean(a, b, 0.5) === 0.0
    @test peuclidean(a, b, 2) === 1.0
    @test cosine_dist(a, b) === 0.0
    @test bhattacharyya(a, b) === bhattacharyya([a], [b]) === -0.0
    @test bhattacharyya(a, b) === bhattacharyya((a,), (b,))
    @test isnan(corr_dist(a, b))
    @test spannorm_dist(a, b) === 0

    bt = [true, false, true]
    bf = [false, true, true]
    @test rogerstanimoto(bt, bf) == 4.0 / 5.0
    @test braycurtis(bt, bf) == 0.5

    for w in (2, (2,))
        @test wsqeuclidean(a, b, w) === 2
        @test weuclidean(a, b, w) === sqrt(2)
        @test wcityblock(a, b, w) === 2
        @test wminkowski(a, b, w, 2) === sqrt(2)
        @test whamming(a, b, w) === 2
    end

    for T in (Float64, F64)
        for (_x, _y) in (([4.0, 5.0, 6.0, 7.0], [3.0, 9.0, 8.0, 1.0]),
                         ([4.0, 5.0, 6.0, 7.0], [3. 8.; 9. 1.0]))
            x, y = T.(_x), T.(_y)
            for (x, y) in ((x, y),
                           (convert(Array{Union{Missing, T}}, x), convert(Array{Union{Missing, T}}, y)),
                           ((Iterators.take(x, 4), Iterators.take(y, 4))), # iterator
                           (((x[i] for i in 1:length(x)), (y[i] for i in 1:length(y)))), # generator
                          )
                xc, yc = collect(x), collect(y)
                @test sqeuclidean(x, y) == 57.0
                @test euclidean(x, y) == sqrt(57.0)
                @test jaccard(x, y) == 13.0 / 28
                @test cityblock(x, y) == 13.0
                @test totalvariation(x, y) == 6.5
                @test chebyshev(x, y) == 6.0
                @test braycurtis(x, y) == 1.0 - (30.0 / 43.0)
                @test minkowski(x, y, 2) == sqrt(57.0)
                @test peuclidean(x, y, fill(10.0, 4)) == sqrt(37)
                @test peuclidean(xc - vec(yc), zero(yc), fill(10.0, 4)) == peuclidean(x, y, fill(10.0, 4))
                @test peuclidean(x, y, [10.0, 10.0, 10.0, Inf]) == sqrt(57)
                @test_throws DimensionMismatch cosine_dist(1.0:2, 1.0:3)
                @test cosine_dist(x, y) ≈ (1.0 - 112. / sqrt(19530.0))
                x_int, y_int = Int64.(x), Int64.(y)
                @test cosine_dist(x_int, y_int) == (1.0 - 112.0 / sqrt(19530.0))
                @test corr_dist(x, y) ≈ cosine_dist(x .- mean(x), vec(yc) .- mean(y))
                @test corr_dist(OffsetVector(xc, -1:length(xc)-2), yc) == corr_dist(x, y)
                @test chisq_dist(x, y) == sum((xc - vec(yc)).^2 ./ (xc + vec(yc)))
                @test spannorm_dist(x, y) == maximum(xc - vec(yc)) - minimum(xc - vec(yc))

                @test gkl_divergence(x, y) ≈ sum(i -> xc[i] * log(xc[i] / yc[i]) - xc[i] + yc[i], 1:length(x))

                @test meanad(x, y) ≈ mean(Float64[abs(xc[i] - yc[i]) for i in 1:length(x)])
                @test msd(x, y) ≈ mean(Float64[abs2(xc[i] - yc[i]) for i in 1:length(x)])
                @test rmsd(x, y) ≈ sqrt(msd(x, y))
                @test nrmsd(x, y) ≈ sqrt(msd(x, y)) / (maximum(x) - minimum(x))

                w = ones(4)
                @test sqeuclidean(x, y) ≈ wsqeuclidean(x, y, w)

                w = rand(Float64, length(x))
                @test wsqeuclidean(x, y, w) ≈ dot((xc - vec(yc)).^2, w)
                @test weuclidean(x, y, w) == sqrt(wsqeuclidean(x, y, w))
                @test wcityblock(x, y, w) ≈ dot(abs.(xc - vec(yc)), w)
                @test wminkowski(x, y, w, 2) ≈ weuclidean(x, y, w)
            end
        end

        # Test ChiSq doesn't give NaN at zero
        @test chisq_dist([0.0], [0.0]) == 0.0

        # Test weighted Hamming distances with even weights
        a = T.([1.0, 2.0, 1.0, 3.0, 2.0, 1.0])
        b = T.([1.0, 3.0, 0.0, 2.0, 2.0, 0.0])
        w = rand(T, size(a))

        @test whamming(a, a, w) === T(0.0)
        @test whamming(a, b, w) === sum((a .!= b) .* w)

        # Minimal test of Jaccard - test return type stability.
        @inferred Jaccard()(rand(T, 3), rand(T, 3))
        @inferred Jaccard()([1, 2, 3], [1, 2, 3])
        @inferred Jaccard()([true, false, true], [false, true, true])

        # Test Bray-Curtis. Should be 1 if no elements are shared, 0 if all are the same
        @test braycurtis([1,0,3],[0,1,0]) == 1.0
        @test braycurtis(rand(10), zeros(10)) == 1.0
        @test braycurtis([1,0],[1,0]) == 0.0

        # Test KL, Renyi and JS divergences
        r = rand(T, 12)
        p = copy(r)
        p[p .< median(p)] .= 0.0
        scale = sum(p) / sum(r)
        r /= sum(r)
        p /= sum(p)
        q = rand(T, 12)
        q /= sum(q)

        klv = 0.0
        for i = 1:length(p)
            if p[i] > 0
                klv += p[i] * log(p[i] / q[i])
            end
        end

        pm = (p + q) / 2
        for (r, p, pm) in ((r, p, pm),
                           (Iterators.take(r, length(r)), Iterators.take(p, length(p)), Iterators.take(pm, length(pm))),
                           ((r[i] for i in 1:length(r)), (p[i] for i in 1:length(p)), (pm[i] for i in 1:length(pm))),
                          )
            @test kl_divergence(p, q) ≈ klv
            @test typeof(kl_divergence(p, q)) == T

            @test renyi_divergence(p, r, 0) ≈ -log(scale)
            @test renyi_divergence(p, r, 1) ≈ -log(scale)
            @test renyi_divergence(p, r, 10) ≈ -log(scale)
            @test renyi_divergence(p, r, rand()) ≈ -log(scale)
            @test renyi_divergence(p, r, Inf) ≈ -log(scale)
            @test isinf(renyi_divergence([0.0, 0.5, 0.5], [0.0, 1.0, 0.0], Inf))
            @test renyi_divergence([0.0, 1.0, 0.0], [0.0, 0.5, 0.5], Inf) ≈ log(2.0)
            @test renyi_divergence(p, q, 1) ≈ kl_divergence(p, q)

            jsv = kl_divergence(p, pm) / 2 + kl_divergence(q, pm) / 2
            @test js_divergence(p, q) ≈ jsv
        end
    end
end # testset

@testset "NaN behavior" begin
    a = [NaN, 0]; b = [0, NaN]
    @test isnan(chebyshev(a, b)) == isnan(maximum(a - b))
    a = [NaN, 0]; b = [0, 1]
    @test isnan(chebyshev(a, b)) == isnan(maximum(a - b))
    @test isnan(renyi_divergence([0.5, 0.0, 0.5], [0.5, 0.5, NaN], 2))
end #testset

@testset "empty vector" begin
    for T in (Float64, F64), (a, b) in ((T[], T[]), (Iterators.take(T[], 0), Iterators.take(T[], 0)))
        @test sqeuclidean(a, b) == 0.0
        @test isa(sqeuclidean(a, b), T)
        @test euclidean(a, b) == 0.0
        @test isa(euclidean(a, b), T)
        @test cityblock(a, b) == 0.0
        @test isa(cityblock(a, b), T)
        @test totalvariation(a, b) == 0.0
        @test isa(totalvariation(a, b), T)
        @test chebyshev(a, b) == 0.0
        @test isa(chebyshev(a, b), T)
        @test braycurtis(a, b) == 0.0
        @test isa(braycurtis(a, b), T)
        @test minkowski(a, b, 2) == 0.0
        @test isa(minkowski(a, b, 2), T)
        @test hamming(a, b) == 0.0
        @test isa(hamming(a, b), Int)
        @test renyi_divergence(a, b, 1.0) == 0.0
        @test isa(renyi_divergence(a, b, 2.0), T)
        @test braycurtis(a, b) == 0.0
        @test isa(braycurtis(a, b), T)

        w = T[]
        @test isa(whamming(a, b, w), T)
        @test peuclidean(a, b, w) == 0.0
        @test isa(peuclidean(a, b, w), T)
    end
end # testset

@testset "DimensionMismatch throwing" begin
    a = [1, 0]; b = [2]
    @test_throws DimensionMismatch sqeuclidean(a, b)
    a = (1, 0); b = (2,)
    @test_throws DimensionMismatch sqeuclidean(a, b)
    a = (1, 0); b = (2.0,); w = (3.0,)
    @test_throws DimensionMismatch wsqeuclidean(a, b, w)
    @test_throws DimensionMismatch peuclidean(a, b, w)
    a = [1, 0]; b = [2.0, 4.0] ; w = [3.0]
    @test_throws DimensionMismatch wsqeuclidean(a, b, w)
    @test_throws DimensionMismatch peuclidean(a, b, w)
    a = (1, 0); b = (2.0, 4.0) ; w = (3.0,)
    @test_throws DimensionMismatch wsqeuclidean(a, b, w)
    @test_throws DimensionMismatch peuclidean(a, b, w)
    p = [0.5, 0.5]; q = [0.3, 0.3, 0.4]
    @test_throws DimensionMismatch bhattacharyya(p, q)
    @test_throws DimensionMismatch hellinger(q, p)
    Q = rand(length(p), length(p))
    Q = Q * Q'  # make sure Q is positive-definite
    @test_throws DimensionMismatch mahalanobis(p, q, Q)
    @test_throws DimensionMismatch mahalanobis(q, q, Q)
    mat23 = [0.3 0.2 0.0; 0.1 0.0 0.4]
    mat22 = [0.3 0.2; 0.1 0.4]
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, mat23)
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, q)
    @test_throws DimensionMismatch colwise!(mat23, Euclidean(), mat23, mat22)
    @test_throws DimensionMismatch colwise!(mat23, Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), mat23, mat22)
    @test_throws DimensionMismatch Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x)([1, 2, 3], [1, 2])
    @test_throws DimensionMismatch Bregman(x -> sqeuclidean(x, zero(x)), x -> [1, 2])([1, 2, 3], [1, 2, 3])
end # testset

@testset "Different input types" begin
    for (x, y) in (([4, 5, 6, 7], [3.0, 9.0, 8.0, 1.0]),
                   ([4, 5, 6, 7], [3//1 8; 9 1]))
        @test (@inferred sqeuclidean(x, y)) == 57
        @test (@inferred euclidean(x, y)) == sqrt(57)
        @test (@inferred jaccard(x, y)) == convert(Base.promote_eltype(x, y), 13 // 28)
        @test (@inferred cityblock(x, y)) == 13
        @test (@inferred totalvariation(x, y)) == 6.5
        @test (@inferred chebyshev(x, y)) == 6
        @test (@inferred braycurtis(x, y)) == convert(Base.promote_eltype(x, y), 13 // 43)
        @test (@inferred minkowski(x, y, 2)) == sqrt(57)
        @test (@inferred peuclidean(x, y, fill(10, 4))) == sqrt(37)
        @test (@inferred peuclidean(x - vec(y), zero(y), fill(10, 4))) == peuclidean(x, y, fill(10, 4))
        @test (@inferred peuclidean(x, y, [10.0, 10.0, 10.0, Inf])) == sqrt(57)
        @test_throws DimensionMismatch cosine_dist(1.0:2, 1.0:3)
        @test (@inferred cosine_dist(x, y)) ≈ (1 - 112 / sqrt(19530))
        @test (@inferred corr_dist(x, y)) ≈ cosine_dist(x .- mean(x), vec(y) .- mean(y))
        @test (@inferred chisq_dist(x, y)) == sum((x - vec(y)).^2 ./ (x + vec(y)))
        @test (@inferred spannorm_dist(x, y)) == maximum(x - vec(y)) - minimum(x - vec(y))

        @test (@inferred gkl_divergence(x, y)) ≈ sum(i -> x[i] * log(x[i] / y[i]) - x[i] + y[i], 1:length(x))

        @test (@inferred meanad(x, y)) ≈ mean(Float64[abs(x[i] - y[i]) for i in 1:length(x)])
        @test (@inferred msd(x, y)) ≈ mean(Float64[abs2(x[i] - y[i]) for i in 1:length(x)])
        @test (@inferred rmsd(x, y)) ≈ sqrt(msd(x, y))
        @test (@inferred nrmsd(x, y)) ≈ sqrt(msd(x, y)) / (maximum(x) - minimum(x))

        w = ones(Int, 4)
        @test sqeuclidean(x, y) ≈ wsqeuclidean(x, y, w)

        w = rand(1:length(x), size(x))
        @test (@inferred wsqeuclidean(x, y, w)) ≈ dot((x - vec(y)).^2, w)
        @test (@inferred weuclidean(x, y, w)) == sqrt(wsqeuclidean(x, y, w))
        @test (@inferred wcityblock(x, y, w)) ≈ dot(abs.(x - vec(y)), w)
        @test (@inferred wminkowski(x, y, w, 2)) ≈ weuclidean(x, y, w)
    end
end

@testset "mahalanobis" begin
    for T in (Float64, F64, ComplexF64)
        x, y = T.([4.0, 5.0, 6.0, 7.0]), T.([3.0, 9.0, 8.0, 1.0])

        Q = rand(T, length(x), length(x))
        Q = Q * Q'  # make sure Q is positive-definite
        for Q in (Q, T <: Complex ? Hermitian(Q) : Symmetric(Q))
            @test sqmahalanobis(x, y, Q) ≈ dot(x - y, Q * (x - y))
            @test eltype(sqmahalanobis(x, y, Q)) == T
            @test mahalanobis(x, y, Q) == sqrt(sqmahalanobis(x, y, Q))
            @test eltype(mahalanobis(x, y, Q)) == T
        end
    end
end #testset

@testset "haversine" begin
    for T in (Float64, F64)
        @test haversine([-180.,0.], [180.,0.], 1.) ≈ 0 atol=1e-10
        @test haversine([0.,-90.],  [0.,90.],  1.) ≈ π atol=1e-10
        @test haversine((-180.,0.), (180.,0.), 1.) ≈ 0 atol=1e-10
        @test haversine((0.,-90.),  (0.,90.),  1.) ≈ π atol=1e-10
        x, y = (1.,-15.625), (-179.,15.625)
        @test haversine(x, y, 6371.) ≈ 20015 atol=1e-1
        @test haversine(x, y) ≈ 20015086 rtol=1e-7
        @test Haversine()(x, y) ≈ 20015086 rtol=1e-7
        @test haversine(Float32.(x), Float32.(y)) isa Float32
        @test haversine(Iterators.take(x, 2), Iterators.take(y, 2), 6371.) ≈ 20015 atol=1e-1
        @test_throws ArgumentError haversine([0.,-90., 0.25], [0.,90.], 1.)
        x, y = (1.0°,-15.625°), (-179.0°,15.625°)
        @test haversine(x, y, 6371.0km) ≈ 20015km atol=1e-1km
    end
end

@testset "spherical angle" begin
    @test spherical_angle([-π, 0],    [π,0]) ≈ 0 atol=1e-10
    @test spherical_angle([0,-π/2], [0,π/2]) ≈ π atol=1e-10
    @test spherical_angle((-π,0),     (π,0)) ≈ 0 atol=1e-10
    @test spherical_angle((0,-π/2), (0,π/2)) ≈ π atol=1e-10
    x, y = rand(2), rand(2)
    x_deg, y_deg = rad2deg.(x), rad2deg.(y)
    @test spherical_angle(x, y) ≈ haversine(x_deg, y_deg, 1) atol=1e-10
    @test spherical_angle(Iterators.take(x, 2), Iterators.take(y, 2)) ≈ haversine(x_deg, y_deg, 1) atol=1e-10
    @test_throws ArgumentError spherical_angle([0.,-π/2, 0.25], [0.,π/2])
end

@testset "bhattacharyya / hellinger" begin
    for T in (Int, Float64, F64)
        x, y = T.([4, 5, 6, 7]), T.([3, 9, 8, 1])
        a = T.([1, 2, 1, 3, 2, 1])
        b = T.([1, 3, 0, 2, 2, 0])
        p = T == Int ? rand(0:10, 12) : rand(T, 12)
        p[p .< median(p)] .= 0
        q = T == Int ? rand(0:10, 12) : rand(T, 12)

        # Bhattacharyya and Hellinger distances are defined for discrete
        # probability distributions so to calculate the expected values
        # we need to normalize vectors.
        px = x ./ sum(x)
        py = y ./ sum(y)
        expected_bc_x_y = sum(sqrt.(px .* py))
        for (x, y) in ((x, y), (Iterators.take(x, 12), Iterators.take(y, 12)))
            @test Distances.bhattacharyya_coeff(x, y) ≈ expected_bc_x_y
            @test bhattacharyya(x, y) ≈ (-log(expected_bc_x_y))
            @test hellinger(x, y) ≈ sqrt(1 - expected_bc_x_y)
        end

        pa = a ./ sum(a)
        pb = b ./ sum(b)
        expected_bc_a_b = sum(sqrt.(pa .* pb))
        @test Distances.bhattacharyya_coeff(a, b) ≈ expected_bc_a_b
        @test bhattacharyya(a, b) ≈ (-log(expected_bc_a_b))
        @test hellinger(a, b) ≈ sqrt(1 - expected_bc_a_b)

        pp = p ./ sum(p)
        pq = q ./ sum(q)
        expected_bc_p_q = sum(sqrt.(pp .* pq))
        @test Distances.bhattacharyya_coeff(p, q) ≈ expected_bc_p_q
        @test bhattacharyya(p, q) ≈ (-log(expected_bc_p_q))
        @test hellinger(p, q) ≈ sqrt(1 - expected_bc_p_q)

        # Ensure it is semimetric
        @test bhattacharyya(x, y) ≈ bhattacharyya(y, x)
    end
end #testset

function test_colwise(dist, x, y, T)
    @testset "Colwise test for $(typeof(dist))" begin
        n = size(x, 2)
        r1 = zeros(T, n)
        r2 = zeros(T, n)
        r3 = zeros(T, n)
        r4 = zeros(T, 1, n)
        for j = 1:n
            r1[j] = dist(x[:, j], y[:, j])
            r2[j] = dist(x[:, 1], y[:, j])
            r3[j] = dist(x[:, j], y[:, 1])
        end
        # ≈ and all( .≈ ) seem to behave slightly differently for F64
        @test all(colwise(dist, x, y) .≈ r1)
        @test all(colwise(dist, (x[:,i] for i in axes(x, 2)), (y[:,i] for i in axes(y, 2))) .≈ r1)
        colwise!(r4, dist, x, y)
        @test all(r4[i] ≈ r1[i] for i in 1:n)
        @test all(colwise(dist, x[:, 1], y) .≈ r2)
        @test all(colwise(dist, x, y[:, 1]) .≈ r3)
    end
end

@testset "column-wise metrics on $T" for T in (Float64, F64)
    m = 5
    n = 8
    X = rand(T, m, n)
    Y = rand(T, m, n)
    A = rand(1:3, m, n)
    B = rand(1:3, m, n)

    P = rand(T, m, n)
    Q = rand(T, m, n)
    # Make sure not to remove all of the non-zeros from any column
    for i in 1:n
        P[P[:, i] .< median(P[:, i]) / 2, i] .= 0.0
    end

    test_colwise(SqEuclidean(), X, Y, T)
    test_colwise(Euclidean(), X, Y, T)
    test_colwise(Cityblock(), X, Y, T)
    test_colwise(TotalVariation(), X, Y, T)
    test_colwise(Chebyshev(), X, Y, T)
    test_colwise(Minkowski(2.5), X, Y, T)
    test_colwise(Hamming(), A, B, T)
    test_colwise(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), X, Y, T);

    test_colwise(CosineDist(), X, Y, T)
    test_colwise(CorrDist(), X, Y, T)

    test_colwise(ChiSqDist(), X, Y, T)
    test_colwise(KLDivergence(), P, Q, T)
    test_colwise(RenyiDivergence(0.0), P, Q, T)
    test_colwise(RenyiDivergence(1.0), P, Q, T)
    test_colwise(RenyiDivergence(Inf), P, Q, T)
    test_colwise(RenyiDivergence(0.5), P, Q, T)
    test_colwise(RenyiDivergence(2), P, Q, T)
    test_colwise(RenyiDivergence(10), P, Q, T)
    test_colwise(JSDivergence(), P, Q, T)
    test_colwise(SpanNormDist(), X, Y, T)

    test_colwise(BhattacharyyaDist(), X, Y, T)
    test_colwise(HellingerDist(), X, Y, T)
    test_colwise(BrayCurtis(), X, Y, T)

    w = rand(T, m)

    test_colwise(WeightedSqEuclidean(w), X, Y, T)
    test_colwise(WeightedEuclidean(w), X, Y, T)
    test_colwise(WeightedCityblock(w), X, Y, T)
    test_colwise(WeightedMinkowski(w, 2.5), X, Y, T)
    test_colwise(WeightedHamming(w), A, B, T)
    test_colwise(PeriodicEuclidean(w), X, Y, T)

    Q = rand(T, m, m)
    Q = Q * Q'  # make sure Q is positive-definite

    test_colwise(SqMahalanobis(Q), X, Y, T)
    test_colwise(Mahalanobis(Q), X, Y, T)
end

function test_pairwise(dist, x, y, T)
    @testset "Pairwise test for $(typeof(dist))" begin
        nx = size(x, 2)
        ny = size(y, 2)
        rxy = zeros(T, nx, ny)
        rxx = zeros(T, nx, nx)
        for j = 1:ny, i = 1:nx
            rxy[i, j] = dist(x[:, i], y[:, j])
        end
        for j = 1:nx, i = 1:nx
            rxx[i, j] = dist(x[:, i], x[:, j])
        end
        # As earlier, we have small rounding errors in accumulations
        @test pairwise(dist, x, y, dims=2) ≈ rxy
        @test pairwise(dist, x, dims=2) ≈ rxx
        @test pairwise(dist, permutedims(x), permutedims(y), dims=1) ≈ rxy
        @test pairwise(dist, permutedims(x), dims=1) ≈ rxx
        vecx = (x[:, i] for i in 1:nx)
        vecy = (y[:, i] for i in 1:ny)
        for (vecx, vecy) in ((vecx, vecy), (collect(vecx), collect(vecy)))
            @test pairwise(dist, vecx, vecy) ≈ rxy
            @test pairwise(dist, vecx) ≈ rxx
            @test pairwise!(similar(rxy), dist, vecx, vecy) ≈ rxy
            @test pairwise!(similar(rxx), dist, vecx) ≈ rxx
        end
    end
end

@testset "pairwise metrics on $T" for T in (Float64, F64)
    m = 5
    n = 8
    nx = 6
    ny = 8

    X = rand(T, m, nx)
    Y = rand(T, m, ny)
    A = rand(1:3, m, nx)
    B = rand(1:3, m, ny)

    P = rand(T, m, nx)
    Q = rand(T, m, ny)

    test_pairwise(SqEuclidean(), X, Y, T)
    test_pairwise(Euclidean(), X, Y, T)
    test_pairwise(Cityblock(), X, Y, T)
    test_pairwise(TotalVariation(), X, Y, T)
    test_pairwise(Chebyshev(), X, Y, T)
    test_pairwise(Minkowski(2.5), X, Y, T)
    test_pairwise(Hamming(), A, B, T)

    test_pairwise(CosineDist(), X, Y, T)
    test_pairwise(CosineDist(), A, B, T)
    test_pairwise(CorrDist(), X, Y, T)

    test_pairwise(ChiSqDist(), X, Y, T)
    test_pairwise(KLDivergence(), P, Q, T)
    test_pairwise(RenyiDivergence(0.0), P, Q, T)
    test_pairwise(RenyiDivergence(1.0), P, Q, T)
    test_pairwise(RenyiDivergence(Inf), P, Q, T)
    test_pairwise(RenyiDivergence(0.5), P, Q, T)
    test_pairwise(RenyiDivergence(2), P, Q, T)
    test_pairwise(JSDivergence(), P, Q, T)

    test_pairwise(BhattacharyyaDist(), X, Y, T)
    test_pairwise(HellingerDist(), X, Y, T)
    test_pairwise(BrayCurtis(), X, Y, T)
    test_pairwise(Bregman(x -> sqeuclidean(x, zero(x)), x -> 2*x), X, Y, T)

    w = rand(m)

    test_pairwise(WeightedSqEuclidean(w), X, Y, T)
    test_pairwise(WeightedEuclidean(w), X, Y, T)
    test_pairwise(WeightedCityblock(w), X, Y, T)
    test_pairwise(WeightedMinkowski(w, 2.5), X, Y, T)
    test_pairwise(WeightedHamming(w), A, B, T)
    test_pairwise(PeriodicEuclidean(w), X, Y, T)

    Q = rand(m, m)
    Q = Q * Q'  # make sure Q is positive-definite

    test_pairwise(SqMahalanobis(Q), X, Y, T)
    test_pairwise(Mahalanobis(Q), X, Y, T)

    m, nx, ny = 2, 8, 6

    X = rand(T, m, nx)
    Y = rand(T, m, ny)
    test_pairwise(Haversine(), X, Y, T)
    test_pairwise(SphericalAngle(), X, Y, T)
end

function test_scalar_pairwise(dist, x, y, T)
    @testset "Scalar pairwise test for $(typeof(dist))" begin
        rxy = dist.(x, permutedims(y))
        rxx = dist.(x, permutedims(x))
        # As earlier, we have small rounding errors in accumulations
        @test pairwise(dist, x, y) ≈ rxy
        @test pairwise(dist, x) ≈ rxx
        @test pairwise(dist, permutedims(x), permutedims(y), dims=2) ≈ rxy
        @test pairwise(dist, permutedims(x), dims=2) ≈ rxx
        @test_throws DimensionMismatch pairwise(dist, permutedims(x), permutedims(y), dims=1)
    end
end

@testset "scalar pairwise metrics on $T" for T in (Float64, F64)
    m = 5
    n = 8
    nx = 6
    ny = 8
    x = rand(T, nx)
    y = rand(T, ny)
    a = rand(1:3, nx)
    b = rand(1:3, ny)
    test_scalar_pairwise(SqEuclidean(), x, y, T)
    test_scalar_pairwise(Euclidean(), x, y, T)
    test_scalar_pairwise(Cityblock(), x, y, T)
    test_scalar_pairwise(TotalVariation(), x, y, T)
    test_scalar_pairwise(Chebyshev(), x, y, T)
    test_scalar_pairwise(Minkowski(2.5), x, y, T)
    test_scalar_pairwise(Hamming(), a, b, T)
    test_scalar_pairwise(CosineDist(), x, y, T)
    test_scalar_pairwise(CosineDist(), a, b, T)
    test_scalar_pairwise(ChiSqDist(), x, y, T)
    test_scalar_pairwise(KLDivergence(), x, y, T)
    test_scalar_pairwise(JSDivergence(), x, y, T)
    test_scalar_pairwise(BrayCurtis(), x, y, T)
    w = rand(1, 1)
    test_scalar_pairwise(WeightedSqEuclidean(w), x, y, T)
    test_scalar_pairwise(WeightedEuclidean(w), x, y, T)
    test_scalar_pairwise(WeightedCityblock(w), x, y, T)
    test_scalar_pairwise(WeightedMinkowski(w, 2.5), x, y, T)
    test_scalar_pairwise(WeightedHamming(w), a, b, T)
    test_scalar_pairwise(PeriodicEuclidean(w), x, y, T)
end

@testset "Euclidean precision" begin
    X = [0.1 0.2; 0.3 0.4; -0.1 -0.1]
    pd = pairwise(Euclidean(1e-12), X, X, dims=2)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(Euclidean(1e-12), X, dims=2)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(SqEuclidean(1e-12), X, X, dims=2)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
    pd = pairwise(SqEuclidean(1e-12), X, dims=2)
    @test pd[1, 1] == 0
    @test pd[2, 2] == 0
end

@testset "non-negativity" begin
    X = [0.3 0.3 + eps()]

    @test all(x -> x >= 0, pairwise(SqEuclidean(), X; dims = 2))
    @test all(x -> x >= 0, pairwise(SqEuclidean(), X, X; dims = 2))
    @test all(x -> x >= 0, pairwise(Euclidean(), X; dims = 2))
    @test all(x -> x >= 0, pairwise(Euclidean(), X, X; dims = 2))
    @test all(x -> x >= 0, pairwise(WeightedSqEuclidean([1.0]), X; dims = 2))
    @test all(x -> x >= 0, pairwise(WeightedSqEuclidean([1.0]), X, X; dims = 2))
    @test all(x -> x >= 0, pairwise(SqMahalanobis(ones(1, 1)), X; dims = 2))
    @test all(x -> x >= 0, pairwise(SqMahalanobis(ones(1, 1)), X, X; dims = 2))
end

@testset "Bregman Divergence" begin
    # Some basic tests.
    @test_throws ArgumentError bregman(x -> x, x -> 2*x, [1, 2, 3], [1, 2, 3])
    # Test if Bregman() correctly implements the gkl divergence between two random vectors.
    F(p) = LinearAlgebra.dot(p, log.(p));
    ∇F(p) = map(x -> log(x) + 1, p)
    testDist = Bregman(F, ∇F)
    p = rand(4)
    q = rand(4)
    p = p/sum(p)
    q = q/sum(q)
    for (p, q) in ((p, q),
                   (Iterators.take(p, 4), Iterators.take(q, 4)),
                   ((p[i] for i in 1:4), (q[i] for i in 1:4)),
                   )
        @test testDist(p, q) ≈ gkl_divergence(p, q)
        # Test if Bregman() correctly implements the squared euclidean dist. between them.
        @test bregman(x -> norm(x)^2, x -> 2 .* x, p, q) ≈ sqeuclidean(p, q)
    end
    # Test if Bregman() correctly implements the IS distance.
    G(p) = -1 * sum(log.(p))
    ∇G(p) = map(x -> -1 * x^(-1), p)
    function ISdist(p::AbstractVector, q::AbstractVector)
        return sum([p[i]/q[i] - log(p[i]/q[i]) - 1 for i in 1:length(p)])
    end
    @test bregman(G, ∇G, p, q) ≈ ISdist(p, q)
end

@testset "Unitful vectors" begin
    x = [1m, 2m, 3m]; y = [2m, 3m, 4m]; w = [1, 1, 1]; p = [2m, 2m, 2m]
    @test @inferred sqeuclidean(x, y) == 3m^2
    @test @inferred euclidean(x, y) == sqrt(3)m
    @test @inferred cityblock(x, y) == 3m
    @test @inferred totalvariation(x, y) == 1.5m
    @test @inferred chebyshev(x, y) == 1m
    @test @inferred minkowski(x, y, 2) == sqrt(3)m
    @test @inferred jaccard(x, y) == 1 - sum(min.(x, y)) / sum(max.(x, y))
    @test @inferred braycurtis(x, y) == sum(abs.(x .- y)) / sum(abs.(x .+ y))
    @test @inferred cosine_dist(x, y) == 1 - dot(x, y) / (norm(x) * norm(y))
    @test @inferred corr_dist(x, y) == cosine_dist(x .- mean(x), y .- mean(y))
    @test @inferred chisq_dist(x, y) == sum((x .- y).^2 ./ (x .+ y))
    @test @inferred spannorm_dist(x, y) == 0m
    @test @inferred hellinger(x, y) == sqrt(1 - sum(sqrt.(x .* y) / sqrt(sum(x) * sum(y))))
    @test @inferred meanad(x, y) == 1m
    @test @inferred msd(x, y) == 1m^2
    @test @inferred rmsd(x, y) == 1m
    @test @inferred nrmsd(x, y) == rmsd(x, y) / (maximum(x) - minimum(x))
    @test @inferred weuclidean(x, y, w) == euclidean(x, y)
    @test @inferred wsqeuclidean(x, y, w) == sqeuclidean(x, y)
    @test @inferred wcityblock(x, y, w) == cityblock(x, y)
    @test @inferred wminkowski(x, y, w, 2) == euclidean(x, y)
    @test @inferred whamming(x, y, w) == hamming(x, y)
    @test @inferred peuclidean(x, y, p) == sqrt(3)m
end

#=
@testset "zero allocation colwise!" begin
    d = Euclidean()
    a = rand(2, 41)
    b = rand(2, 41)
    z = zeros(41)
    colwise!(z, d, a, b)
    # This fails when bounds checking is enforced
    bounds = Base.JLOptions().check_bounds
    if bounds == 0
        @test (@allocated colwise!(z, d, a, b)) == 0
    else
        @test_broken (@allocated colwise!(z, d, a, b)) == 0
    end
end
=#


text:
"""
EXP2 fileformat

See EUROCONTROL NEST Manual for EXP2 fileformat description

Exp2 DataFrame Column Names
#   CellName        Field           Type    Comment
1   :ADEP           origin          String  ICAO code (airport)
2   :ADES           destination     String  ICAO code (airport)
4   :ACTYPE                         String
5   :RFL            RFL             Int64   Requested Flight Level
6   :ZONE_ORIG      zone origin     String  ICAO code (could be the same as airport)
7   :ZONE_DEST      zone destin     String  ICAO code (could be the same as airport)
8   :FLIGHT_ID      flight ID       Int64   SAMAD ID or SAAM ID or ... (must be unique)
9   :ETD_DATETIME                   DateTime Estimated Time and Date of Departure
11  :ETA_DATETIME                   DateTime Estimated Time and Date of Arrival (Not used = !!!!)
12  :CALLSIGN                       String  Char <= 7
13  :COMPANY                        String  Generally the 3 first letters of the callsign
15  :UUID       Universal Unique ID String  CALLSIGN-ADEP-ADES-EOBD-EOBT
16  :FIPS_CLONED    Fips cloned     String  Y/N
18  :FLIGHT_SAAM_ID Flight SAAM ID  Int64   Unique for the day
19  :FLIGHT_SAMAD_ID Flight SAMAD ID Int64  Unique for ever
20  :TACT_ID                        Int64   CFMU ID: Check with TACTID on ALL-FT
21  :SSR_CODE                       String  plan identifiant (12bits mode A), Secondary Surveillance Rader
22  :REGISTRATION                   String
23  :PTD_DATETIME                   DateTime Planned Date and Time Departure
25  :ATFM_DELAY                     Air Traffic Flow Management Delay: Difference between calculated by CASA take off time (CTOT) and estimated take off time (ETOT)
26  :REROUTING_STATE                String
27  :MOST_PEN_REG                   String  Most penelizing regulation. If no regulation 'X' is present: Check with mostPenalizingRegulationID
28  :TYPE_OF_FLIGHT                 String  Letters indicating type of flight (military ?)
29  :EQUIPMENT                      String  Letters indicating equipment of flight
30  :ICAO_EQUIP                     String
31  :COM_EQUIP                      String
32  :NAV_EQUIP                      String
33  :SSR_EQUIP                      String  Secondary Surveillance Radar equipment
34  :SURVIVAL_EQUIP                 String
35  :PERSONS_ON_BOARD               Int64   0 means no information
36  :TOP_FL                         Int64
37  :MAX_RFL                        Int64   Requested Flight Level
38  :FLT_PLN_SOURCE                 String
40  :AOBT                           DateTime Actual Of Block in Time
41  :IFPSID                         String  Prisme (Integrated Initial Flight Plan Processing System)
42  :IOBT                           DateTime Initial Of Block in Time | FLF (planned time departure)
43  :ORIGFLIGHTDATAQUALITY          String  NON, PFD, RPL, FPL
44  :FLIGHTDATAQUALITY              String  FLF (flt plan source), NON, PFD, RPL, FPL
45  :SOURCE                         String  UNK, FPL, RPL, AFI, MFS, FNM, AFP, DIV
46  :EXEMPTREASON                   String  NEXE, EMER, SERE, HEAD, AEAP, empty
47  :EXEMPTREASONDIST               String  NEXE, LONG, empty .Note in NEST the flight exemption is
                                            limited to exempted or not exempted states. Not exempted
                                            corresponds to NEXT in fields 46 and 47. All other values
                                            correspon to exempted. During export, NEST writes empty
                                            fields in 46 and 47 for exempted flights.
48  :LATEFILER                      String  Y/N
49  :LATEUPDATER                    String  Y/N
50  :NORTHATLANTIC                  String  Y/N
51  :COBT                           DateTime Computed Of Block in Time
52  :EOBT                           DateTime Estimate Of Block in Time
53  :FLIGHTSTATE                    String  NE, PL, PS, PR, SR, FI, FS, SI, TA, AA, CA, TE
54  :PREV2ACTIVATIONFLIGHTSTATE     String  NE, PL, PS, PR, SR, FI, FS, SI, TA, AA, CA, TE
55  :SUSPENSIONSTATUS               String  NS, ST, SM, RC, TV, NR, RV
57  :SAMCTOT                        DateTime
58  :SAMSENT                        String  Y/N
59  :SIPCTOT                        DateTime
60  :SIPSENT                        String  Y/N
61  :SLOTFORCED                     String  Y/N
62  :MOSTPENALIZINGREGID            String  Is it the same as most pen reg on FLF
63  :REGAFFECTEDBYNROFINST          Int64
64  :EXCLFROMNROFINST               Int64
65  :LASTRECEIVEDATFMMESSAGETITLE   String  DES, ERR, FCM, FUM, FLS, REA, RFI, RJT, RRN, RRP, SAM, SIP, SLC, SMM, SPA, SRJ, SRM, SWM, UNK
66  :LASTRECEIVEDMESSAGETITLE       String  ABI, ACH, ACT, APL, ARR, CAN, CHG, CNL, DEP, DLA, ERR, EST, FPL, FSA, MFS, PAC, PFD, RPL, UNK
67  :LASTSENTATFMMESSAGETITLE       String  DES, ERR, FCM, FUM, FLS, REA, RFI, RJT, RRN, RRP, SAM, SIP, SLC, SMM, SPA, SRJ, SRM, SWM, UNK
68  :MANUALEXEMPREASON              String  N/S/R
69  :SENSITIVEFLIGHT                String  Y/N
70  :READYFORIMPROVEMENT            String  Y/N
71  :READYFORDEP                    String  Y/N
72  :REVISEDTAXITIME                Int64   0..999999
73  :TIS                            Int64   Time to intert the sequence, 0..999999
74  :TRS                            Int64   Time to remove the sequence, 0..999999
75  :TOBESENTSLOTMESSAGE            String  Related to flight progress, ABI, ACH, ACT, APL, ARR, CAN, CHG, CNL, DEP, DLA, ERR, EST, FPL, FSA, MFS, PAC, PFD, RPL, UNK
76  :TOBESENTPROPMESSAGETITLE       String  Related to flight progress, ABI, ACH, ACT, APL, ARR, CAN, CHG, CNL, DEP, DLA, ERR, EST, FPL, FSA, MFS, PAC, PFD, RPL, UNK
77  :LASTSENTSLOTMESSAGETITLE       String  Related to flight progress, ABI, ACH, ACT, APL, ARR, CAN, CHG, CNL, DEP, DLA, ERR, EST, FPL, FSA, MFS, PAC, PFD, RPL, UNK
78  :LASTSENTPROPMESSAGETITLE       String  Related to flight progress, ABI, ACH, ACT, APL, ARR, CAN, CHG, CNL, DEP, DLA, ERR, EST, FPL, FSA, MFS, PAC, PFD, RPL, UNK
79  :LASTSENTSLOTMESSAGE            DateTime Related to flight progress
80  :LASTSENTPROPMESSAGE            DateTime Related to flight progress
81  :FLIGHTCOUNTOPTION              String  Indicates which flight plan should be / has been ussed when doing flight/count related operations. Used in TACT queries and replies P,N
82  :NORMALFLIGHTTACT_ID            Int64
83  :PROPFLIGHTTACT_ID              Int64
84  :OPERATINGACOPERICAOID          String
85  :REROUTINGWHY                   String  N/M/C/A/O
86  :REROUTINGLFIGHTSTATE           String  P/E/T/R/V/N
87  :RVR                            Int64
88  :FTFMAIRAC                      Int64   Filed traffic flight model (TFM=Profile) airac..
89  :FTFMENVBASELINENUM             Int64   Filed traffic flight model (TFM=Profile) env..
90  :RTFMAIRAC                      Int64   Regulated traffic flight model (TFM=Profile) airac..
91  :RTFMENVBASELINENUM             Int64   Regulated traffic flight model (TFM=Profile) env..
92  :CTFMAIRAC                      Int64   Computer traffic flight model (TFM=Profile) airac..
93  :CTFMENVBASELINENUM             Int64   Computer traffic flight model (TFM=Profile) env..
94  :LASTRECEIVEDPROGRESSMESSAGE    String  DPI, EMPTY, SIZE
"""
module Exp2

export read

#TODO Keep unused fields?
using ..util  # new
using Format
using CSV
using Dates
using DataFrames

const exp2_fileformat = Dict(1=>String, 2=>String, 3=>String, 4=>String,
5=>Int64, 6=>String, 7=>String, 8=>Int64, 9=>String, 10=>String, 11=>String,
12=>String, 13=>String, 14=>Int64, 15=>String, 16=>String, 17=>Int64,
18=>Int64, 19=>Int64, 20=>Int64, 21=>String, 22=>String, 23=>String,
24=>String, 25=>Int64, 26=>String, 27=>String, 28=>String, 29=>String,
30=>String, 31=>String, 32=>String, 33=>String, 34=>String, 35=>Int64,
36=>Int64, 37=>Int64, 38=>String, 39=>Int64, 40=>String, 41=>String, 42=>String,
43=>String, 44=>String, 45=>String, 46=>String, 47=>String, 48=>String,
49=>String, 50=>String, 51=>String, 52=>String, 53=>String, 54=>String, 55=>String,
56=>Int64, 57=>String, 58=>String, 59=>String, 60=>String, 61=>String, 62=>String,
63=>Int64, 64=>Int64, 65=>String, 66=>String, 67=>String, 68=>String, 69=>String,
70=>String, 71=>String, 72=>Int64, 73=>Int64, 74=>Int64, 75=>String, 76=>String,
77=>String, 78=>String, 79=>String, 80=>String, 81=>String, 82=>Int64, 83=>Int64,
84=>String, 85=>String, 86=>String, 87=>Int64, 88=>Int64, 89=>Int64, 90=>Int64,
91=>Int64, 92=>Int64, 93=>Int64, 94=>String, 95=>Int64)

const header_format = ["ADEP", "ADES", "NOT_USED1", "ACTYPE", "RFL",
"ZONE_ORIG", "ZONE_DEST", "FLIGHT_ID", "DEP_DATE", "DEP_TIME",
"ARR_TIME", "CALLSIGN", "COMPANY", "SEPARATOR1", "UUID",
"FIPS_CLONED", "SEPARATOR2", "FLIGHT_SAAM_ID", "FLIGHT_SAMAD_ID", "TACT_ID",
"SSR_CODE", "REGISTRATION", "PLAN_DEP_DATE", "PLAN_DEP_TIME", "ATFM_DELAY",
"REROUTING_STATE", "MOST_PEN_REG", "TYPE_OF_FLIGHT", "EQUIPMENT", "ICAO_EQUIP",
"COM_EQUIP", "NAV_EQUIP", "SSR_EQUIP", "SURVIVAL_EQUIP", "PERSONS_ON_BOARD",
"TOP_FL", "MAX_RFL", "FLT_PLN_SOURCE", "SEPARATOR3", "AOBT",
"IFPSID", "IOBT", "ORIGFLIGHTDATAQUALITY", "FLIGHTDATAQUALITY", "SOURCE",
"EXEMPTREASON", "EXEMPTREASONDIST", "LATEFILER", "LATEUPDATER", "NORTHATLANTIC",
"COBT", "EOBT", "FLIGHTSTATE", "PREV2ACTIVATIONFLIGHTSTATE", "SUSPENSIONSTATUS",
"TACT_ID_DUPLICATE", "SAMCTOT", "SAMSENT", "SIPCTOT", "SIPSENT",
"SLOTFORCED", "MOSTPENALIZINGREGID", "REGAFFECTEDBYNROFINST", "EXCLFROMNROFINST",
"LASTRECEIVEDATFMMESSAGETITLE",
"LASTRECEIVEDMESSAGETITLE", "LASTSENTATFMMESSAGETITLE", "MANUALEXEMPREASON",
"SENSITIVEFLIGHT", "READYFORIMPROVEMENT",
"READYFORDEP", "REVISEDTAXITIME", "TIS", "TRS", "TOBESENTSLOTMESSAGE",
"TOBESENTPROPMESSAGETITLE", "LASTSENTSLOTMESSAGETITLE",
"LASTSENTPROPMESSAGETITLE", "LASTSENTSLOTMESSAGE", "LASTSENTPROPMESSAGE",
"FLIGHTCOUNTOPTION", "NORMALFLIGHTTACT_ID", "PROPFLIGHTTACT_ID",
"OPERATINGACOPERICAOID", "REROUTINGWHY",
"REROUTINGLFIGHTSTATE", "RVR", "FTFMAIRAC", "FTFMENVBASELINENUM", "RTFMAIRAC",
"RTFMENVBASELINENUM", "CTFMAIRAC", "CTFMENVBASELINENUM",
"LASTRECEIVEDPROGRESSMESSAGE", "SEPARATOR4"]

const yymmdd = DateFormat("YYmmdd")
const hhmm = DateFormat("HHMM")
const yyyymmddhhmmss = DateFormat("YYYYmmddHHMMSS")
const yymmddhhmm = DateFormat("YYmmddHHMM")
const year2000 = Year(2000)

function read(file)
    df = CSV.read(file, delim=";", types=exp2_fileformat, header=header_format,
    copycols=true)
    reformat!(df);
    return df
end

function eta_datetime(eta_time::String, etd_datetime::DateTime)
    eta_time_type = format_time(eta_time, hhmm)
    eta_datetime = ""

    if eta_time_type === missing
        return missing
    else
        eta_datetime = Date(etd_datetime) + eta_time_type
    end

    if eta_datetime === missing
        return missing
    elseif eta_datetime < etd_datetime
        return eta_datetime + Dates.Day(1)
    else
        eta_datetime
    end
end

function reformat_exp!(df)
    df[:,:ETD_DATETIME] = format_datetime.(df[:,:DEP_DATE] .*
    df[:,:DEP_TIME], yymmddhhmm, addyear=year2000)

    df[:,:ETA_DATETIME] = eta_datetime.(df[:,:ARR_TIME], df[:,:ETD_DATETIME])
end

function reformat_flf!(df)
    df[:,:PTD_DATETIME] = format_datetime.(df[:,:PLAN_DEP_DATE] .*
    df[:,:PLAN_DEP_TIME], yymmddhhmm, addyear=year2000)
end

function reformat_allft!(df)
    df[:,:TEMP] = format_datetime.(df[:,:AOBT], yyyymmddhhmmss)
    select!(df, Not(:AOBT))
    df[:,:AOBT] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:IOBT], yyyymmddhhmmss)
    select!(df, Not(:IOBT))
    df[:,:IOBT] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:COBT], yyyymmddhhmmss)
    select!(df, Not(:COBT))
    df[:,:COBT] =  df[:,:TEMP]
    select!(df, Not(:TEMP))

    df[:,:TEMP] = format_datetime.(df[:,:EOBT], yyyymmddhhmmss)
    select!(df, Not(:EOBT))
    df[:,:EOBT] =  df[:,:TEMP]
    select!(df, Not(:TEMP))
end

function remove_unused!(df)
    select!(df,[:ADEP, :ADES, :ACTYPE, :RFL, :ZONE_ORIG, :ZONE_DEST,
    :FLIGHT_ID, :ETD_DATETIME, :ETA_DATETIME, :CALLSIGN, :COMPANY, :UUID,
    :FIPS_CLONED, :FLIGHT_SAAM_ID, :FLIGHT_SAMAD_ID, :TACT_ID, :SSR_CODE,
    :REGISTRATION, :PTD_DATETIME, :ATFM_DELAY, :REROUTING_STATE, :MOST_PEN_REG,
    :TYPE_OF_FLIGHT, :EQUIPMENT, :ICAO_EQUIP, :COM_EQUIP, :NAV_EQUIP,
    :SSR_EQUIP, :SURVIVAL_EQUIP, :PERSONS_ON_BOARD, :TOP_FL, :MAX_RFL,
    :FLT_PLN_SOURCE, :AOBT, :IFPSID, :IOBT, :ORIGFLIGHTDATAQUALITY,
    :FLIGHTDATAQUALITY, :SOURCE, :EXEMPTREASON, :EXEMPTREASONDIST, :LATEFILER,
    :LATEUPDATER, :NORTHATLANTIC, :COBT, :EOBT, :FLIGHTSTATE,
    :PREV2ACTIVATIONFLIGHTSTATE, :SUSPENSIONSTATUS, :SAMCTOT, :SAMSENT,
    :SIPCTOT, :SIPSENT, :SLOTFORCED, :MOSTPENALIZINGREGID,
    :REGAFFECTEDBYNROFINST, :EXCLFROMNROFINST, :LASTRECEIVEDATFMMESSAGETITLE,
    :LASTRECEIVEDMESSAGETITLE, :LASTSENTATFMMESSAGETITLE, :MANUALEXEMPREASON,
    :SENSITIVEFLIGHT, :READYFORIMPROVEMENT, :READYFORDEP, :REVISEDTAXITIME,
    :TIS, :TRS, :TOBESENTSLOTMESSAGE, :TOBESENTPROPMESSAGETITLE,
    :LASTSENTSLOTMESSAGETITLE, :LASTSENTPROPMESSAGETITLE, :LASTSENTSLOTMESSAGE,
    :LASTSENTPROPMESSAGE, :FLIGHTCOUNTOPTION, :NORMALFLIGHTTACT_ID,
    :PROPFLIGHTTACT_ID, :OPERATINGACOPERICAOID, :REROUTINGWHY,
    :REROUTINGLFIGHTSTATE, :RVR, :FTFMAIRAC, :FTFMENVBASELINENUM, :RTFMAIRAC,
    :RTFMENVBASELINENUM, :CTFMAIRAC, :CTFMENVBASELINENUM,
    :LASTRECEIVEDPROGRESSMESSAGE])
end

function reformat!(df)
    reformat_exp!(df)
    reformat_flf!(df)
    reformat_allft!(df)
    remove_unused!(df)
end

end # module

ext: jl


text:
__precompile__()

module PlotSVG

#using Rsvg
#using Cairo
using Colors
using ColorSchemes
using LinearAlgebra

#* ===================================

export svg_dot, svg_dots, svg_circle, svg_circles, svg_point, svg_points, svg_line, svg_lines, svg_arrows
include("elementary_svg.jl")

export LineStyle, CircFuncDict

LineStyle =    Dict(  :solid=>(false,false),
                      :dashed=>(true,false),
                      :dotted=>(false,true)  ) # (dashed,dotted)

CircFuncDict = Dict(  :dot=>svg_points,
                      :circle=>svg_circles,
                      :circledot=>svg_dot_circles,
                      :crossdot=>svg_cross_circles  )

export scale_helper,  join_svg,  save_svg,  make_svg_str
include("helper_svg.jl")

export plot_vectors, plot_graph, plot_graph_centre
include("plotters.jl")

export scatter
include("scatter.jl")

export plot
include("plot.jl")

end

max_stars_repo_name: algorithmx/PlotSVG

max_issues_repo_name: algorithmx/PlotSVG

max_forks_repo_name: algorithmx/PlotSVG


text:
"""
    mutate([client::Client], mutation_name::Union{Alias, AbstractString}, args::AbstractDict, output_type::Type=Any; kwargs...)

Perform a mutation on the server, using the global client if `client` not supplied.

If no `output_fields` are supplied, none are returned.

By default `mutate` returns a `GQLReponse{Any}`, where the data for an individual mutation
can be found by `gql_response.data[mutation_name]`.

The mutation uses the `endpoint` field of the `client`.

# Arguments
- `client::Client`: GraphQL client (optional). If not supplied, [`global_graphql_client`](@ref) is used.
- `mutation_name::Union{Alias, AbstractString}`: name of mutation_name.
- `args`: dictionary of mutation argument key value pairs - can be nested with
    dictionaries and vectors.
- `output_type::Type=Any`: output data type for query response object. An object of type
    `GQLResponse{output_type}` will be returned.For further information, see documentation
    for `GQLResponse`.

# Keyword Arguments
- `output_fields=String[]`: output fields to be returned. Can be a string, or
    composed of dictionaries and vectors.
- `direct_write=false`: if `true`, the query is formed by generating a string
    from `args` directly, and the introspected schema is not used. Any ENUMs
    must be wrapped in a `GQLEnum`. See [`directly_write_query_args`](@ref) for
    more information.
- `retries=1`: number of times the mutation will be attempted before erroring.
- `readtimeout=0`: HTTP request timeout length. Set to 0 for no timeout.
- `throw_on_execution_error=false`: set to `true` to throw an exception if the GraphQL server
    response contains errors that occurred during execution.

See also: [`GQLResponse`](@ref)
"""
mutate(mutation_name::Union{Alias, AbstractString}, args::AbstractDict, output_type::Type=Any; kwargs...) = mutate(global_graphql_client(), mutation_name, args, output_type; kwargs...)
function mutate(client::Client, mutation_name::Union{Alias, AbstractString}, args::AbstractDict, output_type::Type=Any; output_fields=String[], kwargs...)
    !in(get_name(mutation_name), get_mutations(client)) && throw(GraphQLClientException("$(get_name(mutation_name)) is not an existing mutation"))

    # If output_fields is empty, return is empty
    output_str = get_output_str(output_fields)

    return generic_gql_query(client, "mutation", mutation_name, args, output_str, output_type; kwargs...)
end


text:
"""
`rotations = grid_rotations(maxradians, rgridsz, SD)` generates
a set of rotations (AffineMap) useful for a gridsearch of
possible rotations to align a pair of images.

`maxradians` is either a single maximum angle (in 2D) or a set of
Euler angles (in 3D and higher). `rgridsz` is one or more integers
specifying the number of gridpoints to search in each of the rotation
axes, corresponding with entries in `maxradians`. `SD` is a matrix
specifying the sample spacing.

e.g. `grid_rotations((pi/8,pi/8,pi/8), (3,3,3), Matrix{Float64}(I,3,3))` would
return an array of 27 rotations with 3 possible angles for each
Euler axis: -pi/8, 0, pi/8. Passing `Matrix{Float64}(I,3,3)` for SD indicates
that the resulting  transforms are meant to be applied to an image with isotropic
pixel spacing.
"""
function grid_rotations(maxradians, rgridsz, SD)
    rgridsz = [rgridsz...]
    maxradians = [maxradians...]
    nd = size(SD,1)
    if !all(isodd, rgridsz)
        @warn("rgridsz should be odd; rounding up to the next odd integer")
    end
    for i = 1:length(rgridsz)
        if !isodd(rgridsz[i])
            rgridsz[i] = max(round(Int, rgridsz[i]) + 1, 1)
        end
    end
    grid_radius = map(x->div(x,2), rgridsz)
    if nd > 2
        gridcoords = [range(-grid_radius[x]*maxradians[x], stop=grid_radius[x]*maxradians[x], length=rgridsz[x]) for x=1:nd]
        rotation_angles = Iterators.product(gridcoords...)
    else
        rotation_angles = range(-grid_radius[1]*maxradians[1], stop=grid_radius[1]*maxradians[1], length=rgridsz[1])
    end
    axmat = Matrix{Float64}(I,nd,nd)
    axs = map(x->axmat[:,x], 1:nd)
    tfeye = tformeye(nd)
    output = typeof(tfeye)[]
    for ra in rotation_angles
        if nd > 2
            euler_rots = map(x->tformrotate(x...), zip(axs, ra))
            rot = foldr(∘, euler_rots, init=tfeye)
        elseif nd == 2
            rot = tformrotate(ra)
        else
            error("Unsupported dimensionality")
        end
        push!(output, AffineMap(SD*rot.linear/SD , zeros(nd))) #account for sample spacing
    end
    return output
end

"""
`best_tform, best_mm = rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD =Matrix{Float64}(I,ndims(fixed),ndims(fixed))))`
Tries a grid of rotations to align `moving` to `fixed`.  Also calculates the best translation (`maxshift` pixels
or less) to align the images after performing the rotation. Returns an AffineMap that captures both the
best rotation and shift out of the values searched, along with the mismatch value after applying that transform (`best_mm`).

For more on how the arguments `maxradians`, `rgridsz`, and `SD` influence the search, see the documentation for
`grid_rotations`.
"""
function rotation_gridsearch(fixed, moving, maxshift, maxradians, rgridsz, SD = Matrix{Float64}(I,ndims(fixed),ndims(fixed)))
    rgridsz = [rgridsz...]
    nd = ndims(moving)
    @assert nd == ndims(fixed)
    rots = grid_rotations(maxradians, rgridsz, SD)
    best_mm = Inf
    best_rot = tformeye(nd)
    best_shift = zeros(nd)
    for rot in rots
        new_moving = transform(moving, rot)
        #calc mismatch
        #mm = mismatch(fixed, new_moving, maxshift; normalization=:pixels)
        mm = mismatch(fixed, new_moving, maxshift)
        thresh = 0.1*maximum(x->x.denom, mm)
        best_i = indmin_mismatch(mm, thresh)
        cur_best =ratio(mm[best_i], 0.0)
        if cur_best < best_mm
            best_mm = cur_best
            best_rot = rot
            best_shift = [best_i.I...]
        end
    end
    return tformtranslate(best_shift) ∘ best_rot, best_mm
end


text:
mutable struct StepAnticipationMethod{F1,F2} <: UnsafeLeapAlgorithm
  rates::Vector{Float64}
  total_rate::Float64
  leap_formula::F1
  validate_leap!::F2
  U::SparseMatrixCSC{Int,Int} # reactant stoichiometry
  V::SparseMatrixCSC{Int,Int} # net stoichiometry
  dxdt::Vector{Float64}
  drdt::Vector{Float64}
end

function generate_events!(f::StepAnticipationMethod, v, rates, s)
  map!((rate, drdt) -> pois_rand(rate * s + 1//2 * drdt * s * s), v, rates, f.drdt)
end

function update!(algorithm::StepAnticipationMethod, state, model)
  # unpack
  rates = algorithm.rates
  total_rate = algorithm.total_rate
  V = algorithm.V
  U = algorithm.U

  dxdt = algorithm.dxdt
  drdt = algorithm.drdt

  # update jump rates
  update_jump_rates!(algorithm, state, model)

  # update auxiliary variables
  update!(algorithm.leap_formula, state, model, V, rates, total_rate)

  return nothing
end


text:
module softmax_regression
export h, J, init_theta, g, predict, accuracy

function h(theta, X)
    mat = exp(theta' * X)
    return (1 ./ sum(mat, 1)) .* mat
end

function J(theta, X, y)
    hh = h(theta, X)
    K, m = size(hh)
    I = sub2ind((K,m), collect(y), 1:m)
    return -sum(log(hh[I]))/m
end

init_theta(n, K) = randn(n, K) * 0.01

function g(K, y)
    idx = map(k -> find(i -> i == k, y), 1:K)
    function grad(theta, X, y)
        hh = h(theta, X)
        K, m = size(hh)
        I = sub2ind((K,m), collect(y), 1:m)
        hh[I] -= 1
        return (X * hh') / m
    end

    return grad
end

function predict(theta, X)
    hh = h(theta, X)
    return map(j -> find(i -> i == maximum(hh[:,j]), hh[:,j])[1], 1:size(hh, 2))'
end

accuracy(y_real, y_pred) = sum(map(i -> i == 0 ? 1 : 0, y_real - y_pred)) / size(y_real, 2)
end

max_stars_repo_name: ilkerkesen/ufldl-tutorial

max_issues_repo_name: ilkerkesen/ufldl-tutorial

max_forks_repo_name: ilkerkesen/ufldl-tutorial


text:
using Plots
path = dirname(@__FILE__)

# Initialize the plotting backend
gr()

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:saveroot, tempdir())
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:use_population_forecast, true)
m <= Setting(:forecast_horizons, 12)
m <= Setting(:impulse_response_horizons, 8)

# Run full-distribution forecast
estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")
output_vars = add_requisite_output_vars([:histobs, :forecastobs, :shockdecobs, :irfobs])

@everywhere using DSGE
m <= Setting(:forecast_block_size, 5)
forecast_one(m, :full, :none, output_vars, verbose = :none)
compute_meansbands(m, :full, :none, output_vars; verbose = :none)
println("The following warning is expected test behavior:")

# Plot history and forecast
plot_history_and_forecast(m, :obs_nominalrate, :obs, :full, :none,
                          bdd_and_unbdd = true,
                          start_date = DSGE.quartertodate("2007-Q1"),
                          verbose = :none)

# Plot forecast comparison
plot_forecast_comparison(m, m, :obs_nominalrate, :obs, :full, :none,
                         bdd_and_unbdd = true,
                         start_date = DSGE.quartertodate("2007-Q1"),
                         verbose = :none)

# Plot shock decomposition
plot_shock_decomposition(m, :obs_nominalrate, :obs, :full, :none, verbose = :none)

# Plot IRF
plot_impulse_response(m, :rm_sh, collect(keys(m.observables)), :obs, :full, :none,
                      verbose = :none)
plot_impulse_response(m, :rm_sh, collect(keys(m.observables)), :obs, :full, :none,
                      flip = true, verbose = :none)

# Plot scenario, zeroing out measurement error
for para in [:e_y, :e_π, :e_R]
    m[para].value = 0
end
alt = Scenario(:altscen, "Test Alternative Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")
forecast_scenario(m, alt, verbose = :none)
scenario_means_bands(m, alt, verbose = :none)
plot_scenario(m, :obs_nominalrate, :obs, alt, untrans = true, verbose = :none)
plot_scenario(m, :obs_nominalrate, :obs, alt, verbose = :none)
plot_scenario(m, :obs_nominalrate, :obs, alt, fourquarter = true, verbose = :none)
@testset "Test scenario exception" begin
    @test_throws ErrorException plot_scenario(m, :obs_nominalrate, :obs, alt,
                                              untrans = true, fourquarter = true)
end

# Hair plot
realized = load_data(m, verbose = :none)
hist_mb = read_mb(m, :full, :none, :histobs)
fcast_mb = read_mb(m, :full, :none, :bddforecastobs)
hair_plot(:obs_nominalrate, realized, [hist_mb], [fcast_mb];
          plotroot = saveroot(m), verbose = :none)


text:
###### type definition ##########

import Base: convert, length, show, getindex, start, next, done, isempty, endof

abstract AbstractTimeSeries

immutable TimeArray{T, N, D<:TimeType, A<:AbstractArray} <: AbstractTimeSeries

    timestamp::Vector{D}
    values::A
    colnames::Vector{String}
    meta::Any

    function TimeArray(timestamp::Vector{D},
                       values::AbstractArray{T,N},
                       colnames::Vector{String},
                       meta::Any)
                           nrow, ncol = size(values, 1), size(values, 2)
                           nrow != size(timestamp, 1) ? error("values must match length of timestamp"):
                           ncol != size(colnames,1) ? error("column names must match width of array"):
                           # length(unique(colnames)) != ncol ? colnames = replace_dupes(colnames):
                           timestamp != unique(timestamp) ? error("there are duplicate dates"):
                           ~(flipdim(timestamp, 1) == sort(timestamp) || timestamp == sort(timestamp)) ? error("dates are mangled"):
                           flipdim(timestamp, 1) == sort(timestamp) ?
                           new(flipdim(timestamp, 1), flipdim(values, 1), replace_dupes(colnames), meta):
                           #new(flipdim(timestamp, 1), flipdim(values, 1), colnames, meta):
                           new(timestamp, values, replace_dupes(colnames), meta)
                           #new(timestamp, values, colnames, meta)
    end
end

TimeArray{T,N,D<:TimeType,S<:AbstractString}(d::Vector{D}, v::AbstractArray{T,N}, c::Vector{S}, m::Any) =
        TimeArray{T,N,D,typeof(v)}(d,v,map(String,c),m)
TimeArray{T,N,D<:TimeType,S<:AbstractString}(d::D, v::AbstractArray{T,N}, c::Vector{S}, m::Any) =
        TimeArray{T,N,D,typeof(v)}([d],v,map(String,c),m)

# when no column names are provided - meta is forced to nothing
TimeArray{T,N,D<:TimeType}(d::Vector{D}, v::AbstractArray{T,N}) = TimeArray(d,v,fill("", size(v,2)),nothing)
TimeArray{T,N,D<:TimeType}(d::D, v::AbstractArray{T,N}) = TimeArray([d],v,fill("", size(v,2)),nothing)

# when no meta is provided
TimeArray{T,N,D<:TimeType}(d::Vector{D}, v::AbstractArray{T,N}, c) = TimeArray(d,v,c,nothing)
TimeArray{T,N,D<:TimeType}(d::D, v::AbstractArray{T,N}, c) = TimeArray([d],v,c,nothing)

###### conversion ###############

convert(::Type{TimeArray{Float64,1}}, x::TimeArray{Bool,1}) = TimeArray(x.timestamp, map(Float64, x.values), x.colnames, x.meta)
convert(::Type{TimeArray{Float64,2}}, x::TimeArray{Bool,2}) = TimeArray(x.timestamp, map(Float64, x.values), x.colnames, x.meta)

convert(x::TimeArray{Bool,1}) = convert(TimeArray{Float64,1}, x::TimeArray{Bool,1})
convert(x::TimeArray{Bool,2}) = convert(TimeArray{Float64,2}, x::TimeArray{Bool,2})

###### length ###################

length(ata::AbstractTimeSeries) = length(ata.timestamp)

###### iterator protocol #########

start(ta::TimeArray)   = 1
next(ta::TimeArray,i)  = ((ta.timestamp[i],ta.values[i,:]),i+1)
done(ta::TimeArray,i)  = (i > length(ta))
isempty(ta::TimeArray) = (length(ta) == 0)

###### show #####################

function show{T,N}(io::IO, ta::TimeArray{T,N})

    # variables
    nrow          = size(ta.values, 1)
    ncol          = size(ta.values, 2)
    intcatcher    = falses(ncol)
    for c in 1:ncol
        rowcheck =  trunc(ta.values[:,c]) - ta.values[:,c] .== 0
        if sum(rowcheck) == length(rowcheck)
            intcatcher[c] = true
        end
    end
    spacetime     = nrow > 0 ? strwidth(string(ta.timestamp[1])) + 3 : 3
    firstcolwidth = strwidth(ta.colnames[1])
    colwidth      = Int[]
    for m in 1:ncol
        (T == Bool) || (nrow == 0) ?
        push!(colwidth, max(strwidth(ta.colnames[m]), 5)) :
        push!(colwidth, max(strwidth(ta.colnames[m]), strwidth(@sprintf("%.2f", maximum(ta.values[:,m]))) + DECIMALS - 2))
    end

    # summary line
    @printf(io, "%dx%d %s", nrow, ncol, typeof(ta))
    nrow > 0 && @printf(io, " %s to %s", string(ta.timestamp[1]), string(ta.timestamp[end]))
    println(io, "")
    println(io, "")

   # row label line

    print(io, ^(" ", spacetime), ta.colnames[1], ^(" ", colwidth[1] + 2 -firstcolwidth))

   for p in 2:length(colwidth)
       print(io, ta.colnames[p], ^(" ", colwidth[p] - strwidth(ta.colnames[p]) + 2))
   end
   println(io, "")

  # timestamp and values line
    if nrow > 7
        for i in 1:4
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end

        println(io, '\u22EE')

        for i in nrow-3:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end
    elseif nrow > 0
        for i in 1:nrow
            print(io, ta.timestamp[i], " | ")
        for j in 1:ncol
            T == Bool ?
            print(io, rpad(ta.values[i,j], colwidth[j] + 2, " ")) :
            intcatcher[j] & SHOWINT ?
            print(io, rpad(round(Integer, ta.values[i,j]), colwidth[j] + 2, " ")) :
            isnan(ta.values[i,j]) ?
            print(io, rpad(MISSING, colwidth[j] + 2, " ")) :
            print(io, rpad(round(ta.values[i,j], DECIMALS), colwidth[j] + 2, " "))
        end
        println(io, "")
        end
    end
end

###### getindex #################

# single row
function getindex{T,N,D}(ta::TimeArray{T,N,D}, n::Int)
    # old code for 0.4
    #TimeArray(ta.timestamp[n], ta.values[n,:], ta.colnames, ta.meta)
    # new code for v0.5 to avoid conversion to column vector
    TimeArray(ta.timestamp[n], ta.values[n:n,:], ta.colnames, ta.meta)
end

# single row 1d
function getindex{T,D}(ta::TimeArray{T,1,D}, n::Int)
    TimeArray(ta.timestamp[n], ta.values[[n]], ta.colnames, ta.meta)
end

# range of rows
function getindex{T,N,D}(ta::TimeArray{T,N,D}, r::UnitRange{Int})
    TimeArray(ta.timestamp[r], ta.values[r,:], ta.colnames, ta.meta)
end

# range of 1d rows
function getindex{T,D}(ta::TimeArray{T,1,D}, r::UnitRange{Int})
    TimeArray(ta.timestamp[r], ta.values[r], ta.colnames, ta.meta)
end

# array of rows
function getindex{T,N,D}(ta::TimeArray{T,N,D}, a::Vector{Int})
    TimeArray(ta.timestamp[a], ta.values[a,:], ta.colnames, ta.meta)
end

# array of 1d rows
function getindex{T,D}(ta::TimeArray{T,1,D}, a::Vector{Int})
    TimeArray(ta.timestamp[a], ta.values[a], ta.colnames, ta.meta)
end

# single column by name
function getindex{T,N,D}(ta::TimeArray{T,N,D}, s::AbstractString)
    n = findfirst(ta.colnames, s)
    TimeArray(ta.timestamp, ta.values[:, n], String[s], ta.meta)
end

# array of columns by name
function getindex{T,N,D}(ta::TimeArray{T,N,D}, args::AbstractString...)
    ns = [findfirst(ta.colnames, a) for a in args]
    TimeArray(ta.timestamp, ta.values[:,ns], String[a for a in args], ta.meta)
end

# single date
function getindex{T,N,D}(ta::TimeArray{T,N,D}, d::D)
    idxs = searchsorted(ta.timestamp, d)
    length(idxs) == 1 ? ta[idxs[1]] : nothing
end

# multiple dates
function getindex{T,N,D}(ta::TimeArray{T,N,D}, dates::Vector{D})
    dates = sort(dates)
    idxs, _ = overlaps(ta.timestamp, dates)
    ta[idxs]
end #getindex

getindex{T,N,D}(ta::TimeArray{T,N,D}, r::StepRange{D}) = ta[collect(r)]

getindex{T,N,D}(ta::TimeArray{T,N,D}, k::TimeArray{Bool,1}) = ta[findwhen(k)]

# day of week
# getindex{T,N}(ta::TimeArray{T,N}, d::DAYOFWEEK) = ta[dayofweek(ta.timestamp) .== d]

# Define end keyword
endof(ta::TimeArray) = length(ta.timestamp)

# helper methods for inner constructor
function find_dupes_index(cnames)
    idx = Int[]
    for c in 1:length(cnames)
        if contains(string(cnames[1:c-1]), cnames[c])
            push!(idx, c)
        end
    end
    idx
end

function replace_dupes(cnames)
    n=1
    while length(unique(cnames)) != length(cnames)
        ds = find_dupes_index(cnames)
        for d in ds
            if n == 1
                cnames[d] = string(cnames[d], "_$n")
            else
                #cnames[d] = string(cnames[d][1:length(cnames[d])-2], "$n")
                cnames[d] = string(cnames[d][1:length(cnames[d])-length(string(n))-1], "_$n")
            end
        end
    n +=1
    end
    cnames
end


text:
#=
simple_module:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-06
=#

module SimpleModule
    # values
    hoge = 2

    # methods
    function hello(name)
        println("Hello $(name)!!")
    end

    function goodbye(name)
        println("Goodbye $(name)!!")
    end

    function main()
        hello("Julia")
        println("SimpleModule.hoge = ", hoge)
        goodbye("Julia")
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .SimpleModule
    SimpleModule.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
@testset "542.01-matrix.jl" begin
    @test update_matrix(hcat([[0, 0, 0], [0, 1, 0], [0, 0, 0]]...)) ==
          hcat([[0, 0, 0], [0, 1, 0], [0, 0, 0]]...)
    @test update_matrix(hcat([[0, 0, 0], [0, 1, 0], [1, 1, 1]]...)) ==
          hcat([[0, 0, 0], [0, 1, 0], [1, 2, 1]]...)
end


text:
module KLDivergences

using Distributions: StatsBase
using Distributions
using LinearAlgebra
using PDMats
using Distances
using SpecialFunctions
using StatsBase: StatsBase, kldivergence

export KL, symmetricKL

"""
    KL(p::Distribution, q::Distribution) -> T
    KL(p::Distribution, q::Distribution, n_samples::Int) -> T

Return the KL divergence of  KL(p||q), either by sampling or analytically
"""
KL

function KLbase(p, q, x)
    # We assume that p(x) > 0 since x is sampled from p
    logpdf(p, x) - logpdf(q, x)
end

## Generic fallback for multivariate Distributions
function KL(p::UnivariateDistribution, q::UnivariateDistribution, n_samples = 1_000)
    return mean(x->KLbase(p, q, x), rand(p, n_samples))
end

function KL(p::MultivariateDistribution, q::MultivariateDistribution, n_samples = 1_000)
    length(p) == length(q) ||
        throw(DimensionMismatch("Dimensions of p and q do not match"))
    return mean(x->KLbase(p, q, x), eachcol(rand(p, n_samples)))
end

function symmetricKL(p, q)
    return KL(p, q) + KL(q, p)
end

include("univariate.jl")
include("multivariate.jl")
include("horrible_ad_workaround.jl")

end


text:
using Base.Test, PyCall, Compat
import Compat.String

if isdefined(Base, :Iterators)
    filter(f, itr) = collect(Iterators.filter(f, itr))
    filter(f, d::Associative) = Base.filter(f, d)
end

PYTHONPATH=get(ENV,"PYTHONPATH","")
PYTHONHOME=get(ENV,"PYTHONHOME","")
PYTHONEXECUTABLE=get(ENV,"PYTHONEXECUTABLE","")
info("Python version $pyversion from $(PyCall.libpython), PYTHONHOME=$(PyCall.PYTHONHOME)\nENV[PYTHONPATH]=$PYTHONPATH\nENV[PYTHONHOME]=$PYTHONHOME\nENV[PYTHONEXECUTABLE]=$PYTHONEXECUTABLE")

roundtrip(T, x) = convert(T, PyObject(x))
roundtrip(x) = roundtrip(PyAny, x)
roundtripeq(T, x) = roundtrip(T, x) == x
roundtripeq(x) = roundtrip(x) == x

# test handling of type-tuple changes in Julia 0.4
import PyCall.pyany_toany
@test pyany_toany(Int) == Int
@test pyany_toany(PyAny) == Any
@test pyany_toany(Tuple{Int,PyAny}) == Tuple{Int,Any}
@test pyany_toany(Tuple{Int,Tuple{PyAny,Int8}}) == Tuple{Int,Tuple{Any,Int8}}
@test pyany_toany(Tuple{PyAny,Int,Vararg{PyAny}}) == Tuple{Any,Int,Vararg{Any}}

@test roundtripeq(17)
@test roundtripeq(0x39)
@test roundtripeq(true) && roundtripeq(false)
@test roundtripeq(3.14159)
@test roundtripeq(1.3+4.5im)
@test roundtripeq(nothing)
@test roundtripeq("Hello world")
@test roundtripeq("Hëllö")
@test roundtripeq("Hello \0\0\0world")
@test roundtripeq("Hël\0\0lö")
@test roundtripeq(Symbol, :Hello)
@test roundtripeq(C_NULL) && roundtripeq(convert(Ptr{Void}, 12345))
@test roundtripeq([1,3,4,5]) && roundtripeq([1,3.2,"hello",true])
@test roundtripeq([1 2 3;4 5 6]) && roundtripeq([1. 2 3;4 5 6])
@test roundtripeq((1,(3.2,"hello"),true)) && roundtripeq(())
@test roundtripeq(Int32)
@test roundtripeq(Dict(1 => "hello", 2 => "goodbye")) && roundtripeq(Dict())
@test roundtripeq(UInt8[1,3,4,5])
@test roundtrip(3 => 4) == (3,4)
@test roundtrip(Pair{Int,Int}, 3 => 4) == Pair(3,4)

@test pycall(PyObject(x -> x + 1), PyAny, 314158) == 314159
@test PyObject(x -> x + 1)(314158) == 314159
@test PyAny(PyObject(3)) == 3
@test roundtrip(x -> x + 1)(314158) == 314159

testkw(x; y=0) = x + 2*y
@test pycall(PyObject(testkw), PyAny, 314157) == 314157
@test pycall(PyObject(testkw), PyAny, 314157, y=1) == 314159
@test roundtrip(testkw)(314157) == 314157
@test roundtrip(testkw)(314157, y=1) == 314159

# check type stability of pycall with an explicit return type
@inferred pycall(PyObject(1)[:__add__], Int, 2)

if PyCall.npy_initialized
    @test PyArray(PyObject([1. 2 3;4 5 6])) == [1. 2 3;4 5 6]
    let A = rand(Int, 2,3,4), B = rand(Bool, 2,3,4)
        @test convert(PyAny, PyReverseDims(A)) == permutedims(A, [3,2,1])
        @test convert(PyAny, PyReverseDims(BitArray(B))) == permutedims(B, [3,2,1])
    end
end
@test PyVector(PyObject([1,3.2,"hello",true])) == [1,3.2,"hello",true]
@test PyDict(PyObject(Dict(1 => "hello", 2 => "goodbye"))) == Dict(1 => "hello", 2 => "goodbye")
@test roundtripeq(BitArray([true, false, true, true]))

let d = PyDict(Dict(1 => "hello", 34 => "yes" ))
    @test get(d.o, 1) == "hello"
    set!(d.o, 34, "goodbye")
    @test d[34] == "goodbye"
    @test sort!(keys(Int, d)) == sort!(collect(d.o[:keys]())) == sort!(collect(keys(d))) == [1, 34]
    @test eltype(d) == eltype(typeof(d)) == Pair{Int, Compat.ASCIIString}
end

let d = Dict(zip(1:1000, 1:1000)), f = (k,v) -> iseven(k)
    @test filter(f, d) == filter(f, PyDict(d)) == filter!(f, PyDict(d)) ==
          Dict(zip(2:2:1000, 2:2:1000))
end

@test roundtripeq(Any[1 2 3; 4 5 6])
@test roundtripeq([])
@test convert(Array{PyAny,1}, PyObject(Any[1 2 3; 4 5 6])) == Any[Any[1,2,3],Any[4,5,6]]
if PyCall.npy_initialized
    @test roundtripeq(begin A = Array{Int}(); A[1] = 3; A; end)
end
@test convert(PyAny, PyObject(begin A = Array{Any}(); A[1] = 3; A; end)) == 3

array2py2arrayeq(x) = PyCall.py2array(Float64,PyCall.array2py(x)) == x
@test array2py2arrayeq(rand(3))
@test array2py2arrayeq(rand(3,4))
@test array2py2arrayeq(rand(3,4,5))

@test roundtripeq(2:10) && roundtripeq(10:-1:2)
@test roundtrip(2:2.0:10) == convert(Vector{Float64}, 2:2.0:10)

@pyimport math
@test math.sin(3) ≈ sin(3)

@test collect(PyObject([1,"hello",5])) == [1,"hello",5]

@test try @eval (@pyimport os.path) catch ex
    isa(ex, ArgumentError)
end

@test PyObject("hello") == PyObject("hello")
@test PyObject("hello") != PyObject("hellö")
@test PyObject(hash) == PyObject(hash)
@test PyObject(hash) != PyObject(println)
@test hash(PyObject("hello")) == hash(PyObject("hello"))
@test hash(PyObject("hello")) != hash(PyObject("hellö"))
@test hash(PyObject("hello")) != hash("hellö")
@test hash(PyObject(hash)) == hash(PyObject(hash))
@test hash(PyObject(hash)) != hash(PyObject(println))
@test hash(PyObject(hash)) != hash(hash)

# issue #92:
let x = PyVector(PyAny[])
    py"lambda x: x.append(\"bar\")"(x)
    @test x == ["bar"]
end

if pyversion >= v"2.7" && isdefined(PyCall, :PyDateTime_CAPI)
    @test roundtripeq(Dates.Date(2012,3,4))
    @test roundtripeq(Dates.DateTime(2012,3,4, 7,8,9,11))
    @test roundtripeq(Dates.Millisecond(typemax(Int32)))
    @test roundtripeq(Dates.Millisecond(typemin(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemax(Int32)))
    @test roundtripeq(Dates.Second, Dates.Second(typemin(Int32)))
    @test roundtripeq(Dates.Day, Dates.Day(999999999)) # max allowed day timedelta
    @test roundtripeq(Dates.Day, Dates.Day(-999999999)) # min allowed day timedelta
end

# fixme: is there any nontrivial mimewritable test we can do?
@test !mimewritable("text/html", PyObject(1))

# in Python 3, we need a specific encoding to write strings or bufferize them
# (http://stackoverflow.com/questions/5471158/typeerror-str-does-not-support-the-buffer-interface)
pyutf8(s::PyObject) = pycall(s["encode"], PyObject, "utf-8")
pyutf8(s::String) = pyutf8(PyObject(s))

# IO (issue #107)
#@test roundtripeq(STDOUT) # No longer true since #250
let buf = IOBuffer(false, true), obuf = PyObject(buf)
    @test !obuf[:isatty]()
    @test obuf[:writable]()
    @test !obuf[:readable]()
    @test obuf[:seekable]()
    obuf[:write](pyutf8("hello"))
    obuf[:flush]()  # should be a no-op, since there's no flushing IOBuffer
    @test position(buf) == obuf[:tell]() == 5
    let p = obuf[:seek](-2, 1)
        @test p == position(buf) == 3
    end
    let p = obuf[:seek](0, 0)
        @test p == position(buf) == 0
    end
    @test String(take!(buf)) == "hello"
    obuf[:writelines](["first\n", "second\n", "third"])
    @test String(take!(buf)) == "first\nsecond\nthird"
    obuf[:write](convert(Vector{UInt8}, "möre stuff"))
    @test String(take!(buf)) == "möre stuff"
    @test isopen(buf) == !obuf[:closed] == true
    obuf[:close]()
    @test isopen(buf) == !obuf[:closed] == false
end
let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
    @test !obuf[:writable]()
    @test obuf[:readable]()
    @test obuf[:readlines]() == readlines(IOBuffer("hello\nagain"))
end
let buf = IOBuffer("hello\nagain"), obuf = PyObject(buf)
    @test obuf[:read](5) == convert(Vector{UInt8}, "hello")
    @test obuf[:readall]() == convert(Vector{UInt8}, "\nagain")
end
let buf = IOBuffer("hello\nagain"), obuf = PyTextIO(buf)
    @test obuf[:encoding] == "UTF-8"
    @test obuf[:read](3) == "hel"
    @test obuf[:readall]() == "lo\nagain"
end
let nm = tempname()
    open(nm, "w") do f
        # @test roundtripeq(f)  # PR #250
        pf = PyObject(f)
        @test pf[:fileno]() == fd(f)
        @test pf[:writable]()
        @test !pf[:readable]()
        pf[:write](pyutf8(nm))
        pf[:flush]()
    end
    @test readstring(nm) == nm
end

# issue #112
@test roundtripeq(Array, [1,2,3,4])
@test roundtripeq(Array{Int8}, [1,2,3,4])

# conversion of numpy scalars
pyanycheck(x::Any) = pyanycheck(typeof(x), PyObject(x))
pyanycheck(T, o::PyObject) = isa(convert(PyAny, o), T)
@test pyanycheck(Int, PyVector{PyObject}(PyObject([1]))[1])
@test pyanycheck(Float64, PyVector{PyObject}(PyObject([1.3]))[1])
@test pyanycheck(Complex128, PyVector{PyObject}(PyObject([1.3+1im]))[1])

pymodule_exists(s::AbstractString) = try
    pyimport(s)
    true
catch
    false
end

# bigfloat conversion
if pymodule_exists("mpmath")
    for x in (big(pi), big(pi) + im/big(pi))
        @test pyanycheck(x)
        # conversion may not be exact since it goes through a decimal string
        @test abs(roundtrip(x) - x) < eps(BigFloat) * 1e3 * abs(x)
    end
end
@test convert(BigInt, PyObject(1234)) == 1234

# buffers
let b = PyCall.PyBuffer(pyutf8("test string"))
    @test ndims(b) == 1
    @test (length(b),) == (length("test string"),) == (size(b, 1),) == size(b)
    @test stride(b, 1) == 1
    @test PyCall.iscontiguous(b) == true
end

let o = PyObject(1+2im)
    @test haskey(o, :real)
    @test :real in keys(o)
    @test o[:real] == 1
end

# []-based sequence access
let a1=[5,8,6], a2=rand(3,4), a3=rand(3,4,5), o1=PyObject(a1), o2=PyObject(a2), o3=PyObject(a3)
    @test [o1[i] for i in eachindex(a1)] == a1
    @test [o1[end-(i-1)] for i in eachindex(a1)] == reverse(a1)
    @test o2[1] == collect(a2[1,:])
    @test length(o1) == length(o2) == length(o3) == 3
    o1[end-1] = 7
    @test o1[2] == 7

    # multiple indices are passed as tuples, but this is apparently
    # only supported by numpy arrays.
    if PyCall.npy_initialized
        @test [o2[i,j] for i=1:3, j=1:4] == a2
        @test [o3[i,j,k] for i=1:3, j=1:4, k=1:5] == a3
        @test o3[2,3] == collect(a3[2,3,:])
        o2[2,3] = 8
        @test o2[2,3] == 8
        o3[2,3,4] = 9
        @test o3[2,3,4] == 9
    end
end

# list operations:
let o = PyObject(Any[8,3])
    @test collect(push!(o, 5)) == [8,3,5]
    @test pop!(o) == 5 && collect(o) == [8,3]
    @test shift!(o) == 8 && collect(o) == [3]
    @test collect(unshift!(o, 9)) == [9,3]
    @test collect(prepend!(o, [5,4,2])) == [5,4,2,9,3]
    @test collect(append!(o, [1,6,8])) == [5,4,2,9,3,1,6,8]
    @test isempty(empty!(o))
end

# issue #216:
@test length(collect(pyimport("itertools")[:combinations]([1,2,3],2))) == 3

# PyNULL and copy!
let x = PyNULL(), y = copy!(x, PyObject(314159))
    @test Int(x) == Int(y) == 314159
end

@test pyimport_conda("inspect", "not a conda package").o != C_NULL
import Conda
if PyCall.conda
    # import pyzmq to test PR #294
    let already_installed = "pyzmq" ∈ Conda._installed_packages()
        @test pyimport_conda("zmq", "pyzmq").o != C_NULL
        @test "pyzmq" ∈ Conda._installed_packages()
        if !already_installed
            Conda.rm("pyzmq")
        end
    end
end

let x = 7
    py"""
    def myfun(x):
        return x + $x
    """
    @test py"1 + 2" == 3
    @test py"1 + $x" == 8
    @test py"1 + $(x^2)" == 50
    @test py"myfun"(10) == 17
end

# issue #352
let x = "1+1"
    @test py"$x" == "1+1"
    @test py"$$x" == py"$$(x)" == 2
    @test py"7 - $$x - 7" == 0 # evaluates "7 - 1 + 1 - 7"
    @test py"7 - ($$x) - 7" == -2 # evaluates "7 - (1 + 1) - 7"
    @test py"1 + $$(x[1:2]) 3" == 5 # evals 1 + 1+ 3
end

# Float16 support:
if PyCall.npy_initialized
    @test roundtripeq(Float16[17 18 Inf -Inf -0.0 0.0])
    @test isa(roundtrip(Float16[17]), Vector{Float16})
end

"""
foobar doc
"""
foobar(x) = x+1

# function attributes
let o = PyObject(foobar)
    @test o[:__name__] == o[:func_name] == string(foobar)
    @test o[:__doc__] == o[:func_doc] == "foobar doc\n"
    @test o[:__module__] == o[:__defaults__] == o[:func_defaults] ==
          o[:__closure__] == o[:func_closure] == nothing
end

# issue #345
let weakdict = pyimport("weakref")["WeakValueDictionary"]
    # (use weakdict for the value, since Python supports
    #  weak references to type objects)
    @test weakdict(Dict(3=>weakdict)) == Dict(3=>weakdict)
end

# Expose python docs to Julia doc system
if isdefined(Docs, :getdoc)
    py"""
    def foo():
        "foo docstring"
        return 0
    """
    foo = py"foo"
    # use 'content' since `Text` objects test equality by object identity
    @test @doc(foo).content == "foo docstring"
end


text:
""" Gibbs sampling for linear regression

    To include an intercept, X must contain ones in first column.
    model: y = Xb + e

    b0 is prior mean for b, iB0 = prior precision matrix for b (inv of cov. matrix)
    a0 and d0 are priors for precision tau ~ Gamma(a0/2, d0/2).

    M = mcmc iterations
    tau is the starting value (= 1/sigma2)

    usage: gsreg(y,X)   # uses default uninformative prior

           gsreg(y,X,M=m,tau=t,b0=priorb,iB0 = invpriorcovb,do=priorgammab,a0=priorgammaa)
           Note: iB0 = prior precision matrix = inv(prior variance matrix)
                 b0 must be a column Vector
"""
function gsreg(y::Array{Float64},X::Array{Float64}; M=100000::Int64, burnin = Int(floor(M/10.0))::Int64,tau=[1.0]::Array{Float64},iB0=[0.0001]::Array{Float64},b0=[0.0]::Array{Float64},d0=0.0001::Float64, a0=3.0::Float64)
    ### check for exception
    M::Int64

    # values outside iterations
    n, k = size(X)
    a1 = a0 + n
    # default uninformative priors
    if b0 == [0.0]
        b0 = zeros(k)
    end
    if iB0 == [0.0001]     # coefficient mean and var.
#        B0 = ones(k)*10000.0
#        mB0 = diagm(B0)
#        iB0 = inv(mB0)
         iB0 = Matrix(1.0I, k, k).*0.0001
    end

    bdraws = zeros(M,k)
    s2draws = zeros(M,1)
# Gibbs algorithm
    for i = 1:(M + burnin)

    # draw betas
        Db = inv(X'*X.*tau[1] + iB0) #\eye(k,k)
        db = X'y*tau[1] + iB0*b0
        H = (cholesky(Hermitian(Db))).U
        betas = Db*db + H'*randn(k,1)

    # draw sigma sq.
    #### N.B. second parameter is inverse of Greenberg/Koop defns.!
        d1 = d0 .+ (y-X*betas)'*(y-X*betas)
        tau = rand(Gamma(a1[1]/2,2/d1[1]),1)
        sig2 = 1.0/tau[1]

    # store draws
        if i > burnin
            j = i - burnin
            bdraws[j,:] = transpose(betas)
            s2draws[j] = sig2::Float64
        end

    end
    return bdraws,s2draws
end


text:
struct ConfirmedStatus <: AbstractStatus end
struct UnconfirmedStatus <: AbstractStatus end


text:
Help = Command(
	"help",
	"pangraph help <command>",
	"prints usage information for a given command",
	"one command name",
	Arg[],
	function(args)
		command = parse(Help, args)
		command === nothing  && return usage(Dispatch)
		length(command) == 0 && return usage(Dispatch)
		length(command)  > 1 && return usage(Help)

		for sub in Dispatch.sub
			if sub.cmd == command[1]
				usage(sub)
				return 0
			end
		end

		return usage(Help)
	end
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
using LibSerialPort
using Test

@testset "hello_libserial" begin
    v_uint8_eoc = [0xff, 0xff, 0xff]  # Array{UInt8,1}
    # v_char_eoc = Char.([0xff, 0xff, 0xff])  # Array{Char,1} ['ÿ', 'ÿ', 'ÿ']
    s_eoc = String(copy(v_uint8_eoc))  # "\xff\xff\xff"

    #=
    function _write_end_of_command(ser)
        for i in 0:2
            write(ser, 0xff)
        end
    end
    =#

    function my_readuntil(sp::SerialPort, delim::Vector{T}, timeout_ms::Real) where {T}
        start_time = time_ns()
        out = Vector{T}()
        lastchars = T[0 for i=1:length(delim)]
        while !eof(sp)
            if (time_ns() - start_time)/1e6 > timeout_ms
                break
            end
            if bytesavailable(sp) > 0
                c = read(sp, T)
                push!(out, c)
                lastchars = circshift(lastchars,-1)
                lastchars[end] = c
                if lastchars == delim
                    break
                end
            end
        end
        return out
    end

    function _format_command(cmd)
        cmd * s_eoc
    end

    function _execute_command(sp, cmd)
        @info "-> $cmd"
        delim = "     "
        s = delim * join(cmd, delim)
        @info "->     $s"
        cmd = _format_command(cmd)
        v_cmd_uint8 = transcode(UInt8, cmd)
        @info "-> $v_cmd_uint8"
        write(sp, cmd)
        #_write_end_of_command(sp)
    end

    function _blockin_read(sp::SerialPort, timeout_ms::Integer)
        nb, msg = sp_blocking_read(sp.ref, bytesavailable(sp), timeout_ms)
        @info "<- $msg"
        println("")
        msg
    end

    function main()
        list_ports()

        portname = "/dev/ttyUSB0"
        baudrate = 9600
        timeout_ms = 3000
        sp = open(portname, baudrate)

        _execute_command(sp, "bkcmd=3")
        _blockin_read(sp, timeout_ms)

        _execute_command(sp, "bkcmd=3")
        _blockin_read(sp, timeout_ms)

        _execute_command(sp, "page 0")
        _blockin_read(sp, timeout_ms)
        #sleep(2)
        _execute_command(sp, "page 1")
        _blockin_read(sp, timeout_ms)
        #sleep(2)
        _execute_command(sp, "t0.txt=\"Hello\"")
        _blockin_read(sp, timeout_ms)
        #sleep(1)

        println("Read")
        _execute_command(sp, "get t0.txt")
        #_blockin_read(sp, timeout_ms)
        r = _blockin_read(sp, timeout_ms)

        #r = my_readuntil(sp, v_uint8_eoc, 1000)
        #println(r)
        s = String(r[2:end-3])
        println(s)
        @test s == "Hello"
        close(sp)
    end

    main()
    @test true

end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
const EXISTING_BOILER_EFFICIENCY = 0.8

struct ExistingBoiler <: AbstractThermalTech  # useful to create AbstractHeatingTech or AbstractThermalTech?
    max_kw::Real
    efficiency::Real
    fuel_cost_series::AbstractVector{<:Real}
end

function ExistingBoiler(;
    max_heat_demand_kw::Real=0,
    production_type::String = "hot_water",
    chp_prime_mover::String = "",
    max_thermal_factor_on_peak_load::Real = 1.25,
    efficiency::Real = 0.0,
    fuel_cost_per_mmbtu::Union{<:Real, AbstractVector{<:Real}} = 0.0,
    time_steps_per_hour::Int = 1
    # fuel_type::String = "natural_gas"  # "restrict_to": ["natural_gas", "landfill_bio_gas", "propane", "diesel_oil"],
    # can_supply_steam_turbine::Bool,
    # emissions_factor_lb_CO2_per_mmbtu::Real,
)
    @assert production_type in ["steam", "hot_water"]

    production_type_by_chp_prime_mover = Dict(
        "recip_engine" => "hot_water",
        "micro_turbine" => "hot_water",
        "combustion_turbine" => "steam",
        "fuel_cell" => "hot_water"
    )

    fuel_cost_per_kwh = fuel_cost_per_mmbtu / MMBTU_TO_KWH
    fuel_cost_series = per_hour_value_to_time_series(fuel_cost_per_kwh, time_steps_per_hour,
                                                     "ExistingBoiler.fuel_cost_per_mmbtu")

    efficiency_defaults = Dict(
        "hot_water" => 0.8,
        "steam" => 0.75
    )

    if efficiency == 0.0
        if !isempty(chp_prime_mover)
            production_type = production_type_by_chp_prime_mover[chp_prime_mover]
        end
        efficiency = efficiency_defaults[production_type]
    end

    max_kw = max_heat_demand_kw * max_thermal_factor_on_peak_load

    ExistingBoiler(
        max_kw,
        efficiency,
        fuel_cost_series
    )
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
# Description of a vibronic system that cannot have inter-surface coupling.

"""
System of `M` coupled harmonic oscillators (modes) across `S` surfaces, with no
coupling between surfaces.
"""
struct DiagonalSystem{S,M} <: System{S,M}
    "Energy offsets (S, S)."
    energy::Matrix{Float64}
    "Frequencies (M, S)."
    freq::Matrix{Float64}
    "Linear prefactors (M, S, S)."
    lin::Array{Float64,3}
    "Quadratic prefactors (M, M, S, S)."
    quad::Array{Float64,4}

    "Energy offsets due to linear terms (S)."
    deltas::Vector{Float64}
    "Position offsets due to linear terms (M, S)."
    ds::Matrix{Float64}

    function DiagonalSystem{S,M}(energy::AbstractMatrix{Float64}, freq::AbstractMatrix{Float64}, lin::AbstractArray{Float64,3}, quad::AbstractArray{Float64,4}) where {S,M}
        check_shape(S, M, energy, freq, lin, quad)
        isdiag(energy, lin, quad) || throw(SurfaceCouplingException())

        deltas = zeros(S)
        ds = zeros(M, S)
        for s in 1:S
            for m in 1:M
                deltas[s] += -0.5 * lin[m, s, s].^2 / freq[m, s]
                ds[m, s] = -lin[m, s, s] / freq[m, s]
            end
        end

        new{S,M}(energy, freq, lin, quad, deltas, ds)
    end

    function DiagonalSystem(energy::AbstractMatrix{Float64}, freq::AbstractMatrix{Float64}, lin::AbstractArray{Float64,3}, quad::AbstractArray{Float64,4})
        M, S = size(freq)
        DiagonalSystem{S,M}(energy, freq, lin, quad)
    end
end

diag(sys::DiagonalSystem) = sys

isdiag(sys::DiagonalSystem) = true


text:
struct GradientDescent{IL, L, T, Tprep<:Union{Function, Nothing}} <: FirstOrderOptimizer
    alphaguess!::IL
    linesearch!::L
    P::T
    precondprep!::Tprep
    manifold::Manifold
end

Base.summary(::GradientDescent) = "Gradient Descent"

"""
# Gradient Descent
## Constructor
```julia
GradientDescent(; alphaguess = LineSearches.InitialHagerZhang(),
linesearch = LineSearches.HagerZhang(),
P = nothing,
precondprep = (P, x) -> nothing)
```
Keywords are used to control choice of line search, and preconditioning.

## Description
The `GradientDescent` method a simple gradient descent algorithm, that is the
search direction is simply the negative gradient at the current iterate, and
then a line search step is used to compute the final step. See Nocedal and
Wright (ch. 2.2, 1999) for an explanation of the approach.

## References
 - Nocedal, J. and Wright, S. J. (1999), Numerical optimization. Springer Science 35.67-68: 7.
"""
function GradientDescent(; alphaguess = LineSearches.InitialPrevious(), # TODO: Investigate good defaults.
                           linesearch = LineSearches.HagerZhang(),      # TODO: Investigate good defaults
                           P = nothing,
                           precondprep = (P, x) -> nothing,
                           manifold::Manifold=Flat())
    GradientDescent(alphaguess, linesearch, P, precondprep, manifold)
end

mutable struct GradientDescentState{Tx, T} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
    f_x_previous::T
    s::Tx
    @add_linesearch_fields()
end

function initial_state(method::GradientDescent, options, d, initial_x::AbstractArray{T}) where T
    initial_x = copy(initial_x)
    retract!(method.manifold, initial_x)

    value_gradient!!(d, initial_x)

    project_tangent!(method.manifold, gradient(d), initial_x)

    GradientDescentState(initial_x, # Maintain current state in state.x
                         copy(initial_x), # Maintain previous state in state.x_previous
                         real(T(NaN)), # Store previous f in state.f_x_previous
                         similar(initial_x), # Maintain current search direction in state.s
                         @initial_linesearch()...)
end

function update_state!(d, state::GradientDescentState{T}, method::GradientDescent) where T
    value_gradient!(d, state.x)
    # Search direction is always the negative preconditioned gradient
    project_tangent!(method.manifold, gradient(d), state.x)
    method.precondprep!(method.P, state.x)
    ldiv!(state.s, method.P, gradient(d))
    rmul!(state.s, eltype(state.s)(-1))
    if method.P != nothing
        project_tangent!(method.manifold, state.s, state.x)
    end

    # Determine the distance of movement along the search line
    lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))

    # Update current position # x = x + alpha * s
    @. state.x = state.x + state.alpha * state.s
    retract!(method.manifold, state.x)
    lssuccess == false # break on linesearch error
end

function trace!(tr, d, state, iteration, method::GradientDescent, options, curr_time=time())
  common_trace!(tr, d, state, iteration, method, options, curr_time)
end


text:
function init_setup()
    param = :complexConjug
    θ₀ = [10.0, -5.0, 5.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(param, θ₀...)
    obs = [[0.0, 0.0], [-0.3, -0.2], [0.4, 0.2]]
    tt = [0.0, 1.0, 1.5]
    P̃ = [BSI.FitzhughDiffusionAux(param, θ₀..., tt[1], obs[1], tt[2], obs[2]),
         BSI.FitzhughDiffusionAux(param, θ₀..., tt[2], obs[2], tt[3], obs[3])]
    setup = BSI.DiffusionSetup(P˟, P̃, BSI.PartObs())
    (setup = setup, θ = θ₀, trgt = P˟, obs = obs, tt = tt, aux = P̃)
end

@testset "mcmc schedule" begin
    num_mcmc_steps = 20
    updt_idx = [[1,2,3],[4,5],[6]]
    actions = (save=10, verbose=3, warm_up=5,
               readjust=(x->x%7==0), fuse=(x->false))
    schedule = BSI.MCMCSchedule(num_mcmc_steps, updt_idx, actions)

    @testset "initialisation" begin
        @test schedule.num_mcmc_steps == num_mcmc_steps
        @test schedule.updt_idx == updt_idx
        @test schedule.actions == actions
    end

    @testset "correct iterations" begin
        i = 0
        for step in schedule
            i += 1
            @test step.iter == i
            @test step.idx == updt_idx[mod1(i, length(updt_idx))]
            @test step.save == (i in [10, 20])
            @test step.verbose == (i % actions.verbose == 0)
            @test step.param_updt == (i > actions.warm_up)
            @test step.readjust == (i % 7 == 0)
            @test step.fuse == false
        end
        @test i == 20
    end
end

@testset "workspace" begin
    out = init_setup()
    setup, θ, trgt, obs, tt = out.setup, out.θ, out.trgt, out.obs, out.tt
    aux = out.aux

    L = [1. 0.; 0. 1.]
    Σ = [0.5 0.0; 0.0 1.0]
    BSI.set_observations!(setup, [L, L], [Σ, 2*Σ], obs, tt)

    dt = 0.01
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    BSI.set_imputation_grid!(setup, dt, τ)

    x0_prior = BSI.KnownStartingPt(obs[1])
    BSI.set_x0_prior!(setup, x0_prior)

    BSI.initialise!(Float64, setup, Vern7(), false, NoChangePt())

    ws = BSI.Workspace(setup).workspace
    @testset "initialisation" begin
        @test typeof(ws.Wnr) == Wiener{Float64}
        @test length(ws.XX) == length(ws.XXᵒ) == 2
        @test length(ws.WWᵒ) == length(ws.WW) == 2
        @test all([ws.XX[i].tt == ws.XXᵒ[i].tt for i in 1:2])
        @test all([ws.XX[i].yy == ws.XXᵒ[i].yy for i in 1:2])
        @test all([ws.WW[i].tt == ws.WWᵒ[i].tt for i in 1:2])
        @test all([ws.WW[i].yy == ws.WWᵒ[i].yy for i in 1:2])
        @test eltype(eltype(ws.WW)) == Float64
        @test eltype(eltype(ws.XX)) == SArray{Tuple{2},Float64,1,2}
        @test ws.WW[1].tt[1] == ws.XX[1].tt[1] == tt[1]
        @test ws.WW[2].tt[1] == ws.XX[2].tt[1] == tt[2]
        @test typeof(ws.P[1].Target) <: BSI.FitzhughDiffusion
        @test typeof(ws.P[1].Pt) <: BSI.FitzhughDiffusionAux
        @test all([ws.P[i].tt == ws.WW[i].tt == ws.XX[i].tt for i in 1:2])
    end
end


text:
###########################
# FunctionCell Interface  #
###########################
"""
	cslice(labels [,idx=1])

Returns a cell that slices data piped to it according
to a vector `labels` provided. If multiple labels are available,
`idx` selects which one to slice on. The returned structed is a
`datasubset` of the original data.

# Examples
```
julia> using j4pr

julia> a = datacell([1,2,3,4,5],["a","b","a","b","a"])
DataCell, 5 obs, 1 vars, 1 target(s)/obs, 2 distinct values: "b"(2),"a"(3)

julia> b = a |> cslice(["a"])
[*]DataCell, 3 obs, 1 vars, 1 target(s)/obs, 1 distinct values: "a"(3)

julia> +b
3-element SubArray{Int64,1,Array{Int64,1},Tuple{Array{Int64,1}},false}:
 1
 3
 5

julia> a = datacell([1,2,3],["a" "b" "a"; "c" "c" "d"])
DataCell, 3 obs, 1 vars, 2 target(s)/obs

julia> b = a |> cslice(["d"],2)
[*]DataCell, 1 obs, 1 vars, 2 target(s)/obs

julia> +b
1-element SubArray{Int64,1,Array{Int64,1},Tuple{Array{Int64,1}},false}:
 3
```
"""
cslice(label_list::T where T<:Vector{<:DataElement}, idx::Int=1) = FunctionCell(cslice, (label_list, idx), "Class slicer")

##########################
# DataCell Interface     #
##########################
"""
	cslice(data, labels [,idx=1])

Slices `data` according to a vector `labels` provided and `idx`.
"""
cslice(x::T where T<:CellData, label_list::U where U<:Vector{<:DataElement}, idx::Int=1) = datacell(cslice(strip(x), label_list, idx))

cslice(x::T where T<:AbstractArray, label_list::U where U<:Vector{<:DataElement}, idx::Int=1) = datasubset(x)

cslice(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractVector, label_list::U where U<:Vector{<:DataElement}, idx::Int=1) =
begin
	@assert !any(isna.(x[2]))
	@assert !any(isnan.(x[2]))
	idxs::Vector{Int} = findin(x[2], label_list)
	return datasubset(x,idxs)
end

cslice(x::Tuple{T,S} where T<:AbstractArray where S<:AbstractMatrix, label_list::U where U<:Vector{<:DataElement}, idx::Int=1) =
begin
	@assert !any(isna.(x[2]))
	@assert !any(isnan.(x[2]))
	idxs::Vector{Int} = findin(x[2][idx,:], label_list)
	return datasubset(x,idxs)
end


text:
using FMRIOpto
using Test

@testset "FMRIOpto.jl" begin
    # Write your own tests here.
end