text:
module TestDeprecated

using Test, DataFrames, Random, Logging, Statistics

const ≅ = isequal

old_logger = global_logger(NullLogger())

@testset "deprecated tuple in sort" begin
    dv1 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv2 = [9, 1, 8, missing, 3, 3, 7, missing]
    dv3 = Vector{Union{Int, Missing}}(1:8)
    cv1 = CategoricalArray(dv1, ordered=true)

    d = DataFrame(dv1 = dv1, dv2 = dv2, dv3 = dv3, cv1 = cv1)
    @test sort(d, (:dv1, :cv1))[!, :dv3] == sortperm(dv1)
    @test sort(d, (:dv1, :dv3))[!, :dv3] == sortperm(dv1)

    df = DataFrame(rank=rand(1:12, 1000),
                   chrom=rand(1:24, 1000),
                   pos=rand(1:100000, 1000))
    ds = sort(df, (order(:rank, rev=true),:chrom,:pos))
    @test issorted(ds, (order(:rank, rev=true),:chrom,:pos))

    ds2 = sort(df, (:rank, :chrom, :pos), rev=(true, false, false))
    @test issorted(ds2, (order(:rank, rev=true), :chrom, :pos))

    sort!(df, (:rank, :chrom, :pos), rev=(true, false, false))
    @test issorted(df, (order(:rank, rev=true), :chrom, :pos))
end

@testset "categorical constructor" begin
    df = DataFrame([Int, String], [:a, :b], [false, true], 3)
    @test !(df[:a] isa CategoricalVector)
    @test df[:b] isa CategoricalVector
    @test_throws DimensionMismatch DataFrame([Int, String], [:a, :b], [true], 3)
end

@testset "DataFrame constructors" begin
    df = DataFrame(Union{Int, Missing}, 10, 3)
    @test size(df, 1) == 10
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Int, Missing}}
    @test typeof(df[3]) == Vector{Union{Int, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 3]) == Vector{Union{Int, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[:, 3]) == Vector{Union{String, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])

    df = DataFrame([Union{Int, Missing}, Union{Float64, Missing}, Union{String, Missing}],
                [:A, :B, :C], [false, false, true], 100)
    @test size(df, 1) == 100
    @test size(df, 2) == 3
    @test typeof(df[1]) == Vector{Union{Int, Missing}}
    @test typeof(df[2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[3]) <: CategoricalVector{Union{String, Missing}}
    @test all(ismissing, df[1])
    @test all(ismissing, df[2])
    @test all(ismissing, df[3])
    @test typeof(df[:, 1]) == Vector{Union{Int, Missing}}
    @test typeof(df[:, 2]) == Vector{Union{Float64, Missing}}
    @test typeof(df[:, 3]) <: CategoricalVector{Union{String, Missing}}
    @test all(ismissing, df[:, 1])
    @test all(ismissing, df[:, 2])
    @test all(ismissing, df[:, 3])
end

df = DataFrame(Union{Int, Missing}, 2, 2)
@test size(df) == (2, 2)
@test eltype.(eachcol(df)) == [Union{Int, Missing}, Union{Int, Missing}]

@test df ≅ DataFrame([Union{Int, Missing}, Union{Float64, Missing}], 2)

@testset "colwise" begin
    Random.seed!(1)
    df = DataFrame(a = repeat(Union{Int, Missing}[1, 3, 2, 4], outer=[2]),
                   b = repeat(Union{Int, Missing}[2, 1], outer=[4]),
                   c = Vector{Union{Float64, Missing}}(randn(8)))

    missingfree = DataFrame([collect(1:10)], [:x1])

    @testset "::Function, ::AbstractDataFrame" begin
        cw = colwise(sum, df)
        answer = [20, 12, -0.4283098098931877]
        @test isa(cw, Vector{Real})
        @test size(cw) == (ncol(df),)
        @test cw == answer

        cw = colwise(sum, missingfree)
        answer = [55]
        @test isa(cw, Array{Int, 1})
        @test size(cw) == (ncol(missingfree),)
        @test cw == answer
    end

    @testset "::Function, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise(length, gd) == [[2,2], [2,2]]
    end

    @testset "::Vector, ::AbstractDataFrame" begin
        cw = colwise([sum], df)
        answer = [20 12 -0.4283098098931877]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length([sum]),ncol(df))
        @test cw == answer

        cw = colwise([sum, minimum], missingfree)
        answer = reshape([55, 1], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length([sum, minimum]), ncol(missingfree))
        @test cw == answer

        cw = colwise([Vector{Union{Int, Missing}}], missingfree)
        answer = reshape([Vector{Union{Int, Missing}}(1:10)], (1,1))
        @test isa(cw, Array{Vector{Union{Int, Missing}},2})
        @test size(cw) == (1, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(["Bob", :Susie], DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Vector, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise([length], gd) == [[2 2], [2 2]]
    end

    @testset "::Tuple, ::AbstractDataFrame" begin
        cw = colwise((sum, length), df)
        answer = Any[20 12 -0.4283098098931877; 8 8 8]
        @test isa(cw, Array{Real, 2})
        @test size(cw) == (length((sum, length)), ncol(df))
        @test cw == answer

        cw = colwise((sum, length), missingfree)
        answer = reshape([55, 10], (2,1))
        @test isa(cw, Array{Int, 2})
        @test size(cw) == (length((sum, length)), ncol(missingfree))
        @test cw == answer

        cw = colwise((CategoricalArray, Vector{Union{Int, Missing}}), missingfree)
        answer = reshape([CategoricalArray(1:10), Vector{Union{Int, Missing}}(1:10)],
                         (2, ncol(missingfree)))
        @test typeof(cw) == Array{AbstractVector,2}
        @test size(cw) == (2, ncol(missingfree))
        @test cw == answer

        @test_throws MethodError colwise(("Bob", :Susie), DataFrame(A = 1:10, B = 11:20))
    end

    @testset "::Tuple, ::GroupedDataFrame" begin
        gd = groupby(DataFrame(A = [:A, :A, :B, :B], B = 1:4), :A)
        @test colwise((length), gd) == [[2,2],[2,2]]
    end
end

@testset "deletecols and deletecols!" begin
    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, :a) == DataFrame(b=[3.0, 4.0])
    deletecols!(df, :a)
    @test df == DataFrame(b=[3.0, 4.0])

    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, Not(2)) == DataFrame(b=[3.0, 4.0])
    deletecols!(df, Not(2))
    @test df == DataFrame(b=[3.0, 4.0])

    df = DataFrame(a=[1,2], b=[3.0, 4.0])
    @test deletecols(df, :a, copycols=false)[1] === df.b
    @test deletecols(df, []) == df
    @test deletecols(df, Not([])) == DataFrame()
end

@testset "df[col] and df[col] for getindex, view, and setindex" begin
    @testset "getindex DataFrame" begin
        df = DataFrame(a=1:3, b=4:6, c=7:9)
        @test df[1] == [1, 2, 3]
        @test df[1] === eachcol(df)[1]
        @test df[1:2] == DataFrame(a=1:3, b=4:6)
        @test df[r"[ab]"] == DataFrame(a=1:3, b=4:6)
        @test df[Not(Not(r"[ab]"))] == DataFrame(a=1:3, b=4:6)
        @test df[Not(3)] == DataFrame(a=1:3, b=4:6)
        @test eachcol(df)[1] === df[1]
        @test eachcol(view(df,1:2))[1] == eachcol(df)[1]
        @test eachcol(df[1:2])[1] == eachcol(df)[1]
        @test eachcol(df[r"[ab]"])[1] == eachcol(df)[1]
        @test eachcol(df[Not(Not(r"[ab]"))])[1] == eachcol(df)[1]
        @test eachcol(df[Not(r"[c]")])[1] == eachcol(df)[1]
        @test eachcol(df[1:2])[1] !== eachcol(df)[1]
        @test df[:] == df
        @test df[r""] == df
        @test df[Not(Not(r""))] == df
        @test df[Not(1:0)] == df
        @test df[:] !== df
        @test df[r""] !== df
        @test df[Not(Not(r""))] !== df
        @test df[Not(1:0)] !== df
        @test eachcol(view(df, :))[1] == eachcol(df)[1]
        @test eachcol(df[:])[1] == eachcol(df)[1]
        @test eachcol(df[r""])[1] == eachcol(df)[1]
        @test eachcol(df[Not(1:0)])[1] == eachcol(df)[1]
        @test eachcol(df[:])[1] !== eachcol(df)[1]
        @test eachcol(df[r""])[1] !== eachcol(df)[1]
        @test eachcol(df[Not(1:0)])[1] !== eachcol(df)[1]
    end
    @testset "getindex df[col] and df[cols]" begin
        x = [1, 2, 3]
        df = DataFrame(x=x, copycols=false)
        @test df[:x] === x
        @test df[[:x]].x !== x
        @test df[:].x !== x
        @test df[r"x"].x !== x
        @test df[r""].x !== x
        @test df[Not(1:0)].x !== x
    end
    @testset "view DataFrame" begin
        df = DataFrame(a=1:3, b=4:6, c=7:9)
        @test view(df, 1) == [1, 2, 3]
        @test view(df, 1) isa SubArray
        @test view(df, 1:2) isa SubDataFrame
        @test view(df, 1:2) == df[1:2]
        @test view(df, r"[ab]") isa SubDataFrame
        @test view(df, r"[ab]") == df[1:2]
        @test view(df, Not(Not(r"[ab]"))) isa SubDataFrame
        @test view(df, Not(Not(r"[ab]"))) == df[1:2]
        @test view(df, :) isa SubDataFrame
        @test view(df, :) == df
        @test parent(view(df, :)) === df
        @test view(df, r"") isa SubDataFrame
        @test view(df, r"") == df
        @test parent(view(df, r"")) === df
        @test view(df, Not(1:0)) isa SubDataFrame
        @test view(df, Not(1:0)) == df
        @test parent(view(df, Not(1:0))) === df
    end
    @testset "getindex SubDataFrame" begin
        df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
        sdf = view(df, 2:4, 2:4)
        @test sdf[1] == [1, 2, 3]
        @test sdf[1] isa SubArray
        @test sdf[1:2] == DataFrame(a=1:3, b=4:6)
        @test sdf[1:2] isa SubDataFrame
        @test sdf[r"[ab]"] == DataFrame(a=1:3, b=4:6)
        @test sdf[r"[ab]"] isa SubDataFrame
        @test sdf[Not(Not(r"[ab]"))] == DataFrame(a=1:3, b=4:6)
        @test sdf[Not(Not(r"[ab]"))] isa SubDataFrame
        @test sdf[:] == df[2:4, 2:4]
        @test sdf[:] isa SubDataFrame
        @test sdf[r""] == df[2:4, 2:4]
        @test sdf[r""] isa SubDataFrame
        @test sdf[Not(1:0)] == df[2:4, 2:4]
        @test sdf[Not(1:0)] isa SubDataFrame
        @test parent(sdf[:]) === parent(sdf)
        @test parent(sdf[r""]) === parent(sdf)
        @test parent(sdf[Not([])]) === parent(sdf)
    end
    @testset "view SubDataFrame" begin
        df = DataFrame(x=-1:3, a=0:4, b=3:7, c=6:10, d=9:13)
        sdf = view(df, 2:4, 2:4)
        @test view(sdf, 1) == [1, 2, 3]
        @test view(sdf, 1) isa SubArray
        @test view(sdf, 1:2) isa SubDataFrame
        @test view(sdf, 2:3) == df[2:4, 3:4]
        @test view(sdf, r"[ab]") isa SubDataFrame
        @test view(sdf, r"[ab]") == df[2:4, r"[ab]"]
        @test view(sdf, Not(Not(r"[ab]"))) isa SubDataFrame
        @test view(sdf, Not(Not(r"[ab]"))) == df[2:4, Not(Not(r"[ab]"))]
        @test view(sdf, :) isa SubDataFrame
        @test view(sdf, :) == df[2:4, 2:4]
        @test view(sdf, r"") isa SubDataFrame
        @test view(sdf, r"") == df[2:4, 2:4]
        @test view(sdf, Not(1:0)) isa SubDataFrame
        @test view(sdf, Not(1:0)) == df[2:4, 2:4]
        @test parent(view(sdf, :)) == parent(sdf)
        @test parent(view(sdf, r"")) == parent(sdf)
        @test parent(view(sdf, Not(1:0))) == parent(sdf)
    end

    @testset "old setindex! tests" begin
        df = DataFrame(reshape(1:12, 4, :))
        df[1, :] = df[1:1, :]

        df = DataFrame(reshape(1:12, 4, :))

        # Scalar broadcasting assignment of rows
        df[1:2, :] = 1
        df[[true,false,false,true], :] = 3

        # Vector broadcasting assignment of rows
        df[1:2, :] = [2,3]
        df[[true,false,false,true], :] = [2,3]

        # Broadcasting assignment of columns
        df[:, 1] = 1
        df[:x3] = 2

        # assignment of subtables
        df[1, 1:2] = df[2:2, 2:3]
        df[1:2, 1:2] = df[2:3, 2:3]
        df[[true,false,false,true], 2:3] = df[1:2,1:2]

        # scalar broadcasting assignment of subtables
        df[1:2, 1:2] = 3
        df[[true,false,false,true], 2:3] = 3

        # vector broadcasting assignment of subtables
        df[1:2, 1:2] = [3,2]
        df[[true,false,false,true], 2:3] = [2,3]

        # test of 1-row DataFrame assignment
        df = DataFrame([1 2 3])
        df[1, 2:3] = DataFrame([11 12])
        @test df == DataFrame([1 11 12])

        df = DataFrame([1 2 3])
        df[1, [false, true, true]] = DataFrame([11 12])
        @test df == DataFrame([1 11 12])
    end

    @testset "old test/dataframes.jl tests" begin
        df = DataFrame(a = Union{Int, Missing}[2, 3],
                    b = Union{DataFrame, Missing}[DataFrame(c = 1), DataFrame(d = 2)])
        df[1, :b][:e] = 5

        x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])

        #test_group("DataFrame assignment")
        # Insert single column
        x0 = x[Int[], :]
        @test_throws ArgumentError x0[:d] = [1]
        @test_throws ArgumentError x0[:d] = 1:3

        # Insert single value
        x[:d] = 3
        @test x[:d] == [3, 3, 3]
    end

    @testset "test getindex using df[col] and df[cols] syntax" begin
        x = [1]
        y = [1]
        df = DataFrame(x=x, y=y, copycols=false)
        @test df.x === x
        @test df[:y] === y
        @test df[1] === x
        @test df[1:1][1] == x
        @test df[r"x"][1] == x
        @test df[1:1][1] !== x
        @test df[r"x"][1] !== x
        @test df[1:2][:y] == y
        @test df[1:2][:y] !== y
        @test df[r""][:y] == y
        @test df[r""][:y] !== y
        @test df[:][:x] == x
        @test df[:][:x] !== x
        @test df[[:y,:x]][:x] == x
        @test df[[:y,:x]][:x] !== x
    end

    @testset "setindex! special cases" begin
        df = DataFrame(rand(3,2), [:x3, :x3_1])
        @test_throws ArgumentError df[3] = [1, 2]
        @test_throws ArgumentError df[4] = [1, 2, 3]
        df[3] = [1,2,3]
        df[4] = [1,2,3]
        @test propertynames(df) == [:x3, :x3_1, :x3_2, :x4]
        df = DataFrame()
        @test_throws MethodError df[true] = 1
        @test_throws MethodError df[true] = [1,2,3]
        @test_throws MethodError df[1:2, true] = [1,2]
        @test_throws MethodError df[1, true] = 1
        @test_throws ArgumentError df[1, 100] = 1
        @test_throws BoundsError df[1:2, 100] = [1,2]
    end

    @testset "handling of end in indexing" begin
        z = DataFrame(rand(4,5))
        for x in [z, view(z, 1:4, :)]
            y = deepcopy(x)
            @test x[end] == x[5]
            @test x[end:end] == x[5:5]
            @test x[end, :] == x[4, :]
            x[end] = 1:4
            y[5] = 1:4
            @test x == y
            x[4:end] = DataFrame([11:14, 21:24])
            y[4] = [11:14;]
            y[5] = [21:24;]
            @test x == y
        end
    end

    @testset "aliasing in indexing" begin
        # columns should not alias if scalar broadcasted
        df = DataFrame(A=[0], B=[0])
        df[1:end] = 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0

        df = DataFrame(A=[0], B=[0])
        df[:, 1:end] = 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0

        # columns should not alias if vector assigned
        df = DataFrame(A=[0], B=[0])
        x = [0.0]
        df[1:end] = x
        x[1] = 1.0
        @test df[1, :A] === 0.0
        @test df[1, :B] === 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0.0

        df = DataFrame(A=[0], B=[0])
        x = [0.0]
        df[:, 1:end] = x
        x[1] = 1.0
        @test df[1, :A] === 0.0
        @test df[1, :B] === 0.0
        df[1, :A] = 1.0
        @test df[1, :B] === 0.0
    end
end

@testset "eltypes" begin
    @test eltypes(DataFrame(x=[1], y=["a"])) == [Int, String]
end

@testset "melt" begin
    mdf = DataFrame(id=[missing,1,2,3], a=1:4, b=1:4)
    @test unstack(melt(mdf, :id), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(melt(mdf, :id), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]
    @test unstack(melt(mdf, Not(Not(:id))), :id, :variable, :value)[1:3,:] == sort(mdf)[1:3,:]
    @test unstack(melt(mdf, Not(Not(:id))), :id, :variable, :value)[:, 2:3] == sort(mdf)[:, 2:3]

    Random.seed!(1234)
    x = DataFrame(rand(100, 50))
    x[!, :id] = [1:99; missing]
    x[!, :id2] = string.("a", x[!, :id])
    x[!, :s] = [i % 2 == 0 ? randstring() : missing for i in 1:100]
    allowmissing!(x, :x1)
    x[1, :x1] = missing
    y = melt(x, [:id, :id2])
    @test y ≅ melt(x, r"id")
    @test y ≅ melt(x, Not(Not(r"id")))

    d1 = DataFrame(a = Array{Union{Int, Missing}}(repeat([1:3;], inner = [4])),
                b = Array{Union{Int, Missing}}(repeat([1:4;], inner = [3])),
                c = Array{Union{Float64, Missing}}(randn(12)),
                d = Array{Union{Float64, Missing}}(randn(12)),
                e = Array{Union{String, Missing}}(map(string, 'a':'l')))
    d1s = stack(d1, [:a, :b])
    d1m = melt(d1, [:c, :d, :e])
    @test d1m == melt(d1, r"[cde]")
    @test d1s == d1m
    d1m = melt(d1[:, [1,3,4]], :a)
    @test propertynames(d1m) == [:a, :variable, :value]
    d1m_named = melt(d1[:, [1,3,4]], :a, variable_name=:letter, value_name=:someval)
    @test propertynames(d1m_named) == [:a, :letter, :someval]
    dx = melt(d1, [], [:a])
    @test dx == melt(d1, r"xxx", r"a")
    @test size(dx) == (12, 2)
    @test propertynames(dx) == [:variable, :value]
    dx = melt(d1, :a, [])
    @test dx == stack(d1, r"xxx", r"a")
    @test size(dx) == (0, 3)
    @test propertynames(dx) == [:a, :variable, :value]
    d1m = melt(d1, [:c, :d, :e], view=true)
    @test d1m == melt(d1, r"[cde]", view=true)
    d1m = melt(d1[:, [1,3,4]], :a, view=true)
    @test propertynames(d1m) == [:a, :variable, :value]
    d1m_named = melt(d1, [:c, :d, :e], variable_name=:letter, value_name=:someval, view=true)
    @test d1m_named == melt(d1, r"[cde]", variable_name=:letter, value_name=:someval, view=true)
    @test propertynames(d1m_named) == [:c, :d, :e, :letter, :someval]
    df1 = melt(DataFrame(rand(10,10)))
    df1[!, :id] = 1:100
    @test size(unstack(df1, :variable, :value)) == (100, 11)
    @test_throws ArgumentError unstack(melt(DataFrame(rand(3,2))), :variable, :value)
end

@testset "insertcols!" begin
    df = DataFrame(x = 1:2)
    @test insertcols!(df, 2, y=2:3) == DataFrame(x=1:2, y=2:3)
end

@testset "join" begin
    name = DataFrame(ID = Union{Int, Missing}[1, 2, 3],
                    Name = Union{String, Missing}["John Doe", "Jane Doe", "Joe Blogs"])
    job = DataFrame(ID = Union{Int, Missing}[1, 2, 2, 4],
                    Job = Union{String, Missing}["Lawyer", "Doctor", "Florist", "Farmer"])
    @test_throws ArgumentError join(name, job)
    @test_throws ArgumentError join(name, job, on=:ID, kind=:other)

    df1 = DataFrame(id=[1,2,3], x=[1,2,3])
    df2 = DataFrame(id=[1,2,4], y=[1,2,4])
    df3 = DataFrame(id=[1,3,4], z=[1,3,4])
    @test_throws ArgumentError join(df1, df2, df3, on=:id, kind=:xxx)
end

@testset "eachcol(df, true)" begin
    df = DataFrame(a=1:3, b=4:6, c=7:9)
    @test eachcol(df)[1] === last(eachcol(df, true)[1])
    @test eachcol(df)[1] === last(eachcol(df, true)[1])

    df = DataFrame(rand(3,4), [:a, :b, :c, :d])
    df2 = DataFrame(eachcol(df, true))
    @test df == df2
    df2 = DataFrame(eachcol(df, true), copycols=false)
    @test df == df2
    @test all(((a,b),) -> a === b, zip(eachcol(df), eachcol(df2)))

    @test Tables.rowtable(df) == Tables.rowtable((;eachcol(df, true)...))
    @test Tables.columntable(df) == Tables.columntable((;eachcol(df, true)...))

    for (a, b, c, d) in zip(Tables.rowtable(df),
                            Tables.namedtupleiterator(eachrow(df)),
                            Tables.namedtupleiterator(eachcol(df)),
                            Tables.namedtupleiterator((;eachcol(df, true)...)))
        @test a isa NamedTuple
        @test a === b === c === d
    end

    @test Tables.getcolumn((;eachcol(df, true)...), 1) == Tables.getcolumn(df, 1)
    @test Tables.getcolumn((;eachcol(df, true)...), :a) == Tables.getcolumn(df, :a)
    @test Tables.columnnames((;eachcol(df, true)...)) == Tuple(Tables.columnnames(df))

    df = DataFrame(A = Vector{Union{Int, Missing}}(1:2), B = Vector{Union{Int, Missing}}(2:3))
    @test size(eachcol(df, true)) == (size(df, 2),)
    @test parent(DataFrame(eachcol(df, true))) == df
    @test names(DataFrame(eachcol(df, true))) == names(df)
    @test IndexStyle(eachcol(df, true)) == IndexLinear()
    @test size(eachcol(df, false)) == (size(df, 2),)
    @test IndexStyle(eachcol(df, false)) == IndexLinear()
    @test length(eachcol(df, true)) == size(df, 2)
    @test length(eachcol(df, false)) == size(df, 2)
    @test eachcol(df, true)[1] == (:A => df[:, 1])
    @test eachcol(df, false)[1] == df[:, 1]
    @test collect(eachcol(df, true)) isa Vector{Pair{Symbol, AbstractVector}}
    @test collect(eachcol(df, true)) == [:A => [1, 2], :B => [2, 3]]
    @test collect(eachcol(df, false)) isa Vector{AbstractVector}
    @test collect(eachcol(df, false)) == [[1, 2], [2, 3]]
    @test eltype(eachcol(df, true)) == Pair{Symbol, AbstractVector}
    @test eltype(eachcol(df, false)) == AbstractVector
    for col in eachcol(df, true)
        @test typeof(col) <: Pair{Symbol, <:AbstractVector}
    end
    for col in eachcol(df, false)
        @test isa(col, AbstractVector)
    end
    @test map(minimum, eachcol(df, false)) == [1, 2]
    @test eltype(map(Vector{Float64}, eachcol(df, false))) == Vector{Float64}

    df = DataFrame([11:16 21:26 31:36 41:46])
    sdf = view(df, [3,1,4], [3,1,4])
    @test eachcol(sdf, true) == eachcol(df[[3,1,4], [3,1,4]], true)
    @test eachcol(sdf, false) == eachcol(df[[3,1,4], [3,1,4]], false)
    @test size(eachcol(sdf, true)) == (3,)
    @test size(eachcol(sdf, false)) == (3,)

    df_base = DataFrame([11:16 21:26 31:36 41:46])
    for df in (df_base, view(df_base, 1:3, 1:3))
        @test df == DataFrame(eachcol(df, true))
    end
end

@testset "deprecated by/combine" begin
    vexp = x -> exp.(x)
    Random.seed!(1)
    df = DataFrame(a = repeat([1, 3, 2, 4], outer=[2]),
                   b = repeat([2, 1], outer=[4]),
                   c = rand(Int, 8))

    @test combine(groupby(df, :a), [:c => sum]) == by(df, :a, c_sum = :c => sum)
    @test combine(groupby(df, :a), [:c => vexp]) == by(df, :a, c_function = :c => vexp)
    @test combine(groupby(df, :a), [:b => sum, :c => sum]) ==
        by(df, :a, b_sum = :b => sum, c_sum = :c => sum)
    @test combine(groupby(df, :a), [:b => vexp, :c => identity]) ==
        by(df, :a, b_function = :b => vexp, c_identity = :c => identity)

    gd = groupby(df, :a)

    @test combine(gd, [:c => sum]) == combine(gd, c_sum = :c => sum)
    @test combine(gd, [:c => vexp]) == combine(gd, c_function = :c => vexp)
    @test combine(gd, [:b => sum, :c => sum]) ==
        combine(gd, b_sum = :b => sum, c_sum = :c => sum) ==
        combine(gd, (:b,) => sum, (:c,) => sum)
    @test combine(gd, [:b => vexp, :c => identity]) ==
        combine(gd, b_function = :b => vexp, c_identity = :c => identity)

    @test map(identity, gd)  == combine(identity, gd, ungroup=false)
    @test map(:b => mean, gd)  == combine(:b => mean, gd, ungroup=false)
    @test map([:b,:c] => x -> x.b+x.c, gd) ==
          combine(AsTable([:b, :c]) => x -> x.b+x.c, gd, ungroup=false)
    @test by(identity, df, :a) == combine(identity, gd)
    @test by(:b => mean, df, :a) == combine(:b => mean, gd)
    @test by([:b,:c] => x -> x.b+x.c, df, :a) ==
          combine(AsTable([:b, :c]) => x -> x.b+x.c, gd)
    @test by(df, :a, identity) == combine(identity, gd)
    @test by(df, :a, :b => mean) == combine(:b => mean, gd)
    @test by(df, :a, [:b,:c] => x -> x.b+x.c) ==
          combine(AsTable([:b, :c]) => x -> x.b+x.c, gd)
    @test by(df, :a, :b => mean, [:b,:c] => x -> x.b+x.c) ==
          combine(gd, :b => mean, AsTable([:b, :c]) => x -> x.b+x.c)
    @test by(df, :a, p = :b => mean, q = [:b,:c] => x -> x.b+x.c) ==
          combine(gd, :b => mean => :p, AsTable([:b, :c]) => (x -> x.b+x.c) => :q)
end

@testset "deprecated aggregate" begin
    df = DataFrame(x1=Int64[1,2,2], x2=Int64[1,1,2], y=Int64[1,2,3])
    @test aggregate(df, sum) == aggregate(df, [], sum) == aggregate(df, 1:0, sum)
    @test aggregate(df, sum) == aggregate(df, [], sum, sort=true, skipmissing=true)

    Random.seed!(1)
    N = 20
    d1 = Vector{Union{Int64, Missing}}(rand(1:2, N))
    d2 = CategoricalArray(["A", "B", missing])[rand(1:3, N)]
    d3 = randn(N)
    d4 = randn(N)
    df7 = DataFrame([d1, d2, d3], [:d1, :d2, :d3])

    @test aggregate(DataFrame(a=1), identity) == DataFrame(a_identity=1)

    df8 = aggregate(df7[:, [1, 3]], sum)
    @test df8[1, :d1_sum] == sum(df7[!, :d1])

    df8 = aggregate(df7, :d2, [sum, length], sort=true)
    @test df8[1:2, :d2] == ["A", "B"]
    @test size(df8, 1) == 3
    @test size(df8, 2) == 5
    @test sum(df8[!, :d1_length]) == N
    @test all(df8[!, :d1_length] .> 0)
    @test df8[!, :d1_length] == [sum(isequal.(d2, "A")), sum(isequal.(d2, "B")), sum(ismissing.(d2))]
    df8′ = aggregate(df7, 2, [sum, length], sort=true)
    @test df8 ≅ df8′
    adf = aggregate(groupby(df7, :d2, sort=true), [sum, length])
    @test df8 ≅ adf
    adf′ = aggregate(groupby(df7, 2, sort=true), [sum, length])
    @test df8 ≅ adf′
    adf = aggregate(groupby(df7, :d2), [sum, length], sort=true)
    @test sort(df8, [:d1_sum, :d3_sum, :d1_length, :d3_length]) ≅ adf
    adf′ = aggregate(groupby(df7, 2), [sum, length], sort=true)
    @test adf ≅ adf′

    df = DataFrame(a = [3, missing, 1], b = [100, 200, 300])
    for dosort in (true, false), doskipmissing in (true, false)
        @test aggregate(df, :a, sum, sort=dosort, skipmissing=doskipmissing) ≅
              aggregate(groupby(df, :a, sort=dosort, skipmissing=doskipmissing), sum)
    end

    # Check column names
    anonf = x -> sum(x)
    adf = aggregate(df7, :d2, [mean, anonf])
    @test propertynames(adf) == [:d2, :d1_mean, :d3_mean,
                                 :d1_function, :d3_function]
    adf = aggregate(df7, :d2, [mean, mean, anonf, anonf])
    @test propertynames(adf) == [:d2, :d1_mean, :d3_mean, :d1_mean_1, :d3_mean_1,
                                 :d1_function, :d3_function, :d1_function_1, :d3_function_1]

    df9 = aggregate(df7, :d2, [sum, length], sort=true)
    @test df9 ≅ df8
    df9′ = aggregate(df7, 2, [sum, length], sort=true)
    @test df9′ ≅ df8
end

@testset "deprecated deleterows!" begin
    @test deleterows!(DataFrame(x=[1, 2]), 1) ==
        deleterows!(DataFrame(x=[1, 2]), [1]) ==
        deleterows!(DataFrame(x=[1, 2]), [true, false]) == DataFrame(x=[2])
end

@testset "by skipmissing and sort" begin
    df = DataFrame(a=[2, 2, missing, missing, 1, 1, 3, 3], b=1:8)
    for dosort in (false, true), doskipmissing in (false, true)
        @test by(df, :a, :b=>sum, sort=dosort, skipmissing=doskipmissing) ≅
            combine(groupby(df, :a, sort=dosort, skipmissing=doskipmissing), :b=>sum)
    end
end

@testset "map skipmissing and sort" begin
    df = DataFrame(a=[2, 2, missing, missing, 1, 1, 3, 3], b=1:8)
    for dosort in (false, true), doskipmissing in (false, true)
        gdf = groupby(df, :a, sort=dosort, skipmissing=doskipmissing)
        @test map(identity, gdf) ≅ combine(identity, gdf, ungroup=false)
        @test map(:b => sum, gdf) ≅ combine(:b => sum, gdf, ungroup=false)
    end
end

@testset "DataFrame!" begin
    x = [1,2,3]
    y = [4,5,6]
    # we allow this to make deprecation message simpler
    @test DataFrame!(x=x, y=y, copycols=true) == DataFrame(x=x,y=y)
    df1 = DataFrame(x=x, y=y)
    df2 = DataFrame!(df1)
    @test df1 == df2
    @test df1.x === df2.x
    @test df1.y === df2.y

    a=[1,2,3]
    df = DataFrame!(:a=>a, :b=>1, :c=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!("a"=>a, "b"=>1, "c"=>1:3)
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!(Dict(:a=>a, :b=>1, :c=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df.a === a

    df = DataFrame!(Dict("a"=>a, "b"=>1, "c"=>1:3))
    @test propertynames(df) == [:a, :b, :c]
    @test df."a" === a

    df = DataFrame!((x, y))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), (:x1, :x2))
    @test propertynames(df) == [:x1, :x2]
    @test df.x1 === x
    @test df.x2 === y

    df = DataFrame!((x, y), ("x1", "x2"))
    @test names(df) == ["x1", "x2"]
    @test df."x1" === x
    @test df."x2" === y

    @test_throws MethodError DataFrame!([1 2; 3 4], copycols=false)
    @test_throws MethodError DataFrame!([1 2; 3 4])
    @test_throws MethodError DataFrame!([Union{Int, Missing}, Union{Float64, Missing}],
                                        [:x1, :x2], 2)
end

global_logger(old_logger)

end # module


text:
## Stokes Examples

f₁ = (x,y)   -> zeros(x)
f₂ = (x,y)   -> zeros(x)
uanalytic = (x,y) -> 20x.*y.^3
vanalytic = (x,y) -> 5x.^4 - 5y.^4
panalytic = (x,y) -> 60x.^2 .*y - 20y.^3# + C
g = (x,y)    -> 0

"Example problem for a homogeneous stationary Stokes equation."
prob_stokes_homogenous = StokesProblem(f₁,f₂,g,uanalytic,vanalytic,panalytic)

f₁ = (x,y)   -> zeros(x)
f₂ = (x,y)   -> zeros(x)
uanalytic = (x,y) -> 0
vanalytic = (x,y) -> 0
panalytic = (x,y) -> 0
g = (x,y)    -> 0
"Example problem for solving the trivial stationary Stokes equation."
prob_stokes_dirichletzero = StokesProblem(f₁,f₂,g,uanalytic,vanalytic,panalytic)

ext: jl


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    MLP(input, spec; hidden_activation = :relu, prefix)

Construct a multi-layer perceptron. A MLP is a multi-layer neural network with
fully connected layers.

# Arguments:
* `input::SymbolicNode`: the input to the mlp.
* `spec`: the mlp specification, a list of hidden dimensions. For example,
          `[128, (512, :sigmoid), 10]`. The number in the list indicate the
          number of hidden units in each layer. A tuple could be used to specify
          the activation of each layer. Otherwise, the default activation will
          be used (except for the last layer).
* `hidden_activation::Symbol`: keyword argument, default `:relu`, indicating
          the default activation for hidden layers. The specification here could be overwritten
          by layer-wise specification in the `spec` argument. Also activation is not
          applied to the last, i.e. the prediction layer. See [`Activation`](@ref) for a
          list of supported activation types.
* `prefix`: keyword argument, default `gensym()`, used as the prefix to
          name the constructed layers.

Returns the constructed MLP.
"""
function MLP(input, spec; hidden_activation::Symbol = :relu, prefix = gensym())
  spec = convert(Vector{Union{Int,Tuple}}, spec)

  n_layer = length(spec)
  for (i, s) in enumerate(spec)
    if isa(s, Tuple)
      n_unit, act_type = s
    else
      n_unit = s
      act_type = hidden_activation
    end
    input = FullyConnected(input, name=Symbol(prefix, "fc$i"), num_hidden=n_unit)
    if i < n_layer || isa(s, Tuple)
      # will not add activation unless the user explicitly specified
      input = Activation(input, name=Symbol(prefix, "$act_type$i"), act_type=act_type)
    end
  end

  return input
end


text:
using StructJuMP

firststage = StructuredModel()
@variable(firststage, x[1:2])
@constraint(firststage, sum(x) == 100)
@NLobjective(firststage, :Min, x[1]^2 + x[2]^2 + x[1]*x[2])

for scen in 1:2
    bl = StructuredModel(parent=firststage, id=scen)
    @variable(bl, y[1:2])
    @constraint(bl, x[3-scen] + sum(y) ≥  0)
    @constraint(bl, x[3-scen] + sum(y) ≤ 50)
    @NLobjective(bl, :Min, y[1]^2 + y[2]^2 + y[1]*y[2])
end


text:
"""
Initialise simulator.
"""
function _initialise(__Problem, state0, __dyn, p, t0, tf)
    tspan = (t0, tf)
    iip = true
    if __Problem == :ODE
        Problem = ODEProblem
    elseif __Problem == :Discrete
        Problem = DiscreteProblem
    else
        error("Not supported Problem: $(__Problem)")
    end
    prob = Problem{iip}(__dyn, state0, tspan, p)  # true: isinplace
    prob
end

"""
Simulator struct.
# NOTICE
- If `p` is not copyable, i.e., `applicable(copy, p) == false`, it would not correctly be logged in log_func.
"""
struct Simulator
    integrator::DEIntegrator
    log_func::Union{Function, Nothing}
    # df::DataFrame
    function Simulator(state0, dyn, p=nothing;
            solver=nothing,
            Problem=:ODE,
            t0=0.0, tf=1.0, kwargs...,
        )
        # DEProblem
        __dyn = (dx, x, p, t) -> dyn(dx, x, p, t)
        prob = _initialise(Problem, state0, __dyn, p, t0, tf)
        log_func = nothing
        # save func (logging func)
        if hasmethod(dyn, Tuple{Any, Any, Any, Any, __LOG_INDICATOR__})
            __log_indicator__ = __LOG_INDICATOR__()  # just an indicator for logging
            log_func = function (x, t, integrator::DEIntegrator; kwargs...)
                x = copy(x)  # `x` merely denotes a "view"; DO NOT copy(integrator.u); if you replace it by copy(integrator.u), we will get step-like responses
                p = applicable(copy, integrator.p) ? copy(integrator.p) : integrator.p
                dyn(zero.(x), x, p, t, __log_indicator__; kwargs...)
            end
        end
        integrator = init(prob, solver; kwargs...)
        integrator
        new(integrator, log_func)
    end
end

"""
Reinitialise simulator.
"""
function DiffEqBase.reinit!(simulator::Simulator, args...; kwargs...)
    DiffEqBase.reinit!(simulator.integrator, args...; kwargs...)
    simulator
end

"""
Step `dt` time.
"""
function DiffEqBase.step!(simulator::Simulator, Δt; stop_at_tdt=true)
    DiffEqBase.step!(simulator.integrator, Δt, stop_at_tdt)
end

"""
Step until `tf`.
"""
function step_until!(simulator::Simulator, tf;
        suppress_termination_warn=true,
        suppress_truncation_warn=false,
    )
    success = true
    integrator = simulator.integrator
    t = integrator.t
    _tf = integrator.sol.prob.tspan[2]
    if (simulator.integrator.tdir > 0 && _tf < tf) || (simulator.integrator.tdir < 0 && _tf > tf)
        if !suppress_truncation_warn
            @warn("step! truncated up to tf = $(_tf)")
        end
        tf = _tf
    end
    if tf ≈ t
        if !suppress_termination_warn
            @warn("step! ignored; simulator seems already terminated")
        end
        success = false
    else
        DiffEqBase.step!(simulator, tf - t)
    end
    success
end

function _push!(df::DataFrame, integrator::DEIntegrator, log_func)
    x = integrator.u
    t = integrator.t
    __logger__ = log_func(x, t, integrator)
    # __log_nt__ = recursive_namedtuple(__logger__)
    push!(df, (; time=t, sol=__logger__))
end

function Base.push!(simulator::Simulator, df::DataFrame, flag=true)
    if flag
        integrator = simulator.integrator
        log_func = simulator.log_func
        _push!(df, integrator, log_func)
    end
end

"""
    solve(simulator; ...)

`solve` for `simulator::Simulator` (similar to that of DifferentialEquations.jl).
This method will automatically reinitialise `simulator`.
"""
function DiffEqBase.solve(simulator::Simulator;
        saveat=nothing,
        savestep = typeof(simulator.integrator.sol.prob) <: ODEProblem ? 0.01 : typeof(simulator.integrator.sol.prob) <: DiscreteProblem ? 1 : nothing,  # ODEProblem -> 0.01, DiscreteProblem -> 1
        callback::DECallback=CallbackSet(),
        kwargs...,
    )
    DiffEqBase.reinit!(simulator)
    prob = simulator.integrator.sol.prob
    solver = simulator.integrator.alg
    log_func = simulator.log_func
    t0, tf = prob.tspan[1], prob.tspan[2]
    if saveat == nothing && savestep == nothing
        saveat = []  # default
    elseif saveat != nothing && savestep == nothing
        # nothing
    elseif saveat == nothing && savestep != nothing
        saveat = t0:savestep:tf
    elseif saveat != nothing && savestep != nothing
        error("Assign values of either `saveat` or `savestep`")
    end
    saved_values = SavedValues(Float64, NamedTuple)
    if log_func != nothing
        cb_save = SavingCallback(log_func, saved_values;
                                 saveat=saveat, tdir=Int(sign(tf-t0)))
        callback = CallbackSet(callback, cb_save)  # save values "after all other callbacks"
    end
    _ = solve(prob, solver;
              callback=callback,
              kwargs...)
    df = DataFrame(
                   time = saved_values.t,
                   sol = saved_values.saveval,
                  )
    df
end

"""
# Notes
- The basic concept is borrowed from [an MRAC example](https://jonniedie.github.io/ComponentArrays.jl/stable/examples/adaptive_control/).
- It is modified to be compatible with [SimulationLogger.jl](https://github.com/JinraeKim/SimulationLogger.jl).
# Limitations
- Conditional method definition is troublesome; see [#16](https://github.com/JinraeKim/FlightSims.jl/issues/16).
Instead of it, I decided to merely define a new method with argument `__log_indicator__::__LOG_INDICATOR__`,
which will provide "empty Dict" in the case of no logging.
"""
maybe_apply(f::Function, x, p, t) = f(x, p, t)
maybe_apply(f, x, p, t) = f
function apply_inputs(func; kwargs...)
    function simfunc(dx, x, p, t)
        func(dx, x, p, t; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...)
    end
    function simfunc(dx, x, p, t, __log_indicator__::__LOG_INDICATOR__)
        if hasmethod(func, Tuple{Any, Any, Any, Any, __LOG_INDICATOR__})
            return func(dx, x, p, t, __log_indicator__; map(f -> maybe_apply(f, x, p, t), (; kwargs...))...)
        else
            return Dict()  # see the above method `NamedTupleTools.namedtuple`
        end
    end
    simfunc
end


text:
module StatGeochemBase

    using NaNStatistics

    # AVX vectorziation tools
    using LoopVectorization
    using VectorizationBase
    using VectorizationBase: Vec

    using SpecialFunctions
    import SpecialFunctions.erf
    erf(x::Vec) = VectorizationBase.verf(x)
    include("Math.jl")
    include("Geochronology.jl")

    using Colors: Color, RGBX, RGB, N0f8
    include("Interpolations.jl")
    include("ArrayStats.jl")

    using IndirectArrays: IndirectArray
    include("Images.jl")
    include("Colormaps.jl")

    import Base.display
    include("Display.jl") # Custom pretty-printing

end

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
immutable Table

    inputs::Vector{Int}
    positives::Vector{Int}
    negatives::Vector{Int}

    function Table(inputs, positives, negatives)

        @assert length(positives) == length(negatives)
        @assert 2^length(inputs) == length(positives)

        new(inputs, positives, negatives)

    end

end

function indexing(vec::Vector{Bool})

    ind = pos = 1

    for bit in vec
        if bit
            ind += pos
        end
        pos <<= 1
    end

    ind

end

function (tab::Table)(x::Vector{Bool}, ties::TieBreak)

    ind = indexing(x[tab.inputs])

    if tab.positives[ind] < tab.negatives[ind]
        false
    elseif tab.positives[ind] > tab.negatives[ind]
        true
    else
        ties == positive
    end

end

function (tab::Table)(x::Vector{Bool}, label::Bool)

    ind = indexing(x[tab.inputs])

    if label
        tab.positives[ind] += 1
    else
        tab.negatives[ind] += 1
    end

end


text:
using Pkg
Pkg.add("GLM")
Pkg.add("StatsModels")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("Lathe")
Pkg.add("Plots")
Pkg.add("MLBase")
Pkg.add("PlotlyJS")
Pkg.add("ORCA")
# ] add https://github.com/diegozea/ROC.jl
# ] add https://github.com/tkf/Kaleido.jl

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
import Base: run

function run(
    policy::AbstractPolicy,
    env::AbstractEnv,
    stop_condition = StopAfterEpisode(1),
    hook = EmptyHook(),
)
    check(policy, env)
    _run(policy, env, stop_condition, hook)
end

"Inject some customized checkings here by overwriting this function"
function check(policy, env) end

function _run(policy::AbstractPolicy, env::AbstractEnv, stop_condition, hook::AbstractHook)

    hook(PRE_EXPERIMENT_STAGE, policy, env)
    policy(PRE_EXPERIMENT_STAGE, env)
    is_stop = false
    while !is_stop
        reset!(env)
        policy(PRE_EPISODE_STAGE, env)
        hook(PRE_EPISODE_STAGE, policy, env)

        while !is_terminated(env) # one episode
            action = policy(env)

            policy(PRE_ACT_STAGE, env, action)
            hook(PRE_ACT_STAGE, policy, env, action)

            env(action)

            policy(POST_ACT_STAGE, env)
            hook(POST_ACT_STAGE, policy, env)

            if stop_condition(policy, env)
                is_stop = true
                break
            end
        end # end of an episode

        if is_terminated(env)
            policy(POST_EPISODE_STAGE, env)  # let the policy see the last observation
            hook(POST_EPISODE_STAGE, policy, env)
        end
    end
    hook(POST_EXPERIMENT_STAGE, policy, env)
    hook
end


text:
# Mutual Information

function _mutualinfo(A::AbstractMatrix{<:Integer}, normed::Bool)
    N = sum(A)
    (N == 0.0) && return 0.0

    rows = sum(A, dims=2)
    cols = sum(A, dims=1)
    entA = entropy(A)
    entArows = entropy(rows)
    entAcols = entropy(cols)

    hck = (entA - entAcols)/N
    hc = entArows/N + log(N)
    hk = entAcols/N + log(N)

    mi = hc - hck
    return if normed
        2*mi/(hc+hk)
    else
        mi
    end
end

"""
    mutualinfo(a, b; normed=true) -> Float64

Compute the *mutual information* between the two clusterings of the same
data points.

`a` and `b` can be either [`ClusteringResult`](@ref) instances or
assignments vectors (`AbstractVector{<:Integer}`).

If `normed` parameter is `true` the return value is the normalized mutual information (symmetric uncertainty),
see "Data Mining Practical Machine Tools and Techniques", Witten & Frank 2005.

# References
> Vinh, Epps, and Bailey, (2009). “Information theoretic measures for clusterings comparison”.
Proceedings of the 26th Annual International Conference on Machine Learning - ICML ‘09.
"""
mutualinfo(a, b; normed::Bool=true) = _mutualinfo(counts(a, b), normed)


text:
import Base: vec

export
ruleSPNonlinearEOutNG,
ruleSPNonlinearEOutNGX,
ruleSPNonlinearEIn1GG,
ruleSPNonlinearEInGX,
ruleMNonlinearEInGX

"""
Concatenate a vector (of vectors and floats) and return with original dimensions (for splitting)
"""
function concatenate(xs::Vector)
    ds = size.(xs) # Extract dimensions
    x = vcat(vec.(xs)...)

    return (x, ds)
end

ForneyLab.vec(d::Float64) = [d] # Extend vectorization to Float

"""
Return local linearization of g around expansion point x_hat
for Nonlinear node with single input interface
"""
function localLinearizationSingleIn(g::Function, x_hat::Float64)
    a = ForwardDiff.derivative(g, x_hat)
    b = g(x_hat) - a*x_hat

    return (a, b)
end

function localLinearizationSingleIn(g::Function, x_hat::Vector{Float64})
    A = ForwardDiff.jacobian(g, x_hat)
    b = g(x_hat) - A*x_hat

    return (A, b)
end

"""
Return local linearization of g around expansion point x_hat
for Nonlinear node with multiple input interfaces
"""
function localLinearizationMultiIn(g::Function, x_hat::Vector{Float64})
    g_unpacked(x::Vector) = g(x...)
    A = ForwardDiff.gradient(g_unpacked, x_hat)'
    b = g(x_hat...) - A*x_hat

    return (A, b)
end

function localLinearizationMultiIn(g::Function, x_hat::Vector{Vector{Float64}})
    (x_cat, ds) = concatenate(x_hat)
    g_unpacked(x::Vector) = g(split(x, ds)...)
    A = ForwardDiff.jacobian(g_unpacked, x_cat)
    b = g(x_hat...) - A*x_cat

    return (A, b)
end

#-----------------------
# Extended Update Rules
#-----------------------

# Forward rule
function ruleSPNonlinearEOutNG(g::Function,
                               msg_out::Nothing,
                               msg_in1::Message{<:Gaussian})

    (m_in1, V_in1) = unsafeMeanCov(msg_in1.dist)
    (A, b) = localLinearizationSingleIn(g, m_in1)
    m = A*m_in1 + b
    V = A*V_in1*A'

    return Message(variateType(m), GaussianMeanVariance, m=m, v=V)
end

# Multi-argument forward rule
function ruleSPNonlinearEOutNGX(g::Function, # Needs to be in front of Vararg
                                msg_out::Nothing,
                                msgs_in::Vararg{Message{<:Gaussian}})

    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearizationMultiIn(g, ms_fw_in)
    (m_fw_in, V_fw_in, _) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)
    m = A*m_fw_in + b
    V = A*V_fw_in*A'

    return Message(variateType(m), GaussianMeanVariance, m=m, v=V)
end

# Backward rule with given inverse
function ruleSPNonlinearEIn1GG(g::Function,
                               g_inv::Function,
                               msg_out::Message{<:Gaussian},
                               msg_in1::Nothing)

    (m_out, V_out) = unsafeMeanCov(msg_out.dist)
    (A, b) = localLinearizationSingleIn(g_inv, m_out)
    m = A*m_out + b
    V = A*V_out*A'

    return Message(variateType(m), GaussianMeanVariance, m=m, v=V)
end

# Multi-argument backward rule with given inverse
function ruleSPNonlinearEInGX(g::Function, # Needs to be in front of Vararg
                              g_inv::Function,
                              msg_out::Message{<:Gaussian},
                              msgs_in::Vararg{Union{Message{<:Gaussian}, Nothing}})

    (ms, Vs) = collectStatistics(msg_out, msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearizationMultiIn(g_inv, ms)
    (mc, Vc) = concatenateGaussianMV(ms, Vs)
    m = A*mc + b
    V = A*Vc*A'

    return Message(variateType(m), GaussianMeanVariance, m=m, v=V)
end

# Backward rule with unknown inverse
function ruleSPNonlinearEIn1GG(g::Function,
                               msg_out::Message{<:Gaussian},
                               msg_in1::Message{<:Gaussian})

    m_in1 = unsafeMean(msg_in1.dist)
    (m_out, W_out) = unsafeMeanPrecision(msg_out.dist)
    (A, b) = localLinearizationSingleIn(g, m_in1)
    xi = A'*W_out*(m_out - b)
    W = A'*W_out*A

    return Message(variateType(xi), GaussianWeightedMeanPrecision, xi=xi, w=W)
end

# Multi-argument backward rule with unknown inverse
function ruleSPNonlinearEInGX(g::Function,
                              inx::Int64, # Index of inbound interface inx
                              msg_out::Message{<:Gaussian},
                              msgs_in::Vararg{Message{<:Gaussian}})

    # Approximate joint inbounds
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearizationMultiIn(g, ms_fw_in)

    (m_fw_in, V_fw_in, ds) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)
    m_fw_out = A*m_fw_in + b
    V_fw_out = A*V_fw_in*A'
    C_fw = V_fw_in*A'

    # RTS Smoother
    (m_bw_out, V_bw_out) = unsafeMeanCov(msg_out.dist)
    (m_in, V_in) = smoothRTS(m_fw_out, V_fw_out, C_fw, m_fw_in, V_fw_in, m_bw_out, V_bw_out)

    # Marginalize joint belief on in's
    (m_inx, V_inx) = marginalizeGaussianMV(m_in, V_in, ds, inx)
    W_inx = cholinv(V_inx) # Convert to canonical statistics
    xi_inx = W_inx*m_inx

    # Divide marginal on inx by forward message
    (xi_fw_inx, W_fw_inx) = unsafeWeightedMeanPrecision(msgs_in[inx].dist)
    xi_bw_inx = xi_inx - xi_fw_inx
    W_bw_inx = W_inx - W_fw_inx # Note: subtraction might lead to posdef violations

    return Message(variateType(xi_bw_inx), GaussianWeightedMeanPrecision, xi=xi_bw_inx, w=W_bw_inx)
end

function ruleMNonlinearEInGX(g::Function,
                             msg_out::Message{<:Gaussian},
                             msgs_in::Vararg{Message{<:Gaussian}})

    # Approximate joint inbounds
    (ms_fw_in, Vs_fw_in) = collectStatistics(msgs_in...) # Returns arrays with individual means and covariances
    (A, b) = localLinearizationMultiIn(g, ms_fw_in)

    (m_fw_in, V_fw_in, _) = concatenateGaussianMV(ms_fw_in, Vs_fw_in)
    m_fw_out = A*m_fw_in + b
    V_fw_out = A*V_fw_in*A'
    C_fw = V_fw_in*A'

    # RTS Smoother
    (m_bw_out, V_bw_out) = unsafeMeanCov(msg_out.dist)
    (m_in, V_in) = smoothRTS(m_fw_out, V_fw_out, C_fw, m_fw_in, V_fw_in, m_bw_out, V_bw_out)

    return ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=m_in, v=V_in)
end


text:
using PSUtil
using Test, Suppressor

using PSUtil: current_system, Linux, MacOS

@testset "PSUtil.jl" begin
    include("cpu.jl")
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Classify features ###
function classify(features::Dict{Int,Feature})
    feats = Dict{Int,Int}()

    for (key, feature) in features
        if haskey(FEATURE_CLASSES, feature.class)
            feats[key] = FEATURE_CLASSES[feature.class]
        end
    end

    return feats
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule rr_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("rr")
JLLWrappers.@generate_main_file("rr", UUID("e86bdf43-55f7-5ea2-9fd0-e7daa2c0f2b4"))
end  # module rr_jll


text:
# Copyright 2014 Mitchell Kember. Subject to the MIT License.
# Google Code Jam Round 1A 2008
# Problem A: Minimum Scalar Product

# Calculates the minimum dot product of the permutations of u and v.
function mindot(u, v)
	sum = 0
	while !isempty(u)
		i = indmin(u)
		j = indmax(v)
		sum += splice!(u, i) * splice!(v, j)
	end
	sum
end

function main()
	n = int(readline())
	for i = 1:n
		readline() # vector length
		u = int(split(readline()))
		v = int(split(readline()))
		println("Case #$i: ", mindot(u, v))
	end
end

main()


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Depth-first visit / traversal

#################################################
#
#  Depth-first visit
#
#################################################
"""
    DepthFirst
## Conventions in Breadth First Search and Depth First Search
### VertexColorMap
- color == 0    => unseen
- color < 0     => examined but not closed
- color > 0     => examined and closed

### EdgeColorMap
- color == 0    => unseen
- color == 1     => examined
"""
mutable struct DepthFirst <: AbstractGraphVisitAlgorithm end

function depth_first_visit_impl!(
    g::AbstractGraph,      # the graph
    stack,                          # an (initialized) stack of vertex
    vertexcolormap::AbstractVertexMap,    # an (initialized) color-map to indicate status of vertices
    edgecolormap::AbstractEdgeMap,      # an (initialized) color-map to indicate status of edges
    visitor::AbstractGraphVisitor)  # the visitor

    while !isempty(stack)
        u, udsts, tstate = pop!(stack)
        found_new_vertex = false

        while !done(udsts, tstate) && !found_new_vertex
            v, tstate = next(udsts, tstate)
            u_color = get(vertexcolormap, u, 0)
            v_color = get(vertexcolormap, v, 0)
            v_edge = Edge(u,v)
            e_color = get(edgecolormap, v_edge, 0)
            examine_neighbor!(visitor, u, v, u_color, v_color, e_color) #no return here

            edgecolormap[v_edge] = 1

            if v_color == 0
                found_new_vertex = true
                vertexcolormap[v] = vertexcolormap[u] - 1 #negative numbers
                discover_vertex!(visitor, v) || return
                push!(stack, (u, udsts, tstate))

                open_vertex!(visitor, v)
                vdsts = out_neighbors(g, v)
                push!(stack, (v, vdsts, start(vdsts)))
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vertexcolormap[u] *= -1
        end
    end
end

function traverse_graph!(
    g::AbstractGraph,
    alg::DepthFirst,
    s::Integer,
    visitor::AbstractGraphVisitor;
    vertexcolormap = Dict{eltype(g), Int}(),
    edgecolormap = DummyEdgeMap())

    T = eltype(g)
    vertexcolormap[s] = -1
    discover_vertex!(visitor, s) || return

    sdsts = out_neighbors(g, s)
    sstate = start(sdsts)
    stack = [(T(s), sdsts, sstate)]

    depth_first_visit_impl!(g, stack, vertexcolormap, edgecolormap, visitor)
end

#################################################
#
#  Useful applications
#
#################################################

# Test whether a graph is cyclic

mutable struct DFSCyclicTestVisitor <: AbstractGraphVisitor
    found_cycle::Bool
    DFSCyclicTestVisitor() = new(false)
end

function examine_neighbor!(
    vis::DFSCyclicTestVisitor,
    u::Integer,
    v::Integer,
    ucolor::Int,
    vcolor::Int,
    ecolor::Int)

    if vcolor < 0 && ecolor == 0
        vis.found_cycle = true
    end
end

discover_vertex!(vis::DFSCyclicTestVisitor, v) = !vis.found_cycle

"""
    is_cyclic(g)

Return `true` if graph `g` contains a cycle.

### Implementation Notes
Uses DFS.
"""
function is_cyclic end
@traitfn is_cyclic(g::::(!IsDirected)) = ne(g) > 0
@traitfn function is_cyclic(g::::IsDirected)
    cmap = zeros(Int, nv(g))
    visitor = DFSCyclicTestVisitor()

    for s in vertices(g)
        if cmap[s] == 0
            traverse_graph!(g, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end
        visitor.found_cycle && return true
    end
    return false
end

# Topological sort using DFS

mutable struct TopologicalSortVisitor{T} <: AbstractGraphVisitor
    vertices::Vector{T}
end

function TopologicalSortVisitor(n::T) where T<:Integer
    vs = Vector{T}()
    sizehint!(vs, n)
    return TopologicalSortVisitor(vs)
end

function examine_neighbor!(visitor::TopologicalSortVisitor, u::Integer, v::Integer, ucolor::Int, vcolor::Int, ecolor::Int)
    (vcolor < 0 && ecolor == 0) && error("The input graph contains at least one loop.")
end

function close_vertex!(visitor::TopologicalSortVisitor, v::Integer)
    push!(visitor.vertices, v)
end

function topological_sort_by_dfs(g::AbstractGraph)
    nvg = nv(g)
    cmap = zeros(Int, nvg)
    visitor = TopologicalSortVisitor(nvg)

    for s in vertices(g)
        if cmap[s] == 0
            traverse_graph!(g, DepthFirst(), s, visitor, vertexcolormap=cmap)
        end
    end

    reverse(visitor.vertices)
end

mutable struct TreeDFSVisitor{T} <:AbstractGraphVisitor
    tree::DiGraph
    predecessor::Vector{T}
end

TreeDFSVisitor(n::T) where T<:Integer = TreeDFSVisitor(DiGraph(n), zeros(T,n))

function examine_neighbor!(visitor::TreeDFSVisitor, u::Integer, v::Integer, ucolor::Int, vcolor::Int, ecolor::Int)
    if (vcolor == 0)
        visitor.predecessor[v] = u
    end
    return true
end

"""
    dfs_tree(g, s)

Return an ordered vector of vertices representing a directed acylic graph based on
depth-first traversal of the graph `g` starting with source vertex `s`.
"""
function dfs_tree(g::AbstractGraph, s::Integer)
    nvg = nv(g)
    visitor = TreeDFSVisitor(nvg)
    traverse_graph!(g, DepthFirst(), s, visitor)
    # visitor = traverse_dfs(g, s, TreeDFSVisitor(nvg))
    h = DiGraph(nvg)
    for (v, u) in enumerate(visitor.predecessor)
        if u != 0
            add_edge!(h, u, v)
        end
    end
    return h
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
export web_chanspec, seis_www, track_on!, track_off!

# Generic handler for getting data by HTTP
# Returns:
#   R::Array{UInt8,1}, either request body or error data
#   parsable::Bool, whether or not R is parsable
function get_http_req(url::String, req_info_str::String, to::Int; status_exception::Bool=false)
  (R::Array{UInt8,1}, parsable::Bool) = try
    req = request(  "GET", url, webhdr,
                    readtimeout = to,
                    status_exception = status_exception  )
    if req.status == 200
      (req.body, true)

    else
      @warn(string("Request failed", req_info_str,
      "\nRESPONSE = ", req.status, " (", statustext(req.status), ")",
      "\n\nHTTP response is in misc[\"data\"]"))
      (Array{UInt8,1}(string(req)), false)
    end

  catch err
    T = typeof(err)
    @warn(  string( "Error thrown", req_info_str,
                    "\n\nERROR = ", T,
                    "\n\nTrying to store error message in misc[\"data\"]"
                    )
          )
    msg_data::Array{UInt8,1} = Array{UInt8,1}( try; string(getfield(err, :response)); catch; try; string(getfield(err, :msg));  catch; ""; end; end )
    (msg_data, false)
  end

  return R, parsable
end

function get_http_post(url::String, body::String, to::Int; status_exception::Bool=false)
  try
    req = request(  "POST", url, webhdr, body,
                    readtimeout = to,
                    status_exception = status_exception)
    if req.status == 200
      return (req.body, true)
    else
      @warn(string( "Request failed!\nURL: ", url, "\nPOST BODY: \n", body, "\n",
                    "RESPONSE: ", req.status, " (", statustext(req.status), ")\n" ) )
      return (Array{UInt8,1}(string(req)), false)
    end

  catch err
    @warn(string( "Error thrown:\nURL: ", url, "\nPOST BODY: \n", body, "\n",
                  "ERROR TYPE: ", typeof(err), "\n" ) )
    msg_data::Array{UInt8,1} = Array{UInt8,1}( try; string(getfield(err, :response)); catch; try; string(getfield(err, :msg));  catch; ""; end; end )
    return (msg_data, false)
  end
end

datareq_summ(src::String, ID::String, d0::String, d1::String) = ("\n" * src *
  " query:\nID = " * ID * "\nSTART = " * d0 * "\nEND = " * d1)

# ============================================================================
# Utility functions not for export
hashfname(str::Array{String,1}, ext::String) = string(hash(str), ".", ext)

# Start tracking channel IDs and data lengths
"""
    track_on!(S::SeisData)

Track changes to S.id, changes to channel structure of S, and the sizes of data
vectors in S.x. Does not track data processing operations to any channel i
unless length(S.x[i]) changes for channel i.

**Warning**: If you have or suspect gapped data in any channel, do not use
ungap! while tracking is active.
"""
function track_on!(S::SeisData)
  if S.n > 0
    ids = unique(getfield(S, :id))
    nxs = zeros(Int64, S.n)
    for i = 1:S.n
      nxs[i] = length(S.x[i])
    end
    S.misc[1]["track"] = (ids, nxs)
  end
  return nothing
end

# Stop tracking channel IDs and data lengths; report which have changed
"""
    u = track_off!(S::SeisData)

Turn off tracking in S and return a boolean vector of which channels have
been added or altered significantly.
"""
function track_off!(S::SeisData)
  k = findfirst([haskey(S.misc[i],"track") for i = 1:S.n])
  if S.n == 0
    return nothing
  elseif k == nothing
    return trues(S.n)
  end
  u = falses(S.n)
  (ids, nxs) = S.misc[k]["track"]
  for (n, id) in enumerate(S.id)
    j = findfirst(x -> x == id, ids)
    if j == nothing
      u[n] = true
    else
      if nxs[j] != length(S.x[n])
        u[n] = true
      end
    end
  end
  delete!(S.misc[k], "track")
  return u
end

function savereq(D::Array{UInt8,1}, ext::String, id::String, s::String)
  if ext == "miniseed"
    ext = "mseed"
  elseif occursin("sac", ext)
    ext = "SAC"
  end
  s_str = int2tstr(tstr2int(s))
  yy = s_str[1:4]
  mm = s_str[6:7]
  dd = s_str[9:10]
  HH = s_str[12:13]
  MM = s_str[15:16]
  SS = s_str[18:19]
  nn = lpad(div(parse(Int64, s_str[21:26]), 1000), 3, '0')
  jj = lpad(md2j(yy, mm, dd), 3, '0')
  fname = join([yy, jj, HH, MM, SS, nn, id, "R", ext], '.')
  safe_isfile(fname) && @warn(string("File ", fname, " contains an identical request; overwriting."))
  f = open(fname, "w")
  write(f, D)
  close(f)
  return fname
end

"""
| String | Source             |
|:------:|:-------------------|
|BGR    | http://eida.bgr.de |
|EMSC   | http://www.seismicportal.eu |
|ETH    | http://eida.ethz.ch |
|GEONET | http://service.geonet.org.nz |
|GFZ    | http://geofon.gfz-potsdam.de |
|ICGC   | http://ws.icgc.cat |
|INGV   | http://webservices.ingv.it |
|IPGP   | http://eida.ipgp.fr |
|IRIS   | http://service.iris.edu |
|IRISPH5| http://service.iris.edu/ph5ws/ |
|ISC    | http://isc-mirror.iris.washington.edu |
|KOERI  | http://eida.koeri.boun.edu.tr |
|LMU    | http://erde.geophysik.uni-muenchen.de |
|NCEDC  | http://service.ncedc.org |
|NIEP   | http://eida-sc3.infp.ro |
|NOA    | http://eida.gein.noa.gr |
|ORFEUS | http://www.orfeus-eu.org |
|RESIF  | http://ws.resif.fr |
|SCEDC  | http://service.scedc.caltech.edu |
|TEXNET | http://rtserve.beg.utexas.edu |
|USGS   | http://earthquake.usgs.gov |
|USP    | http://sismo.iag.usp.br |
"""
seis_www = Dict("BGR" => "http://eida.bgr.de",
                "EMSC" => "http://www.seismicportal.eu",
                "ETH" => "http://eida.ethz.ch",
                "GEONET" => "http://service.geonet.org.nz",
                "GFZ" => "http://geofon.gfz-potsdam.de",
                "ICGC" => "http://ws.icgc.cat",
                "INGV" => "http://webservices.ingv.it",
                "IPGP" => "http://eida.ipgp.fr",
                "IRIS" => "http://service.iris.edu",
                "ISC" => "http://isc-mirror.iris.washington.edu",
                "KOERI" => "http://eida.koeri.boun.edu.tr",
                "LMU" => "http://erde.geophysik.uni-muenchen.de",
                "NCEDC" => "http://service.ncedc.org",
                "NIEP" => "http://eida-sc3.infp.ro",
                "NOA" => "http://eida.gein.noa.gr",
                "ODC" => "http://www.orfeus-eu.org",
                "ORFEUS" => "http://www.orfeus-eu.org",
                "RESIF" => "http://ws.resif.fr",
                "SCEDC" => "http://service.scedc.caltech.edu",
                "TEXNET" => "http://rtserve.beg.utexas.edu",
                "USGS" => "http://earthquake.usgs.gov",
                "USP" => "http://sismo.iag.usp.br")
ph5_www = Dict("IRISPH5" => "http://service.iris.edu")

function fdsn_uhead(src::String)
  if !occursin("PH5",src)
    return haskey(seis_www, src) ? seis_www[src] * "/fdsnws/" : src
  else
    return haskey(ph5_www, src) ? ph5_www[src] * "/ph5ws/" : src
  end
end

@doc """
    web_chanspec

## Specifying Channel IDs in Web Requests
| Str | L  | Meaning              | Example |
|:--- |:---|:-----                |:-----   |
| NET | 2  | Network code         | "IU"    |
| STA | 5  | Station code         | "ANMO"  |
| LOC | 2  | Location identifier  | "00"    |
| CHA | 3  | Channel code         | "BHZ"   |

A channel is uniquely specified by four substrings (NET, STA, LOC, CHA), which
can be formatted as a String or a String array. Each substring has a maximum
safe length of `L` characters (column 2 in the table).

### Acceptable Channel ID Formats
| Type              | Example                                           |
|:-----             |:-----                                             |
| String            | "PB.B004.01.BS1, PB.B004.01.BS2"                  |
| Array{String, 1}  | ["PB.B004.01.BS1","PB.B004.01.BS2"]               |
| Array{String, 2}  | ["PB" "B004" "01" "BS?"; "PB" "B001" "01" "BS?"]  |

The `LOC` field can be blank in FDSN requests with get_data; for example,
`chans="UW.ELK..EHZ"; get_data("FDSN", chans)`.

#### SeedLink only
For SeedLink functions (`seedlink!`, `has_stream`, etc.), channel IDs can
include a fifth field (i.e. NET.STA.LOC.CHA.T) to set the "type" flag (one of
DECOTL, for Data, Event, Calibration, blOckette, Timing, or Logs). Note that
SeedLink calibration, timing, and logs are not in the scope of SeisIO.

See also: `get_data`, `seedlink`
"""
function web_chanspec()
  return nothing
end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

"Abstract type for messages to actors."
abstract type Message end

"""
    Dispatch

Depending on its `Dispatch` mode an actor composes the arguments
to the behavior function:

- `full`: from the [`Become`](@ref) `args...` and the `msg.x...`.
    This is the default dispatch mode.
- `state`: from the actor state and the `msg.x...`. In this case the
    actor updates its state with the result of the behavior
    function. The result is saved in a `Tuple`.
"""
@enum Dispatch full state

"""
    Link{C}(chn::C, pid::Int, type::Symbol)

A mailbox for communicating with actors.

# Fields/Parameters
- `chn::C`:
- `pid::Int`:
- `type::Symbol`:
"""
struct Link{C}
    chn::C
    pid::Int
    type::Symbol
end


text:
# ---
# title: 888. Fair Candy Swap
# id: problem888
# author: Indigo
# date: 2021-02-01
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/fair-candy-swap/description/>
# hidden: true
# ---
#
# Alice and Bob have candy bars of different sizes: `A[i]` is the size of the
# `i`-th bar of candy that Alice has, and `B[j]` is the size of the `j`-th bar
# of candy that Bob has.
#
# Since they are friends, they would like to exchange one candy bar each so that
# after the exchange, they both have the same total amount of candy.  ( _The
# total amount of candy  a person has is the sum of the sizes of candy bars they
# have._)
#
# Return an integer array `ans` where `ans[0]` is the size of the candy bar that
# Alice must exchange, and `ans[1]` is the size of the candy bar that Bob must
# exchange.
#
# If there are multiple answers, you may return any one of them.  It is
# guaranteed an answer exists.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,1], B = [2,2]
#     Output: [1,2]
#
#
# **Example 2:**
#
#
#
#     Input: A = [1,2], B = [2,3]
#     Output: [1,2]
#
#
# **Example 3:**
#
#
#
#     Input: A = [2], B = [1,3]
#     Output: [2,3]
#
#
# **Example 4:**
#
#
#
#     Input: A = [1,2,5], B = [2,4]
#     Output: [5,4]
#
#
#
#
# **Note:**
#
#   * `1 <= A.length <= 10000`
#   * `1 <= B.length <= 10000`
#   * `1 <= A[i] <= 100000`
#   * `1 <= B[i] <= 100000`
#   * It is guaranteed that Alice and Bob have different total amounts of candy.
#   * It is guaranteed there exists an answer.
#
#
## @lc code=start
using LeetCode

function fair_candy_swap(A::Vector{Int}, B::Vector{Int})
    set_a, set_b = Set(A), Set(B)
    dif = (sum(A) - sum(B)) ÷ 2
    for num_a in set_a
        ((num_a - dif) in set_b) && return [num_a, num_a - dif]
    end
end
## @lc code=end


text:
function GaussNewtonOptimizer()
    return GaussNewtonOptimizer(())
end

function GaussNewtonOptimizer(arg0::MatrixDecomposer, arg1::jboolean)
    return GaussNewtonOptimizer((MatrixDecomposer, jboolean), arg0, arg1)
end

function get_decomposer(obj::GaussNewtonOptimizer)
    return jcall(obj, "getDecomposer", MatrixDecomposer, ())
end

function is_form_normal_equations(obj::GaussNewtonOptimizer)
    return jcall(obj, "isFormNormalEquations", jboolean, ())
end

function optimize(obj::GaussNewtonOptimizer, arg0::LeastSquaresProblem)
    return jcall(obj, "optimize", LeastSquaresOptimizer_Optimum, (LeastSquaresProblem,), arg0)
end

function to_string(obj::GaussNewtonOptimizer)
    return jcall(obj, "toString", JString, ())
end

function with_decomposer(obj::GaussNewtonOptimizer, arg0::MatrixDecomposer)
    return jcall(obj, "withDecomposer", GaussNewtonOptimizer, (MatrixDecomposer,), arg0)
end

function with_form_normal_equations(obj::GaussNewtonOptimizer, arg0::jboolean)
    return jcall(obj, "withFormNormalEquations", GaussNewtonOptimizer, (jboolean,), arg0)
end


text:
#=
# Normal shocks
In the last notebook, we saw that normal shocks develop in a nozzle when it is
not perfectly expanded. This notebook demonstrates the use of the compressible tools
for dealing with **normal shock waves**.
=#

# ### Set up the module
using Gasdynamics1D
#-
using Plots

#=
### General use of normal shock relations
Normal shock relations are designated with the argument `NormalShock`. For example,
let's look at the ratio of stagnation pressures $p_{02}/p_{01}$ across a normal shock
 with the Mach number $M_1$ entering the shock equal to 2:
=#
p02_over_p01 = StagnationPressureRatio(MachNumber(2),NormalShock)

# and the density $\rho_2/\rho_1$ increases by the factor:
ρ2_over_ρ1 = DensityRatio(MachNumber(2),NormalShock)

# What is the entropy change across this shock?
s2_minus_s1 = Entropy(Entropy(0),MachNumber(2),NormalShock)

# Let's plot the entropy increase as a function of Mach number:
M1 = range(1,3,length=101)
s2 = []
pratio = []
for M in M1
    push!(s2,ustrip(Entropy(Entropy(0),MachNumber(M),NormalShock)))
    push!(pratio,ustrip(PressureRatio(MachNumber(M),NormalShock)))
end
#-
plot(M1,s2,xlim=(1,3),ylim=(0,400),xlabel="Mach number",ylabel="Entropy increase (J/kgK)")
# Not surprisingly, the entropy jump gets bigger as the Mach number increases

#=
### Example 1
For air entering a shock at Mach number 3, 1 atm, and 50 degrees F, what are the
Mach number, pressure, and temperature on the other side of the shock?
=#
p1 = Pressure(1u"atm")
M1 = MachNumber(3)
T1 = Temperature(50u"°F")

# The Mach number exiting the shock is
M2 = MachNumber(M1,NormalShock)

# Pressure exiting the shock (in atm) is
p2 = Pressure(p1*PressureRatio(M1,NormalShock))
value(p2,u"atm")

# and the temperature exiting the shock (in F) is
T2 = Temperature(T1*TemperatureRatio(M1,NormalShock))
value(T2,u"°F")

#=
### Example 2
A blow-down supersonic windtunnel is supplied with air from a large reservoir.
A Pitot tube is placed at the exit plane of a converging-diverging nozzle. The test
section lies at the end of the diverging section. The Mach number in the test
section (M2) is 2 and the pressure is below atmospheric, so a shock is formed at
the exit. The pressure just after the shock (p3) is 14.7 psi. Find the pressure in
the reservoir (p01), the pressure in the throat (pth), the Mach number just after
the shock (M3), the pressure at the Pitot tube (p4), and the temperature at the Pitot
 tube (T4).
=#
M2 = MachNumber(2)
p3 = Pressure(14.7u"psi")

#=
First, let's find the pressure just before the shock and the Mach number
just after the shock
=#
# p3/p2
p3_over_p2 = PressureRatio(M2,NormalShock)

# p03/p02
p03_over_p02 = StagnationPressureRatio(M2,NormalShock)

# M3
M3 = MachNumber(M2,NormalShock)

# We can now calculate the pressure just before the shock, using $p_2 = p_3/(p_3/p_2)$
p2 = Pressure(p3/p3_over_p2)

#=
Now, let's find the conditions in the throat and reservoir, based on what we have
 been given here. We can immediately find the stagnation pressure upstream of the
 shock using the isentropic relation. This stagnation pressure is the same as the
 reservoir pressure.
=#
p02 = StagnationPressure(p2,M2,Isentropic)
p01 = p02

#=
We know that the throat is choked. It must be, because the flow is subsonic before
and supersonic after the throat. Therefore, $M_{th}$ is 1, and we can get the local
pressure ($p_{th}$) from knowing the stagnation pressure and Mach number there:
=#
Mth = MachNumber(1)
p0th = p02
pth = Pressure(p0th,Mth,Isentropic)

#=
Now we can calculate the conditions at the Pitot tube at the exit. At the nose of
the Pitot tube, the flow stagnates. Thus, the pressure and temperature are equal
to the stagnation values. The stagnation pressure at the exit is the same as the
stagnation pressure at point 3 (just after the shock):
=#
p03 = Pressure(p03_over_p02*p02)
#-
value(p03,u"psi")


text:
using PolynomialTestSystems
using Test

@testset "getters" begin
    F = chandra(5)
    @test length(equations(F)) == 5
    @test nvariables(F) == 5
    @test bezout_number(F) == 32
    @test multi_bezout_number(F) !== nothing
    @test mixed_volume(F) == 16

    @test nsolutions(heart()) == 4
    @test nreal_solutions(heart()) == 2
end

@testset "systems" begin
    for n=3:11
        @test chandra(n) isa TestSystem
        @test cyclic(n) isa TestSystem
        @test katsura(n) isa TestSystem
    end

    @test ipp() isa TestSystem
    @test ipp2() isa TestSystem
    @test fourbar() isa TestSystem
    @test rps10() isa TestSystem

    @test boon() isa TestSystem
    @test heart() isa TestSystem
    @test d1() isa TestSystem
    @test bacillus_subtilis() isa TestSystem
    @test griewank_osborne() isa TestSystem
    @test tritangents() isa TestSystem
    @test cyclooctane() isa TestSystem
end


text:
@recipe function f(model::Union{UnobservedComponents,
                                ExponentialSmoothing},
                  kf::FilterOutput{<:AbstractFloat})
    if !isunivariate(model)
        error("This plot recipe currently works for univariate models only.")
    end
    components = dict_components(model)
    att = get_filtered_state(kf)
    # Plot the series
    layout := (length(components) + 1, 1)
    sp = 1 # subplot
    @series begin
        seriestype := :path
        label := ""
        seriescolor := "black"
        title := "Observed"
        subplot := sp
        model.system.y
    end
    sp += 1
    for (component, idx) in components
        @series begin
            seriestype := :path
            label := ""
            seriescolor := "black"
            title := "Filtered " * component
            subplot := sp
            att[:, idx]
        end
        sp += 1
    end
end

@recipe function f(model::Union{UnobservedComponents,
                                ExponentialSmoothing},
                   ks::SmootherOutput{<:AbstractFloat})
    if !isunivariate(model)
        error("This plot recipe currently works for univariate models only.")
    end
    components = dict_components(model)
    alpha = get_smoothed_state(ks)
    # Plot the series
    layout := (length(components) + 1, 1)
    sp = 1 # subplot
    @series begin
        seriestype := :path
        label := ""
        seriescolor := "black"
        title := "Observed"
        subplot := sp
        model.system.y
    end
    sp += 1
    for (component, idx) in components
        @series begin
            seriestype := :path
            label := ""
            seriescolor := "black"
            title := "Smoothed " * component
            subplot := sp
            alpha[:, idx]
        end
        sp += 1
    end
end


text:
#!/usr/bin/env julia

using Klio
using Sockets

Klio.settings.server_host = @ip_str "0.0.0.0"
Klio.settings.expl_sqlite_file = joinpath(@__DIR__, "db/expl.sqlite")

Klio.run()


text:
module TestOptional
using StaticArrays
using Rotations
using CoordinateTransformations
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using Setfield
using PhaseSpaceIO: position

@testset "position, direction" begin
    dir = StaticArrays.normalize(@SVector(randn(3)))
    pos_iaea = @SVector(randn(3))
    pos_egs = @SVector(randn(2))
    for (p, pos, dir) in [
                          (arbitrary(EGSParticle{Float32}), pos_egs, dir),
                          (arbitrary(IAEAParticle{2,3}), pos_iaea, dir),
           ]
        q = @set position(p) = pos
        @test typeof(q) == typeof(p)
        @test q.E === p.E
        @test position(q) ≈ pos

        q = @set direction(p) = dir
        @test typeof(q) == typeof(p)
        @test q.E === p.E

        @test direction(q) ≈ dir
        @test dir ≈ [q.u, q.v, q.w]
    end

    p = arbitrary(EGSParticle{Nothing})
    @test position(p) == [p.x, p.y]
    @test position(p,z=3) == [p.x, p.y, 3]

    p = arbitrary(IAEAParticle{1,2})
    @test position(p) == [p.x, p.y, p.z]
end

@testset "coordinate transformations" begin
    p = arbitrary(IAEAParticle{1,2})

    t = Translation(randn(3))
    @test position(t(p)) ≈ t(position(p))
    @test direction(p) ≈ direction(t(p))

    l = LinearMap(rand(AngleAxis))
    @test position(l(p)) ≈ l(position(p))
    @test l(direction(p)) ≈ direction(l(p))

    tl = t ∘ l
    @test position(tl(p)) ≈ position(t(l(p)))
    @test direction(tl(p)) ≈ direction(l(p))
end

end#module


text:
# Wrap a Factorization in a AbstractSparseMatrix
struct OptimCholWrapper{T} <: SparseArrays.AbstractSparseMatrix{Float64, Int}
    data::T
end
function Base.similar(::OptimCholWrapper, ::Type{T}, n::Int, m::Int) where T
    @assert n == 0 == m
    zeros(0,0)
end
Base.:\(x::OptimCholWrapper, y::Vector{T}) where T = x.data\y
Base.copy(x::OptimCholWrapper) = x #OptimCholWrapper(copy(x.data))
Base.show(io::IO, X::OptimCholWrapper) = show(io, X.data)
LinearAlgebra.logdet(x::OptimCholWrapper) = logdet(x.data)

function lambertwexp(x)
    if x < 700
        LambertW.lambertw(exp(x))
    else
        x - log(x)
    end
end
function weightstonodes(locationweight)
    weightstonodes!(zeros(size(locationweight)), locationweight)
end
function weightstonodes!(A, locationweight, b=1)
    for (i, w) in locationweight
        # w is the sum of weights for the grid point.
        # The contributions of different datapoints may be combined or kept separate.
        A[i] += w*b
    end
    A
end

# Solution for the problem assuming that Q = qI
function initial_guess(fixed_field, locationweight, q)
    aa = weightstonodes!(zero(fixed_field), locationweight)

    @. aa = aa/q - lambertwexp(fixed_field - log(q) + aa/q)
end


text:
###################################################
## selection.jl
##      Build graph for scenario-based network flow method to solve full problem
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

"""
    Represents a node in scenario graph
    Required attributes:
        - id::Int           : unique reference of node
"""
abstract type ScenarioNode end

"""
    Node in scenario graph, represents buses arriving at the beginning of a route for a school
"""
struct ArrivalNode <: ScenarioNode
    id::Int

    "School associated with this node"
    school::Int
    "Route corresponding to the node"
    route::Route
    "Total service time of route"
    serviceTime::Float64
    "Scenario number"
    scenario::Int
end

"""
    Node in scenario graph, represents buses departing from a school
"""
struct DepartureNode <: ScenarioNode
    id::Int

    "School associated with this node"
    school::Int
    "Capacity, i.e. number of buses that can leave"
    capacity::Int
    "Scenario number"
    scenario::Int
end

"""
    Node in scenario graph, represents a bus yard
"""
struct YardNode <: ScenarioNode
    id::Int

    "Yard associated with this node"
    yard::Int
end

"""
    represents graph of scenarios for schools (master problem)
"""
mutable struct ScenarioGraph
    "graph object"
    graph::DiGraph
    "list of nodes"
    nodes::Vector{ScenarioNode}
    "dictionary mapping node type to node ids"
    nodeTypes::Dict{Symbol,Vector{Int}}
    "dictionary mapping pairs of nodes to costs"
    costs::Dict{Tuple{Int,Int}, Float64}
    "number of scenarios per school"
    numScenarios::Vector{Int}
end

"""
    Assumes the scenarios have already been computed in the data
    Builds ScenarioGraph using this input
"""
function buildScenarioGraph(data::SchoolBusData)
    if !data.withRoutingScenarios
        error("Cannot compute scenario graph without scenarios")
    end
    nodes = ScenarioNode[]
    nodeTypes = Dict(elt => Int[] for elt in [:YardNode, :ArrivalNode, :DepartureNode])
    # get the yard nodes
    yardNodes, currentNodeId = getYardNodes!(data, nodeTypes)
    append!(nodes, yardNodes)
    numScenarios = [length(data.scenarios[sID]) for sID=eachindex(data.schools)]
    for schoolID = eachindex(data.schools)
        for scenarioNum = 1:numScenarios[schoolID]
            schoolNodes, currentNodeId = getSchoolScenarioNodes!(data, schoolID,
                                                                 currentNodeId, scenarioNum,
                                                                 nodeTypes)
            append!(nodes, schoolNodes)
        end
    end
    graph = DiGraph(length(nodes))
    costs = Dict{Tuple{Int,Int}, Float64}()
    possibleEdgeTypes = [(:DepartureNode, :ArrivalNode), (:DepartureNode, :YardNode),
                         (:ArrivalNode, :DepartureNode), (:YardNode, :ArrivalNode)]
    for (type1, type2) in possibleEdgeTypes
        for idx1 in nodeTypes[type1], idx2 in nodeTypes[type2]
            node1 = nodes[idx1]
            node2 = nodes[idx2]
            edgeWasAdded, cost = createEdge!(data, graph, node1, node2)
            if edgeWasAdded
                costs[node1.id, node2.id] = cost
            end
        end
    end
    return ScenarioGraph(graph, nodes, nodeTypes, costs, numScenarios)
end

"""
    Get nodes for one scenario for a given school
    Args:
        data::SchoolBusData    : the data
        schoolId::Int          : the id of the school
        currentNodeId          : id of node, incremented for each new node
        scenarioNum            : id of scenario
        nodeTypes              : the attributes of the ScenarioGraph that we are updating
    Returns:
        the nodes for that scenario
        the id of the current node (to be passed along to the next call)
"""
function getSchoolScenarioNodes!(data::SchoolBusData, schoolID::Int,
                                 currentNodeId::Int, scenarioNum::Int,
                                 nodeTypes::Dict{Symbol, Vector{Int}})
    schoolNodes = ScenarioNode[]
    scenario = data.scenarios[schoolID][scenarioNum]
    for routeID in scenario.routeIDs
        route = data.routes[schoolID][routeID]
        # create arrival nodes
        push!(schoolNodes, ArrivalNode(currentNodeId, schoolID, route,
                                       serviceTime(data, schoolID, route),
                                       scenario.id))
        push!(nodeTypes[:ArrivalNode], currentNodeId)
        currentNodeId += 1
    end
    push!(schoolNodes, DepartureNode(currentNodeId,schoolID,length(scenario.routeIDs),scenarioNum))
    push!(nodeTypes[:DepartureNode], currentNodeId)
    currentNodeId += 1
    return schoolNodes, currentNodeId
end

"""
    Get nodes for yards
"""
function getYardNodes!(data::SchoolBusData, nodeTypes::Dict{Symbol,Vector{Int}})
    yardNodes = ScenarioNode[]
    nodeId = 1
    for yard in data.yards
        push!(yardNodes, YardNode(nodeId, yard.id))
        push!(nodeTypes[:YardNode], nodeId)
        nodeId += 1
    end
    return yardNodes, nodeId
end

"""
    Decide whether a bus from one school can serve another school (time feasible)
    By convention, school 1 is the school from which the bus would depart in the morning, and
        the route belings to school 2
"""
function isFeasibleInTime(data::SchoolBusData, school1::Int, school2::Int,
                          route::Route, routeTime::Float64)
    return (data.schools[school1].starttime +
            traveltime(data, data.schools[school1], data.stops[school2][route.stops[1]]) +
            routeTime + data.schools[school2].dwelltime <= data.schools[school2].starttime)
end

"""
    Add edge between two nodes and compute cost of the edge
"""
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::ScenarioNode, node2::ScenarioNode)
    edgeWasAdded = false
    cost = nothing
    return edgeWasAdded, cost
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::DepartureNode, node2::ArrivalNode)
    edgeWasAdded = false
    if node1.school != node2.school
        # can only go to another school and route must be feasible
        if isFeasibleInTime(data, node1.school, node2.school,
                            node2.route, node2.serviceTime)
            edgeWasAdded = add_edge!(graph, node1.id, node2.id)
        end
    end
    return edgeWasAdded, 0.
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::DepartureNode, node2::YardNode)
    return add_edge!(graph, node1.id, node2.id), 0.
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::ArrivalNode, node2::DepartureNode)
    edgeWasAdded = false
    if node1.school == node2.school && node1.scenario == node2.scenario
        edgeWasAdded = add_edge!(graph, node1.id, node2.id)
    end
    return edgeWasAdded, 0.
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::YardNode, node2::ArrivalNode)
    return add_edge!(graph, node1.id, node2.id), 1.
end

"""
    Wrapper function replacing LightGraphs.out_edges since it was deprecated
    Returns a generator with the desired edges
"""
function edges_out(graph::DiGraph, node::Int)
    return (Edge(node, j) for j in outneighbors(graph, node))
end

"""
    Wrapper function replacing LightGraphs.in_edges since it was deprecated
    Returns a generator with the desired edges
"""
function edges_in(graph::DiGraph, node::Int)
    return (Edge(j, node) for j in inneighbors(graph, node))
end

"""
    Flow method for scenario selection
    Args:
        - the data
        - the scenario graph we've constructed
    Returns:
        - a vector of length the number of schools, containing the index of
                the scenario that is used for each school
"""
function selectScenario(data::SchoolBusData, sg::ScenarioGraph; args...)
    edgeCost(edge::Edge) = sg.costs[src(edge), dst(edge)]

    model = Model(solver = GurobiSolver(Threads=getthreads(), MIPGap=0.01; args...))

    @variable(model, useScenario[i=eachindex(data.schools), j=1:sg.numScenarios[i]], Bin)
    @variable(model, busFlow[edge=edges(sg.graph)] >= 0, Int)

    @constraint(model, capacityDeparture[i=sg.nodeTypes[:DepartureNode]],
                sum(busFlow[edge] for edge=edges_out(sg.graph, i)) <= sg.nodes[i].capacity)
    @constraint(model, capacityArrival[i=sg.nodeTypes[:ArrivalNode]],
                sum(busFlow[edge] for edge=edges_out(sg.graph, i)) ==
                useScenario[sg.nodes[i].school, sg.nodes[i].scenario])

    @constraint(model, oneScenario[i=eachindex(data.schools)],
                sum(useScenario[i,j] for j=1:sg.numScenarios[i]) == 1)

    @constraint(model, flowConservation[i=vertices(sg.graph)],
                sum(busFlow[edge] for edge=edges_in(sg.graph, i)) ==
                sum(busFlow[edge] for edge=edges_out(sg.graph, i)))

    @objective(model, Min, sum(busFlow[edge] * edgeCost(edge) for edge=edges(sg.graph)))

    status = solve(model)
    scenarioUsed = [indmax(getvalue(useScenario[i,j])
                           for j=1:sg.numScenarios[i]) for i=eachindex(sg.numScenarios)]
    return scenarioUsed
end

"""
    Represents node in FullRoutingGraph
    Mandatory attributes:
    - id::Int        : unique identifier
    - yardId::Int    : yard number
"""
abstract type FullRoutingNode end

"""
    In full routing flow graph, represents a bus with a certain capacity
"""
struct BusNode <: FullRoutingNode
    id::Int
    yardId::Int
    "Route associated with node"
    route::Route
    "School associated with node"
    school::Int
    "Service time of route"
    serviceTime::Float64
end

"""
    In full routing flow graph, represents a yard
"""
struct FullYardNode <: FullRoutingNode
    id::Int
    yardId::Int
end

"""
    Represents the situation after the scenario has been picked: one bus per route per school
    Similar to ScenarioGraph
"""
mutable struct FullRoutingGraph
    "The graph itself"
    graph::DiGraph
    "list of node objects"
    nodes::Vector{FullRoutingNode}
    "dictionary mapping pairs of nodes to costs"
    costs::Dict{Tuple{Int,Int},Float64}
end

"""
    Given routes object (such as that output by flow master problem), set up full routing graph
"""
function buildFullRoutingGraph(data::SchoolBusData, usedScenario::Vector{Int})
    nodes = FullRoutingNode[]
    costs = Dict{Tuple{Int,Int},Float64}()

    currentId = 1
    # yard nodes
    for yard in data.yards
        push!(nodes, FullYardNode(currentId, yard.id))
        currentId += 1
    end
    # school bus nodes
    for schoolID in eachindex(data.schools)
        scenario = data.scenarios[schoolID][usedScenario[schoolID]]
        for routeID in scenario.routeIDs
            route = data.routes[schoolID][routeID]
            for yard in data.yards
                push!(nodes, BusNode(currentId, yard.id, route, schoolID,
                                     serviceTime(data, schoolID, route)))
                currentId += 1
            end
        end
    end
    graph = DiGraph(length(nodes))
    # edges
    for node1 in nodes, node2 in nodes
        edgeWasAdded, cost = createEdge!(data, graph, node1, node2)
        if edgeWasAdded
            costs[node1.id, node2.id] = cost
        end
    end
    return FullRoutingGraph(graph, nodes, costs)
end

"""
    Method creates an edge between two nodes, with one function per type of node pair
"""
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::FullRoutingNode, node2::FullRoutingNode)
    return false, 0.
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::FullYardNode, node2::BusNode)
    edgeWasAdded, cost = false, 0.
    if node1.yardId == node2.yardId
        edgeWasAdded = add_edge!(graph, node1.id, node2.id)
        cost = traveltime(data, data.yards[node1.yardId],
                          data.stops[node2.school][node2.route.stops[1]])
    end
    return edgeWasAdded, cost
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::BusNode, node2::FullYardNode)
    edgeWasAdded, cost = false, 0.
    if node1.yardId == node2.yardId
        edgeWasAdded = add_edge!(graph, node1.id, node2.id)
        cost = traveltime(data, data.schools[node1.school], data.yards[node2.yardId])
    end
    return edgeWasAdded, cost
end
function createEdge!(data::SchoolBusData, graph::DiGraph,
                     node1::BusNode, node2::BusNode)
    edgeWasAdded, cost = false, 0.
    if node1.yardId == node2.yardId && node1.school != node2.school
        if isFeasibleInTime(data, node1.school,
                            node2.school, node2.route, node2.serviceTime)
            edgeWasAdded = add_edge!(graph, node1.id, node2.id)
            cost = traveltime(data, data.schools[node1.school],
                              data.stops[node2.school][node2.route.stops[1]])
        end
    end
    return edgeWasAdded, cost
end

"""
    Given a full routing graph, constructs a YardId => node dictionary
"""
function yardToNodeDict(frg::FullRoutingGraph)
    d = Dict{Int,Int}()
    for node in frg.nodes
        if typeof(node) == FullYardNode
            d[node.yardId] = node.id
        end
    end
    return d
end

"""
    Bus scheduling problem
"""
function solveFullRouting(data::SchoolBusData, frg::FullRoutingGraph; args...)
    model = Model(solver = GurobiSolver(Threads=getthreads(), MIPGap=0.01; args...))
    # variables - morning
    @variable(model, busFlow[edge=edges(frg.graph)], Bin)
    length([node for node in frg.nodes if typeof(node) == FullYardNode]) > 1 &&
        error("Too many yards")
    # variables - yard capacity
    @variable(model, yardCapacity[y=eachindex(frg.nodes);
                                  typeof(frg.nodes[y]) == FullYardNode] >= 0, Int)
    # a route can be served by at most one bus from a given yard - morning
    @constraint(model, oneBusPerRoute[nodeId=eachindex(frg.nodes);
                                      typeof(frg.nodes[nodeId]) == BusNode],
                sum(busFlow[edge] for edge = edges_in(frg.graph, nodeId)) == 1)
    # flow conservation constraints - morning
    @constraint(model, flowConservation[nodeId=vertices(frg.graph)],
                sum(busFlow[edge] for edge=edges_in(frg.graph, nodeId)) ==
                sum(busFlow[edge] for edge=edges_out(frg.graph, nodeId)))
    @constraint(model, flowBound[y=eachindex(frg.nodes);
                                 typeof(frg.nodes[y]) == FullYardNode],
                sum(busFlow[edge] for edge=edges_out(frg.graph, y)) <= yardCapacity[y])
    # minimize the total cost of buses (First Order = total number of bus)
    @objective(model, Min, sum(yardCapacity[y]
                               for (y,node)=enumerate(frg.nodes) if typeof(node) == FullYardNode))
    status = solve(model)
    flows = Dict(edge => getvalue(busFlow[edge]) for edge=edges(frg.graph))
    return interpretFlows!(data, frg, flows)
end

"""
    Given flows, get final buses
"""
function interpretFlows!(data::SchoolBusData, frg::FullRoutingGraph, flows::Dict)
    finalBuses = Bus[]
    currentBusId = 1
    dict = yardToNodeDict(frg)
    for yard in data.yards
        yardNode = dict[yard.id]
        yardIsEmpty = false
        # keep following bus routes until yard is empty
        while !yardIsEmpty
            schools, routes = followBusAlongRoute!(frg, flows, yardNode)
            if isempty(schools)
                yardIsEmpty = true
            else
                push!(finalBuses, Bus(currentBusId, yard.id, schools, routes))
                currentBusId += 1
            end
        end
    end
    return finalBuses
end

"""
    Extract a bus route from a flows object, removing flows as they are found
"""
function followBusAlongRoute!(frg::FullRoutingGraph, flows::Dict, yardNode::Int)
    schools, routes = Int[], Int[]
    if yardNode == -1 && length(flows) == 0
        return schools, routes
    end
    currentNode = yardNode
    backToYard = false
    while !backToYard
        for edge in edges_out(frg.graph, currentNode)
            if flows[edge] > 0.5
                if typeof(frg.nodes[dst(edge)]) == BusNode
                    push!(schools, frg.nodes[dst(edge)].school)
                    push!(routes, frg.nodes[dst(edge)].route.id)
                    currentNode = dst(edge)
                else
                    backToYard = true
                end
                flows[edge] = 0.
                break
            end
        end
        if isempty(schools) # couldn't find anything leading out of the yard
            backToYard = true
        end
    end
    return schools, routes
end

function routeBuses!(data::SchoolBusData; args...)
    sg = buildScenarioGraph(data)
    data.usedScenario = selectScenario(data, sg; args...)
    frg = buildFullRoutingGraph(data, data.usedScenario)
    data.buses = solveFullRouting(data, frg; args...)
    data.withFinalBuses = true
    return data
end

max_stars_repo_name: jihoyeo/SchoolBusRouting

max_issues_repo_name: jihoyeo/SchoolBusRouting

max_forks_repo_name: jihoyeo/SchoolBusRouting


text:
import Base.@kwdef

@kwdef mutable struct Person
    name::String
    address::String
end

const people = Person[]

Jus.present(people;
            Person = (;
                      # specify what to gen
            ))